We present a novel framework for high-fidelity simulations of inert and reacting sprays at transcritical conditions with highly accurate and computationally efficient models for complex real-gas effects in high-pressure environments, especially for the hybrid subcritical/supercritical mode of evaporation during the mixing of fuel and oxidizer. The high-pressure jet disintegration is modeled using a diffuse interface method with multiphase thermodynamics, which combines multi-component real-fluid volumetric and caloric state equations with vapor–liquid equilibrium calculations for the computation of thermodynamic properties of mixtures at transcritical pressures. 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 is used as a physically consistent turbulence model for large eddy simulation (LES). The proposed method represents multiphase turbulent fluid flows at transcritical pressures without relying on any semi-empirical breakup and evaporation models. All multiphase thermodynamic model equations are presented for general cubic state equations coupled with a rapid phase-equilibrium calculation method that is formulated in a reduced space based on the molar specific volume function. LES results show a very good agreement with available experimental data for the reacting and non-reacting engine combustion network benchmark spray A at transcritical operating conditions.

## I. 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 require the rapid injection of cold liquid or liquid-like supercritical fuels into a chamber environment 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. The term *transcritical* 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.

Although the operating conditions are often nominally supercritical with respect to the fuel stream and a direct transition from a liquid to a gas-like supercritical state is expected, aligned with Gibbsian thermodynamics and experimental reports,^{1–4} the mixture can locally enter the two-phase regime and interfaces between liquid-like and gas-like phases may form. Transcritical jets and sprays, therefore, resemble a 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–7} or by Eulerian single-phase dense-gas (DG) approaches with diffuse vapor–liquid interfaces.^{8–10} LES-LPT and LES-DG approaches are efficient for many flows, but their inherent model assumptions impose significant limitations at transcritical conditions. The standard LPT method is very sensitive to empirical tuning parameters and was developed for dilute mixtures, neglecting real-fluid effects. The LES-DG approach, on the other hand, excludes transcritical phase separation and may lead to nonphysical or ill-defined states when a part of the flow passes the metastable boundaries, in particular, at lower transcritical pressures.^{11} Furthermore, some important transcritical effects such as the high solubility of the saturated liquids or the different evaporation rates of the components of the fuel are not captured by these models; see also the discussion in Refs. 12–16, for example.

Using multiphase thermodynamics (MT) in the context of a diffuse-interface method has been demonstrated recently to be a promising technique to overcome the aforementioned limitations.^{14,17} In this formulation, the fully conservative Navier–Stokes equations (NSE) are solved for a hypothetical multi-component fluid mixture with thermo-transport 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 flow physics including the subcritical region of coexisting multi-component vapor and liquid as well as the non-ideal fluid behavior such as dissolution of the ambient gas in the compressed liquid.

LES-MT studies show excellent agreement with experimental data for the non-reacting transcritical sprays;^{14,18,19} however, the applicability of the method to reacting flows remained as an open question mainly due to the high computational cost of the VLE calculations, which increases with the number of components in the mixture, and due to the need for chemistry models that can capture the departure from ideal-gas (IG) behavior in regions with high pressures and low temperatures. The latter can be solved using an appropriate reduced reaction mechanism^{20} and utilizing the fugacity values of species in the mixture for the evaluations of real fluid effects on the reaction rates.^{21}

The non-ideal fluid behavior becomes apparent when using an appropriate EOS. Two-parameter cubic EOS, such as the Soave–Redlich–Kwong (SRK)^{22} and Peng–Robinson (PR)^{23} models, are very popular because of their computational efficiency.^{17,24–27} The intrinsic drawback of these two-parameter cubic EOS is the use 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.^{28} Another option is utilizing a cubic EOS with three parameters. The Redlich–Kwong–Peng–Robinson (RKPR) model of Cismondi and Mollerup^{29} 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.^{30}

