Analyzing single and multicomponent supercritical jets using volume-based and mass-based numerical approaches

Trans- and supercritical fluids have attracted the interest of many researchers, due to their broad range of applications from food processing and energy storage to modern combustion engines. For example, supercritical CO₂ is widely used in many applications such as dry cleaning, dyeing, energy storage, and power cycles.¹ In modern combustion systems, increasing the combustion chamber pressure to enhance engine performance has resulted in the fuel being injected into the chamber at pressure and temperature conditions that are significantly above its critical point. Common jet fuels such as Jet A, JP, and their variants have critical temperatures and pressures in the range of 650–700 K and 2–3 MPa, respectively,² whereas in modern gas turbine engines, the compressors are capable of pressure ratios over 40:1,³ which results in combustor operating at pressures greater than 4 MPa (40 bar) and temperatures up to 2000 K. Similarly, automotive diesel engines operate above 3 MPa (30 bar) and can attain up to 20 MPa (200 bar) cylinder pressure, which is usually trans-/supercritical for the diesel fuel.⁴ Furthermore, cryogenic fluids are being utilized in various high-end modern applications and usually find themselves in a high-pressure environment exceeding their critical pressures. Cryogenic rocket engine is a well-known example whose combustion chamber pressure can reach up to 20 MPa, which is supercritical for the O₂ injected as a cryogenic liquid. The “Cryopower” variant of high-efficiency Recuperated Split Cycle Engine (RSCE) is another example of use of cryogenic fluids at supercritical conditions, where cryogenic nitrogen is injected into the compression chamber at supercritical pressures to achieve near-isothermal compression.^5–7 

Despite continuous efforts in the scientific community to understand the multi-physics of the trans- and supercritical fluid flows, there are still many aspects that have not been explored, especially toward the understanding of multicomponent mixing characteristics. When a liquid jet is injected in a gaseous environment close to or beyond the critical point several thermodynamic and physical processes take place. First, unlike sub-critical injection conditions, phenomena like ligament formation or droplet evaporation are less likely to occur. The vanishing latent heat and surface tension effects allow for turbulent mixing and diffusion to become the dominant mechanisms of jet disintegration.⁸ As a result, diffusion layers develop around the jets. A second characteristic commonly attributed to the supercritical regimes is that there is no clear separation of vapor and liquid phases. It has been suggested from previous works that the transition of the liquid phase to a supercritical state is very fast, and this is the reason why the two-phase behavior cannot be sustained.^9,10 However, this conclusion has been recently revised. The concept of “Widom line,” the extension of the coexistence line of gas and liquid above the critical point in the (p,T) plane, has been introduced, and it is considered as an indication of the separation of the liquid-like behavior and the gas-like behavior in the supercritical states.¹¹ Various thermodynamic response functions can be used to locate the Widom line, but approaches do not always coincide, as explained by Banuti et al.¹² The most common definition is the locus of heat capacity maxima, which will be used throughout this study.

Initial experimental and numerical studies on fuel injection at high-pressure supercritical conditions were mainly focused on space propulsion systems. Mayer and Tamura¹³ and Oefelein and Yang¹⁴ are among the pioneers. The past few decades of research toward the same direction has been extended to conditions relevant to transportation and power generation applications (gas turbines and diesel engines). The first numerical studies on supercritical jet mixing focused on single-species configuration corresponding to Mayer et al.'s experiments where cryogenic LN₂ was injected into warm N₂ maintained at high supercritical pressures and temperatures. Following the studies conducted by Chehroudi et al.¹⁵ and Mayer et al.,¹⁶ various complementary numerical investigations on the flow behavior and mixing process characteristics in single-species configuration have followed, as for example.^17–23 Though not exactly related to practical conditions, these studies provided significant insight into the fundamental thermophysics of jets where the injected liquid transitions into a supercritical fluid. Notably, it was shown that the density gradients arising due to the continuous decrease in fluid density as the jet heats up were identified to have a stabilizing effect by suppressing radial fluctuations and hence delaying the jet breakup and mixing.^17,18 This was also confirmed by various other researchers through their numerical simulations.^19,22–27 Recently, Gopal et al.²⁸ have explored further the three dimensional aspect of such density gradients and their influence on both axial and radial fluctuations.

Bellan and co-workers have also provided pioneering contributions toward numerical studies of the mixing process in cold flows at high pressure. With many of their studies focusing on the temporal mixing layers,^29–31 the authors give insight into the influence of the thermodynamic characteristics of a binary system on the evolution of the mixing layer. However, these works were limited to conditions above the pseudo-boiling point, and hence, the rapid transition of thermophysical properties of the fluid around the pseudo-boiling was not explored. In their more recent work on temporal mixing layers, the dissimilar diffusion characteristics existing between single-species (N₂) and binary-species (CH₄–N₂) mixing at supercritical pressures was revealed with the aid of DNS simulations.³² 

