Supercritical Injection Modeling by an Incompressible but Variable Density Approach

Abstract

Supercritical nitrogen jet behavior is modeled using an incompressible but variable density approach developed for variable density jets. Following mechanical and thermal breakup concepts, several injection conditions relevant to liquid rocket propulsion are analyzed, considering heat transfer in the injector. Regarding axial density distributions, different levels of agreement with experimental data are encountered for potential core, subsided core, and plateau formations. Further comparisons with compressible formulations from the literature are a good indicator of the proposed methodology’s suitability for the simulation of supercritical injection behavior.

1. Introduction

While pressure and temperature increase in the combustion (or thrust) chambers of LREs—liquid rocket engines—increase the rocket’s specific impulse, fuel, and oxidizer enter a supercritical state characterized by both pressure and temperature above their critical point values [1]. Substantially different from their subcritical counterparts [2,3,4], supercritical fluids exhibit peculiar behaviors not encountered at subcritical conditions. While mass diffusivity, surface tension, and latent heat are zero at the critical point, isentropic compressibility, specific heat, and thermal conductivity diverge to infinity. The critical point also marks the end of the discontinuity between liquid and gas phases, being where the single-term fluid is widespread. Consequently, fluids at these conditions have liquid-like density and gas-like properties governed by mass diffusion [5,6]. Recently [7,8], it was demonstrated that a phenomenological difference exists inside the supercritical regime when the Widom line is crossed. This division line, commonly labeled in the literature as the pseudo-boiling line, is identified by a peak in the isobaric specific heat, which separates liquid-like (high-density) states from gas-like (low-density) states. The effect of pseudo-boiling is more pronounced the closer one gets to the critical point. In order to have a clearer view over this “pseudo-state” transition, in Figure 1, the isobaric specific heat of nitrogen is depicted for several pressure levels (1–10 MPa), with data available from the NIST (National Institute of Standards) database [9].

As shown in Figure 1, the isobaric specific heat of nitrogen, $c_p$, is represented as a function of the temperature, T. With a critical pressure of 3.4 MPa and a critical temperature of 126.2 K, pseudo-boiling behavior is more evident in proximity to the critical point of nitrogen, namely for the pressure level of 4 MPa, which dissipates as the pressure is further increased. Essentially dividing the supercritical regime in two, the Widom line establishes a transition [10,11] from liquid- to gas-like conditions, resembling subcritical boiling. The transition across the Widom line is only physically relevant for $p_c<p<3p_c$, where the locus of the maxima in isobaric-specific heat can track the phenomenon of pseudo-boiling. The reason for approximating the Widom line to the locus of the maxima in isobaric-specific heat is addressed and justified by the lack of correlation lengths to estimate the microscopic interface [12]. If this state is heated to $T>2.5T_c$, the supercritical fluid behaves as an ideal gas. The Widom line and the pseudo-boiling phenomena have been experimentally verified for water [10]. However, a quantitative theory of pseudo-boiling has been proposed [7], where it is explained that, at supercritical conditions and in contrast to the subcritical case, the energy supplied to the fluid is used to overcome molecular attraction and to raise its temperature, corresponding to structural and thermal contributions. This differentiation will ultimately impact the jet evolution, owing to the crossing of the Widom line and the associated pseudo-boiling temperature. This temperature then replaces the critical temperature as the transition criterion. Through the analysis of Figure 1 and Figure 2, it can be concluded pseudo-boiling leads to the storage of a large amount of heat with only a marginal increase in temperature. The heat transfer contributes to the fluid’s expansion through its volume. In this sense, a variation of a couple of Kelvin leads to a significant density variation that needs to be considered in any successful numerical effort attempting to replicate the behavior of supercritical fluid flows. High levels of pseudo-boiling, such as those registered close to the critical point, contribute to jet stabilization, hindering the development of Kelvin–Helmholtz instabilities, playing a similar role to surface tension at subcritical conditions [13], with the terminology “thermal shield” being used. In addition, a DNS-direct numerical simulation study [14] highlighted the difficulties associated with the higher levels of grid refinement needed to capture the pseudo-boiling phenomena. While the continuous development of numerical tools [15,16,17,18,19,20,21,22] capable of predicting supercritical fluid flow behavior increasingly incorporates previously unknown or undefined physical phenomena, there is also a need for broad data sets of experiments against which the models themselves can be validated.

