Large eddy simulations of reacting and non-reacting transcritical fuel sprays using multiphase thermodynamics

Accurate simulations of high-pressure transcritical fuel sprays are essential for the design and optimization of next-generation gas turbines, internal combustion engines, and liquid propellant rocket engines. Most important and challenging is the accurate modelling of complex real-gas effects in high-pressure environments, especially the hybrid subcritical-to-supercritical mode of evaporation during the mixing of fuel and oxidizer. In this paper, we present a novel modeling framework for high-fidelity simulations of reacting and non-reacting transcritical fuel sprays. In this method, the high-pressure jet disintegration is modeled using a diffuse interface method with multiphase thermodynamics, which combines multi-component real-fluid kinetic and caloric state equations with vapor-liquid equilibrium calculations in order to compute thermodynamic properties of the mixture at transcritical pressures. All multiphase thermodynamic formulations are presented for general cubic state equations coupled with a rapid phase-equilibrium calculation method. The proposed method represents multiphase turbulent fluid flows at transcritical pressures without relying on any semi-empirical break-up and evaporation models. Combustion source terms are evaluated using a finite-rate chemistry model, including real-gas effects based on the fugacity of the species in the mixture. The adaptive local deconvolution method (ALDM) is used as a physically consistent turbulence model for large-eddy simulation (LES). LES results show a very good agreement with available experimental data for the reacting and non-reacting ECN Spray A at transcritical operating conditions.