One of the main challenges with respect to multicomponent jets and their mixing characteristics is the estimation of the critical point of the mixture and the calculation of the mixture properties of the fluid, which depend on the mixture fraction and are different from those of its constituent fluids. As a result even if the environment is supercritical for the injected fluid, subcritical phenomena can be observed especially when multi-species cases are under investigation. For example, as far back as two decades ago, Mayer et al.¹⁶ observed that when LN₂ is injected into Helium (and not gaseous N₂ as his baseline experiments) at supercritical pressures around twice the critical pressure of N₂, the jet presented both liquid-like and gas-like behavior. Moreover, droplet/ligament formation and surface tension were widely observed in fuel injection experiments by various researchers at high supercritical temperatures and low supercritical pressures (with respect to the respective fuel's critical pressure).^33,34 In other words, a sub-critical regime may be observed prior to the transition to a full diffusion mixing. These fuel injection experiments always involve multiple species.

The aforementioned complexities of trans- and supercritical injection makes the numerical simulations at such conditions very challenging. Not only the choice of the thermodynamics and transport properties models are essential but also the numerical approach plays an important role. In addition, depending on the initial thermodynamics conditions of the injected fluid and environment, condensation may occur. Therefore, some studies have also considered including phase separation in their approach.^35–37 

More recent numerical works of Poblador-Ibanez et al.³⁷ and Davis et al.³⁸ have focused on the laminar mixing layers between a liquid stream and a gas stream under supercritical environments in the near field of a splitter plate. In the work of Davis et al.,³⁸ validation of boundary layer approximation as well as similarity was investigated, and a similarity transformation combined with real-fluid thermodynamic model and a local phase equilibrium approach was also developed by Poblador-Ibanez et al.³⁷ for such laminar two-phase mixing layers at supercritical pressures. Although they provided useful insight, it is not clear if this simplified configuration is representative of realistic three-dimensional jets, also because of the relatively low velocities and the short time covered by the simulations.

The influence of the thermodynamics and transport properties to the physics of these flows characteristics has been stressed by Bellan.³⁹ In a recent work, Toki and Bellan³² with the aid of DNS simulations revealed that the specific mass diffusion is considerably different between non isothermal single-species (N₂) and binary-species (CH₄–N₂) cases at supercritical pressures. In the binary-species case, uphill diffusion was shown to exist due to the Soret effect. On the other hand, Ma et al. emphasized that the numerical dissipation from the chosen conservative scheme (fully conservative or quasiconservative) influences the predicted mixing behavior for under-resolved cases.⁴⁰ An extensive review of the recent contributions in the context of the numerical modeling of trans- and supercritical injection is given by Ries and Sadiki⁴¹ who showed that a consensus on a standard approach has not yet been reached.

One of the popular choices for fuel injection simulations, aside from the VoF methods, is the Eulerian homogeneous mixture model due to its ease to include chemical reactions. Hence, most of the supercritical fuel and even supercritical cryogenic numerical studies are performed using the Eulerian homogeneous mixture model combined with real-fluid thermodynamic and transport properties, for instance, in Refs. 19, 35, 37, and 42. This approach is based on the transport of species mass fractions and can be cast in terms of single-fluid mixing or multi-fluid mixing, in the sense of Ref. 11. In contrast, the VoF model is traditionally utilized for non-reacting immiscible fluids, as it is able to include the surface tension effects and numerically reconstruct the interface between these fluids. In Ref. 6, it was shown that a pressure-based compressible VoF model (with diffused interface) by Miller et al.⁴³ is able to simulate supercritical cryogenic fluids when combined with a suitable real fluid thermophysical model. In the subsequent work,²⁸ the model was enhanced with temperature adaptive surface tension and mass diffusion (both molecular and Soret) although the simulations corresponded to entirely supercritical conditions alone. This model applies surface tension effects where the fluid is locally subcritical and mass diffusion effects, where it is supercritical. Recently, Poblador-Ibanez and Sirignano^44,45 have also proposed an interesting and suitable VoF framework for variable density two-phase supercritical fluid flows using a low-Mach-number assumption, which however treats the thermodynamic pressure as constant, while only the dynamic pressure varies.

While various dissimilar modeling approaches have been utilized by previous research works in the numerical simulations of cryogenic and supercritical jets,^{17–25,27,35,37,42} the choice of the modeling approach and associated numerical methods for trans-critical conditions is seldom justified in the literature. The majority of these approaches predict with high accuracy single-species cases or binary-species cases of selected fluids such as liquid and supercritical/gas phases of cryogenic fluids, atmospheric gases, and light hydrocarbons. However, there has been quite varied results and discrepancies observed in the experiments of heavy hydrocarbon injection into supercritical environment where multiple fluids/species and their respective mixtures are involved.^{10,34,46–48} For any mixture, its critical point lies between those of its constituent fluids/species.⁴⁹ However, the exact value depends on the mixture fraction, which varies throughout the jet flow domain. As a result, the local critical point also varies throughout the domain which in turn dictates the local conditions for two-phase behavior inside the jet domain. This study comparatively explores the influence of two different solver algorithms (Eulerian homogeneous and Eulerian VoF) along with respective mixture models (volume fraction weighted average and non-linear Van der Waals mixing rules) on the thermophysical properties and thermodynamic (state) evolution of the injected fluid in supercritical binary-species environment. This is deemed necessary to determine the numerical model factors that affect the injected fluid's transition across the critical point during its evolution into multi-species mixture of injected and ambient fluid, as a precursor for future numerical studies to implement appropriate two-phase behavior that fundamentally determine the spray/jet characteristics such as surface tension and diffusion effects.

In our work, a comprehensive comparison of the two methods mentioned above (Eulerian homogeneous mixture model realFluidReactingFoam and Eulerian VoF model coolFoam) is presented. The implementation of the models is the ones proposed by Gopal et al.^6,28 and by Ningegowda et al.,^50,51 both using OpenFOAM⁵² as basis. The aim of this work is twofold:

1. To complement the previous studies by providing further insight specifically targeting the understanding of the similarities and differences between the single and two components supercritical mixing with focus on fluid mechanics, heat transfer, and mass transfer processes;

2. To highlight the effect that volume-based vs. mass-based modeling approaches and their mixture models have on the prediction of the mixing dynamics corresponding to the different thermodynamic behavior and with respect to the two-phase regions.

To this end, after models descriptions, two different cases are investigated in a two-dimensional computational domain. The first configuration is cryogenic nitrogen injected into warm nitrogen environment (similar to earlier Mayer's experiments¹⁶) and the second is n-dodecane fuel injected into gaseous nitrogen at conditions typical of diesel engine combustion chambers (Engine Combustion Network, ECN, spray A conditions⁵³). The two cases are designed to share the same reduced pressure and temperature and the same Reynolds numbers in order to emphasize the role of multi-species mixing.

A comprehensive description of the two computational fluid dynamics (CFD) frameworks is given in Sec. II, followed by a detailed presentation of the computational domain and the simulations setup in Sec. III. Finally, in Sec. IV, the numerical results obtained through the two approaches are compared for an analysis of the physical aspects and mechanisms that drive the mixing process under trans- and supercritical flow conditions.

A comparison of two multi-species solvers originating from different numerical frameworks based on OpenFOAM is used in order to provide insight into the complex dynamics of single- and multicomponent supercritical mixing. One approach, coolFoam (cF) uses a compressible VoF framework, which has been improved with adaptive surface tension and interface diffusion capabilities in order to be able to capture the fluid transitioning from subcritical liquid state to gas-like phases at supercritical state. A rather simplified approach for modeling thermophysical properties is also included based on polynomial fitting to NIST data.^6, coolFoam was originally designed to simulate cryogenic fluids; however, here we test its capabilities for non-cryogenic fluids as well.⁶ The other approach, realFluidReactingFoam (rFRF) is a compressible multi-species solver accounting for the rigorous calculation of real-fluid thermophysical properties through a cubic equation of state for the one-fluid mixture with appropriate thermodynamic relations.^50,51 The code was originally designed for supercritical combustion simulations, although here only non-reactive cases are considered.

One major difference between the two frameworks is that the evolution of species/fluid distribution is simulated through transport of volume fraction in coolFoam and mass fraction in realFluidReactingFoam. An overview of both solvers with the details of their respective frameworks and capabilities are presented in Table I.

The in-house solver coolFoam is designed for compressible non-isothermal two-fluid simulations where diffusive transport of species is not negligible. The equations are based on volume fraction quantities. The injected fluid is treated through a hybrid approach based on a combination of previously developed VoF methodologies⁵⁴ and through a conservative species approach that allows for different fluids to diffuse into each other for the regions with miscible phases (gas and/or supercritical). The volume fraction equation includes terms for diffusive transport of species, rearranged in volume fraction quantities.

The governing conservation equations for volume fraction, momentum, and energy are solved for a single effective fluid weighted by the local volume fraction. This weighted averaging eliminates the generation of negative pressures and drastic pressure fluctuations, sometimes observed when cubic EoS are utilized to estimate thermophysical properties of mixture inside the two-phase region.⁴⁰ 

When dealing with bounded scalars such as α, it is important to respect the boundedness of the variable. For the α transport equation, Eq. (7), the convective term can be source of unboundedness. Among the different strategies to mitigate the issue of the boundness, in this study the Multidimensional Universal Limiter with Explicit Solution (MULES), proposed by Weller,⁵⁵ is adopted. The MULES limiter consists of an explicit solver, which is based on the flux-corrected transport in order to maintain the boundedness of the scalar. A detailed description of the algorithm is beyond the scope of this study, and more details about its implementation can be found in the literature.⁵⁶ Despite being initially derived in an incompressible framework, the MULES scheme can be directly applied to the compressible case. To be suitable to this scheme, the volume fraction equation is rearranged and the contribution from the compressibility is accounted numerically as a source term. This allows to hold a similar numerical methodology for the transport of the volume fraction between the compressible and incompressible solvers.

In coolFoam, the thermophysical properties of N₂, i.e., density, heat capacity, viscosity, and thermal conductivity, are specified as polynomial functions obtained by polynomial fitting of the NIST data for the temperature range encountered in the simulations. The concentration-driven diffusion coefficient D of nitrogen at atmospheric pressure is estimated through models specified by Fuller,^57,58 which are then translated for high pressures using the generalized chart and method described by Takahashi.⁵⁹ For more details on the polynomial fitting as real fluid thermophysical model and diffusion coefficients employed in coolFoam, the reader is referred to Gopal et al.^28,60

For the C₁₂H₂₆ simulations, the transport properties are modeled differently. The NIST data are available up to 700 K. The viscosity and thermal conductivity are modeled through the polynomial fitting as for the N₂. With respect to C_p, using a polynomial fitting for this range will lead to negative C_p for $T > 800$ K. For this reason, the heat capacity is modeled as a function of temperature T from a set of coefficients taken from JANAF tables of thermodynamics, while the diffusion coefficients are calculated by the same method utilized for nitrogen, by Fuller et al.⁵⁷ and Takahashi.,⁵⁹ For both fluids, liquid phase surface tension is assumed negligible and not modeled close to the critical temperature. This is appropriate as the fluid's surface tension decreases linearly with increase in temperature and completely vanishes at the critical temperature beyond which the fluid is supercritical and does not exhibit surface tension.

The parameter $F c , m$ is used to account for the shape and the polarities of the fluid, which allows to obtain a good prediction for a broad range of fluid states. The parameter E₆ belongs to a group of parameters E_i (⁠ $i = 1 , 2 … 10$⁠) that are function of the acentric factor ω and constants given in Ref. 62 and are also used to evaluate the other terms. For the sake of brevity, they are not developed here as well as the combination rules used to evaluate the mixture properties indicated by the subscript “m” but are detailed in Ref. 62.

Furthermore, details regarding the solver and the validations of the thermophysical quantities against reference data can be found in Refs. 50 and 63.

The geometry used for the investigation is a 2D simplification of the Mayer's configuration.¹⁶ The computational domain consists of $150 d × 71 d$⁠, where d is the hole diameter, with the liquid injector exit located at (0, 0). Two different liquid injection scenarios into supercritical environments are investigated:

• Single species: the injection of liquid nitrogen into gaseous nitrogen (LN₂–N₂), and

• Two species: the injection of n-dodecane into nitrogen (C₁₂H₂₆–N₂).

For both the system configurations, different operating points investigated are presented on the phase diagrams of C₁₂H₂₆ and N₂, in absolute (in Fig. 1) and reduced (in Fig. 2) values of pressure and temperature. It is important to underline that the operating conditions for C₁₂H₂₆–N₂ are investigated at the same reduced pressure p_r, reduced temperature T_r, and injected Reynolds number Re_inj with the corresponding LN₂–N₂ cases in order to be able to reveal any potential flow similarity arising in the case of single- and multi-species supercritical mixing.

The case summary is provided in Table II. The boundary conditions and the schematic of the computational domain are shown in Fig. 3. It consists of a 2D uniform and structured mesh of which each cell has dimensions of $0.05 d × 0.05 d$⁠. The gray dashed area of $75 d × 50 d$ in Fig. 3 shows the region where time averages are performed for analyzing the results, while the red one of $30 d × 10 d$ represents the area of focus for the analysis.

It is worth pointing out that turbulence was not included in the simulations as a choice in 2D. Comparisons between DNS vs LES have already been performed (see  Appendix A), and results showed that the turbulence model did not affect significantly the conclusions on the thermodynamic aspects. Various authors have also avoided including a turbulence model, in favor of a consistent and direct analysis of the thermodynamic state.⁴⁰ 

Figures 4 and 5 present the instantaneous and the time-averaged contour plots of volume fraction for the two reduced pressures for the liquid nitrogen jet case injected into nitrogen (LN₂–N₂). The dark region in the time-averaged contour plots represents the core of the jet occupied by the injected LN₂ (α₁ = 0.5–1). Both the solvers result in a similar intact core region for these LN₂–N₂ cases and predict a reduction of the core with increase in supercritical pressure. This is expected due to the smoother transition of thermophysical properties at higher supercritical pressures, which leads to weaker density gradients that in turn are unable to dampen the turbulent fluctuations as effectively and result in earlier breakup when compared with lower supercritical pressures.^28,50 The instantaneous contour plots of volume fraction reveal the diffusive nature of LN₂ when injected into the same (species) fluid at supercritical pressures and temperatures.

More evident differences between the two frameworks are manifested in the time-averaged contour plots of temperature and density in Figs. 6 and 7, respectively. Particularly, the transitioning (from liquid-to-supercritical) injected LN₂ occupies a deeper and wider part of the mixing region (beyond the intact dark blue core) as it heats up in rFRF, as opposed to cF. The transitioning fluid can be observed by the yellowish-green (approximately 150–250 K) region in the figures. The temperature distributions reveal that the injected fluid transitions and heats up more rapidly in cF compared to rFRF. Here, the cooler bluish regions are more prevalent in the jet domain of rFRF. Ultimately, this is the reason for the transitioning region occupied by comparatively cooler LN₂ in time-averaged contours of rFRF. These differences observed in the temperature distribution are also reflected in the density distribution. The density field better highlights the regions where the injected liquid nitrogen (high-density regions) transitions to gas-like fluid (low-density regions). We observe that the transitioning fluid occupies a larger region in rFRF as opposed to cF.

As for the effect of the supercriticality level, it is confirmed by both codes that an increase in the chamber gas density, probed through the variation of the pressure at constant temperature, induces a broader jet, as it is visible in the density plots. This is in agreement with numerous literature works,^15,50,64 which explain that the entrainment and the spreading angle increase with a decrease in the density ratio (between injected and chamber densities). Furthermore, details and quantitative analyses supporting this argument will be provided later in the manuscript.

The strong density gradients, which are a unique feature of transitioning jets due to the fluid's continuous change in density with respect to temperature, are presented in Fig. 8 for the near nozzle region. The tubular structure of the strong density gradients near the nozzle exit is identical across the solvers. The tubular structure of the density gradients and their collapse in such three-dimensional supercritical jets have been discussed in-depth in Ref. 28.

With increase in supercritical pressure, the extent of the strong density gradients decreases. One difference observed between the solvers is that the collapsed and dispersed tubular density gradients²⁸ in the mixing region, after 10d, are much stronger in rFRF than cF.

Despite the small differences in the fluid transition/heat-up predicted by the two solvers, the flow dynamics close to injection seems unaffected, as the core region up to about 10d in time averaged contour plots of volume fraction, temperature, and density are very similar. Furthermore, both the solvers result in identical intact liquid core length, establishing the independence of flow dynamics under these conditions.

As the pressure increases, the density gradients reduce due to the smoother transition of the thermophysical properties from liquid to supercritical gas-like at higher supercritical pressures. As these density gradients are responsible for suppression of instabilities^17,19,25,26 the reduction of these gradients with higher pressures allows for more instabilities to be generated around the liquid core of the injected nitrogen, which in turn results in earlier breakup explaining the shorter liquid core at lower pressures.

Figures 9 and 10 show the instantaneous and time-averaged contour plots of the volume fraction for the two reduced pressure cases with the n-dodecane injected into hot nitrogen. The dark colors here represent the core region of the jet dominated by injected C₁₂H₂₆ (α₁ = 0.5–1). With increase in supercritical pressure, the core region reduces for the two-species (C₁₂H₂₆–N₂) cases, as discussed already for the single-species cases. However, in contrast to the single-species (LN₂–N₂) cases, the instantaneous contour plots of volume fraction (see Fig. 9) present the injected C₁₂H₂₆ behaving like a liquid into gas. In fact, it is interesting to note that the single-species cases present miscible fluid behavior (species diffusing into each other), whereas the two-species cases present characteristics of fluid column breaking up and shedding dense packets which then diffuse, solely based on the thermophysical properties.

Larger differences between the results of the two solvers are present in the contour plots of temperature for the two-species cases (see Fig. 11), compared to the observations for single-species cases. The injected C₁₂H₂₆ heats up quicker in cF compared to rFRF. This can be observed in the above figures where the mixing region of the jet in rFRF is dominated by lower temperatures (620–700 K), whereas the mixing region in cF is dominated by higher temperatures (700–800 K).The normalized density contours shown in Fig. 12 do not seem to be influenced by the temperature difference across the results of cF and rFRF. It is important to note that in this binary mixture configuration the thermophysical properties of the injected fluid (C₁₂H₂₆) at the chamber state are rather different than the ambient chamber fluid, unlike the single-species cases.

Finally, similar to the single-species cases, identical density gradients across the solvers are also observed in the two-species cases C₁₂H₂₆–N₂ and their extent decreases with increase in supercritical pressures. As with the single-species cases, the rFRF predicts stronger density gradients along the centerline beyond about 10d, seen as reddish hue in Fig. 13.

In summary, the mechanical features are predicted in a very similar way across the two solvers, but the mixture thermodynamics are rather different, suggesting the role that the different mixing models employed by the two solvers play in the prediction of the mixing behavior. In other words, for two-species mixing, the injected fluid has drastically different thermophysical properties from that of the chamber fluid; therefore, the thermodynamic state (controlled by the enthalpy in isobaric cases) has a marginal effect on the density, as major differences exist mainly due to composition variations, but has a profound impact on the temperature predictions.

In Secs. IV B 1, IV B 2, and IV B 3, normalized mean profiles of density, temperature, and velocity fields in the axial and the radial directions at various locations are discussed for single species (LN₂ injected into N₂) as well as two species (C₁₂H₂₆ injected into N₂). A normalized profile correlates the fluid's properties predicted across the two thermophysical models at injection for different supercritical pressures.

Figure 14 shows the $⟨ ρ + ⟩ = ⟨ ρ ⟩ / ρ inj$ profiles for the LN₂-N₂ cases. Figure 15 shows the corresponding $⟨ T + ⟩ = ⟨ T ⟩ / T inj$ profiles. Overall, axial and radial profiles predicted by the two solvers demonstrate good agreement with each other. The density profiles of the two solvers agree well for the entire range (i.e., up to $x / d = 30$⁠) with minor differences are observed beyond $x / d = 20$⁠. For both solvers, the centerline density is constant near the injection region representing the liquid potential core. Then, a steep fall in centerline density is observed between $10 ≤ x / d ≤ 20$⁠, and after that, the density tends to the warm chamber nitrogen density. For a fully developed jet, disintegration and the mixing process can be characterized through three regions: (a) potential core, (b) transitional, and (c) self-similar. In the single-species cases, this categorization also corresponds with the fluid phases: (a) liquid, (b) transitional, and (c) ambient/supercritical gas-like.

Increase of the supercritical pressure from $p r = 2.95$ to $p r = 5.05$ shifts the breakup length from about $x / d = 8$ to about $x / d = 6$ and can be associated with the reduced extent of density gradient formations, which are known to suppress instabilities and delay breakup.²⁸ Similarly, with increase in supercritical pressure, the centerline temperature rapidly increases toward ambient (chamber) temperature suggesting a more enhanced mixing process. Although the density profiles across the solvers are in good agreement, the temperature profiles $⟨ T + ⟩ = ⟨ T ⟩ / T inj$ in Fig. 15 indicate that the injected N₂ heats up faster in cF when compared to rFRF, particularly beyond $x / d = 25$⁠, that is after the transitional region. This observation is consistent with the observations of the contour plots presented in Sec. IV A.

Figure 16 shows the density $⟨ ρ + ⟩$ profiles. As with the single-species cases, the centerline density of the two-species cases reveals the distinct transition regions established in the fully developed jet. The centerline density is constant in the near nozzle region; then, it drops rapidly as the injected C₁₂H₂₆ transitions (heats up) while also simultaneously mixes with the warm ambient N₂. With increase in supercritical pressure, the breakup occurs earlier across both the solvers, from $x / d > 10$ to about $x / d = 5$ for the two reduced pressures. Although the centerline and radial density profiles from both solvers are in good agreement for the higher reduced pressure, the density drop in cF starts slightly further downstream for the lower reduced pressure. Unlike the density profiles of the two solvers that are in good agreement for the entire range of consideration, different behaviors are observed for the temperature profiles. The mean temperature profiles are presented in Fig. 17, and it can be observed that immediately after the constant temperature region (corresponding to intact core distance), they deviate and the discrepancies between the two solvers here are higher than in the previous single-species cases. The centerline temperature rises faster in cF when compared to rFRF. This is consistent with what has been observed earlier in the temperature contours in Fig. 12, where in cF, the incoming C₁₂H₂₆ jet quickly rises its temperature once in contact with the surrounding N₂, while for rFRF, the jet presents a slower temperature rise once in contact with the surrounding N₂.

To investigate this difference, the thermophysical properties at both the pure fluid state and mixed state predicted by each solver need also to be considered. We stress again that, while in the single-species case, the state of the mixed fluid is controlled by its local temperature, this is not the case for the two-species situation. Here, since the injected and chamber fluids are of different species, the critical properties (and Widom line) of the mixed fluid depends also on the local mass (or volume) fractions. Hence, as such, the state of the mixed fluid in the transitional region cannot be calculated without utilizing a model for predicting the mixture properties starting from the pure species data.

Figure 18 shows the comparison of $⟨ m ̇ + ⟩$ for both the LN₂–N₂ and C₁₂H₂₆–N₂ cases. Independently of the solver, overall the entertainment is higher for the single-species cases compared to the binary mixture results. Even though same reduced pressure and same reduced temperature are used, the density ratios between the injected fluid and the ambient environment are different for the single- and multi-species cases. At reduced pressure $p r = 2.95$ for instance, the density ratio is on the order of 20 for C₁₂H₂₆–N_2, which is approximately four times of that $LN 2 − N 2$ case. For the single-species case as the pressure increases, the mixing rate also rises. Moreover, the differences between the two solvers are more evident for this case than the two-species case. In other words, for the single-species cases, not only the mixing process is intensified because of the similarity at molecular level, but also because of the denser ambient gas. On the solvers side, rFRF records higher values of the entrainment growth rate for the cryogenic nitrogen cases, while curves agree for the binary mixture cases.

One of the main differences in the codes lies in the modeling of thermophysical properties of the fluids and in particular in how the mixture thermophysical properties are derived from the constituent fluids. Previous studies have highlighted the importance of the accurate modeling of these properties although until now there are no direct comparisons available between different approaches in order to identify their effect on the physical predictions. In the next paragraphs, we will try to close this gap by discussing in more details the predicted thermophysical properties and transport properties profiles of the fluids by the two solvers and investigate how they are linked with the differences observed previously (density and temperature profiles as well as mixing rates).

Figure 19–21 show the scatterplot distribution of the dynamic viscosity, the constant pressure heat capacity, and the thermal conductivity vs the mixture temperature, obtained from each CFD control volume (cell). In addition, Fig. 22 reports the CFD scatter data of the temperature with respect to the volume fraction. All the plots are given for both $LN 2 − N 2$ and C₁₂H₂₆–N₂ cases. The viscosity is a crucial property affecting the mechanics of the flow, whereas heat capacity and thermal conductivity influence the thermodynamic behavior of the flow. Overall, for the $LN 2 − N 2$ cases, the two solvers show similar trend and the results are in excellent agreement with the NIST data for all the thermophysical properties, while for the C₁₂H₂₆–N₂ differences are observed, as expected.

For the single-species and two-species cases, as aforementioned, the cF utilizes polynomial fitting to NIST data for N₂ and JANAF for $C 12 H 26 − N 2$⁠, whereas the rFRF utilizes the NASA polynomials with real-fluid corrections derived from PR EOS. The figures present not only the injected fluid's (LN₂ or C₁₂H₂₆) thermophysical properties but also the properties of the injected fluid as it transitions and mixes with the chamber fluid (N₂). The transition can be identified through the temperature, which translates to the fluid's state as the cases are isobaric and the mixing can be identified through the α₁ volume fraction. The properties of the individual fluids (N₂ or $C 12 H 26 − N 2$⁠) from the NIST data for the pure species are also plotted within each profile as limits, to give a perspective on the thermophysical properties of the mixture predicted by each solver.

For the single-species LN₂–N₂ cases [Figs. 19(a), 20(a), and 21(a)], the thermophysical property profiles of the injected N₂ as it transitions and mixes with the surrounding gas are in very good agreement with the NIST data for both solvers (cF and rFRF). While this is not surprising for cF since it uses NIST data based polynomials, the PR-EOS combined with the Chung methods for high-pressure conditions also matches very closely with NIST data with minor deviations on the pure liquid side, which are expected because of the selected EOS intrinsic accuracy. Since both the injected and chamber fluid are the same species (N₂), the mixture properties also lie on the curve. Looking closely, it can be observed that the volume fraction against temperature is consistent across both the solvers, indicating that the predicted flow evolution and mixing processes are also similar between the solvers, as observed in Secs. IV A and IV B.

In contrast, for the binary mixture $C 12 H 26 − N 2$ cases [Figs. 19(b), 20(b), and 21(b)], differences are observed between the solvers. In the case of cF since the mixture properties are calculated through volume fraction weighted averages, it can be observed that the thermophysical properties transition gradually from the injected fluid state to the chamber fluid state as the injected C₁₂H₂₆ heats up and mixes with the chamber N₂. On the other hand, rFRF assumes equilibrium without phase change and phase splitting, the calculated temperature corresponds to the frozen adiabatic mixing temperature, non-linearity emerges, and as a consequence, a slightly lower temperature (under the injection temperature) is recorded for $0.5 < α 1 < 0.8$⁠. Moreover, discrepancies are observed for the pure n-dodecane computed properties and NIST data for rFRF due to the known deviation behavior of the chosen EOS on the region below the critical point. As for the mixture properties, it is important to recall that in rFRF, both the PR-EOS and the Chung correlations for the transport properties make use of non-linear mixing rules for real fluid mixtures, which could explain the significant differences between the solvers observed particularly on the mixing region.

These profiles also help us explain the differences in temperature contours predicted by the two codes, especially for the binary mixture $C 12 H 26 − N 2$ cases. Since the thermal diffusivity quantifies the rate of heat transfer through the relation $α = κ / ρ C p$⁠, differences in the mixture heat capacity and the mixture thermal conductivity translate into difference in temperature. Although the thermal conductivity is generally higher in rFRF that enhances the rate of heat transfer, the heat capacity also has higher values particularly for the higher volume fractions (⁠ $α 1 > 0.5$⁠), which in contrast reduces the rate of heat transfer. As a result, although these aspects identify key thermophysical property differences of the mixture predicted by the solvers, which translate into the observed difference in temperature contours, a straightforward comparison is not possible. Enthalpy-to-temperature relationships are also very different in the two models, being one linear and one highly non-linear. To decouple the effects, it would be necessary to model the heat flux in an identical way, so that the remaining differences can be attributed to the mixing rules. However, this is complicated by constantly evolving fluid properties and the transient nature of the flow; therefore, it is beyond the scope of the current study.

Figure 22 depicts the scatter data obtained through the numerical simulations in the form of temperature-composition phase diagram. For the single-species cases, the trends are very similar between the two solvers [Fig. 22(a)]. For the $C 12 H 26 − N 2$ cases [Fig. 22(b)], the solvers display different mixing behavior. In an effort to explain this behavior, the bubble point and dew point lines are also plotted for the $C 12 H 26 − N 2$ cases along with the isobaric–adiabatic temperature and the isobaric–isochoric temperature calculated following Refs. 36, 66, and 67. It can be observed that the rFRF results follow closely the adiabatic temperature (which is here the frozen temperature) as expected, since it assumes equilibrium with no phase separation. On the other hand, cF shows a different behavior and predicts a temperature that falls within the two limits. For the latter, mixing occurs at higher temperature and hence the injected liquid temperature rises faster. In addition, since the injected temperature is below the critical temperature of n-dodecane, both approaches show that the mixture goes through the two-phase region when crossing the vapor liquid equilibrium (VLE) curve, but for rFRF, it is for a substantially larger amount of cells compared to cF. This is consistent with the previous observations of mean profiles, where in the centerline, the density profiles do not show much differences but as the jet spreads out, rFRF cases demonstrate lower and wider temperature distribution and hence higher density distribution in the radial direction. Studies in Refs. 66 and 67 show that the numerical solution of fully conservative formulation of the energy equation follows the adiabatic mixing, whereas quasi conservative approach follows the isochoric mixing. Nevertheless, in this work, the different behavior is potentially associated with the two different mixing models adopted by the two present approaches to compute the thermophysical properties since both of the approaches adopt fully conservative schemes and simulations are run on a common mesh. Mixing rules play also an important role in the mixing process prediction for these binary mixture cases. Temperature values predicted by the volume based method with linear mixing rules are higher than those obtained by the mass-based method with non-linear mixing rules. This causes a different extent of the two-phase region as seen in Fig. 22(b).

The aim of this work is to gain further insight into the mechanism that drives jet mixing at trans- and supercritical conditions. To this end, a thorough comparative study of the results obtained by two different and independently developed numerical solutions was conducted. First, injection of cryogenic nitrogen into a nitrogen warm surrounding at high pressure was considered. Then, using identical flow parameters (reduced pressure, reduced temperature, and Reynolds numbers) liquid n-dodecane injection into nitrogen was also analyzed. The two solvers used the same 2D computational domain and grid resolutions to perform the simulations. The two numerical approaches are pressure based but use two different thermodynamic and mixture models. While the VoF based coolFoam uses polynomial fitting, mass-based realFluidReactingFoam accounts for real-gas thermodynamic effects by making use of PR EOS and Chung's methods for transport properties. Contour, axial, and radial profiles at different locations of the mean density, velocity, and temperature are provided for the different conditions and the two different injected fluids. The main conclusions can be summed up in the following points:

• Similar trends in the flow dynamics were recorded for both solvers with discrepancies in the heated up fluid temperature along the transitional region. The increase of the supercritical pressure results in a shorter breakup length.

• The simulation of single-species cases by both solvers showed that the simulated dynamics of the jet are in excellent agreement with each other, and the small differences are largely due to the different thermophysical property models employed by the solvers. The influence of the difference between volume-based vs. mass-based transport equations is minor on the flow characteristics in case of single-species cases.

• For the binary mixture (⁠ $C 12 H 26 − N 2$⁠) cases, significant differences in the profiles were observed mainly for the temperature, where lower temperature were recorded for realFluidReactingFoam, namely, in the mixing region, in consistency with the fact that this approach considers adiabatic mixing. This leads eventually to different ranges of the two-phase region and ultimately is bound to result in a different behavior in case of reacting scenarios. In addition, it should be noted that the state of the mixed fluid in the transitional region cannot be calculated without utilizing a model for predicting the mixture properties starting from the pure species data. Hence, the difference in mixing rules has a greater influence for the multi-species case. Despite the many studies in the literature, none is currently useful to discern whether one approach or the other is closer to the experimental mixture temperature, mainly because of lack of spatial and temporal accuracy as it would be requested for these validations. Therefore, further investigations are still required in this regard.

• Despite the different thermodynamic closures, the density profiles are nearly identical. A similar prediction of the flow field dynamics by the two numerical approaches leads to the same density distribution.

• On the other hand, for the entrainment rate, the single-species cases show more significant differences between the two frameworks. In addition, despite the fact that the same flow parameters have been adopted for both the single- and binary-species cases, the mixing rate is found to be not only influenced by the similarity at molecular level but also by the density ratio between the injected fluid and the chamber.

Overall, this study highlights how high-pressure devices, operating in supercritical and trans-critical conditions (such as rockets, diesel engines, and gas turbine burners) still need further model development to accurately predict thermo- and fluid-dynamic processes of the injected fluid, to enhance their efficiencies. For instance, auto-ignition and combustion phenomena depend significantly on local temperatures.

Future research directions could explore the inclusion of surface tension under sub-critical and trans-critical regimes, with progressively diminishing effects. In addition, to deal with finite rate phase change processes, the current formulations of the real-fluid models based on the local thermodynamic equilibrium assumption should be modified to allow modeling of such phenomena. In addition, also with the current models, further investigations could be carried out with a wider range of fluids and considering fully 3D configurations.

The University of Perugia authors acknowledge the support from King Abdullah University of Science and Technology, Saudi Arabia, under the CRG Grant No. OSR-2017-CRG6-3409.03. The first author also acknowledges the support for student mobility from the University of Perugia through the EU Erasmus traineeship program. The University of Brighton and the University of Oxford authors would like to acknowledge funding by the UK Engineering and Physical Science Research Council support through the Grant No. EP/S001824/1.

The authors have no conflicts to disclose.

Faniry Nadia Zazaravaka Rahantamialisoa: Conceptualization (lead); Data curation (equal); Formal analysis (lead); Investigation (equal); Methodology (lead); Software (lead); Validation (lead); Visualization (equal); Writing – original draft (lead); Writing – review & editing (equal). Jaya Vignesh Madana Gopal: Conceptualization (lead); Data curation (equal); Formal analysis (lead); Investigation (equal); Methodology (lead); Software (lead); Validation (lead); Visualization (equal); Writing – original draft (lead); Writing – review & editing (equal). Giovanni Tretola: Conceptualization (supporting); Investigation (supporting); Methodology (supporting); Software (supporting); Writing – original draft (supporting). Nasrin Sahranavardfard: Writing – original draft (supporting); Writing – review & editing (supporting). Konstantina Vogiatzaki: Conceptualization (equal); Formal analysis (supporting); Funding acquisition (lead); Writing – original draft (supporting); Writing – review & editing (equal). Michele Battistoni: Conceptualization (equal); Formal analysis (supporting); Funding acquisition (lead); Writing – original draft (supporting); Writing – review & editing (equal).

Data sharing is not applicable to this article as no new data were created or analyzed in this study.

To validate the solver, experimental measurements from the well-known Mayer et al.’s work¹⁶ are used. Specifically, the density values measured in case 9 of the cryogenic dataset are utilized, where liquid N₂ is injected into a pressurized N₂ chamber maintained at supercritical pressure. This specific case is summarized in the Table III. It is necessary here to clarify the discrepancies in the injected temperature. Although the target temperature was 120 K, when the thermocouple to measure injection temperature was set such that it extends 1 mm into the flow (perpendicular to the flow axis), the flow was heavily disturbed and it measured 135 K as the thermocouple was not wet properly. When it was set such that it extends 1 mm into the injector at the axis (along the centerline axis), it measured 122.8 K. Although this 122.8 K is a reliable measurement for the core temperature, it is impossible to achieve a perfectly adiabatic injector, and hence, the temperature near the walls is supposed to be higher. The injection velocity, temperature, and other fluid properties in reality possess a parabolic profile whose characteristics are unknown. Hence, it is impossible for any simulation to match precisely the density measurements close to the injector, as already revealed by past simulations by various researchers.^19,24,25 It should be noted here that due to the lack of more accurate experimental data on cryogenic and/or supercritical fluids, we are forced to utilize this case. Regardless, the past simulations by various researchers have identified that for simulating this case with uniform injection temperature, 135 K is the most appropriate value.

A 3D, axi-symmetric and cylindrical mesh is used, in which the LN₂ is injected through a circular 2.2 mm inlet (Fig. 23). The mesh resolution increases both toward the injection in axial direction and toward the centerline in the radial direction, ensuring that sufficient resolution is present along the jet boundary. The 3D mesh is obtained by extruding a Discretized plane domain around the centerline. More details on the expansion ratios associated with the mesh and grading for mesh1 are available in our previous publication.²⁸ 

Previous numerical simulations by various researchers using RANS and LES have revealed that the common RANS models reproduce average quantities with almost the same accuracy as that of the LES models, and the simulation accuracy depends much more on the chosen thermophysical model.^19,24,25 Moreover, the LES models used commonly were not developed for cryogenic fluids especially in conditions where thermophysical properties of the fluid vary drastically. While Selle and Ribert⁶⁸ based on their a priori study expected existing low-pressure LES models to perform well, their research in-fact showed that SGS (Sub Grid Scale) dissipation was too great. In addition, our previous work⁶³ in aiming to simulate high-Reynolds-number supercritical flow with large density gradients using realFluidReactingFoam has also shown that the filtering approach does not have major impact on the mean flows. Similarly, simulations performed using coolFoam with LES-SGS models have also revealed the same thing where the time-averaged results were largely unchanged, but the large SGS dissipation diffused the interface.

Furthermore, since experimental results of centerline density in cryogenic jets by Mayer et al.¹⁶ are compared against time averaged simulation results and due to the already identified discrepancies in the measurement of injection temperature,¹⁶ it is impossible to identify from previous numerical simulation studies, which of the LES models are more suitable for such supercritical and/or cryogenic fluids. Therefore, for all the above reasons any turbulence model is not used in this present work.

To validate the solvers and the simulation results, a grid independence study and experimental validation are presented together in this section. As we have mentioned earlier, the temperature discrepancies at the injection adds some uncertainty to the experimental results. In the Mayer's experiment, the temperature profile of the injected fluid is actually parabolic, due to inevitable heat exchange taking place between the injector pipe and cryogenic fluid.¹⁶ However, the simulation utilizes a uniform injection temperature (in tune with previous numerical simulations by other researchers^19,24,25) again due to the uncertainties in the injection temperature profile. As a result of uniform injection temperature assumption in the simulation, heat transfer to the injected fluid starts to take place after the fluid is injected into the chamber and the centerline density stays constant until the heat transfer reaches the centerline. In contrast, heat transfer into the injected fluid begins even before injection in Mayer's experiments,¹⁶ resulting in the observed discrepancies up to 5d.

For the grid independency study, three meshes with identical domain and boundary conditions but varying level of refinement are employed. The expansion ratios in all the three meshes are kept similar such that, although the number of cells varies, similar grading exists across the three meshes. Brief details of these three meshes are presented in Table IV. The time-averaged axial centerline density distribution simulation results of this case performed using coolFoam and realFluidReactingFoam with the aforementioned three meshes are presented in Figs. 24 and 25 along with the experimental results by Mayeret al.'s¹⁶ case 9. It can be seen that the results of the simulation in mesh1 improve significantly compared to mesh0, whereas there is only minor improvement downstream (from 7.5d) with further discretization of the mesh, even though the mesh2 has twice the number of cells than mesh1 and the cell size is almost halved in the grid domain near the injection. Thus, mesh1 provides sufficient accuracy to explore the thermophysics, in relation to the computational cost.