Quantitative cold flow data provide insight into jet mixing at supercritical conditions without introducing the complexities of combustion. Raman scattering measurements on a cryogenic nitrogen jet [23] injected into a quiescent nitrogen environment with pressure above nitrogen’s critical point have become a canonical test case for validating numerical solvers. It remains the most comprehensive source of information for single-species injection at supercritical conditions. Albeit, being a highly sought-after validation tool by many researchers [14,17,21,24,25,26,27,28,29,30,31,32,33,34,35], the focus has been to demonstrate the behavior of supercritical fluid flows in proximity to the critical point for a couple of test cases. Recently, the need to account for heat transfer inside the injector in the simulations for nitrogen jets has been demonstrated [18], to obtain the characteristic shape of a disintegrated core observed at gas-like conditions. This is sustained on thermal disintegration [18] in opposition to a purely mechanical view of supercritical jet disintegration. Through the concept of thermal disintegration, it was possible to evaluate [7] if the amount of energy received by the jet in the injector would be sufficient to overcome the Widom line and transition to gas-like conditions. This analysis would later be extended to a broader range of conditions [36]. However, we also need to consider uncertainties inherent to the experimental data, namely in the initial conditions [37], which lead to a search not for a quantitative matching between proposed numerical tools and experimental results but a qualitative replication of experiments.

Experimental studies on the injection of nitrogen/carbon dioxide mixtures at critical conditions using shadowgraphy [38] concluded that the jet could be treated as a variable density jet at pressures ranging from sub- to supercritical conditions, as long as the temperature remained above the critical point. Quantitative evidence in favor of the incompressible variable-density hypothesis was provided through the comparison of jet growth rates [39] at and above the critical point conditions with the theory of incompressible variable-density gaseous mixing layers. Therefore, the jet spreading rate became the first quantitative parameter that indicated the similarity between both regimes went beyond a qualitative physical appearance. Lastly, the quantitative measurements of density and decay rates that serve as a validation for numerical solvers [23,40,41] provide further evidence in support of the hypothesis. Using real gas relationships for transport and thermodynamic properties allows the physical model to capture the weak compressibility effects when using an incompressible variable-density approach.

As reviewed by Banuti [42], the seemingly overwhelming variety of computational fluid dynamics (CFD) codes is misleading, because successful approaches and methods are used as almost canonical procedures. Given the amount of experimental evidence, the variable-density behavior of supercritical jets can then be modeled through the usage of real gas relationships for density (in the form of an EoS) and transport properties such as dynamic viscosity and thermal conductivity, detailing ideal gas behavior and departure functions to account for high-pressure effects. Incompressible but variable density formulations have also started to be used recently in DNS simulations, such as in the study of Ries et al. [28], where jet disintegration was studied with a focus on the pseudo-boiling effects. Moreover, the different stages of jet disintegration were identified and related to mixing efficiency. However, the Reynolds numbers considered were very small and of little interest in the present practical context due to computational limitations. In contrast, the numerical modeling of supercritical fluid flows has been traditionally conducted by resorting to compressible solvers. However, as pointed out by Lapenna and Creta [43], flows at such conditions fall below the significant Mach number threshold of 0.3 and, in conjunction with the low injection velocities, will impose severe restrictions on the numerical solver. The use of the incompressible but variable density jet approach resulted from a hypothesis from the similarity of their visual appearance with the supercritical ones. The assumption was tested and evaluated by Barata et al. [44] and confirmed in other, broader situations given the previous encouraging results.

In the present manuscript, based on the thermal breakup mechanism, we intend to evaluate the suitability of the incompressible, but variable density approach to model supercritical fluid flows, comparing the obtained results with those of compressible formulations available from the literature. Several configurations [18]—potential core, subsided core, and plateau—are considered and evaluated. Focus is particularly given to injection conditions where the plateau [18] is experimentally observed, extending the work carried out by Banuti and Hannemann [18].

