A tabulated real-fluid modeling approach applied to cryogenic LN2-H2 jets evaporation and mixing at transcritical regime

Accurate, robust, and fast fully compressible real-fluid simulation of fuel jets is today’s one of the highly debated topics in various research laboratories and industries. Indeed, the use of real-fluid equations of state has proved to be computationally very expensive and is one of the current challenges. In this paper, a thermodynamic equilibrium tabulation approach is proposed to overcome most of the limitations and to make real-fluid simulations affordable in the industry. This tabulation approach is implemented in the Converge software as a closure to the fully compressible two-phase and multi-component real-fluid model (RFM). This modeling approach has been applied to the simulation of a classical cryogenic injection of liquid nitrogen (LN 2 ) coaxially with warm hydrogen (H 2 ) jet, in a transcritical regime using thermodynamic tables generated by two different equations of state: Peng-Robinson (PR), and Soave-Redlich-Kwong (SRK). The numerical results are compared to available experiments and published numerical studies. The computational efficiency, accuracy, and robustness of the proposed RFM model are thereby confirmed.


1-Introduction
Generally, using a real fluid equation of state (EoS) has proved to be computationally timeconsuming, especially for three-dimensional simulations [1]. Indeed, it is found that the complexity of cubic EoS and solving iteratively Vapor-Liquid-Equilibrium (VLE) using the isoenergetic-isochoric flash set of equations are computationally expensive [1][2][3]. Hence, an accurate, robust, and most importantly a faster real-fluid simulation is required. One of the remedies to elevate this problem is using a VLE tabulation approach before the beginning of the simulation and stores the required thermodynamic properties values into a table. Then, interpolation is carried out in this table during the simulations. This interpolation process is usually called "a look-up table" in the literature. Indeed, there are some studies focused on thermodynamic property tables [4][5][6][7][8][9][10] proposing various methods for tabulation, interpolation, and look-up of data. Azimian et al. [5] introduced an artificial neural network model for generating the water/steam thermodynamic tables. Lorenzo et al. [6][7] suggested a 2D look-up table method for water-steam simulation by transforming the irregular shape of the thermodynamic domain of − to the regular shape so that the nonlinear behavior of thermodynamic properties can be predicted using bicubic interpolation. De Föll et al. [8] proposed a tabulation mechanism based on a projection approach for single-component fluids using a quadtree data structure. Praneeth and Hickey [9] introduced a systematic error quantification and computational cost estimation of EoS for compressible single-component tabulation. They reported that, due to the inherent thermodynamic non-linearities around pseudo-boiling points, small errors in the thermodynamic property estimation can lead to big variations in the computed pressure and temperature, which directly affects the thermal and transport properties calculations in the flow. Most of the above-mentioned studies relating to the tabulated thermodynamic properties have focused on single-component fluids. However, few studies have been carried with multi-component fluid mixture tabulation. One of the notable studies has been done by Brown et al. [4]. They suggested a multi-component twophase homogeneous equilibrium model (HEM) and flow solver for binary and quaternary constant and uniform concentration mixtures, allowing a two-dimensional (P,T) based adaptive inverse interpolation approach. Another notable study has recently been done by Koukouvinis et al. [10]. In that study, the authors proposed a uniform tabulated thermodynamic approach based on NIST-database and VLE using PC-SAFT EoS for the simulation of the ECN spray A [10]. They reported that by use of the PC-SAFT and log10(P)-T tabulation approach, the results are independent of table size when bilinear interpolation is adopted, even with linear mixing rules. As a result, to sum up, most of the available studies show that the tabulation should be highly accurate, with low-storage requirements, and a fast lookup algorithm. However, it seems not clear how to make real-fluid simulations affordable in the industry. In this study, a fully compressible multi-component two-phase real fluid model (RFM) [1,2,11] has been closed using a generalized three-dimensional (3D) tabulation method. In this method, an in-house thermodynamic library IFPEN-Carnot is used to generate the 3D- Table (with T-P-Y as the axis for a binary mixture) based on various real fluid EoS. This thermodynamic library involves many EoS such as Peng-Robinson (PR), and Soave-Redlich-Kwong (SRK), which represent also one of its main advantages. The tabulated thermodynamic properties of the mixture are calculated as a function of temperature (T), pressure (P), and species composition (Y). Noteworthy, the uniform look-up table is constructed based on an isothermal-isobaric flash (TP flash) [12]. This tabulation method also improves the efficiency of the tabulation approach previously developed by Yi et al. [11] based only on PR-EoS [1,2]. The RFM model along with the generalized 3D tabulation method has been implemented in the CONVERGE v3 software [13]. The present work attempts to investigate the accuracy of the proposed tabulation methodology compared to Matheis et al. [14] numerical results obtained using VLE directly solved during their simulations. Hence, section 2 is dedicated to a brief explanation of the numerical and theoretical method including governing equations, thermodynamic tabulation, and look-up methods. Section 3 describes the simulation setup for the cryogenic co-injection of liquid nitrogen (N2) and gaseous hydrogen (H2) at transcritical conditions while two different EoS: PR and SRK in the cryogenic conditions. Then, the simulation results are examined and compared with Matheis et al. [14] numerical results. Finally, the conclusions are presented in section 4.

