Theoretical and Numerical Modeling of Multicomponent Transcritical Diffuse Interfaces Under LRE Conditions

In this study, a theoretical and numerical framework for simulating transcritical flows under a variety of conditions of interest for aerospace propulsion applications is presented. A real-fluid multicomponent and multiphase thermodynamic model, based on a cubic equation of state (EoS) and vapor–liquid equilibrium (VLE) assumptions, is presented to describe transcritical mixtures properties. The versatility of this thermodynamic model is reported since it can represent at the same time the supercritical states as well as subcritical stable two-phase states at equilibrium, via a homogeneous mixture approach. The effect this model has on the evaluation of the thermophysical variables will be emphasized. From the Computational Fluid Dynamics (CFD) point of view, the well-known numerical challenges that arise with the coupling between real-fluid thermodynamics and governing equations under transcritical conditions, are addressed by comparing a fully conservative (FC) to a quasi-conservative (QC) numerical schemes, in the context of the advection problem of a transcritical contact discontinuity.


Introduction
Liquid rocket engines (LREs) combustion chambers often operate by injecting cryogenic propellants at elevated pressure, typically higher than the critical values of pure components and most combustion products. The thermodynamic conditions, which involve high pressure and cryogenic temperature, are usually referred to as transcritical injections; where one or more propellants undergo a transition from a liquid-like state to a gas-like state. These particular injection trajectories induce nonlinearities due to the real-gas effect, as a consequence, the fluid's thermophysical properties are strongly dependent on both pressure as well as temperature. From a physical point of view, the propellants in a supercritical state exhibit no clear distinction between the liquid and gas phases. Consequently, surface tension and the enthalpy of vaporization are negligible. Approaching the critical point, fluid thermodynamic and transport properties become highly non-linear, and slight variations in temperature and pressure produce significant changes in fluid thermophysical properties. Considering a pure substance at supercritical pressure, during the injection process, the jet does not break up into droplets or with the formation of ligament structures, but rather, the mixing takes place with the disintegration of the jet through a turbulent process. At that regime, diffuse mixing occurs governed by heat and diffusion phenomena [1,6,18,19,23]. It is traditionally accepted that the key factor in transition phenomena is the surface tension force. In a supercritical regime, it vanishes, not allowing the existence of the so-called phase separation process governed by the atomization mechanism. The same concept cannot be extended to the mixtures, the effects induced by non-ideal mixing cause the rise of the mixture critical point w.r.t the values of the pure propellants. Consequently even at supercritical nominal pressures, during the injection process, in the jet locally subcritical zones may form, which allow the phenomenon of phase separation to occur. As a result, the fluids are subject to surface tension forces despite the high pressure. At the local subcritical regime, the mixture can present itself in thermodynamic states within the coexistence region; droplets and ligaments can form, separated from the vapor phase by a sharp interface, several experimental studies [1,3,15,17,26] have reported evidence of this phenomenon.
From the computational fluid dynamics (CFD) standpoint, transcritical flow simulation has been traditionally performed through Diffuse-Interface methods, mainly with density-based solvers with fully conservative (FC) or quasi-conservative (QC) numerical schemes. In these methods, the fluid is considered a homogeneous mixture, and the surface tension is neglected. The FC schemes coupled with the real thermodynamics have highlighted numerical problems due to the generation of spurious pressure oscillations [8,9,16,28,29], especially in transcritical conditions where the thermodynamic properties exhibit strongly non-linear behavior. Schmitt et al. [28] and Ruiz et al. [27] have developed stabilization methods using artificial viscous terms in the governing equations and setting the pressure differential to the total energy to zero. Alternatively, QC conservative methods were used to overcome the problem of spurious oscillations. QC models solve a pressure-transport equation instead of the total energy equation [9,16,29] or the Double-Flux (DF) method [13,14], where the relation between pressure and internal energy is imposed frozen in time and space. This leads to a quasi-conservative formulation, where the energy is not exactly conserved and the pressure equilibrium across the interfaces is maintained. In this context pressure-based with segregated approaches, despite having shown promising results in a single species context [10], or with moderate stratified flows [5], have found minor applications so far. The solution of a pressure-equation is a promising alternative to avoid the mentioned pressure oscillations, as long as consistent thermodynamic modeling of the involved flow variables is employed. In this contribution are presented: 1. An extensive validation for the evaluation of multicomponent and multiphase fluid mixtures in a wide variety of extreme thermodynamic conditions, ranging from slightly below the critical pressure up to largely supercritical pressure levels. 2. A critical assessment of the numerical issues ensuing from the multicomponent transcritical diffuse interfaces.