The mathematical model used is described in the Section 2, followed by the Section 3 and accurate modeling of Section 4 and Section 5. Next, a brief description of the numerical treatment is undertaken in Section 6, followed by the analysis of initial and boundary conditions and their relation to the experimental test matrix in Section 7. Finally, the obtained results are analyzed for low and high pseudo-boiling conditions in Section 8. The main conclusions are summarized in Section 9.

2. Governing Equations

The Favre-averaged conservation equations for mass, momentum, and energy are given in Equations (1)–(3), respectively, with $\rho$ being the density, $u_i$ the velocity vector, $x_j$ the coordinate vector, p the pressure, h the enthalpy, ${\tau }_{ij}$ the viscous stress tensor, and $q_j$ the heat flux, written following an incompressible but variable density approach [44].

$$\frac{\partial \overline{\rho }}{\partial t}+\frac{\partial \overline{\rho }{\tilde{u}}_i}{\partial x_i}=0$$

(1)

$$\frac{\partial \overline{\rho }{\tilde{u}}_i}{\partial t}+\frac{\partial \overline{\rho }{\tilde{u}}_i{\tilde{u}}_j}{\partial x_j}=-\frac{\partial \overline{p}}{\partial x_i}+\frac{\partial \overline{{\tau }_{ij}}}{\partial x_j}-\frac{\partial \overline{\rho }\widetilde{u_i^{\prime \prime }{u}_j^{\prime \prime }}}{\partial x_j}$$

(2)

$$\frac{\partial \overline{\rho }{\tilde{u}}_j\tilde{H}}{\partial x_j}=\frac{\partial \overline{{\tau }_{ij}{u}_i}}{\partial x_j}-\frac{\partial {\tilde{u}}_i\overline{\rho }\widetilde{u_i^{\prime \prime }{u}_j^{\prime \prime }}}{\partial x_j}-\frac{\partial \left({\overline{q}}_j+\overline{\rho }\widetilde{u_j^{\prime \prime }{h}^{\prime \prime }}\right)}{\partial x_j}$$

(3)

The Boussinesq hypothesis of Equation (4) is used to model the Reynolds stress tensor, with $\kappa$ representing the turbulence kinetic energy, ${\delta }_{ij}$ the Kronecker’s delta function, and ${\mu }_t$ the eddy viscosity.

$$-\overline{\rho }\widetilde{u_i^{\prime \prime }{u}_j^{\prime \prime }}=-\frac{2}{3}\rho k{\delta }_{ij}+{\mu }_t\left(\frac{\partial {\tilde{u}}_i}{\partial x_j}+\frac{\partial {\tilde{u}}_j}{\partial x_i}-\frac{2}{3}\frac{\partial {\tilde{u}}_k}{\partial x_k}{\delta }_{ij}\right)$$

(4)

The turbulence heat flux is given in Equation (5), with ${\mathrm{Pr}}_t$ representing the turbulence Prandtl number, framed analogously to its laminar counterpart, given by Fourier’s law.

$$\overline{\rho }\widetilde{u_j^{\prime \prime }{h}^{\prime \prime }}=-\frac{c_p{\mu }_t}{{\mathrm{Pr}}_t}\frac{\partial \tilde{T}}{\partial x_j}=-\frac{{\mu }_t}{{\mathrm{Pr}}_t}\frac{\partial \tilde{h}}{\partial x_j}$$

(5)

Turbulent stresses $\widetilde{u_i^{\prime \prime }{u}_j^{\prime \prime }}$ are modeled by the $\kappa -ϵ$ turbulence model [45]. Eddy viscosity is given according to Equation (6), while transport equations for the turbulence kinetic energy ($\kappa$) and its dissipation ($ϵ$) are modeled according to Equations (7) and (8), respectively. The turbulence quantities are written in terms of mean velocity components, effective viscosity, and density. So, any changes made to these variables by other models are updated continuously. The model constants are defined as ${\sigma }_{\kappa }=1.0$, ${\sigma }_{ϵ}=1.3$, $C_{\mu }=0.09$, $C_{ϵ1}=1.44$, and $C_{ϵ2}=1.92$. In the past, several studies were conducted on the performance of several RANS turbulence closures [46,47] where no correlation was found between model complexity and accuracy.