2-1 Governing equations
The real-fluid RFM model used in the current work is presented in this section. It is a fully compressible two-phase and multicomponent flow model which consists of the four balance equations (1)(2)(3)(4) written under the assumptions of full thermodynamic equilibrium, as follows.
Where ( , , , , ) are the mixture's density, velocity, pressure, temperature, and specific internal energy, respectively. ( , ℎ )are mass fraction and enthalpy of species respectively. ( ) is the total number of species. The thermal conductivity ( ) and the dynamic viscosity ( ) cover laminar and turbulent contributions. The laminar contribution of ( , ) is computed by Chung et al. [15] correlations. The turbulent conductivity is calculated using a given turbulent Prandtl number, and the turbulent viscosity is computed by the adopted turbulence model. In this work, the standard turbulent subgrid-scale Smagorinsky model is used in the large-eddy simulation (LES) framework [16]. Also, the laminar and turbulent diffusion coefficients are estimated using a given Schmidt number ( ) as = / and = / , respectively.

2-2 Thermodynamic closure, tabulation look-up, and its validation
The current work has adopted a tabulation approach where the thermal and transport properties as well as the phase state and composition are tabulated before the CFD simulation. During the simulation, the various tabulated parameters are robustly interpolated over the entire range of the thermodynamic states based on three inputs for the tables which are (T, P, Ym). The coupling between the flow solver and the thermodynamic table is implemented in CONVERGE v3 software [13]. The flow solver offers pressure-velocityenergy coupling using a modified Pressure Implicit with Splitting of Operator (PISO) [17]. A second-order central difference scheme with flux limiter is used for the spatial discretization of each equation. The time discretization is achieved by the first-order Euler scheme. The automatic mesh refinement (AMR) feature is also used in this work, as described in section 3. The inverse-distance weighting (IDW) is implemented for the table interpolation for binary mixtures [11]. The 3D-table is generated using an in-house thermodynamic library named IFPEN-Carnot. This library uses robust isothermal-isobaric (TP) flash [1,2,12] algorithms coupled to various real-fluid EoS to compute the Vapor-liquid equilibrium (VLE) as well as the thermal, transport properties, and phase composition. The interpolation in the table during the simulation is mainly carried out using two main functions: 1-A Table look-up function: compute the thermal, transport properties as well as the phase state using the inputs (T, P, Ym). 2-A Reverse look-up function: compute/update the temperature using the inputs (e, Ym, P). Also, the Inverse-Distance Weighting method is used for the interpolation of the tabulated quantities of the 3D table. A general form of finding an interpolated value at a given point based on samples = ( ), = , … using IDW can be expressed as equation 5, where ( ) = ( , ) , denotes an interpolated (arbitrary) point, is an interpolating (known) point, is the given distance from the known point to the unknown point , is the total number of known points used in each local interpolation. In equations 6-7, the SRK-EoS and PR-EoS are presented along with the various thermodynamic properties that are calculated. A general form of a cubic EoS can be written as equation 6. Two commonly used values for calculating the pure component parameters for two EoS are listed in Table 1. Also, when the SRK-EoS or the PR-EoS is used for mixtures, Van-der Waals mixing rules are applied as equation 7. Figure 1 shows the prediction of density for H2 and N2 from both EoS in comparison with NIST.  a) Nitrogen density distribution b) Hydrogen density distribution Figure 1 Computed density for N2 and H2 with SRK-EoS and PR-EoS using IFPEN-Carnot Library compared to NIST data at 4MPa [20].
Noteworthy, this Figure has shown that for cryogenic conditions the SRK-EoS gives better predictions of density. However, PR-EoS has been shown to have a better representation of the thermodynamic properties close to the critical point [18][19].