To alleviate the high computational cost of VLE calculations of reacting fluids, which typically have a very large number of components, Fathi and Hickel^{31} have recently introduced a new phase-splitting method, which is formulated in a reduced space and based on the molar specific values of the volume function. In comparison with the conventional method (e.g., Refs. 14 and 32), 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 robust and quick convergence near the critical point and phase boundaries. In addition, the novel method is directly applicable to isoenergetic–isochoric 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.^{31}

In this paper, we present several novel physical and numerical models that improve the accuracy and computational efficiency of high-fidelity 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 with multiphase thermodynamics models and the adaptive local deconvolution method (ALDM) for LES turbulence modeling.^{33} Improved accuracy compared to previous LES-MT simulations is provided by the RKPR EOS^{29} for transcritical vaporization and real-fluid properties, and by fugacity-based finite rate chemistry for combustion modeling. The computational cost is significantly reduced by a rapid VLE solver,^{31} which uses reduction methods for the phase-splitting calculations, and by utilizing a highly optimized reaction mechanism,^{20} which reduces the required number of transported species. 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.

## II. PHYSICAL AND NUMERICAL MODELS

We use the LES-MT method for solving reacting multiphase compressible Navier–Stokes equations 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.

### A. Governing equations

The three-dimensional compressible reacting Navier–Stokes equations (NSE) 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 $et=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,\u2026,N$ with *N* being the total number of components comprising the mixture, *Y _{i}* is the mass fraction, $ji$ is the diffusive mass-flux vector, and $\omega \u0307i$ is the net mass production rate. $\tau $ and

**q**denote the viscous stress tensor and the vector of heat fluxes.

**I**is the unit tensor and $\u2207$ is the nabla operator.

The viscous stress tensor is modeled by assuming a Newtonian fluid and Stokes' hypothesis,

The molecular viscosity coefficient *μ* of the multi-component mixture is estimated using Chung's correlations.^{34} 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.

^{34}

By neglecting Soret and Dufour effects, the mass diffusion $ji$ 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 $Xi=YiM/Mi$ from the mass fractions and the mean molecular weight $M=1/(\u2211iNYi/Mi)$ of the mixture, with

*M*being the molar mass of the species

_{i}*i*. The product of density and the binary mass diffusion coefficient of species

*i*and

*j*, $\rho Dij$, can be approximated accurately with Chapman and Enskog theory

^{34}without high-pressure corrections if the system pressure is under 100 bar.

^{35}

### B. Multiphase thermodynamics

Simulations of multi-component fluid flows at elevated pressure require proper volumetric 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 equations required for the computation of temperature, pressure, and partial mass enthalpies of a general fluid in single phase or two phases in equilibrium.

#### 1. Volumetric equation of state

A volumetric EOS represents the thermodynamic relation 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 volumetric EOS is the ideal gas (IG); however, its applicability is limited to gaseous media with a compressibility factor close to unity, that is, it is 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^{22} and PR.^{23} Both are formulated based on two model parameters. They have the intrinsic limitation of predicting a unique and universal compressibility factor at the critical point, which therefore 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^{29} proposed a three-parameter cubic EOS. This so-called RKPR EOS considers the effect of the actual critical compressibility factor as the third EOS parameter. In Fig. 1, we compare the accuracy of various cubic EOS against the highly accurate reference EOS of Lemmon and Huber^{36} for pure n-dodecane at the ECN spray A operating pressure $60\u2009bar$. The data show that RKPR EOS significantly improves on the known deficiencies of two-parameter EOS (SRK and PR). The RKPR EOS can also yield more accurate predictions of caloric properties than the perturbed chain-statistical associating fluid theory (PC-SAFT) EOS,^{37} such as shown for the heat capacity in the temperature range $450\u2264T\u2264750\u2009K$ in Fig. 1. The advantages of using the RKPR EOS for the prediction of the fluid properties at the transcritical and supercritical pressures have also been highlighted previously.^{14,30,38}

The RKPR EOS is employed for all application examples discussed in the present paper. However, as the general cubic EOS,