However, the development of tailored turbulence models for supercritical conditions was also attempted [48,49,50,51], while the quantification of turbulence model uncertainty has also been gaining relevance [52]. Since the turbulence model is not valid up to the solid walls, wall functions set the $y^{+}$ value in the first grid cell to 11.63 to ensure the boundary layer is adequately resolved. Additionally, 20 points are taken inside the boundary layer, and turbulence intensity is set to 5%.

$${\mu }_t=C_{\mu }\frac{\overline{\rho }{\kappa }^2}{ϵ}$$

(6)

$$\frac{\partial }{\partial x_i}\left(\overline{\rho }{\tilde{u}}_i\kappa \right)={\tau }_{ij}\frac{\partial {\tilde{u}}_j}{\partial x_i}-\overline{\rho }ϵ+\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_{\kappa }}\right)\frac{\partial \kappa }{\partial x_j}\right]$$

(7)

$$\frac{\partial }{\partial x_i}\left(\overline{\rho }{\tilde{u}}_iϵ\right)=C_{ϵ1}\frac{ϵ}{\kappa }{\tau }_{ij}\frac{\partial {\tilde{u}}_j}{\partial x_i}-C_{ϵ2}\overline{\rho }\frac{{ϵ}^2}{\kappa }+\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_{ϵ}}\right)\frac{\partial ϵ}{\partial x_j}\right]$$

(8)

3. Equation of State

For the compromise cubic equations of state represent between accuracy and computational cost, they are an attractive choice for describing supercritical fluid flows [53]. In the present work, the Peng–Robinson [54] equation of state, whose formulation is represented in Equation (9), is used. Here, $\mathcal{R}=8.314{\mathrm{Jmol}}^{-1}{\mathrm{K}}^{-1}$ corresponds to the gas constant and $v_m$ is the molar volume, with a (Equation (10)) and b (Equation (11)) accounting for the effect of intermolecular attractive and repulsive forces. Accordingly, $\alpha$ is defined following Equation (12) and $\kappa$, dependent upon the acentric factor $\omega$, is defined in Equation (13).

$$p=\frac{\mathcal{R}T}{v_m-b}-\frac{a\alpha }{v_m^2+2v_{m}b-b^2}$$

(9)

$$a=\frac{0.45724{\mathcal{R}}^2{T}_c^2}{p_c}$$

(10)

$$b=\frac{0.0778\mathcal{R}{T}_c}{p_c}$$

(11)

$$\alpha ={\left[1+\kappa \left(T_r^{1/2}\right)\right]}^2$$

(12)

$$\kappa =0.37464+1.54226\omega -0.26992{\omega }^2\omega \le 0.49$$

(13)

4. Thermodynamic Properties

Thermodynamic properties, such as enthalpy in Equation (14), are determined by combining the ideal gas property denoted by the 0 subscript added with a dense fluid correction for high-pressure effects dependent upon the equation of state used. The ideal component is evaluated by the seven-coefficient NASA polynomial [55] whose coefficients $a_1$ to $a_6$ are given in

Table 1. Nitrogen coefficients of the polynomial for the determination of ideal gas enthalpy [55].

Coefficient	Value
$a_1$	0.24159429 × 101
$a_2$	0.17489065 × 10−3
$a_3$	−0.11902369 × 106
$a_4$	0.30226245 × 10−10
$a_5$	−0.20360982 × 10−14
$a_6$	0.56133773 × 105

. In addition, the departure function for enthalpy is given in Equation (16) [56], while the variation of $\alpha$ with temperature is evaluated according to Equation (17).

$$h(p,T)=h_0\left(T\right)+{\int }_{p_0}^p{\left[\frac{1}{\rho }+\frac{T}{{\rho }^2}{\left(\frac{\partial p}{\partial T}\right)}_{\rho }\right]}_{T}dp$$