Introduction
In-depth understanding of turbulent reacting multiphase flows at transcritical pressures is essential for the design and optimization of efficient energy conversion systems, such as liquid rocket engines, or modern diesel engines and gas turbines. Such systems typically work based on rapid injection of cold liquid or liquid-like supercritical fuels into a chamber with an elevated pressure and temperature. Combustion occurs after transcritical evaporation during mixing of the fuel with a hot and pressurized gas or gas-like supercritical oxidizer. Transcritical pressure refers to an operating pressure higher than the critical pressure of the pure fuel or oxidizer streams, but lower than the cricondenbar pressures of their possible mixtures. Since the cricondenbar point of hydrocarbon/air mixtures is unreachable even at elevated pressures of several hundred bars, the operating pressure of most energy conversion systems is transcritical in practice.
Whereas the operating condition of the transcritical chamber is nominally supercritical with respect to the fuel stream and a direct phase change from liquid to supercritical is expected, aligned with Gibbsian thermodynamics and experimental reports [1][2][3][4], the transcritical mixture can locally enter the two-phase regime and interfaces between liquid-like and gas-like phases may form. Transcritical injectors, therefore, resemble a hybrid combination of the classical two-phase disintegration with breakup and evaporation of droplets and the supercritical turbulent mixing of two dense fluids. This hybrid behavior complicates their numerical simulation; despite comprehensive experimental campaigns, it is not well understood which type of the jet disintegration dominates under exactly which conditions. Traditional large-eddy simulation (LES) of the high-speed transcritical injection have either modeled the transcritical multiphase fluid flows by Lagrangian particle tracking (LPT) methods with sharp vapor-liquid interfaces [5][6][7], or by Eulerian single-phase dense-gas (DG) approaches with diffuse vapor-liquid interfaces [8][9][10]. Both LES-LPT and LES-DG approaches can be justified, but they have some inherent limitations at the transcritical conditions. The standard LPT method is very sensitive to empirical tuning parameters and was developed for the dilute mixtures, neglecting the real-fluid thermodynamics. The LES-DG approach, on the other hand, excludes the effect of the transcritical phase separation and may lead to nonphysical or ill-defined states when a part of the flow passes the meta-stable boundaries, specially at lower transcritical pressures [11]. Furthermore, some important transcritical effects such as the high solubility of the saturated liquids or the different components' evaporation rates in case of surrogate fuels are not captured by these models [12][13][14].
Using multiphase thermodynamics (MT) in the context of the diffuse-interface method has been demonstrated recently to be a promising technique to overcome the aforementioned limitations [15,16]. In this formulation, the fully conservative Navier-Stokes equations (NSE) are solved for a hypothetical multi-component fluid mixture with thermotransport properties computed using a suitable equation of state (EOS) coupled with vapor-liquid equilibrium (VLE) calculations. Although the weak surface tension force is neglected, the method can accurately capture the physics of the problem including the subcritical region of coexisting multi-component vapor and liquid as well as the real fluid behavior such as dissolution of the ambient gas in the compressed-liquid.
Whereas LES-MT studies show excellent agreement with experimental data for the non-reacting transcritical sprays [16][17][18], the applicability of the method to reacting flows remained as an open question mainly due to the high computational cost of the vaporliquid equilibrium (VLE) calculations, which increases with the number of components in the mixture, and due to the need for chemistry models that can also capture the departure from ideal-gas behavior in regions with high pressures and low temperatures. The latter can be solved using an appropriate reduced reaction mechanism [19] and utilizing the fugacity values of species in the mixture for the evaluations of real fluid effects on the reaction rates [20].
The real fluid behavior becomes apparent when using an appropriate EOS. Twoparameter cubic EOS, such as the Soave-Redlich-Kwong (SRK) [21] and Peng-Robinson (PR) [22] models, are very popular because of their computational efficiency [15,[23][24][25][26]. The intrinsic drawback of two-parameter cubic EOS is the usage of a universal critical compressibility factor, which can severely limit their accuracy close to the critical point. To overcome this limitation, volume translation methods can be employed to improve the density predictions, but the calculation of consistent caloric properties is computationally very expensive and, thus, simplified approximations are typically used in practice [27]. Another option is utilizing a cubic EOS with three parameters. The Redlich-Kwong-Peng-Robinson (RKPR) model of Cismondi and Mollerup [28] calibrates the cubic EOS using the actual value of the compressibility at the critical point, and provides a consistent framework for real-fluid thermodynamic modeling at tractable computational cost [29].
To alleviate the high computational cost of VLE calculations of reacting fluids, which typically have a very large number of components, Fathi and Hickel [30] have recently introduced a new phase-splitting method that is formulated in a reduced space and based on the molar specific values of the volume function. In comparison to the conventional method (e.g., Refs. [16,31]), the reduction method prevents the quadratic growth of the computational time with the number of components. The formulation based on the volume function leads to a better convergence behavior near the critical point and phase boundaries. In addition, the novel method directly applicable to isoenergeticisochoric conditions aligned with the transported variables in conservative compressible Navier-Stokes flow solvers. Combined, this yields a very considerable speed-up for the simulation of transcritical flows with many components [30].
In this paper, we present several novel physical and numerical models for the highfidelity simulation of turbulent reacting and non-reacting multiphase flows at transcritical pressures. The framework is provided by a fully conservative formulation of the multi-component compressible Navier-Stokes equations, using the adaptive local deconvolution method (ALDM) for LES turbulence modelling [32], the RKPR EOS [28] coupled with the newly-introduced rapid VLE calculator [30] for transcritical vaporization and real-fluid properties, and fugacity-based finite rate chemistry for combustion modeling. The proposed LES-MT method is validated by comparing computational results with experimental data reported for the transcritical reacting and non-reacting Engine Combustion Network (ECN) Spray-A benchmark test cases.

Physical and Numerical Models
We use the LES-MT method for solving reacting multiphase compressible NSEs in a fully conservative form with real-fluid thermo-transport properties and fugacity-based finite-rate chemistry. In this section, the required physical and numerical models for the MT-based simulations are presented.

Governing Equations
The three-dimensional compressible reacting NSEs describe the conservation of mass, species, momentum, and total (absolute) internal energy: where ρ is the mixture mass density, u is the velocity vector, p is the thermodynamic pressure of the mixture, and e t = e + |u| 2 /2 is the total absolute specific internal energy of the mixture and e is the corresponding absolute specific internal energy. For species i = 1, 2, ..., N with N being the total number of components comprising the mixture, Y i is the mass fraction, j i is the diffusive mass-flux vector, andω i is the net mass production rate. τ and q denote the viscous stress tensor and the vector of heat fluxes. I is the unit tensor and ∇ is the gradient operator. The viscous stress tensor is modeled by assuming a Newtonian fluid and Stokes' hypothesis τ = µ ∇u + (∇u) T − 2/3µ(∇ · u)I.
The molecular viscosity coefficient µ of the multi-component mixture is estimated using Chung's correlations [33]. We neglect bulk viscosity effects in our computations because accurate models to be used at pressures close to the critical point are unavailable. For multi-component fluid flows, the total heat flux vector consists of heat conduction and interspecies enthalpy diffusion, and is a function of the thermal conductivity of the mixture λ, temperature T , and the partial mass absolute enthalpy h i of component i. Similar to the viscosity, the thermal conductivity of the mixture is modeled by Chung's correlations [33]. By neglecting the Soret and Dufour effects, the mass diffusion j i is modeled using Fick's law through a simplified correlation based on the mixture averaged diffusion approximation along with an extra term to ensure zero total mass diffusion The effective binary diffusion coefficient between species i and the bulk mixture is approximated by where X i is the mole fraction of species i, which can be computed via X i = Y i M/M i from the mass fractions and the mean molecular weight M = 1/( N i Y i /M i ) of the mixture, with M i being the molar mass of the species i. The product of density and the binary mass diffusion coefficient of species i and j, ρD ij , can be approximated accurately with Chapman and Enskog theory [33] without high-pressure corrections if the system pressure is under 100 bar [34].

Multiphase Thermodynamics
Simulations of multi-component fluid flows at elevated pressure require proper kinetic and caloric EOS that account for the volume of molecules and interactive forces between them in order to accurately calculate pressure and temperature from the density, internal energy, and mass composition of the mixture. The partial mass enthalpies of all components in the mixture are additionally required for the evaluation of the interspecies enthalpy diffusion flux. At transcritical pressures, additional phase-splitting calculation should be carried out to account for the coexistence of vapor and liquid phases. In this section, we present the MT relationships required for the computations of temperature, pressure, and partial mass enthalpies of a general real fluid that can be single phase or two phases in equilibrium.

Kinetic Equation of State
A kinetic EOS represents the thermodynamic relationship between the pressure, specific volume (inverse of the density), and temperature of a single-phase substance that can be pure or a multi-component mixture. The most popular kinetic EOS is the ideal gas (IG); however, its applicability is limited to gaseous media with a compressibility factor close to unity, that is, it only appropriate at relatively high temperatures and low pressures. For the typical pressures and temperatures of transcritical fuel sprays, cubic EOS are an attractive compromise between accuracy, model complexity, and computational cost.
The most widely used cubic EOS are SRK [21] and PR [22]. Both are formulated based on two model parameters and therefore have the intrinsic limitation of predicting a unique and universal compressibility factor at the critical point. This results in a systematic error for the specific volume (or density) at conditions close to the critical point. To overcome this limitation, Cismondi and Mollerup [28] proposed three-parameter cubic EOS. This so-called RKPR EOS considers the effect of the actual critical compressibility factor as the third EOS parameter. The advantages of using the RKPR EOS for the prediction of the fluid properties at the transcritical and supercritical pressures has been highlighted previously [16,29,35].
The RKPR EOS is employed for all examples discussed in the present paper, and the real-fluid thermodynamic calculations are presented in form of the general cubic EOS where R is the universal gas constant, ϑ = M/ρ is the molar specific volume, a is the attractive energy parameter, and b is the co-volume parameter. Whereas δ 1 and δ 2 are two constant parameters in case of two-parameter cubic EOS, they are variables in case of three-parameter cubic EOS and are computed with one extra constraint. In the RKPR EOS, the energy parameter is , and m = (−2.4407Z c + 0.0017)ω 2 + (7.4513Z c + 1.9681)ω + (12.5040Z c − 2.7238), (11) with Z c = 1.168Z c being the tuned critical compressibility factor. The co-volume parameter of the RKPR EOS is The third parameter is a function of only Z c and δ 2 follows from the constraint In these equations, p c , T c , Z c , and ω are the critical pressure, critical temperature, critical compressibility factor, and acentric factor of the pure fluid. The conventional approach to extend this pure-fluid cubic EOS to multi-component mixtures, is based on considering the mixture as a pure hypothetical substance with the parameters estimated via the van der Waals mixing rule: where a ij is obtained through a combination rule that can include any possible binary interaction effects among the components. We use the classical combination rule with ij being the binary interaction coefficient between components i and j in the mixture. In contrast to the pseudo-critical combination rule, in which a ij is estimated by the pure-fluid formula of the energy parameter with critical quantities estimated based on those for components i and j, the classical combination rule (16) provides the possibility of applying reduction methods for the phase equilibrium calculations, which can significantly reduce the computational costs for mixtures with many components in the mixture. Using the reduction theory [30], the energy parameter of the mixture can be computed through Here, λ k and s ki are significant (non-zero) eigenvalues and their corresponding eigenvectors for a symmetric matrix with entries ij = 1 − ij . N m is the size of the significant eigenvalues vector. More often than not, binary interaction coefficients are set to zero due to lack of accurate data or just to simplify the computations. In that case, the matrix of ij becomes a diagonal matrix, and N m = 1 regardless of the number of components in the mixture.
Single-Phase Pressure: The RKPR EOS explicitly provides the thermodynamic pressure of a stable mixture whose molar composition, specific volume, and temperature are specified. The required temperature is computed via the caloric EOS, which is explained in the next subsection for non-ideal multi-component fluids.

Caloric Equation of State
A caloric EOS provides a thermodynamic relationship between a caloric property like specific internal energy, and two other properties such as temperature and specific volume for a mixture with specified molar or mass composition. Real-fluid caloric EOS can be derived by the departure function formalism using the kinetic EOS. For the general cubic EOS (9), the molar specific internal energy e = M e of a multi-component real fluid is obtained as where the first two terms account for the absolute internal energy of the mixture at the actual temperature but at the standard pressure, and the last term accounts for the internal energy change via an isothermal thermodynamic path from the standard reference pressure to the actual value. In Eq. 18, the molar specific enthalpy at standard pressure can be computed using the NASA polynomials [36] where A i,1...8 are the polynomial coefficients for each component i, which include the formation enthalpy.
Single-Phase Temperature: The temperature of a stable mixture with specified molar composition, molar specific volume, and molar specific internal energy is determined implicitly via the caloric Eq. 18. It can be computed numerically by an initial guess and Newton iterations: where L is the line search parameter in order to ensure global convergence, e is the target energy, and e * is computed via Eq. 18 using T * . Here, c * V is the molar specific heat capacity at constant volume and computed at temperature T * . The thermodynamic relation required for the evaluation of c V using the departure function formalism for a general cubic EOS is The specific molar heat capacity of the component i at the standard pressure is computed by taking the derivative of Eq. 19 with respect to the temperature, which yields

Vapor-Liquid Equilibrium
When vapor and liquid phases coexist, phase-splitting calculations are necessary in order to correctly determine the thermodynamic properties of the multi-component mixture, of which overall values of molar specific internal energy e, molar specific volume ϑ, and molar composition X i are known. The required phase-splitting or flash calculations are briefly explained in this section. The two-phase equilibrium state of a fluid with N components is described by N equations for species mass conservation and 2 + N equations for temperature, pressure, and species chemical potentials of the two phases. This set of 2N + 2 equations must be supplemented with two more constraints in order to uniquely determine the molar specific volume, temperature, and molar composition of the multi-component liquid and vapor phases. These two constraints define the type of the flash problem. As we solve the conservative form of the governing Navier-Stokes equations, that is, transport equations for the mixture energy and density, isoenergetic-isochoric phase-splitting calculations also known as UV-flash calculations must be carried out. The two additional constraints for UV-flash calculations are and Subscripts l and v refer to liquid and vapor values, and α is the vapor mole fraction, defined as the ratio of the mole of the vapor phase to the total mole of all phases. In order to evaluate the required specific internal energies of the liquid and vapor phases, we use Eq. 18 for each phase separately, and to compute α, we use total mass balance rewritten as with M as overall average molecular weight computed using overall mole fractions X i . M l and M v are average molecular weights for the liquid and vapor mixtures that are computed from the liquid mole fractions X l,i and vapor mole fractions X v,i , respectively. The molar compositions of liquid and vapor phases represent the solution of the phasesplitting calculations.
In order to solve the flash equations more efficiently, Fathi and Hickel [30] recently introduced a new method that performs UV-flash calculations very fast and robust via Newton iterations with the exact Jacobian of the equilibrium temperature used for VTflash calculations. The VT-flash here refers to isochoric-isothermal phase-splitting calculations, which this method formulates in an effectively reduced space in terms of the molar specific value of the volume function of Mikyška and Firoozabadi [37] and reduced parameters similar to those used by Nichita and Graciaa [38]. In the following, we briefly explain how this method can be used to determine the equilibrium temperature and pressure in a general cubic EOS framework. For a comprehensive review and practical implementation guidelines, interested readers are referred to the original paper [30].
Equilibrium Temperature: The temperature in the two-phase region is obtained through an iterative method similar to that used for the single-phase case, but using the specific vapor and liquid internal energies and heat capacities to determine the overall values Vapor and liquid quantities require the molar composition and specific volume of that phase, which are results of the rapid VT-flash calculations together with the vapor mole fraction. The iterative method can be terminated when |e − e * |/|e| < 10 −6 . According to Fathi and Hickel [30], the VT-flash problem can be formulated based on effective reduced parameters ) derived from Helmholtz free energy using species molar specific volume functions. The reduced parameters are determined iteratively by applying the Newton-Raphson method for N m + 2 error equations for k = 1, . . . , N m + 2 with the Jacobian matrix for k, j = 1, . . . , N m + 2. In order to compute − → H v and − → H l as well as the required partial derivatives analytically, we fist compute the K-factors with K i ≡ y i /x i being the K-factor of the component i in the mixture. Afterwards, the vapor mole fraction is initially obtained by the solving the classic Rachford-Rice equation Next, the vapor and liquid molar compositions are obtainable from the material balances and the K-factors through for i = 1, 2 . . . N . From the molar compositions for both phases, we compute the required parameters of the RKPR EOS, including a, q, b, δ 1 , δ 2 through the relations given in the previous section for a stable single-phase mixture. The specific molar volume of the vapor is then evaluated through the isochoric constraint Requiring equality of pressures between two phases, the liquid specific molar volume for the general cubic EOS is the solution of a fifth-order polynomial The polynomial coefficients ς 0..5 are listed in Ref. [30]. They are functions of the parameters that we computed in the previous steps. Finally, for k = N m + 1, and for k = N m + 2. Note that the subscripts l and v are dropped for simplicity as these equations apply to both phases. The required equations are listed in Ref. [30] including the analytical expression of the components of the Jacobian matrix. Starting the procedure requires an initial guess for the vector ). This can be conducted using K-factor values by following the steps from the solution of the Rachford-Rice until the evaluation of H k,l and H k,v . Then, the vector of reduced parameters can be computed using its definition via for k = 1, . . . , N m + 2. Typically the K-factors are available from previous time steps in computational fluid dynamics simulations. In case of blind flashes where no previous equilibrium information is available, one can use Wilson's correlation as an initial guess.
Equilibrium Pressure Neglecting the surface tension force, the pressure is equal for both phases in equilibrium. As we used this equality for the determination of specific volumes, the equilibrium pressure is the saturation pressure of vapor or liquid obtained by the flash calculations at the equilibrium temperature.