includes IG, SRK, PR, RKPR, and other EOS as special cases, all thermodynamic model relations are presented in form of this general cubic EOS, in which $R$ is the universal gas constant, $\u03d1\xaf=M/\rho $ is the molar specific volume, *a* is the attractive energy parameter, and *b* is the co-volume parameter. *δ*_{1} and *δ*_{2} are two constant parameters in the case of two-parameter cubic EOS, whereas they are variables computed from one extra constraint in the case of three-parameter cubic EOS. With the RKPR EOS, the energy parameter is

where $d=(1+\delta 12)/(1+\delta 1),\u2009\xi =1+[2(1+\delta 1)]13+[4/(1+\delta 1)]13$, and

with $Zc=1.168Zc$ being the tuned critical compressibility factor. The co-volume parameter of the RKPR EOS is

The third parameter,

is an empirical function of $Zc$, and *δ*_{2} follows from the constraint,

In these equations, *p _{c}*,

*T*,

_{c}*Z*, and

_{c}*ω*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 and estimating the EOS parameters 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 $\u03f1ij$ 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. Using the reduction theory,

^{31}the energy parameter of the mixture can be computed through

Here, *λ _{k}* and

*s*are significant (non-zero) eigenvalues and their corresponding eigenvectors for a symmetric matrix with entries $\u03f1ij\u2032=1\u2212\u03f1ij$.

_{ki}*N*is the number of significant eigenvalues. 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 $\u03f1ij\u2032$ becomes a diagonal matrix, and

_{m}*N*= 1 regardless of the number of components in the mixture.

_{m}##### a. Single-phase pressure

The RKPR EOS explicitly provides the thermodynamic pressure of a thermodynamically stable mixture whose molar composition, specific volume, and temperature are specified. The required temperature is computed via the caloric EOS, which is explained in Sec. II B 2 for non-ideal multi-component fluids.

#### 2. Caloric equation of state

A caloric EOS provides a thermodynamic relation 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 volumetric EOS. For the general cubic EOS (9), the molar specific internal energy $e\xaf=Me$ of a multi-component real fluid is obtained as follows:

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.

The molar specific enthalpy at standard pressure is usually computed using polynomials,^{39}

where $Ai,1\u20268$ are the calibrated polynomial coefficients for each component *i*, which include the formation enthalpy.

##### a. 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\xaf$ is the target energy, and $e\xaf*$ is computed via Eq. (18) using $T*$. Here, $c\xafV*$ is the molar specific heat capacity at constant volume and computed at temperature $T*$. The thermodynamic relation required for the evaluation of $c\xafV$ 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 the polynomial (19) with respect to the temperature, which yields

#### 3. 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\xaf$, molar specific volume $\u03d1\xaf$, 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 the total mass balance rewritten as follows:

with *M* as overall average molecular weight computed using overall mole fractions *X _{i}*.

*M*and

_{l}*M*are the average molecular weights for the liquid and vapor mixtures that are computed from the liquid mole fractions $Xl,i$ and vapor mole fractions $Xv,i$, respectively. The molar compositions of liquid and vapor phases represent the solution of the phase-splitting calculations.

_{v}In order to solve the flash equations more efficiently, Fathi and Hickel^{31} 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 VT-flash 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^{40} and reduced parameters similar to those used by Nichita and Graciaa.^{41} 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.^{31}

##### a. 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 the results of the rapid VT-flash calculations together with the vapor mole fraction. The iterative method can be terminated when $|e\xaf\u2212e\xaf*|/|e\xaf|<10\u22126$.

According to Fathi and Hickel,^{31} the VT-flash problem can be formulated based on effective reduced parameters $H\u2192=(H1\Delta ,\u2026,HNm+2\Delta )$ derived from Helmholtz free energy using species molar specific volume functions. The reduced parameters are determined iteratively by applying the Newton–Raphson method for $Nm+2$ error equations,

for $k=1,\u2026,Nm+2$ with the Jacobian matrix,

for $k,j=1,\u2026,Nm+2$. In order to compute $H\u2192v$ and $H\u2192l$ as well as the required partial derivatives analytically, we first compute the K-factors,