(14)

$$h_0\left(T\right)=\left(a_1+a_2\frac{T}{2}+a_3\frac{T^2}{3}+a_4\frac{T^3}{4}+a_5\frac{T^4}{5}+\frac{a_6}{T}\right)\mathcal{R}T$$

(15)

$$\frac{h(p,T)-h_0\left(T\right)}{\mathcal{R}}=\frac{a\alpha \left(T\right)\frac{\partial \alpha }{\partial T}}{\mathcal{R}T\sqrt{8b^2}}ln\left(\frac{2v_m+2b-\sqrt{8b^2}}{2v_m+2b+\sqrt{8b^2}}\right)-1+Z$$

(16)

$$\frac{\partial \alpha }{\partial T}=-\frac{\kappa \alpha \left(T\right)}{T_c{T}_r^{1/2}}$$

(17)

5. Transport Properties

Dynamic viscosity and thermal conductivity are evaluated following the same methodology used for the thermodynamic properties [57]. The values at ideal gas conditions are added contributions from departure functions that account for the high-pressure effects. Furthermore, a third contribution is considered for the particular case of thermal conductivity, resulting from the critical divergence; it reaches an infinite value at the critical point. Viscosity [57] is expressed according to Equation (18), evaluated from a dilute gas contribution, ${\mu }^0\left(T\right)$ (Equation (19)), and a residual component, ${\mu }^r(\tau ,\delta )$, (Equation (20)). The $\sigma$ represents the Lennard–Jones parameter and $\Omega$ the collision integral, while the remaining coefficients are constants given in

Table 2. Nitrogen coefficients for the collision integral equation [57].

i	${\mathit{b}}_{\mathit{i}}$
0	0.431
1	−0.4623
2	0.08406
3	0.005341
4	−0.00331

and

Table 3. Nitrogen coefficients and exponents of the residual fluid viscosity equations [57].

i	${\mathit{N}}_{\mathit{i}}$	${\mathit{t}}_{\mathit{i}}$	${\mathit{d}}_{\mathit{i}}$	${\mathit{l}}_{\mathit{i}}$
1	10.72	0.1	2	0
2	0.03989	0.25	10	1
3	0.001208	3.2	12	1
4	−7.402	0.9	2	2
5	4.620	0.3	1	3

. Reduced density and temperature are given by $\delta$ and $\tau$, respectively. Coefficients $b_i$, pertaining to the collision integral, are reproduced for the sake of completeness in Table 2, while the coefficients and exponents of Equation (20) are given in Table 3.

$$\mu ={\mu }^0\left(T\right)+{\mu }^r(\tau ,\delta )$$

(18)

$${\mu }^0\left(T\right)=\frac{0.0266958\sqrt{MT}}{{\sigma }^2\Omega \left(T^{*}\right)}$$

(19)

$${\mu }^r(\tau ,\delta )=\sum _{i=1}^n{N}_i{\tau }^{t_i}{\delta }^{d_i}exp(-\gamma {\delta }^{l_i})$$

(20)

Similarly, thermal conductivity is defined following Equation (21), with dilute and residual gas contributions given by Equations (22) and (23), respectively. Model constants are reproduced in

Table 4. Nitrogen coefficients and exponents of the residual fluid thermal conductivity equations [57].

i	${\mathit{N}}_{\mathit{i}}$	${\mathit{t}}_{\mathit{i}}$	${\mathit{d}}_{\mathit{i}}$	${\mathit{l}}_{\mathit{i}}$
1	1.511	-	-	-
2	2.117	−1.0	-	-
3	−3.332	−0.7	-	-
4	8.862	0.0	1	0
5	31.11	0.03	2	0
6	−73.13	0.2	3	1
7	20.03	0.8	4	2
8	−0.7096	0.6	8	2
9	0.2672	1.9	10	2