Governing Equations
The governing equations for the diffuse interface approach considered in this study are the compressible, inviscid multispecies Navier-Stokes equations; written in a Fully Conservative (FC) form as where is the density, U is the velocity vector, E = e + 1 2 U ⋅ U is the specific total energy, P the pressure, and Y i the mass fraction of species i. The system of Eqs.
(1)-(4) is closed with a real-fluid equation of state. Further details are given in the next section.

Single-Phase Real-Fluid Modeling
At high pressures and in rocket engine relevant transcritical injection conditions, the ideal gas equation of state leads to erroneous results; the attractive and repulsive interactions between fluid molecules gain importance. In these extreme thermodynamic conditions, cubic Equations of State (EoS) are typically used to compute a fluid's thermophysical properties. Among the various cubic EoS available in the literature in this work, the Peng-Robinson (PR) equation of state was chosen due to its relative simplicity and common usage in other works [6,8,9,14,27]. It can be written as [24]: Here, T is the temperature, and R is the universal gas constant. The parameter a = a c (T) can be regarded as a measure of the intermolecular attraction force, and b takes into account real-gas effects related to finite packing volume. Their mathematical expressions are reported in Table 1, where the subscript c denotes the critical value, and is the acentric factor. To extend the PR-EoS to a multicomponent mixture, conventional mixing rules [25] are used for calculating the EoS parameters: where X i is the mole fraction of species i. N s is the number of the species. The off-diagonal a ij terms are evaluated using the pseudo-critical combination rules [25]. The other thermophysical properties are calculated using the departure function formalism [2,25]. Figure 1 shows the comparison between the PR-EoS results against the NIST data [12] for pure Methane, while Fig. 2 shows the results for a mixture of Methane-Oxygen at P = 150bar.

Two-Phase Real-Fluid Equilibrium Modeling
The thermodynamic model presented in Sect. 3.1 is valid under the assumptions of a single-phase state; as previously mentioned, during the mixing and combustion processes, locally subcritical states can form, leading to the formation of the phase separation phenomenon. In such cases, to capture the correct behavior of the mixture's thermophysical properties, a specific thermodynamic model is required. To meet this need, in this section a multiphase thermodynamic framework based on the concept of vapor-liquid equilibrium (VLE) and the cubic PR-EoS is reported. For a homogeneous mixture, it can extend the single-phase approach for the equilibrium coexistence states.