with $Ki\u2261yi/xi$ being the K-factor of the component *i* in the mixture. Afterward, the vapor mole fraction is initially obtained by the solving the classic Rachford–Rice equation,

Next, the vapor and liquid molar compositions are obtained from the material balances and the K-factors through

for $i=1,2,\u2026,N$. From the molar compositions for both phases, we compute the required parameters of the RKPR EOS, including $a,q,b,\delta 1,\delta 2$ through the relations given in Sec. II B 1 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 $\u03c20..5$ are listed in Ref. 31. They are functions of the parameters that we computed in the previous steps. Finally, $Hk,l$ or $Hk,v$ is computed from

for $k\u2264Nm$,

for $k=Nm+1$ and

for $k=Nm+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. 31 including the analytical expression of the components of the Jacobian matrix.

Starting the procedure requires an initial guess for the vector $H\u2192=(H1\Delta ,\u2026,HNm+2\Delta )$. This can be calculated from K-factor values by following the steps from the solution of the Rachford–Rice until the evaluation of $Hk,l$ and $Hk,v$. Then, the vector of reduced parameters can be computed using its definition via

for $k=1,\u2026,Nm+2$. Typically, the K-factors are available from previous time steps in computational fluid dynamics simulations. In the case of blind flashes, where no previous equilibrium information is available, one can use Wilson's correlation as an initial guess.

##### b. 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.

#### 4. 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 relations between 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: $hi=h\xafi/Mi$. 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 $\u03d1\xafi\u2261(\u2202\u03d1\xaf/\u2202Xk)|T,P,Xj\u2260k$ is the partial molar volume, which can be obtained from its definition for the considered volumetric EOS (9) as follows:

where $ai\u2261(\u2202a/\u2202Xi)|T,P,Xj\u2260i$ can be computed using reduced parameters,^{31}

and

The partial molar internal energy $e\xafi\u2261(\u2202e\xaf/\u2202Xk)|T,P,Xj\u2260k$ can be computed from the caloric EOS as follows:

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.*^{42} for SRK and PR EOS.

##### a. Equilibrium partial enthalpy

If the mixture is unstable in a single phase, then phase-equilibrium 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 follows:

### C. Real finite-rate chemistry

Chemical reactions change the composition of a fluid mixture via the formation and destruction of species, which is expressed in the species mass balance equations by net mass production rates $\omega \u0307i$. At elevated pressures, finite-rate chemistry models have to account for the non-ideal thermodynamic state. Differences from the usually used ideal-gas definitions are addressed in this section.

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 follows:

where *N _{r}* is the total number of reactions, and $\nu irP$ and $\nu irR$ 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 $Qr$ of the reaction

*r*depends on the concentrations of the reactants through

The forward and backward rate coefficients *k _{fr}* and

*k*are usually expressed by the Arrhenius law,

_{br}as an exponential function of the temperature, with $A$ being the pre-exponential factor and *E _{a}* the reaction activation energy. In addition, $njr\u2032$ and $njr\u2033$ are reaction orders with respect to species

*j*in the forward and backward reactions. Reaction orders 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 $Cj$ of the species *j* in Eq. (44) should be computed in a thermodynamically consistent way for dense gaseous mixtures at high pressure. We propose to use the species' fugacity for this purpose, that is,

with *f _{j}* being the fugacity of species

*j*in the mixture. The fugacity is computed from the fugacity coefficient $\phi j=fj/(Xjp)$, 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. 31). 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 non-ideal multi-component gases:

where *ψ _{j}* can be computed using the effective reduced parameters,

with $H1,\u2026,HNm+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 toward unity and the ideal-gas mixture formula $Cj=Xjp/(RT)$ is recovered consistently.

The fugacity-based evaluation of concentrations via 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 $Kc=kf/kb$. 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 Ref. 43 from

with $\nu r=\u2211i=1N(\nu irP\u2212\nu irR)$ being the net change in the number of species in the reaction, and $\Delta Hr\xb0$ and $\Delta Sr\xb0$ being the reaction enthalpy and entropy net change at the standard pressure $p\xb0=1\u2009atm$.