Partial Enthalpy
By definition, the partial derivative of a quantity with respect to the mole fraction of a species while keeping temperature, pressure, and mole fractions of the other species unchanged is called the partial molar value of that quantity. It can be shown that thermodynamic relationships among extensive quantities are also valid for corresponding partial values. In this section, the calculation of the partial (mass) enthalpy in the multiphase thermodynamics framework is briefly explained, as it is required for the heat flux (6) and for some discretization schemes. The (mass-basis) partial enthalpies can be computed easily from the mole-basis ones: The partial molar enthalpy for a single-phase mixture with known temperature, specific volume and molar composition can be computed via the thermodynamic relation where ϑ i ≡ (∂ϑ/∂X k )| T,P,X j =k is the partial molar volume, which can be obtained from its definition for the considered kinetic EOS (9) as where a i ≡ (∂a/∂X i )| T,P,X j =i can be computed using reduced parameters [30] with i being species index, and The partial molar internal energy e i ≡ (∂e/∂X k )| T,P,X j =k can be computed from the caloric EOS as Note that we have neglected the dependence of δ 1 and δ 2 on the molar composition in the above proposed equations following the approach used by Gmehling et al. [39] for SRK and PR EOS.
Equilibrium Partial Enthalpy: If the mixture is unstable in a single phase, then phaseequilibrium calculations are necessary in order to determine the thermodynamic properties of a stable mixture of saturated vapor and liquid phases. The overall partial enthalpy of the component i can then be estimated as