. Here, a third component [58] is needed for evaluating thermal conductivity in the critical region due to the critical divergence of thermal conductivity at the critical point, as defined in Equation (24). The ${\Omega }_i$ coefficients are calculated from the definition of specific heat at constant pressure and volume [59].

$$\lambda ={\lambda }^0\left(T\right)+{\lambda }^r(\tau ,\delta )+{\lambda }^c(\tau ,\delta )$$

(21)

$${\lambda }^0=N_1\left[\frac{{\mu }^0\left(T\right)}{1^{-}\mathrm{Pa}.\mathrm{s}}\right]+N_2{\tau }^{t_2}+N_3{\tau }^{t_3}$$

(22)

$${\lambda }^r=\sum _{i=4}^n{N}_i{\tau }^{t_i}{\delta }^{d_i}exp(-{\gamma }_i{\delta }^{l_i})$$

(23)

$${\lambda }^c=\rho C_p\frac{K{R}_{0}T}{6\pi \xi \eta (T,\rho )}(\tilde{\Omega }-{\tilde{\Omega }}_0)$$

(24)

Figure 3, Figure 4 and Figure 5 detail the Peng and Robinson [54] EoS and Lemmon and Jacobsen [57] models for dynamic viscosity and thermal conductivity validation with the NIST [9] database.

6. Computational Methods

The governing equations are discretized algebraically through a finite volume/finite difference method. Conservation at the discrete level is automatically ensured in the finite volume method due to using the conservation equations directly in the integral formulation. Finite differences are used at the faces of the control volumes.

A finite volume method is employed, with diffusion modeled by the second-order central scheme, while advective terms follow the QUICK scheme [60]. In this way, the appearance of non-physical pressure oscillations related to the formulation of the equation of state and the transition across the Widom line is mitigated [20] by dampening the stable central scheme with the upwinding of the QUICK scheme [44].

A staggered grid formulation prevents the decoupling between pressure and velocity. Velocity and pressure values are stored in different positions for which the control volumes are no longer the same. Ultimately, the pressure values are calculated directly for the cell face, and no interpolation is needed. This eliminates the decoupling of the pressure and velocity fields and any possible oscillations, assuming particular relevance for the use of a pressure-velocity formulation, in which the conservation equation for mass (Equation (1)) is forfeited in favor of the development of a conservation equation for pressure [26], constructed from the divergence of the momentum conservation equation. A study comparing density- and pressure-based formulations [61] indicated that no foreseeable differences were found between the two formulations. The highly coupled behavior between pressure and velocity is then translated into a need to retrieve both values simultaneously. After this, the conservation equation for energy and the turbulence-related transport variables are solved until convergence is reached.

Lastly, high levels of under-relaxation are used to accelerate the solution’s convergence. Furthermore, it is practical since it can help avoid the divergence of the iterative solution in strongly non-linear equations.

7. Experimental Conditions Analysis

The experimental conditions [23] for which the present numerical computations are validated are given in

Table 5. Experimental test matrix [23].

Case	${\mathit{v}}_{\mathit{inj}}\left[{\mathbf{ms}}^{-1}\right]$	${\mathit{p}}_{\mathit{ch}}$ [MPa]	${\mathit{T}}_{\mathit{inj}}$ [K]	$\mathbf{\Delta }{\dot{\mathit{H}}}_{\mathit{pb}}$ [W] [18]
5	2.0	4.90	131.0	103.59
6	1.9	4.90	140.0	−113.30
7	4.5	5.01	126.2	448.72
8	4.9	5.00	135.7	−51.88
9	2.0	5.85	135.0	81.22
10	1.9	5.88	140.5	−34.27

. The pseudo-vaporization powers [18] are also given ($\Delta {\dot{H}}_{pb}$) based on the experimental injection temperature. High pseudo-vaporization power is related to forming a dense potential core [18]. In contrast, low values result in a subsided core formation, while plateau formation occurs at intermediate levels. In Figure 6, the boundary conditions (BCs) are represented. As depicted in the figure, the injector and combustion chamber are coupled in the same simulation. While a localized refinement of the meshes is needed for an as-smooth-as-possible transition, this prevents the need to perform two different computations, one for the injector, another to the chamber, and the subsequent matching at $x/d=0$, x being the axial distance measured from the injector exit plane and d its diameter.