### D. Large eddy simulation

In this work, LES is utilized for turbulence modeling. LES provides 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 afterward using standard approximation theory. For implicit LES, however, the coarse-grained discrete numerical model is directly derived in a single approximation step that includes models for the effects of unresolved interactions consistent with a specifically tailored high-order numerical discretization.^{44}

The LES presented in this paper have been performed with the adaptive local deconvolution method (ALDM) of Hickel *et al.*,^{33} 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 with common explicit LES models has been reported by Müller *et al.*^{45} for several transcritical and supercritical nitrogen jets. As discussed by Matheis and Hickel,^{14} 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 toward modeling 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.*^{28} for Tani's cryogenic coaxial injector, by Müller *et al.*^{46} for Oschwald's coaxial injector, and by Matheis and Hickel^{14} for the non-reacting ECN spray-A.

### E. 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.^{14} 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^{47} 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^{48} for this purpose, which is a variable-coefficient implicit solver based on fifth-order backward differentiation formulas. Second, the obtained solution is set as initial condition for a full advection and diffusion time step with the explicit third-order strong stability preserving Runge–Kutta scheme of Gottlieb and Shu.^{49} 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.

## III. CASE STUDY: TRANSCRITICAL ECN SPRAY-A

### A. Experimental setups

The standard reacting and non-reacting test cases of ECN spray-A^{50} are selected for validation and demonstration. The fuel is pure n-dodecane ($C12H26$) 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 undergoes a transcritical vaporization process, since the chamber pressure is much higher than the critical pressure (18.2 bar) of n-dodecane. The injection rail pressure is 1500 bar, and the injector consists of a hydro-eroded single-hole nozzle with the nominal diameter of $D=0.09\u2009mm$, leading to the injection of about 3.5 mg of n-dodecane fuel with a relatively constant jet velocity of about 600 ms^{−1} during the injection time of about 1.5 ms. 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. Table I summarized the initial molar composition of the gas in the chamber.

### B. Computational setup

#### 1. Grid specifications

All simulations have been carried out using a mesh generated for a three-dimensional rectangular domain with a length of $934\u2009D$ in the axial direction and $467\u2009D$ in the lateral directions, with $D=0.09\u2009mm$ being the nominal diameter of the injector nozzle.

To minimize the total number of computational cells, we instruct the flow solver INCA to perform static zonal mesh refinement within a user-define region of interest. This region of interest is defined as a 10° cone that encloses the injected jet. The target level of refinement inside the cone is defined based on the distance from the injector nozzle, with resolution steps at axial locations 665, 411, 254, 157, 97 and 60 D from the injector nozzle. The resulting mesh is shown in Fig. 2. The multi-block structured grid consists of 2864 blocks and 12.7 × 10^{6} cells with seven resolution levels, which are marked by L1–L7 in Fig. 2. Around 40% of the cells belong to the finest level with $\Delta ymin=\Delta zmin\u224310.25\u2009\mu m$ and $\Delta xmin=2\Delta ymin$ located in the near-nozzle region ($x/D<60$); see the zoomed area on the left side in the figure.

#### 2. Boundary conditions

The injector nozzle is not resolved; instead, a suitable transient inflow condition 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. Figure 3 shows the generated transient injection velocity. Subsonic inlet conditions are used as the Mach number of the pure liquid is approximately 0.5 at the nozzle exit. Mass densities and the velocity vector are imposed as Dirichlet boundary condition and a homogeneous Neumann condition is used for the static pressure. The transient profile is directly used as boundary condition without adding any artificial turbulent fluctuations at the boundary. The induced shear and hydrodynamic pressure fluctuations at 600 ms^{−1} are expected to be strong enough to create turbulence almost instantaneously.

Subsonic pressure outflow conditions are imposed at the exit face, with a constant static pressure of 60 bar as Dirichlet condition and extrapolation of the other independent flow variables from the domain inside. Adiabatic no-slip conditions are applied for all other boundaries of the rectangular domain.