Phase Split Calculation: Isothermal-Isobaric Flash Problem
Given a generic mixture identified from the thermodynamic state (T, P, X) , the thermodynamic equilibrium is reached if,  where f i is the fugacity of component i, the superscript V and L refer to the vapor and liquid phases, respectively. The fugacities equality of Eq. (7) can be rearranged in a more convenient form: y i P are the fugacity coefficients, x i and y i are the species liquid and vapor phase molar fractions, respectively, and their ratio K i is often called K-ratio. Eq. (8) is the starting point for VLE calculation via an equation of state; as for the other thermodynamic variables, the Peng-Robinson EoS is used for the calculation of fugacity. For a generic mixture, the fugacity for component i in both liquid and vapor phases can be specialized as [2] where A and B are the adimensional EoS parameters: and Z is the compressibility factor. The VLE model implies the equality of fugacities for each component in each phase, at the same time, the mass conservation between the mole fraction for each component i, in each phase ( x i and y i ), and the overall mole fraction ( X i ) of the same components, must be satisfied [2]. The mass conservation constraint, and the minimization of Gibbs energy through the fugacities equalities, lead to a non-linear system in the independent variables (T, P, X i ) , and the unknown variables ( , x i , y i ) . Where is an additional thermodynamic property, which indicates the ratio between the mixture number of moles in the vapor (or liquid) phase over the total number of moles. This problem is also known as the isothermal-isobaric flash calculation TPn . The overall mass balance can be written as the initial number o moles denoted by F (Feed composition), and it is separated into V moles of vapor and L moles of liquid: The balance for the i-component, can be written as Combining Eqs. (12) and (13), the phase fractions explicit form can be derived as Since the two phase fractions sum to one, hereafter only the vapor phase fraction V = is considered. Eq. (13) can be recasted as where x i and y i satisfy the condition ∑ N s i y i = 1 and ∑ N s i x i = 1 . This two constraint can be manipulated using the definition of the K i ratio and Eq. (15) to lead the Rachford-Rice equation [2], a unique relation which provides the link between x i , y i , and X i . therefore, the TPn problem is a non-linear system of N s + 1 equations in the unknowns x i ,y i and for a given thermodynamic state (T, P, X i ): System 17 is solved for a given temperature T, pressure P and composition X, by implementing a dedicated algorithm that combines the Successive-Substitution-Iterative (SSI) and Newton-Raphson methods (see [4,11,21,22] for practical implementation details). It is very important to keep in mind the temperature, pressure and composition dependency of the unknowns: = (T, P, X) , x = x(T, P, X) , y = y(T, P, X) . Additionally is important to underline the assumption of mechanical equilibrium P L = P V = P has been used.Considering mechanical non-equilibrium would mean providing a model for the pressure ratio, which involves capillarity effects. Figure 3 shows the validation of the VLE solver, in the pressure-composition phase diagram, against experimental data for two binary mixtures, Methane-Nitrogen 3a, and

3
Ethane-Heptane 3b. The symbols denote the experimental data, taken from [7,20]. The dashed lines are the bubble points, while the dew points are plotted as solid lines. At a fixed temperature, the pressure along the bubble-line is the pressure where an infinitesimal bubble of vapor coexists with liquid, while the pressure along the dew-line is the pressure where an infinitesimal drop of liquid coexists with vapor. The region enclosed by these curves is the twophase region. The comparison with the experimental data shows reasonable agreement for all conditions. In Fig. 4, an example of the VLE's feature is reported. Here the phase fraction (T, X) for a hydrogen-oxygen mixture, and (P, X) for a Nitrogen-Hexane mixture, fields are shown. The figures represent the multi-phase mapping encountered in the thermodynamic space. In particular, the blue color ( = 0) denotes the state in which the mixture is in a purely liquid phase, the red color ( = 1) in a purely gaseous one, while the color changes (inside the coexistence VLE dome) denote the multiphase states.

Real-Fluid VLE Thermodynamics
Once the phase fraction and the composition of the mixture in the respective phases are known, a specific formulation to define the properties of the mixture in phase equilibrium is necessary. Additional blend rules are required to define the coexistence of bulk properties; as suggested in other works [30,31], two properties are required to be blended, the mixture density and the internal energy: The phase-related quantities are obtained through the singlephase model, but with the appropriate phase composition ( x i or y i ) calculated with the VLE solver. To have a direct interpretation of VLE effects on the thermophysical properties, an equimolar X = [0.5, 0.5] Methane-Oxygen mixture is analyzed. To better visualize this effect, three deeply subcritical cases are considered for the properties evaluation, respectively, at 20, 30, and 40 bar. The pressure-temperature phase diagram is shown in Fig. 5 from which it can be seen, how the selected cases are all below the mixture critical pressure. Figure 6 shows the comparison between the density and the isobaric specific heat of the mixture, obtained with the single-phase assumption (dashed lines) against the VLE model (solid line). As can be seen for the temperature values which fall inside the vapor dome, the density calculated with the single-phase assumption exhibits a non-physical discontinuity, this is since the minimum Gibbs criterion suddenly jumps between the phases and selects a completely different compressibility root. In contrast, the VLE model captures the transition between the liquid and vapor phase, exhibiting a smooth function.