Constant velocity and temperature profiles are imposed at the inlet according to the experimental data in Table 5. Values for the pressure in the chamber are also given in Table 5, while chamber temperature is kept constant between the computations at 298 K. The face-plate of the chamber is set as an adiabatic wall, while the injector and chamber walls are considered to be isothermal, with imposed temperatures of 298 K [18]. Due to the experimental apparatus’s configuration, only half of the domain is considered, with the symmetry axis where normal velocity components and radial derivatives of the other variables are zero. Furthermore, gradients of dependent variables are set to zero at the outlet, ensuring global mass conservation and a pressure conservation equation solution. The injector has a diameter of 2.2 mm and a length of 90 mm, with the combustion chamber measuring 250 mm in length and 122 mm in diameter.

Figure 7 depicts the experimental test cases considered in the present work, with pressure levels approximate to those in the chamber and injection temperatures.

8. Results

Figure 8 depicts the axial density distribution obtained for experimental test case 7 as a function of the normalized distance from the injector. A dense potential core can be observed in the results over a distance of approximately nine injector diameters, in close agreement with the experimental data [23]. From this point forward, the large densities present start to have a preponderant role in the jet dynamics, leading to the break of nitrogen from the core and nitrogen entrainment from the chamber into the jet. Overall the obtained numerical results can replicate the experimental data. The results are compared with the LES computations of Müller et al. [61], who compared pressure- and density-based approaches following a compressible formulation.

A grid independence study was performed (Figure 8) with three refinement levels to ensure that the grid resolution does not affect results. First, a coarse mesh with 180,000 points, a standard mesh with 280,000 points, and a refined mesh with 495,000 points are tested by considering the $y^{+}$ restrictions. Since no slope variation is observed between the three mesh iterations, the standard mesh is used in the computations.

In Figure 9, we can observe how the error decreases as a function of the grid convergence index (GCI) [62] from the coarse to the refined mesh. From the slope, it is possible to retrieve the achieved order of convergence, 2.512, versus a nominal convergence rate of 3, being first- and second-order slopes represented for comparative purposes. The difference between the nominal and achieved orders of convergence is explained due to initial and boundary conditions and the discrete representation of the physical domain.

Moving on to the remaining injection conditions, it is possible to retrieve the axial density distribution for the injection configuration of case 8 in Figure 10. Here, a subsided core [18] is expected ($\Delta {\dot{H}}_{pb}$ = −51.88 W), meaning that the jet needs to receive less energy in the injector to trigger the pseudo-phase transition. Several striking differences are observed concerning the case 7 configuration. First, a small region of constant density is observed over a distance of two injector diameters, followed by a slight decrease in density until an $x/d$ of 8. Second, a faster rate of density decay is observed for its experimental counterpart, owing to the thermal stratification resulting from the heat exchange between the injector and gas-like nitrogen. Finally, the entrainment of chamber nitrogen into the jet is responsible for the sharper decrease in density after a distance of 10 injector diameters until the end of the domain, corresponding to 30 diameters. Thus, while an over-prediction of the axial density values is observed, the proposed numerical solver can replicate the qualitative behavior of the jet, namely in terms of predicting different density rates of decay. From all the considered experimental cases, case 8 is more susceptible to variations in the injection temperature due to its proximity to the peak in isobaric-specific heat. Moreover, no other formulation was found in the literature against which the present computations could be compared.

So far, two different configurations have been considered: the high pseudo-vaporization condition of case 7, characterized by the appearance of a dense potential core, and the low pseudo-vaporization configuration of case 8, influenced by the thermal stratification occurring in the injector and leading to a decay in density at the beginning of the chamber portion and the consequent formation of a subsided core. However, it may be the case that the pseudo-vaporization power may be intermediate between these two cases leading to the formation of an intermediate plateau as in cases 5, 6, 9, and 10 (Table 5).