#### 3. Reaction mechanism

The finite-rate chemistry model formulated for the real-fluid effects in Sec. II C is used for the calculation of chemical source terms in the species mass conservation equations. To alleviate the already high computational load of the LES-MT method, the highly reduced global two-step reaction mechanism of Hakim *et al.*^{51} is selected. This mechanism was 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.^{51} The reduced mechanism has six 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 forward 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}–*θ*_{6} found by matching the ignition delay time with the skeletal mechanism of Narayanaswamy *et al.*^{52} are listed in Table II. The local fresh gas temperature *T*_{0} and equivalence ratio $\varphi 0$ can be estimated using the local value of mixture fraction.^{9} The mixture fraction is computed here by tracking the inert nitrogen mass fraction,

with $YN2ox=0.72$ in the oxidizer stream of reacting case.

θ_{0}
. | θ_{1}
. | θ_{2}
. | θ_{3}
. | θ_{4}
. | θ_{5}
. | θ_{6}
. |
---|---|---|---|---|---|---|

27.38 | −2.13 | −2.05 | 1.89 | −0.01 | 2.87 × 10^{−4} | 8.43 |

θ_{0}
. | θ_{1}
. | θ_{2}
. | θ_{3}
. | θ_{4}
. | θ_{5}
. | θ_{6}
. |
---|---|---|---|---|---|---|

27.38 | −2.13 | −2.05 | 1.89 | −0.01 | 2.87 × 10^{−4} | 8.43 |

We evaluate the performance of the reduced two-step six-species reaction mechanism of Hakim *et al.*^{20} by calculating the ignition delay time in a homogeneous constant-volume reactor for a stoichiometric mixture of the fuel and oxidizer 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 the RKPR EOS. Table III lists the critical values and the acentric factor of the species required for the simulations of ECN spray-A with this reduced reaction mechanism. Figure 4 shows the time elapsed until the temperature of the reactor has increased by 400 K above its initial value. We observe an excellent agreement between the ignition delay time predicted by the reduced two-step mechanism used in this study and the skeletal mechanism of Narayanaswamy *et al.*^{52} with 876 reactions and 164 species.

Species . | $Tc\u2009(K)$ . | $pc\u2009(bar)$ . | $Zc\u2009(\u2212)$ . | $\omega \u2009(\u2212)$ . |
---|---|---|---|---|

$C12H26$ | 658.0 | 18.20 | 0.251 | 0.576 |

$O2$ | 154.6 | 50.43 | 0.288 | 0.022 |

$N2$ | 126.2 | 34.00 | 0.289 | 0.038 |

$CO2$ | 304.2 | 73.83 | 0.274 | 0.224 |

$H2O$ | 647.1 | 220.6 | 0.229 | 0.345 |

$CO$ | 132.9 | 34.99 | 0.299 | 0.048 |

Species . | $Tc\u2009(K)$ . | $pc\u2009(bar)$ . | $Zc\u2009(\u2212)$ . | $\omega \u2009(\u2212)$ . |
---|---|---|---|---|

$C12H26$ | 658.0 | 18.20 | 0.251 | 0.576 |

$O2$ | 154.6 | 50.43 | 0.288 | 0.022 |

$N2$ | 126.2 | 34.00 | 0.289 | 0.038 |

$CO2$ | 304.2 | 73.83 | 0.274 | 0.224 |

$H2O$ | 647.1 | 220.6 | 0.229 | 0.345 |

$CO$ | 132.9 | 34.99 | 0.299 | 0.048 |

## IV. NUMERICAL RESULTS

In this section, we use the 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.

### A. Non-reacting case

We compare the experimental and numerical flow visualizations of the inert spray-A at seven specific instants after the start of injection (ASOI) in Fig. 5. The left column shows experimental schlieren images,^{53} and the right column shows predictions of our LES-MT model. 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 (LPL) is reported to be about 10 mm in the experiment, which is in excellent agreement with our numerical results (see also Fig. 6).