Diffuse Interface Methods
Two approaches are used in the context of a densitybased solver, namely the standard (FC) scheme and the (QC) one, where the latter features the Double-Flux (DF) method extended by Ma et al. [13] for transcritical flows. The system of Eqs. (1)-(4), is discretized using the Finite-Volume-Method (FVM), using an in-house solver written in MATLAB environment. A first-order spatial integration is employed for both approaches, and the Harten-Lax-van-Leer-contact (HLLC) approximate Riemann-solver is used for the flux calculation. The temporal integration is performed using the low-dissipative third-order Runge-Kutta (RK) scheme. Two cases of onedimensional (1D) advection problems are simulated, which are common benchmark problems used in the literature to compare the performances of the numerical schemes against the exact solution. Before analyzing the test cases, here a brief outline of the differences between the two methods is presented, and details of the DF-model can be found in the Ma et al. work [13]. The DF method uses an effective heatspecific ratio ( * ) , and an effective reference energy (e * 0 ) , to relate the pressure and internal energy [13]: where c is the speed of sound. In Algorithm 1 is summarized the time advancement procedure from time t n to time t n+1 for the QC double-flux scheme. The idea behind the double flux model is to compute and store * and e * 0 for each computational cell and for each time-step advancement. Then the conservative Euler flux reconstruction and computation at the faces of the cell i are evaluated, particularly the total energy is computed from * and e * 0 for cell i through [13]: where j denotes the spatial stencil of cell i. After that, the primitive variables for each cell are updated, specifically, the pressure is computed using [13]: It is very important to underline the fact that, in this step * and e * 0 are frozen for cell i in time. When the pressure is updated, the temperature can be calculated analytically for , cubic equations of state like the Peng-Robinson one; from the pressure, density, and species composition. After the RK sub-stages, the total energy is updated using the primitive variable, this step ensures that the final thermodynamic state is consistent between primitive and conservative variables.
To resume, with the DF treatment, the link between pressure and internal energy is frozen in time and space, converting the local system to a calorically perfect gas system with a constant specific heat ratio [13], making the spurious pressure oscillation stable. The main drawback with the DF formulation is the loss of total energy conservation, due to the discrepancies of energy fluxes at faces, induced by different * , e * 0 values calculated for each cell. The procedure for the time advancement of the FC scheme is reported in the Algorithm 2. For the FC scheme, the conservative variables are updated without using the DF formulation for the total energy (Eq. (21)), then the primitive variables are calculated. First, the temperature is calculated given the density, internal energy, and species with an iterative process. After that, the pressure is evaluated from temperature, density, and species with the EoS, Eq. (5). Due to the strongly endemic non-linearities present with transcritical flows, non-physical pressure oscillations are generated with that scheme, leading to the failure of the solver. Contrary to the QC scheme, the conservative fluxes calculated from the face of two adjacent cells are the same, ensuring the full conservation of the flow variables.