Figure 11 depicts the obtained axial density evolution of cases 5 and 9, which have relatively pseudo-vaporization powers: 103.59 W for case 5 and 81.22 W for case 9. Accordingly, the experimental data [23] indicate a similar jet development, forming intermediate plateaus. However, the numerical predictions cannot replicate the experimental plateau formation, instead describing a behavior similar to the one encountered for low pseudo-vaporization power (case 8 in Figure 10). Comparing these results with the compressible formulation of Gopal et al. [63] in Figure 12, where, in addition to the experimental injection temperature of 135 K, the authors considered another injection temperature of 142.5 K, some more insights can be gathered. These do not, however, consider the injector heat transfer, which leads to the appearance of dense potential cores, as indicated in the figure, extending up to a distance of 10 injector diameters downstream of the injector exit plane. An incompressible but variable density approach with an adiabatic injector configuration (not shown) yielded similar results to those of Gopal et al. [63].

All in all, no definitive conclusion can be drawn in relation to the numerical capture of the plateau. Moreover, employing an incompressible but variable density approach or a compressible formulation leads to a similar prediction. To the best of the authors’ knowledge, no further attempt to replicate the experimental plateau has been published in the literature.

Lastly, it remains to consider the intermediate pseudo-vaporization powers of cases 6 and 10, to the right of the peak in isobaric-specific heat (Figure 7). These are depicted in Figure 13, where trends similar to those of Figure 11 are encountered. Comparing the qualitative behavior of the numerical predictions and the experimental data, a low decay rate is present over a distance of 8 injector diameters. At this point, a portion of supercritical nitrogen starts to break from the jet, and chamber fluid entrains into it, increasing the axial decay of density for the remaining domain.

Generally, it is possible to observe the formation of the experimental intermediate plateau at a distance between 5 and 10 injector diameters, a feature that the numerical solver is incapable of reproducing. The comparison of cases 6 and 10 in Figure 13 highlights a precise plateau formation in case 10. Of all the considered experimental cases, this is the one further away from the critical point, where the effects of pseudo-boiling are the weakest and where the most well-defined plateau is observed.

While the intermediate plateau could not be recovered numerically, the results highlight the similarity between the experimental conditions to the left and right of the peak in isobaric-specific heat (cases 9 and 10), as indicated by the values of the pseudo-vaporization powers [18] replicated in Table 5. Figure 14 depicts the density fields of cases 9 (top) and 10 (bottom). Here, it is clear that density stratification occurs in both cases inside the injector due to the heat transfer between the injector and the flow.

There are still open questions to be answered in the future, such as under which conditions the intermediate plateau formation takes place. From a quantitative point of view, further experimental data would also be of value to develop generalized models and establish a precise definition between liquid-like, gas-like, and transitional states. However, one point is clear: the influence of injector heat transfer in supercritical jet development.

9. Conclusions

Employing an incompressible but variable density approach, the simulation of several supercritical injection configurations is undertaken based on the concepts of thermal breakup and the pseudo-vaporization powers, both available in the literature. The proposed solver is tested in several typical configurations: potential core, subsided core, and plateau formation, for which different degrees of agreement with the experimental data are observed.

In the case of high pseudo-vaporization powers, leading to the formation of dense potential cores, a higher level of agreement with the experiments is reached, highlighting the proposed solver’s ability to replicate the jet’s mean axial density distribution. In contrast, in the case of the lowest pseudo-vaporization power, a qualitative agreement is found between experiments and numerical predictions. Lastly, for intermediate pseudo-vaporization powers characterized by the formation of an intermediate plateau, it was impossible to retrieve its formation, obtaining a similar axial density evolution as the one observed in the subsided core cases. Further work will be conducted toward replicating the plateau formation since it is currently absent from the literature.

Comparing the obtained numerical results with compressible formulations from the literature indicates similar performance between both formulations in terms of the mean axial density decay rate predictions. Moreover, the choice of the equation of state and the thermophysical properties modeling is more preponderant than considering the flow as compressible or incompressible but with a variable density. As such, an incompressible but variable density approach could help alleviate the computational load of numerical computations at such conditions.