Figure 6 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.^{54,55} 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 to in previous studies.^{17} 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 than those measured in the experiment. This could be due 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 the ensemble averaged profiles of the mixture fraction on the centerline and at two downstream locations in Fig. 7. The statistics have been computed by ensemble averaging LES data collected at eight circumferential sections every 0.5 *μ*s during the time interval highlighted in Fig. 3. 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 8 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 non-reacting 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/mixture-fraction 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 a 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. 14 and 17. 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. 8. The diagrams indicate that the saturated liquid contains only about 90% n-dodecane and the rest is mostly nitrogen. This is a consequence of the high solubility of nitrogen in n-dodecane at transcritical pressures and questions the standard pure-fuel assumption made for liquid droplets in traditional LPT methods.

Figure 9 shows contours of the molar composition of the overall mixture as well as the compositions of the liquid and vapor phases for all species. The figure indicates the spatial evolution of mixing process controlled by multiphase thermodynamics. While the heavy n-dodecane molecules represent about 90% of the total moles of the liquid core, they are in minority in the vapor phase, which consists mainly of light nitrogen molecules corresponding to the initial composition of the chamber gas. The fuel-rich liquid core extends up to about 3 mm from the injector nozzle. Far from the nozzle and close to the tip of the liquid core, where the vapor mole fraction is close to unity, the mixture composition is close to the saturated vapor phase and, accordingly, the mixture fraction is close to 0.34. Water has roughly the same concentration in the liquid and vapor phases, which is in contrast to nitrogen, n-dodecane, and carbon dioxide. This means that water in contact with the cold n-dodecane dissolves quickly into the liquid phase, while the other species in the environment mostly stay in the saturated vapor phase.

### B. Reacting case

The temporal evolution of the reacting ECN spray-A is shown in Fig. 10. This figure shows experimental schlieren images^{53} on the left side and corresponding snapshots of the 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, and the jet develops very similar to inert spray up to the auto-ignition. The ignition delay time is approximately 400 *μ*s for both experiment and simulation. Around 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 occurs after the low-temperature reactions created a sufficient concentration of 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. An abrupt radial expansion of the reacting jet, compared to the inert case, is observed near the location of the flame front. We observe good agreement of the experimental flame liftoff length (LOL) with our LES-MT result which is about 15 mm as shown in the last snapshot.

The reacting transcritical jet can be characterized by three important lengths scales: LPL, VPL, and LOL. Figure 11 shows the time evolution of these lengths for our LES-MT results and the experimental measurements.^{54,55} 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. 11. 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 afterward 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 same time instant on the right side of Fig. 11. The numerical and experimental flame snapshots are in good qualitative agreement.