Real Finite-Rate Chemistry
Chemical reactions change the composition of the real-fluid mixture via the formation and destruction of species, which are expressed in the species mass balance equations by net mass production ratesω i . At elevated pressures, the evaluation ofω i requires some real-gas considerations concerning the species concentrations, which are different from the usually used ideal-gas definitions. This is addressed in this section for finite-rate chemistry modeling. Similar to the mass action law of the elementary reactions, the species net mass production rates of reactions with non-integer stoichiometric coefficients can also be expressed asω where N r is the total number of reactions, and ν P ir and ν R ir are the positive molar stoichiometric coefficients of species i on the right (product) and left (reactant) side of the reaction r, which might be non-integer numbers in general. The reaction rate Q r of the reaction r depends on the concentrations of the reactants through The forward and backward rate coefficients k f r and k br are are usually expressed by the Arrhenius law as an exponential function of the temperature, with A being the pre-exponential factor and E a the reaction activation energy. In addition, n jr and n jr are reaction orders with respect to species j in the forward and backward reactions. They might be listed separately for global reactions, otherwise they are considered equal to the molar stoichiometric coefficients of that species in the reactants and products.
It is important to note that the concentrations C j of the species j in Eq. 44 should be computed in a thermodynamically consistent way for the dense gaseous mixture. We propose to use the species' fugacity for this purpose, using the simple relation with f j being the fugacity of species j in the mixture. The fugacity is computed from the fugacity coefficient ϕ j = f j /(X j p), which can be obtained using the molar specific volume function with ψ j being the molar specific value of volume function for species j, see Ref. [30].
Using the definition of the fugacity coefficient and the molar specific volume function, the required concentration can therefore be computed via the following simple relation for real multi-component gases where ψ j can be computed using the effective reduced parameters with H 1 , . . . , H Nm+2 computed using the gas or supercritical mixture temperature, composition and specific volume via Eqs. 34-36. Note that in the limit of very high temperature, the fugacity coefficient tends towards unity and the ideal-gas mixture formula C j = X j p/(RT ) is recovered consistently. Using Eq. 47 simplifies the computation of the backward reaction coefficient whenever it is not explicitly provided. In such cases, the backward rate coefficient should be computed using the equilibrium constant K c = k f /k b . If the reaction rates are expressed as proposed in terms of the fugacity, the equilibrium constant can be computed similarly to the ideal-gas formula in [40] from being the net change in the number of species in the reaction, and ∆H • r and ∆S • r being the reaction enthalpy and entropy net change at the standard pressure p • = 1atm.