3-The coaxial cryogenic LN2/H2 test case 3-1 Simulations setup
The study of coaxial cryogenic nitrogen with a co-flow of warm hydrogen into supercritical nitrogen was done employing Smagorinsky turbulence model in the LES framework. The thermodynamic models are based on two different EoS, PR as well as SRK in order to compare theirs performances. In modern Liquid Rocket Engines, the chamber operating pressure (4MPa) lies above the critical pressure of the propellants, (Pc=3.35MPa for N2 and Pc=1.296 MPa for H2). At these conditions, known as transcritical conditions, the fluid properties significantly differ from an ideal gas. Transcritical jet behavior differs significantly from injection at lower pressures. In the last decades, several research has led to a better understanding of the cryogenic rocket engine processes [18-19, 14, 21]. In this paper, we will study a selected operating point of a series of experiments from [21] in which quantitative density measurements were obtained in a coaxial LN2 and GH2 jets at supercritical conditions (compared to the pure N2 and pure H2 critical pressures). Figure 2(left) shows a schematic of the experimental setup that has been investigated [21], along with the computational domain highlighting the AMR region near the nozzle. The Gaseous GH2 (outer stream) at T=270 K and LN2 (inner stream) at T=118 K are injected through a coaxial injector into a cylindrical tank (D = 10 cm) initially filled with GN2 at 4 MPa and 298.15 K. The inner and outer diameter of the GH2 annulus is DH2,i = 2.4 mm and DH2,o = 3.4 mm, respectively. The inner LN2 injector diameter is Di = 1.9 mm. The inlet velocity of LN2 and GH2 are 5 m/s and 120 m/s, respectively. To assess the VLE calculation based on the EoS, the VLE of a binary mixture of N2-H2 has been computed, and will be compared with the data from [14] in next section. A 3D simulation setup was employed for which the computational grid has been generated in a rectangular domain with the dimensions Lx= 100 mm in the streamwise and Ly = Lz = 40 mm in the lateral directions, similar to [14]. The base mesh cell size is set to 0.5 mm and refined to 0.1 mm in the flow field using adaptive mesh refinement, AMR. The subgrid criterion for adaptive mesh refinement has been chosen for the cells based on the minimum velocity of 0.1 m/s. The time step is automatically adjusted based on a maximum CFL number of 0.5 reaching a value in the range of (10 -8 s-10 -9 s). Typical CPU time for the simulation discussed below is 240 hours which is much more efficient than our previous simulations [1,2] using VLE thermodynamic flashes during the simulations. Figure 3 depicts an instantaneous snapshot of the H2 density, temperature, and phase indicator distribution (PHI=0, 1, 2 denotes liquid, gas, and two-phase state, respectively) using PR, and SRK EoS, respectively. In this Figure, the cryogenic LN2 at 118K may be identified as the dark core in the temperature contour (Figure 3(a , d)) which is surrounded by a co-flow warmer GH2 injected at 270K into the GN2 at 298.15K and 4MPa. Because the latter pressure is higher than the critical pressure of both propellants, (Pc=3.35MPa for N2 and Pc=1.296 MPa for H2) and the cryogenic LN2 is injected at 118K a temperature much smaller than its critical temperature (Tc,N2= 126.2 K), the jet flow belongs to transcritical flow regime, It is worth mentioning that the critical line of H2-N2 mixture reaches higher pressures than 4MPa. In this regime, Crua et al. [23] experiments varying the ambient pressure above pure fuels critical pressure have revealed the existence of a subcritical relatively sharp liquidgas interface that is progressively thickened [22] and replaced by a diffusion-dominated mixing regime at a sufficiently high ambient pressure. These experiments have been corroborated recently by Yang et al. [1] LES simulations of the Spray-A injector in the transcritical regime. However, when and how this transition to a diffusion-dominated mixing regime occurs is not well understood. According to Yang et al. [1], this transition should take place when the ambient pressure exceeds the mixture critical pressure. Nevertheless, it is worth recalling that the latter is a local variable that depends on the local composition. Similar to Diesel injection [1], a two-phase region around the liquid core is identified by the phase indicator (PHI=2, yellow color in Figure 3(c), (f)). It is interesting to see also in Figure  (3(b), (e)) that the temperature within the mixing layer in the two-phase region drops below its inflow value of 118 K to approximately 114 K. Another remarkable effect also discussed by Matheis et al. [14], can be observed for the H2 partial density in Figure (3(a), (d)) for both EoS. Indeed, H2 is injected with a density of 3.56 kg/m 3 and reaches a much higher partial density of almost 4.7 kg/m 3 within the two-phase region. As the first milestone for this study, the binary phase diagram for N2-H2 has been compared with the study of Matheis et al. [14], and very similar numerical results are obtained, as shown in the T-zH2 diagram is plotted in Figure 4. Pure LN2 (zH2 = 0) at 118 K as injection temperature, and GN2 at 298 K as reservoir temperature as well as zH2 = 1 for pure H2 at 270 K, on the right-hand side, can be identified in Figure 4. In between, either cryogenic LN2 from the main injection temperature or warm GN2 from the reservoir is mixed with the warm GH2. Also, adiabatic mixing temperature (TAM) as an analytical approach for calculating the ideal mixture temperature is plotted in Figure 4. This (TAM) is obtained by solving the enthalpy balance equation as  Figure 4. On the second hand, the highest deviation of the scatter plots is obtained in the fully turbulent mixing region between the injected GH2 layer and ambient GN2 as shown around the dashed line in the upper side part of Figure 4. Figure 5 also shows the temporal evolution (at 2 and 3ms) of mixture density, and H2 mass fraction in the central cut section of the geometry using SRK EoS. In the following, to validate these LES simulations, the comparison between the experimental data [21] and the numerical study of Matheis et al. [14] as well as the current study are shown in Figure 6. The N2 radial and axial partial density profiles extracted at 4 mm from the nozzle exit and the centerline, respectively are depicted in Figures 6(a, b), using PR and SRK EoS. With PR-EoS, very similar radial profiles are obtained when compared to Matheis et al. [14]. Both studies have overestimated the maximum LN2 density value of around 608 kg/m 3 , compared to the experimental value of 390 kg/m 3 [21]. However, the current results using SRK-EoS show a smaller deviation from the experiments with a maximum of LN2 density value around 530 kg/m 3 . This better prediction of SRK compared to the PR for cryogenic fuels has also been reported by other researchers in recent years [18][19]. Due to the H2 partial density high value (~2-4 kg/m 3 ) close to the injector axis in Figure 6(c), the injection conditions seem to be at the origin of these large deviations with the experiments. Figure (6(c), (d)) also illustrates the radial and axial partial density profile of H2. Figure 6(c) displays maximal values of the radial H2 density distribution. This Figure shows similar profiles as in [14] and the same trends as experiments [21]. Also, Figure 6(d) shows a fairly good agreement with small deviations between the two different EoS and the experimental results. Matheis et al. [14] explained that due to the conformity of EoS with the NIST [20], as shown in Figure 1, the deviation with the experiments cannot be attributed only to the inaccuracy of the EoS, and it may be the result of uncertainties in the experimental measurements.