1D Transcritical LN 2 -GN 2 Interface Advection
The transcritical interface is a contact discontinuity between liquid-and gas-like nitrogen. The operating supercritical pressure is 5 MPa, and the advection velocity is 50 m/s. The computation domain has a length of L = 1 m. The initial conditions are reported in Table 2. A one-point jump condition is initially applied, in which the interface exhibits the transcritical pseudo-phase change across it. The initial discontinuity thickness is 1 mm, which is a consistent value in supercritical pressure conditions [6]. Three different uniform grids were used for the calculation, with N x = 100, 1 3 500 and 1000 points, corresponding with a grid spacing of Δx = 10, 2, and 1 mm respectively; both left-and rightside boundaries are imposed as a periodic. In this kind of test the initial pressure, velocity, temperature, and density will be maintained after one periodic advection cycle, which is t = L∕u = 0.02s . For all simulations, the Courant-Friedrichs-Lewy (CFL) number is set to 0.5. In Fig. 7, the results of the FC scheme are reported, with that method the pressure and temperature are directly calculated from the internal energy and the density. Due to the non-linearities present in the proximity of the critical point, small variations of the calculated flow-field density and internal energy, due to numerical diffusion, cause the non-maintaining of the pressure and velocity equilibriums. This leads to several numerical errors in the estimate of the temperature field, which can be partially mitigated by increasing the grid resolution. The QC method, results are shown in Fig. 8, maintains the pressure and velocity equilibriums, without the generation of spurious oscillation for each grid resolution. Density and temperature are also captured with good accuracy for medium to high-resolution grids. As mentioned earlier, the main drawback of quasi-conservative methods is the non-perfect energy conservation, these numerical errors are caused by the forcing of the pressure equilibrium during the interface advection. Figure 9a, b shows the conservation loss of the conservative variables with the two schemes; the error as a function of time has been quantified as where the conservative variables are W = [ , u, E] T . As can be seen the conservation losses on , u have the same order of magnitude for both approaches. While on E , as expected, there are several orders of magnitude of difference (see Fig. 9c). At first glance, this would lead to thinking of greater deterioration in the accuracy of the temperature field capture by the QC method; however, this behavior is not highlighted in the results in Figs. 7 and 8. From a more careful analysis reported in Fig. 9d, it emerges that, after an advection cycle, the percentage error on the specific energy (and as a direct consequence on the estimate of the temperature flow field) is always lower for each grid resolution. With a decreasing discrepancy as to the number of points increases. This is due to a greater influence of the pressure equilibrium on the estimation of the thermophysical properties of the flow.

1D Transcritical LOx-GH 2 Interface Advection
The second test case under analysis is in the context of a multispecies flow, considering the advection of a contact discontinuity between liquid Oxygen LO x and gaseous Hydrogen GH 2 . The pressure is set to 100 Bar and, the advection velocity is 50 m/s. The computational domain has a length Fig. 9 Errors analysis. a conservation losses for the FC scheme; b conservation losses for the QC scheme; c comparison of total energy losses between FC and QC schemes; d specific total energy errors percentage after one advection cycle as a function of grid resolution Fig. 10 Initial conditions of the LOx-GH 2 advection test of L = 1 nm. A representation of the initial conditions is reported in Fig. 10. The computation is conducted with an N=1000 points uniform grid. The simulation is run for t = L∕u = 20 ps. Figure 11 shows the comparison between the two methods. As in the previous case, the QC schemes maintain the pressure and velocity at the same values as the initial conditions. The density field is in good agreement w.r.t the analytic solution, while the temperature exhibits a deviation in the proximity of the discontinuity center, due to the non-perfect energy conservation. On the other hand, the FC scheme does not maintain pressure and velocity equilibrium. The density is well captured, and the overall underestimate of the pressure field from its uniform value, affects the temperature field, particularly in correspondence to the Hydrogen side.

Conclusions
A description of a comprehensive theoretical and numerical framework for simulating transcritical flows under LRE conditions has been presented. A high-fidelity thermodynamic model based on a cubic equation of state and vapor-liquid equilibrium assumptions that can effectively describe transcritical mixtures in a wide variety of conditions has been validated. This modeling framework has allowed the discussion of the impact of VLE assumptions on the fluid's thermophysical properties calculation. In the context of transcritical interfaces also the well-known numerical issue caused by the coupling between thermodynamics and governing equations has been discussed. A fully conservative and a quasi-conservative scheme have been compared showing the superior performance of the second approach. The problems arising from the non-maintenance of pressurevelocity equilibria for the FC schemes and the non-conservation of energy for the QC schemes have been addressed. So far, although the consequences that these behaviors have on the mixing regimes for jets and flames have been addressed, there are no numerical tools available that guarantee stability and high fidelity in the capture of the physical regimes that are exhibited in transcritical conditions. The experience acquired and the results obtained with the QC method, will be used for future development of a pressure-based numerical framework able to simulate high-fidelity, transcritical flows in extreme-stratified conditions.