Large Eddy Simulation
In this work, LES is utilized for turbulence modeling. The basic idea is to focus the computational effort on the time-accurate evolution of the large scales of fluid flows that encompass almost all of the mechanical energy of the turbulent fluid motion. The evolution of the large scales is governed by a coarse-grained or low-pass filtered form of the Navier-Stokes equations. The effect of nonlinear interactions between the large resolved scales and the small truncated ones is taken into account via an appropriate modeling of the subgrid-scale (SGS) stress tensor.
In explicit LES, the SGS tensor is approximated based on the continuous filtered differential equations and model expression depending on, e.g., the filtered strain rate and possibly other quantities for which additional transport equations are needed, and the discretization of the continuous filtered transport and model equations in space and time is carried out afterwards using standard approximation theory. In implicit LES, however, the coarse-grained discrete numerical model equations are directly derived in a single approximation step that includes modeling of the effects of unresolved interactions consistent with a specifically tailored high-order numerical discretization [32,41].
The LES presented in this paper have been performed with the adaptive local deconvolution method (ALDM) of Hickel et al. [32], which is a physically consistent implicit LES method based on a nonlinear, solution adaptive finite-volume discretization scheme and spectral turbulence theory. ALDM is appropriate for LES of the full Mach number range without any user-defined model parameters. Its superior performance in comparison to common explicit LES models has been reported by Müller et al. [42] for the transcritical and supercritical injections of pure nitrogen. As discussed in [16], ALDM does not account for all additional SGS quantities that appear after low-pass filtering the non-linear equations used for the evaluation of the thermo-transport properties and chemistry sources in transcritical flows. Although novel ideas towards modelling some of these terms have been proposed, their performance for transcritical conditions at high Reynolds number remains unclear due to a lack of experimental data for high-order statistics and spatial details. ALDM has been extensively verified and validated for transcritical injector flows; for example, by Matheis et al. [27] for LES of Tani's cryogenic coaxial injector, by Müller et al. [43] for LES of Oschwald's coaxial injector, and by Matheis and Hickel [16] for LES of the non-reacting ECN Spray-A.