4-Conclusions
To gain computational efficiency and robustness, the fully compressible multicomponent real-fluid (RFM) model has been coupled to the IFPEN-Carnot thermodynamic library using a generalized tabulation method. The phase change quantity, thermal and transport properties, and the thermodynamic derivatives such as sound speed, heat capacity, volume fraction, etc are computed efficiently using a bijective look-up tabulation method. This modeling approach has been applied to the simulation of a classical cryogenic injection of LN2 coaxially with a warm GH2 jet, in a transcritical regime using two different equations of state: PR, and SRK. The model predictions have been compared with experimental observations [21] and available numerical simulation results [14]. The main achievements of this study are as follows. 1-The computational efficiency, accuracy, and robustness of this proposed tabulated RFM model as a remedy to the direct evaluation of costly phase equilibrium solver have been confirmed. 2-The simulation results demonstrated some interesting thermodynamic phenomena such as condensation because of the increasing of the H2 partial density within the turbulent mixing layer in which the temperature has been found to decrease as shown in figure 3. 3-The current study has confirmed the importance of using a powerful EoS in the modeling of fluid behavior. These investigations have well-illustrated that SRK EoS has a better prediction of fluid density compared to the PR EoS, in particular for the H2, corroborating the results reported by [14]. 4-Multiple realizations and their ensemble averaging will be performed in future work in order to analyse the turbulent behaviour of such flows. a) Radial N2 density profile at 4mm b) Axial N2 density profile at the centerline c) Radial H2 density profile at 4mm d) Axial H2 density profile at ( ) Figure 6 Axial and radial density profiles compared to the experimental [21] and the LES study Matheis et al. [14]