Figure 12 shows species mole fraction contours on a plane normal to the axial direction at $x=18\u2009mm$ and $t=590\u2009\mu s$, that is, at a distance larger than LOL of the developed flame. 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. 13, where we show global scatter plots 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 the temperature is slowly increasing around the stoichiometric mixture fraction. The spatial distribution of these low-temperature kernels is shown in Fig. 10. The second snapshot is taken at 420 *μ*s, that is, shortly after the auto-ignition time ($t\u2248400\u2009\mu 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.

## V. DISCUSSION

The above-presented results of ECN spray-A simulations demonstrate the potential of LES-MT modeling 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 arbitrary inaccurate in the two-phase region (where they may predict negative pressures, e.g.). Moreover, high-pressure combustion typically initiates in low-temperature reacting zones in which the fluid 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 mechanisms in the liquid and vapor phases. Classical DG methods cannot make this distinction and have to apply the reaction mechanism on the overall composition, including possible condensates. For these reasons, we emphasize that real-fluid and phase-equilibrium effects should be considered in transcritical 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 liftoff length, such as shown in Fig. 14 and discussed in Secs. I–IV.

The LES-MT method is accurate and calibration-free. To put its performance into perspective, we compare the results for the vapor and liquid penetration trajectories of spray-A computed with the present LES-MT method, LES-DG results of Hakim *et al.*,^{9} results of the highly calibrated LES-LPT method of Gadalla *et al.*,^{7} and the reference experiments^{54,55} in Fig. 15. The figure shows the initial 0.2 ms of the transient startup during which we expect a similar behavior for the reacting and inert cases. Predictions of the current method are very close to the measurements for both liquid and vapor length. Gadalla *et al.*^{7} have calibrated the initial droplet size for their LES-LPT method in such a way that the liquid length is close to the experimental data. This calibration, however, apparently affected the vapor length prediction unfavorably. The single-phase DG method yields the least accurate prediction for this case. As the liquid phase boundaries are not directly predicted by the single-phase method, Hakim *et al.*^{9} defined the boundary of the liquid regions via the mixture fraction value at which the density gradient peaks.

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. The wall-clock time for simulating 10 *μ*s of the reacting spray, including time needed for writing snapshots every 0.5 *μ*s, is approximately 8.5 h on 16 nodes with two Intel Xeon Gold 6248R CPU (24 cores at 3.00 GHz), which corresponds to approximately 6500 CPU h. On average, about 40% of the CPU time is used for evaluating the EOS and transport properties, and the majority of the rest is used for evaluating the reaction rates, the spatial derivatives of the Navier–Stokes equations, communication, and snapshot output.

Last but not least, it should be also mentioned that we have applied LES in a rather straightforward way without specifically addressing turbulence–combustion interactions other than by providing reasonable spatial and temporal resolution. The unresolved thermodynamic microstructure (SGS variations of pressure, temperature, and composition) has without doubt nonlinear effects on the resolved-scale solution,^{57} and it seems reasonable to assume that transcritical multiphase flows are much stronger affected than ideal-gas or single-phase real-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.

## VI. CONCLUSIONS

We have presented a real-fluid finite-rate chemistry multiphase thermodynamics model for the numerical simulation of transcritical vaporization, auto-ignition, and combustion 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 volumetric 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 be based on the fugacity or the molar specific volume function 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 both 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 has been reported experimentally.

The very good agreement of auto-ignition time, flame liftoff 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, real-gas vapor–liquid equilibrium calculations accurately determine 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.

## ACKNOWLEDGMENTS

This paper is submitted as part of the themed issue “Development and Validation of Models for Turbulent Reacting Flows” in honor of Professor M. Pfitzner on the occasion of his retirement. The second author would like to thank Michael for motivating him to start working on this topic and for the very productive collaboration over many years.

This work is funded by The Netherlands Organization for Scientific Research (NWO) under Contract No. 680.91.082. Simulations have been performed on DelftBlue^{58} at the Delft High Performance Computing Center (DHPC).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Mohamad Fathi:** Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (lead); Writing – review and editing (supporting). **Stefan Hickel:** Conceptualization (lead); Funding acquisition (lead); Methodology (equal); Project administration (lead); Resources (lead); Software (equal); Supervision (lead); Writing – review and editing (lead). **Dirk Roekaerts:** Supervision (supporting); Writing – review and editing (supporting).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX: SENSITIVITY TO GRID RESOLUTIONS

To evaluate the grid sensitivity of the presented numerical results, a grid convergence study is conducted by refining and coarsening the base grid. The finest grid has the minimum spacing of 6.84 *μ*m, while the coarsest has 20.50 *μ*m. By applying the same local grid refinement procedure as used for the base grid (with 12.7 × 10^{6} cells), we obtain a total cell number of about 1.2 × 10^{6} for the coarse grid and 36.1 × 10^{6} for finest grid.

Figure 16 visualizes the temporal evolution of the jet on the three grids. The results obtained on the medium and fine grid are almost indistinguishable, whereas a too rapid penetration and a different jet breakup topology are observed for the coarsest grid. Such a rapid nonphysical penetration was also reported by authors for too coarse meshes.^{14} The main reason for this error is an insufficient number of cells in the radial direction resulting in extra momentum in the axial direction. The corresponding liquid and vapor penetration trajectories are shown in Fig. 17. Only negligible differences can be observed between the results obtained on the medium grid and on the fine grid.