Numerical Implementation
Real-fluid multiphase thermodynamics and the fugacity-based finite rate chemistry model are implemented in our fluid flow solver INCA (https://www.inca-cfd.com). We employ the same discretization schemes for the transport equations as Matheis and Hickel [16]. In this subsection, these schemes are briefly reviewed.
The governing equations are discretized in space by a conservative finite-volume scheme that uses a second-order central difference method for the diffusion terms, and ALDM for the inviscid fluxes. The van Albada limiter [44] is utilized for the mass and energy flux reconstruction to avoid spurious oscillations with possible under-or overshoots at sharp density gradients.
A second-order Strang-splitting scheme is utilized to separate the stiff chemical reactions from the advection and diffusion mechanisms. The method consists of three major steps: First, the solution is updated by considering only the reaction source terms in half a time step by means of a stiff ODE solver. We used VODE [45] for this purpose, which is a variable-coefficient implicit solver based on 5th-order backward differentiation formulas. Second, the obtained solution is set as initial condition for a full advection and diffusion time step with the strong stability preserving third-order explicit Runge-Kutta scheme of Gottlieb and Shu [46]. Finally, the solution is updated with the second half reaction step. The time step size is adapted dynamically according to the Courant-Friedrichs-Lewy stability condition with CFL=1.

Experimental Setups
The standard reacting and non-reacting test cases of ECN Spray-A [47] are selected for validation and demonstration. The fuel is pure n-dodecane (C 12 H 26 ) with a low temperature of 363 K, and is injected into a chamber with a pressure of 60 bar and a temperature of 900 K. At these conditions, the cold liquid n-dodecane experiences a transcritical vaporization process, since the chamber pressure is much higher than the critical pressure (18.  initial molar compositions of the gas in the chamber are listed in Table 1 for the reacting and non-reacting cases. The chamber gas is mostly nitrogen, to which oxygen is added for the reacting case. Elevated concentrations of carbon dioxide and water vapor result from burning acetylene along with hydrogen in the preparation stage.

Grid Specifications
All simulations have been carried out in a three-dimensional rectangular domain with a length of 934 D in the axial direction and 467 D in the lateral directions, with D = 0.09 mm being the nominal diameter of the injector nozzle. To minimize the total number of computational cells, we applied adaptive mesh refinement with several predefined local refinement regions towards the injector nozzle and within a hypothetical cone with spreading angle of 10°; the resulting mesh is shown in Fig. 1. The multi-block structured grid consists of 2864 blocks and 12.7 × 10 6 cells with seven resolution levels, which are marked by L1 to L7 in Fig. 1. Around 40 % of the cells belong to the finest level with ∆y min = ∆z min 10.25 µm and ∆x min = 2∆y min located in the near-nozzle region (x/D < 60), see the zoomed area on the left side in the figure.

Boundary Conditions
Subsonic pressure outflow conditions are imposed at the exit face, using a constant static pressure of 60 bar and by extrapolating all conservative flow variables from the domain inside. The injector nozzle is not resolved, instead a transient velocity inflow is used at the injector's nozzle exit plane. This inflow is pure n-dodecane at a temperature of 363 K and a pressure of 60 bar with a transient axial velocity that provides the same amount of momentum as the experiment. We calculated the mass flow rate with the CMT virtual injection rate generator (http://www.cmt.upv.es) with input parameters according to the experimental conditions and a fuel density of 687.24 kg m −3 consistent with the RKPR EOS at the given pressure and temperature. Fig. 2

Reaction Mechanism
The finite-rate chemistry model formulated for the real-fluid effects in section 2.3 is used for the calculation of chemical source terms in the species mass conservation equations. Due to the already high computational load of the LES-MT method, a highly-reduced reaction mechanism is selected. It is a global two-step reaction mechanism calibrated using Bayesian inference for the precise prediction of the ignition delay time of n-dodecane combustion at the experimental condition of the ECN Spray-A test case [48]. The reduced mechanism has 6 species. Similar to the approach of Westbrook and Dryer for the oxidation of paraffins, the incomplete oxidation of n-dodecane is accounted for by two-step reactions as given below: The reaction rates (in cgs units) are expressed in Arrhenius form The logarithm of the Arrhenius pre-exponential factor of the first reaction varies according to the local fresh gas condition dynamically through The best possible values of θ 0 to θ 6 found by matching the ignition delay time with the Table 2: The optimal parameters for evaluation of the pre-exponential factor of ndodecane reaction in the reduced mechanism of Hakim et al. [48]. θ 0 θ 1 θ 2 θ 3 θ 4 θ 5 θ 6 27.38 -2.13 -2.05 1.89 -0.01 2.87e-4 8.43 skeletal mechanism of Narayanaswamy et al. [49] are listed in Table 2. The local fresh gas temperature T 0 and equivalence ratio φ 0 can be estimated using the local value of mixture fraction [9].   [49] and the current work using the reduced mechanism of Hakim et al. [19] for the stoichiometric mixture of the fuel and oxidizer of ECN Spray-A. Fig. 3 shows an excellent agreement between the ignition delay time predicted by the reduced two-step 6-species reaction mechanism of Hakim et al. [19] used in this study, and the skeletal mechanism of Narayanaswamy et al. [49] with 876 reactions and 164 species. The ignition delay time is reported as the elapsed time until the initial temperature passes 400 K raise in a homogeneous constant-volume reactor. The fluid is a stoichiometric mixture between the fuel and oxidizer with the composition of the reacting Spray-A case. The volume of the reactor is set according to the initial temperature, stoichiometric composition, and pressure of 60 bar using RKPR EOS. Table 3 lists the critical values and the acentric factor of the species required for the simulations of ECN Spray-A with this reduced reaction mechanism.

Numerical Results
In this section, we use experimental reference data of the reacting and non-reacting ECN Spray-A in order to evaluate and validate the proposed multiphase thermodynamics and real-fluid finite rate chemistry model for transcritical fuel sprays.

Non-reacting case
In Fig. 4, snapshots of inert Spray-A at 7 specific times after the start of injection (ASOI) predicted by the LES-MT model (right column) are compared with experimental Schlieren images [50]. The non-reacting jet first penetrates in axial direction and spreads in radial direction while evaporating and mixing with the hot chamber gas, and then maintains an approximately constant radius after some downstream distance from the injector nozzle. In the experimental Schlieren images, the saturated dark regions represent the liquid phase. For LES-MT snapshots, we show colored contours of the liquid volume fraction (LVF), which also indicate the amount of the liquid phase, and density gradients in the single-phase gaseous regions, where LVF=0. The liquid penetration length is reported to be about 10 mm in the experiment, which is in excellent agreement with our numerical results, see also Fig. 5. Figure 5 shows the temporal evolution of the liquid penetration length (LPL) and the vapor penetration length (VPL) for the LES-MT numerical simulation and experimental measurements [51,52]. For the LES-MT, LPL and VPL are defined as the maximum axial locations with a LVF of 15 % and a mixture fraction of 0.01 %, respectively, similar as in previous studies [15]. The LPL signals of LES and experiment are in excellent agreement. For the VPL, we observe excellent match up to about 0.6 ms. At later times, the simulation predicts slightly larger values that measured in the experiment. This could be do to the coarsened mesh at those locations far from the nozzle, or due to measurement uncertainties. The estimated uncertainty of the experimental measurements is indicated by the gray-shaded area and significantly increases with time. On the right side of the same figure, we also compare an experimental snapshot with highlighted liquid and vapor boundaries with a numerical visualization. The agreement is good, considering that both snapshots are instantaneous samples of independent realizations of highly turbulent flows.
We present ensemble averaged profiles of the mixture fraction on the center line and Figure 4: Time sequence snapshots of non-reacting Spray-A using Schlieren images [50] and predictions by the LES-MT method. The background contour of LES data shows the gradient of density, which is overprinted in the two-phase region by liquid volume fraction contours.
at two downstream locations in Fig. 6. The statistics have been computed by ensemble averaging LES data collected at 8 circumferential sections every 0.5 µs during the time interval highlighted in Fig. 2. LES-MT results follow the measured axial profiles very well. At the first station, x =18 mm, the n-dodecane mass fraction on the jet axis is overestimated by the LES compared to the experiment; however, the LES data fully agrees with the experimental data further downstream, which can be also seen from the radial profiles at x =36 mm. Figure 7 shows global scatter plots of the temperature and the mole fractions of the major species n-dodecane and nitrogen as a function of the mixture fraction. Each point represents an instantaneous local state of the resolved LES-MT flow field for the nonreacting case at 680 µs. The data points are colored with the molar vapor fraction in a way that green is completely vapor and purple is completely liquid. The two-phase region at the nominal operating pressure of 60 bar is indicated in the temperature/mixturefraction diagram. A hypothetical temperature profile predicted by multiphase thermodynamics assuming isenthalpic mixing at 60 bar is also shown. Two-phase boundary   and isenthalpic curve cross each other at the VLE mixture fraction of about 0.34 for this case using RKPR EOS. The LES data points follow the isenthalpic line strikingly well. This is explained in more detail in Refs. [15,16]. In the two phase region, liquid and vapor phases have different molar composition, as determined by phase-equilibrium computations. The vapor and liquid compositions are shown for the two major species, n-dodecane and nitrogen, in the mole-fraction/mixture-fraction diagrams in Fig. 7. The diagrams indicate that the saturated liquid contains only about 90 % n-dodecane and the rest is mostly nitrogen. This is a consequence of high solubility at transcritical pressures, and questions the standard pure-fuel assumption made for liquid droplet in traditional LPT methods.

Reacting case
A fist impression of the time evolution of the reacting ECN Spray-A case is provided in Fig. 9 based on experimental Schlieren images [50] on the left side and corresponding snapshots of LES-MT solution on the right side. Similar to the inert test case, the liquid phase appears as a saturated black region in the experimental Schlieren images and contours of the LVF show the predicted amount of liquid for the LES-MT. There is a very good agreement between the experiment and simulation as before. The jet develops very similar to non-reacting spray up to the auto-ignition. The ignition delay time is approximately 400 µs for both experiment and simulation. Around to the ignition time, low-temperature reactions are activated in a significant portion of the vaporized fuel. This is detectable between 20 and 25 mm after 314 µs via brighter regions in the Schlieren images, and visualized by blue iso-temperature surfaces at 920 K for the LES-MT results. The high-temperature ignition after the low-temperature reactions created enough radicals and intermediate species. The highly reacting region is highlighted by the red iso-temperature surface at 2000 K for the LES results in the last snapshot at 680 µs. The high-temperature ignition is a rapid volumetric process. The abrupt radial expansion of the reacting jet, compared to the inert case, indicates the location of quasi- Figure 8: Instantaneous snapshots for the vapor-liquid molar composition for the nonreacting case. The background contour represents the temperature field from dark to light shades. In the two-phase region, contours of the species mole fraction of saturated liquid, saturated vapor and overall mixture are presented.
steady the lift-off length (LOL), which is about 17 mm. We observe good agreement of the experimental LOL with our LES-MT result. The reacting transcritical jet can be characterised by three important lengths scales: LPL, VPL, and LOL. Figure 10 shows the time evolution of these lengths for our LES-MT results and the experimental measurements [51,52]. For the LES-MT, LPL and VPL are computed by the method explained above for the non-reacting case. The LOL is computed as the minimum axial location where the temperature exceeds 1800 K. The uncertainty of the experimental measurements is indicated by the gray area around the dotted lines in Fig. 10. The numerically predicted and experimentally measured LPL evolution is in excellent agreement. The VPL evolution is in very good agreement up to about 40 mm and afterwards the values predicted by the simulation are slightly higher than the experimental data. This is consistent with the results for the non-reacting case. The evolution of the flame LOL is again in excellent agreement with the experiment, even though we use a highly reduced reaction mechanism for the simulation. In order to visualize the flame shape and volume, an experimental snapshot with highlighted high temperature boundaries is compared with the temperature contour plot for an LES-MT snapshot at the the same time instant on the right side of Fig. 10. The numerical and experimental flame snapshots are strikingly similar. Figure 11 shows species mole fraction contours on a plane normal to the axial direction at x = 18 mm and t = 590 µs, that is, at distance and time greater than LOL and IDT. Accordingly, there exists no liquid phase at this location. The molar composition indicates partially premixed conditions of the fuel and environment even in the core of the reacting jet. The molar composition at the core is more close to the saturated vapor phase.
The auto-ignition process is illustrated in Fig. 12   of the temperature in mixture fraction space at three different time instants. The first snapshot taken at 380 µs is representative for the very late pre-ignition state, where significant low-temperature reactions are occurring and temperature is slowly increasing around the stoichiometric mixture fraction. The spatial distribution of these lowtemperature kernels is shown in Fig. 9. The second snapshot is taken at 420 µs, that is, shortly after the auto-ignition time (t ≈ 400 µs). Several high-temperature flame kernels have been formed in regions with a stoichiometric mixture fraction, where enough radicals and intermediate species coexist. In the third snapshot at 590 µs, the ignition process is completed and has provided enough thermal energy for the propagation of the flame. We see that the high-temperature reactions are spreading into the fuel-rich regime.

Discussion
The above presented results of ECN Spray-A simulations demonstrate the potential of LES-MT modelling for the accurate prediction of reacting and non-reacting turbulent multiphase fluid flows at transcritical pressures. The proposed LES-MT method is free from the main limitations of classical Lagrangian methods, which do not account for the solubility of real-fluid mixtures at elevated pressures and additionally may suffer from inaccurate density predictions if ideal-gas models are used, and Eulerian dense-gas (DG) models, which can become arbitrarily inaccurate in the two-phase region (where they predict negative pressures, e.g.). Moreover, high-pressure combustion typically initiates in low-temperature reacting zones in which the fluid's state and transport properties deviate strongly from ideal-gas laws. Making ideal-gas assumptions, which is very common in combustion simulations, or the assumption of a dense-gas without accounting for local two-phase regions, can lead to large errors and uncertainties in transcritical combustion simulations. Furthermore, the LES-MT method provides the possibility of applying different reaction mechanism in the liquid and vapor phases, whereas classical DG methods have to apply the reaction mechanism on the overall composition, including possible condensates. For these reasons, we emphasise that real-fluid and phase-equilibrium effects should be considered in transcritial jets and flames. Doing so can yield very accurate predictions for the spreading angle, liquid penetration lengths, vapor penetration lengths, ignition delay time, and flame lift-off length, such as shown in Fig. 13 and discussed in the previous sections. The major drawback of the LES-MT method might be the highly non-uniform computation costs, which limits the parallel scalability of flow solvers with domain decomposition methods. It should be also mentioned that we have applied LES in a rather straightforward way without specifically addressing turbulence-combustion interactions Non-reacting jet Reacting jet Figure 13: LES data of the transcritical injection of n-dodecane (C 12 H 26 at 363 K) into a preheated chamber at 900 K with and without oxygen. Operating pressure is 60 bar which exceeds the critical pressure of n-dodecane.
other than by providing reasonable spatial and temporal resolution. The unresolved thermodynamic microstructure has without doubt nonlinear effects on the resolved-scale solution, and it seems reasonable to assume that transcritical high-pressure flows are much stronger affected than ideal-gas flows. Hence, the quantification of uncertainties resulting from the unresolved microstructure on the resolved pressure, temperature, and subgrid-scale terms are subject of our future studies.

Conclusions
We have presented a real-fluid finite-rate chemistry multiphase thermodynamics model for the numerical simulation of transcritical vaporization and auto-ignition of a cold or cryogenic fuel injected into a hot high-pressure environment. The methodology is based on solving the fully conservative form of the compressible multi-species Navier-Stokes equations along with real-fluid kinetic and caloric state equations. These state equations provide accurate real-fluid thermodynamic properties for multi-component fluids that can exist either in a single phase or undergo phase transitions during vaporization and condensation. Multiphase and real-gas effects are also considered in the finite-rate chemistry model, which we propose to base on the fugacity of the species. The methodology has been demonstrated and validated for the transcritical reacting and non-reacting ECN Spray-A. LES results obtained with the proposed thermodynamics models are in excellent agreement with experimental reference data for the non-reacting and the reacting cases. The time evolution of temperature in the mixture fraction space proves the existence of the low-temperature and high-temperature combustion stages in the auto-ignition process of Spray-A that have been reported experimentally.
The very good agreement of auto-ignition time, flame lift-off length, and flame structure might be surprising because only a simple two-step reaction mechanism has been applied. However, one should note that Hakim's mechanism has been specifically calibrated for the considered flow conditions including low-and high-temperature ignition stages. In our LES, the real-gas vapor-liquid equilibrium accurately determines the saturated vapor composition, and this composition and the rate with which it mixes with the oxidizer determines when and where auto-ignition takes place.