Songzhi Yang

Ping Yi

Chaouki Habchi
email: chaouki.habchi@ifp.fr

Real-fluid injection modelling and LES simulation of the ECN Spray A injector using a fully compressible two-phase flow approach

Keywords: two-phase flow, real fluid, Vapor-liquid equilibrium, Spray A, in-nozzle flow

A fully compressible two-phase flow model consisting of four balance equations including two mass, one momentum, and one internal energy equation, formulated with the mechanical and thermal equilibrium assumptions is developed in this article. This model is closed with a real fluid equation of state (EoS) and has been applied to the simulation of different 1D academic cases, in addition to the 3D Large-Eddy Simulation (LES) of the Engine Combustion Network (ECN) Spray A injector including the needle to target part with and without the phase change (i.e. frozen) assumptions. The obtained numerical results from the model with phase change have proven to be able to accurately predict the liquid, vapor penetrations and rate of injection compared to experimental data. However, the frozen model has presented some uncertainties and deviations in predicting the penetration length as with different measure criterions, even though an excellent agreement can be achieved in the estimation of rate of injection, nearnozzle mass and velocity distribution. Several conclusions are drawn from the simulations: (1) the initial in-nozzle flow has a strong effect on the early jet development; (2) considering phase change is still essential in the high temperature, high pressure (HTHP) injection modelling since it strongly affects the temperature distribution, turbulence intensity and thereby the jet development; (3) significant variations of liquid compressibility factor and density, as well as the cooling effect through the nozzle are highlighted. Overall, the detailed analysis of the numerical results reported in this article may complement the Engine Combustion Network (ECN) experimental database.

Introduction

The widespread applications of fuel injection at high temperature, high pressure (HTHP) conditions in compression ignition engines, gas turbines and rocket engines have stimulated great interest in studying the liquid injection experimentally and numerically. Researches related to real-fluid injection were previously concentrated on the liquid rocket engine field, and then extended recently to the diesel engine industry. As confirmed in abundant studies of liquid rocket fuel injections [START_REF] Chehroudi | Recent experimental efforts on high-pressure supercritical injection for liquid rockets and their implications[END_REF]- [START_REF] Hickey | Large Eddy Simulation of Supercritical Mixing and Combustion for Rocket Applications[END_REF], the spray has been through an evident transition from two-phase atomization, breakup and droplets evaporation dominated physical processes to continuous diffusion and mixing phenomena as the pressure increases from subcritical condition towards supercritical condition. One primary reason for such transition is the gradual diminishing surface tension and latent heat as the ambient condition is above the critical point of the injected fuel. In fact, similar transition phenomena also occur during the injection in diesel engines, as investigated by several researchers [START_REF] Crua | On the transcritical mixing of fuels at diesel engine conditions[END_REF], [START_REF] Manin | Microscopic investigation of the atomization and mixing processes of diesel sprays injected into high pressure and temperature environments[END_REF]. For example, Crua et al. [START_REF] Crua | On the transcritical mixing of fuels at diesel engine conditions[END_REF] proposed a criterion for the mixing transition based on their recent experimental observations. They observed the droplets undergo a gradual transition from subcritical evaporation to mixing regime at different pressure and temperature above the pure fuel's critical point. Thereby, they deduced that the fuel still stays in the subcritical twophase state for some time before fully entering the diffusion mixing regime and the transition time varies with fuel types and droplets size. Thus, it is reasonable to conclude that in the whole injection process, the combination of classical evaporation regime for the main liquid core and transition to the dense gas mixing state is possible at high ambient temperature especially for the droplets formed by possible primary atomization near the nozzle or at the end of the injection events, the same as the droplets tracked experimentally by Crua et al. [START_REF] Crua | On the transcritical mixing of fuels at diesel engine conditions[END_REF]. This conclusion indirectly justifies the development of the proposed real-fluid diffuse interface model (DIM) model (see next Section), which aims to model the subcritical regime with the consideration of phase change, supercritical regime, as well as the transition from one to the other simultaneously. All the above experimental observations [START_REF] Crua | On the transcritical mixing of fuels at diesel engine conditions[END_REF], [START_REF] Manin | Microscopic investigation of the atomization and mixing processes of diesel sprays injected into high pressure and temperature environments[END_REF] have provided valuable references and guidance for modelling. One classical benchmark case corresponding to the mixing transition regime is the Engine Combustion Network (ECN) Spray A case operated under the high pressure and temperature evaporating conditions [START_REF]ECN-Engine Combustion Network (Spray A&B)[END_REF]. Many researchers have contributed to the modelling of the spray issuing from this typical Diesel injector. Generally, the involved models are based on the Eulerian-Lagrangian (EL) approach [START_REF] Xue | Large eddy simulation of fuel-spray under nonreacting IC engine conditions[END_REF], [START_REF] Wehrfritz | Large eddy simulation of high-velocity fuel sprays: studying mesh resolution and breakup model effects for spray A[END_REF] and the liquid phase is treated as dispersed particles with various diameters, smaller than the grid size, while the gas phase is regarded as the continuous Eulerian carrier fluid. As known that fuel injection using the EL approach is highly dependent on some tunable coefficients, and cannot capture the spray physics accurately, especially for the near-nozzle region, an Eulerian-Eulerian (EE) approach coupled with an Eulerian-Lagrangian (EL) Spray Atomization (ELSA) method has been proposed [START_REF] Desantes | A comparison of diesel sprays CFD modeling approaches: DDM versus Σ-Y Eulerian atomization model[END_REF], [START_REF] Desantes | Simulation of the First Millimeters of the Diesel Spray by an Eulerian Spray Atomization Model Applied on ECN Spray A Injector[END_REF]. In such model, autonomous transition from EE incompressible liquid jet to EL spray is implemented using an additional surface density transport equation, as the jet develops from the near nozzle dense zone to the downstream dilute zone. Recent work from Xue et al. [START_REF] Xue | Eulerian CFD Modeling of Coupled Nozzle Flow and Spray with Validation Against X-Ray Radiography Data[END_REF] has confirmed that the turbulent-mixing based EE model can predict better physics in the near nozzle region than the EL method for the simulation of ECN Spray A injector. However, the absence of specific modelling for the primary atomization in their model has limited the accuracy of the downstream dispersed spray results. Based on this, Devassy et al. [START_REF] Devassy | Atomization Modelling of Liquid Jets Using a Two-Surface-Density Approach[END_REF] have developed an Eulerian-Eulerian atomization model similar to the ELSA approach using a fully compressible model including a two surfaces density balance equations (TwoSD) model for the liquid core atomization and the droplets secondary breakup. This model has shown the huge effect of in-nozzle cavitation on primary atomization [START_REF] Habchi | A comprehensive Two-Fluid Model for Cavitation and Primary Atomization Modelling of liquid jets -Application to a large marine Diesel injector[END_REF], thereby demonstrating the necessity of considering physical primary atomization models directly linked to internal nozzle flows for more accurate injection simulations. In these conditions, the two successive Eulerian and then Lagrangian simulations coupled through collecting the flow information at the orifices exit as boundary conditions for the downstream Lagrangian spray simulation cannot provide better accuracy for the primary atomization of the liquid jet, especially for multiple short injections usually performed in advanced calibration of engines. In contrast, the fully compressible Eulerian method based on the continuum fluid theory can easily realize the accurate modelling from the in-nozzle flow to the dilute downstream spray continuously. While the computational cost is much higher than the Lagrangian strategy, it is still tractable if applied with varying mesh resolution based on different regions or using automatic mesh refinement technique.

Currently, the compressible Eulerian based model has been successfully applied to the ECN Spray A modelling mainly focusing on the downstream spray as demonstrated in the earlier work of Lacaze et al. [START_REF] Lacaze | Analysis of high-pressure Diesel fuel injection processes using LES with real-fluid thermodynamics and transport[END_REF], as well as in the recent work of Matheis and Hickel [START_REF] Matheis | Multi-component vapor-liquid equilibrium model for LES of highpressure fuel injection and application to ECN Spray A[END_REF], and Ma and Ihme [START_REF] Ma | Modeling and Simulation of Diesel Injection at Transcritical Conditions[END_REF].

Since the physical properties of the fuel at HTHP condition show an obvious deviation from the ideal gas state, the real-fluid EoS especially the Peng Robinson (PR) equation of state (EoS) is widely employed for the modelling of transcritical flows and the transition from subcritical regime to supercritical mixing regime, because of its ideal compromise between computational efficiency and accuracy. The involved thermodynamics models vary according to whether considering the phase change or not. Indeed, in the case of the flow with phase change, the vapor-liquid equilibrium (VLE) theory is considered. Otherwise, if no evaporation and condensation are expected in the whole computation, a direct solving of PR EoS for the multicomponent single-phase mixture can readily provide the solution. However, the latter singlephase simulation strategy may be risky if the flow thermodynamics state enters the vapor dome. As a consequence, the simulation may crash due to the appearance of negative pressures, as reported in [START_REF] Ma | On the numerical behavior of diffuse-interface methods for transcritical real-fluids simulations[END_REF]. According to Castier [START_REF] Castier | Solution of the isochoric-isoenergetic flash problem by direct entropy maximization[END_REF], negative pressure usually indicates that the fluid is in two-phase state. Since the cost of solving VLE equations is high especially for the full injector simulation, a compromised solution proposed by Knudsen et al. [START_REF] Knudsen | Compressible Eulerian needle-totarget large eddy simulations of a diesel fuel injector[END_REF] is to find an approximate saturation pressure corresponding to the point of (( ∂P ∂ρ ) 𝑇 = 0) instead of solving the real saturation state. Even if this method is not theoretically valid and the phase change process is neglected, fairly good numerical results such as jet penetration length, mass and momentum flux can still be achieved in the simulation of ECN Spray A injector [START_REF] Knudsen | Compressible Eulerian needle-totarget large eddy simulations of a diesel fuel injector[END_REF]. The undeniable fact based on the experimental results is that two-phase subcritical regime indicators like droplets observation with relatively significant surface tension indeed exists when injecting fuel like n-dodecane or hexadecane at HTHP conditions, as discussed by Crua et al. [START_REF] Crua | On the transcritical mixing of fuels at diesel engine conditions[END_REF]. However, the numerical results of Knudsen et al. [START_REF] Knudsen | Compressible Eulerian needle-totarget large eddy simulations of a diesel fuel injector[END_REF] discussed above have stimulated us to wonder whether the evaporation process really play an important role in HTHP injection modelling. This is among the topics the current paper would like to investigate and discuss.

All the previous work experimentally or numerically has enlightened us to explore more about the transcritical modelling in HTHP diesel engine. Since most of the previous studies about Spray A modelling are concentrated on the downstream spray, the integral continuous simulation from the innozzle flow to the downstream spray is very rare especially for the injection at HTHP evaporating condition. Therefore, to obtain more understanding about the effect of the in-nozzle flow on the ensuing spray development, as well as the effect of phase change on the jet evolving, an entire Spray A injector containing the needle to target part is simulated for the first time through the fully compressible Eulerian approach with phase change (VLE) [START_REF] Yi | A multicomponent real-fluid fully compressible fourequation model for two-phase flow with phase change[END_REF], but also without phase change assumption. The flow solver employed in the current work has been implemented in the in-house code IFP-C3D [START_REF] Bohbot | IFP-C3D: an Unstructured Parallel Solver for Reactive Compressible Gas Flow with Spray[END_REF], in which a fully compressible non-equilibrium two-phase flow 7-Equation model is resolved as presented in previous work of Habchi and Devassy [START_REF] Devassy | Atomization Modelling of Liquid Jets Using a Two-Surface-Density Approach[END_REF], [START_REF] Habchi | a Gibbs Energy Relaxation (Germ) Model for Cavitation Simulation[END_REF]. Since the involved EoS in the original system is employing the Stiffed Gas (SG) EoS for liquid phase and ideal gas EoS for gas phase, it is not applicable for the HTHP injection modelling. Therefore, a real fluid phase equilibrium solver based on PR EoS has been developed and implemented into this 7-Equation two-phase flow solver, and has been applied to analyze a series of 1D academic test cases and 2D cavitation modelling [START_REF] Yang | Cavitation Modelling Using Real-Fluid Equation of State[END_REF]. Following this, a further reduced 4-Equation model with the assumption of mechanical and thermal equilibrium has been proposed and combined with the real fluid phase equilibrium solver [START_REF] Yi | A multicomponent real-fluid fully compressible fourequation model for two-phase flow with phase change[END_REF]. This 4-Equation model has been successfully applied to 2D supercritical and transcritical injection modelling in [START_REF] Yang | Towards a multicomponent real-fluid fully compressible two-phase flow model[END_REF], and achieved a good agreement with experimental results in the 1D flash boiling conditions, as reported in [START_REF] Yi | A multicomponent real-fluid fully compressible fourequation model for two-phase flow with phase change[END_REF]. Based on these previous studies, the real fluid DIM solver will be further utilized to explore its potential in solving multi-scale injection problems at HTHP conditions. This article is organized as follows. First, a concise description of the mathematical model and thermodynamics theory is presented in the following Section 2. In Section 3, two 1D academic test cases including the classical advection and shock tube are simulated at conditions relevant to the ECN Spray A condition. Next, the numerical results of the Spray A 3D simulations assuming equilibrium thermodynamic as well as the frozen non-equilibrium thermodynamics are reported and discussed along with a detailed thermodynamics analysis. Finally, the conclusions and future work are presented.

Physical and Numerical models 2.1 Governing equations

The governing equations adopted in current study is the 4-Equation model as described in detail in our previous work [START_REF] Yi | A multicomponent real-fluid fully compressible fourequation model for two-phase flow with phase change[END_REF]. For completeness, we recall here its main features. It is a fully compressible multicomponent two-phase flow model obtained from the classical two-phase flow non-equilibrium 7-Equation model [START_REF] Saurel | A Multiphase Godunov Method for Compressible Multifluid and Multiphase Flows[END_REF], [START_REF] Baer | A two-phase mixture theory for the deflagration to detonation (DDT) transition in reactive granular materials[END_REF] with the assumption of mechanical and thermal equilibrium. As formulated in the following Equations ((1.1)-(1.4)), the 4-Equation model includes the balance equations for different species in gas and liquid phases, mixture momentum and mixture specific internal energy, respectively. The legacy from the 7-Equation model is maintained by transporting the specific density (𝛼 𝑝 𝜌 𝑝,𝑘 ) of each species (𝑘) in each phase (p) separately, as follows, In these equations, index 𝑙 and 𝑔 denote the liquid and gas phase respectively; 𝑘 represents different species; the right hand side (RHS) terms 𝑚̇𝑙 ,𝑘 𝑎𝑛𝑑 𝑚̇𝑔 ,𝑘 are the phase transition mass of each species in the liquid and vapor phases, respectively, restricted by 𝑚̇𝑙 ,𝑘 + 𝑚̇𝑔 ,𝑘 = 0; 𝜏 𝑖𝑗 𝐿,𝑇 is the shear stress tensor covering the laminar (L) and turbulent (T) contributions, which can be written as 𝜏 𝑖𝑗 𝐿,𝑇 = 𝜏 𝑖𝑗 𝐿 + 𝐾 0 𝜏 𝑖𝑗 𝑇 with 𝐾 0 = 1 for turbulent flow. As described in our previous work [START_REF] Habchi | a Gibbs Energy Relaxation (Germ) Model for Cavitation Simulation[END_REF], a standard Boussinesq approximation is used for the modeling of 𝜏 𝑖𝑗 𝐿,𝑇 in which a subgrid-scale turbulent viscosity given by the Smagorinsky LES model is used for the turbulent viscosity. However, the laminar viscosity is computed from Chung's correlation [START_REF] Horng | Generalized Multiparameter Correlation for Nonpolar and Polar Fluid Transport Properties[END_REF]. In Eq. (1.4), 𝑒 represents the specific internal energy; 𝑞 𝑖 𝐿,𝑇 is the heat conduction flux and usually modelled as 𝑞 𝑖 𝐿,𝑇 = -𝜆 𝜕𝑇 𝜕𝑥 𝑖 based on Fourier's law. The heat conduction coefficient 𝜆 contains the laminar and turbulent contributions. The laminar contribution is computed from Chung's correlation [START_REF] Horng | Generalized Multiparameter Correlation for Nonpolar and Polar Fluid Transport Properties[END_REF] and the turbulent one is estimated by a constant Prandtl number, 𝑃𝑟 𝑡 = 0.9. α 𝑝 denotes the volume fraction of phase p. It is computed in the phase equilibrium solver (see next Section) along with 𝑚̇𝑝 ,𝑘 . It is important to note that the species diffusion inside the liquid, Eq. (1.1) and inside the gas, Eq. (1.2) have not been considered in order to avoid additional diffusion of the liquid-gas interface in the DIM approach. However, the effects of such assumption on the results accuracy should be investigated in future work as the primary atomization could be considered using a surface density approach, for instance.

Compared with previous researches [START_REF] Chiapolino | A simple and fast phase transition relaxation solver for compressible multicomponent two-phase flows[END_REF], [START_REF] Chiapolino | A simple phase transition relaxation solver for liquidvapor flows[END_REF], the novelty of the current research lies in the application of the real fluid phase equilibrium theory to the two-phase system. In addition, each phase is considered multi-components. The latest research from Chiapolino [START_REF] Chiapolino | A simple phase transition relaxation solver for liquidvapor flows[END_REF] has applied the VLE model to single component system. Nevertheless, the current 4-Equation model with phase equilibrium model has proved to share high similarity with the very recent research from Matheis and Hickel [START_REF] Matheis | Multi-component vapor-liquid equilibrium model for LES of highpressure fuel injection and application to ECN Spray A[END_REF]. The first difference resides in their fully conservative formulation using the total energy for the energy transport equation; while a non-conservative internal energy equation is used in current work, as expressed in Equation (1.4).

The second difference lies in the mass transport equations. As noted above, the mass conservation is carried out by considering each phase and component separately instead of a classical homogeneous equilibrium mixture, in which species diffusion between liquid and gas is permitted in addition to phase change [START_REF] Matheis | Multi-component vapor-liquid equilibrium model for LES of highpressure fuel injection and application to ECN Spray A[END_REF].

Thermodynamics model

Equation of state

The deficiency of ideal gas EoS and Stiffened gas EoS for liquid in describing high pressure and high temperature (HTHP) properties has drawn researchers to real fluid EoS. Among them, the cubic equation series have been widely employed for its good compromise between computational efficiency and accuracy [START_REF] Matheis | Multi-component vapor-liquid equilibrium model for LES of highpressure fuel injection and application to ECN Spray A[END_REF], [START_REF] Ma | On the numerical behavior of diffuse-interface methods for transcritical real-fluids simulations[END_REF], [START_REF] Qiu | Simulation of supercritical fuel injection with condensation[END_REF]. In this work, Peng-Robinson EoS is selected as the thermal closure equation, formulated as Eq.(1.5). Since the interest is to solve the multicomponent system, van der Waals mixing rule (Eq. (1.6)) is adopted for the computation of mixture thermal properties. 𝑥 𝑖 is the molar fraction for each species 𝑖, 𝑗. 𝑘 𝑖𝑗 is the binary interaction parameter between different components which is usually fitted based on the experimental data.

𝑃 = 𝑅𝑇 𝑣 -b - a(𝑇) 𝑣(𝑣 + b) + b(𝑣 -b) (1.

Description of thermodynamics solver

To close the 4-Equation system (Eq.(1.1)-(1.4)), the equation of state is essential to construct the relation between specific internal energy (𝑒 * ), density (𝜌 * ), species composition 𝑌 𝑖 * and temperature (𝑇) and pressure (𝑃 ). All the variables marked with (*) exponent denote the output value of transport equations. After solving the transport equation, the new (𝑃, 𝑇) need to be computed from the (𝑒 * , 𝜌 * , 𝑌 𝑖 * ) to start a new cycle. The transformation of (𝑒 * , 𝜌 * ) to molar quantities (𝑢 * , 𝑣 * ) is through molar weight (𝑀 𝑤 ) (Eq.(1.7)). Unlike the explicit relation between internal energy and temperature as in SG EoS [START_REF] Saurel | A general formulation for cavitating, boiling and evaporating flows[END_REF] or ideal gas EoS, an iterative root searching algorithm is necessary for cubic EoS to find the right 𝑃, 𝑇 from 𝑢, 𝑣. Hence the final objective functions are formulated as Eq. (1.8). 𝑢 and 𝑣 are the mixture amount of gas phase and liquid phase, formulated as Eq.(1.9). The internal energy (𝑢 𝑝 ) in each phase is divided into the ideal gas part (𝑢 0 ) and the departure part (𝑢 𝑑,𝑝 ) (Eq. (1.10)). The ideal gas internal energy (𝑢 0 ) is computed with empirical coefficient equation referred to [START_REF] Aly | Self-consistent equations for calculating the ideal gas heat capacity, enthalpy, and entropy[END_REF]. The departure function is derived according to the PR EoS, formulated as Eq. (1.11). The coefficients 𝑎, 𝑏 in Eq. (1.11) are calculated with the van der Walls mixing rule (Eq. (1.6)) which differs significantly with the widely applied ideal gas mixing rule [START_REF] Chiapolino | A simple and fast phase transition relaxation solver for compressible multicomponent two-phase flows[END_REF], [START_REF] Wang | Simulating cavitating liquid jets using a compressible and equilibrium two-phase flow solver[END_REF]. In addition, the molecular mutual effect between different components is controlled by binary interaction parameter (𝑘 𝑖𝑗 or BIP).The density (𝜌 𝑝 ) or specific volume (𝑣 𝑝 ) in each phase is computed directly by solving the cubic equation [START_REF] Wilczek-Vera | Understanding cubic equations of state: A search for the hidden clues of their success[END_REF]. Phase composition ( 𝑥 𝑘 (𝑙𝑖𝑞𝑢𝑖𝑑 𝑝ℎ𝑎𝑠𝑒), 𝑦 𝑘 (𝑣𝑎𝑝𝑜𝑟 𝑝ℎ𝑎𝑠𝑒) : molar fraction of species) is needed to compute the thermal properties of each phase (𝑢 𝑝 , 𝑣 𝑝 ). Generally, the calculation of phase composition is through the phase equilibrium computation [START_REF] Yi | A multicomponent real-fluid fully compressible fourequation model for two-phase flow with phase change[END_REF]. However, in current model, with the transporting of each individual specific density (𝛼 𝑝

𝑢 * =
1) The first method is to assume the local flow field reaching the thermodynamics equilibrium instantaneously which implies that the characteristic time used to reach thermodynamics equilibrium is smaller than the time step. According to the phase equilibrium theory, the final equilibrium state requires the equaling of pressure, temperature and chemical potential between the two phases, as listed in Eq. (1.14), which corresponds to the maximum point of entropy and the minimum of Gibbs free energy. The procedure of computing 𝑃 and 𝑇 from 𝑢 * ,𝑣 * and 𝑧 𝑘 * or solving Eq. (1.8) is realized by the Isoenergetic-Isochoric (UV) flash. The detailed steps about this flash implemented in current work can be found in our recent publication [START_REF] Yi | A multicomponent real-fluid fully compressible fourequation model for two-phase flow with phase change[END_REF]. The real phase composition (𝑥 𝑘 , 𝑦 𝑘 ) and vapor mole fraction (𝜓 𝑣 ) are obtained through the phase equilibrium computation at given 𝑇, 𝑃 and molar fraction 𝑧 𝑘 . An integral phase equilibrium computation is started with the stability test to determine the phase state. If the systems turns out to be unstable, a follow-up phase split computation or Isothermal-Isobaric (TP) flash is proceeded to obtain the new phase composition (𝑥 𝑘 , 𝑦 𝑘 ). The validation of the phase equilibrium solver implemented in this study has been elaborated in Appendix A in which the VLE computation is applied to n-C 12 H 26 -N 2 and n-C 7 H 16 -N 2 systems. The results have correlated with experimental results well in most regions except the critical point zone in which PR EoS has some deficiency. The phase equilibrium computation is embedded into the internal iteration of the (UV) flash. A simplified implementation flowchart of this (UV) flash is depicted in Figure 1. The 4-Equation model assuming thermodynamics equilibrium is referred to 4EQ-PR-EQ model. 2) The second assumption is constructed by the fact that the time scale of relaxing chemical potential at the interface (𝑡 𝑟 ) is much longer than the mechanical and thermal relaxation process [START_REF] Zein | Modeling phase transition for compressible two-phase flows applied to metastable liquids[END_REF]. Here, we consider a limiting case for which (𝑡 𝑟 ) tends to infinity and no phase change occur. Therefore, the final relaxation will only realize the relaxation of pressure and temperature. Accordingly, the final phase composition (𝑥 𝑘 , 𝑦 𝑘 ) and vapor molar fraction (𝜓 𝑣 ) will sustain the initial nonequilibrium value ( 𝑥 𝑘 * , 𝑦 𝑘 * , 𝜓 𝑣 * ) without resorting to VLE calculation. The Newton iterative algorithm is adopted to resolve the objective functions Eq. (1.8) to obtain the final 𝑃,𝑇. The computation steps are concisely presented in the following flowchart, as shown in Figure 2. The 4-Equation model assuming no phase change is named as 4EQ-PR-Wo-EQ model or called frozen model. Both aforementioned 4-Equation models can deal with pure single-phase liquid or gas without adding any trifle impurities which are usually essential in 6-or 7-Equation models [START_REF] Habchi | a Gibbs Energy Relaxation (Germ) Model for Cavitation Simulation[END_REF], [START_REF] Zein | Modeling phase transition for compressible two-phase flows applied to metastable liquids[END_REF]. However, it is important to note that from the thermodynamics point of view, though single-phase, the injected liquid fuel contains a small amount of dissolved air corresponding to the vehicle reservoir pressure ( 𝑌 𝑎𝑖𝑟 ≈ 1𝐸 -5 at atmospheric pressure). This amount of dissolved air may have a significant effect on the inception and collapse of cavitation [START_REF] Duke | X-ray fluorescence measurements of dissolved gas and cavitation[END_REF], [START_REF] Battistoni | Comparison of Mixture and Multifluid Models for In-Nozzle Cavitation Prediction[END_REF], and thereby on the ensuing fuel injection process. Being able to simulate dissolved air in the fuel is one of the main features made by the suggested 4-Equation model. This can make the simulation more approach reality in conditions where numerous nucleus may exist (at wall heterogeneous nucleation for instance) in the sub-grid scale. In this condition, bubbles grow in low pressure regions and finally collapse and the formed shock wave may bring erosions to the wall around [START_REF] Koukouvinis | Large Eddy Simulation of Diesel injector including cavitation effects and correlation to erosion damage[END_REF], [START_REF] Örley | Large-Eddy Simulation of turbulent, cavitating fuel flow inside a 9-hole Diesel injector including needle movement[END_REF].

Coupling of flow solver with thermal solver

In this Section, the coupling of the thermodynamics solver with the flow solver is summarized in order to explain the main implementation stages and recall key numerical features. The detailed introduction about the IFP-C3D solver can be found in [START_REF] Bohbot | IFP-C3D: an Unstructured Parallel Solver for Reactive Compressible Gas Flow with Spray[END_REF]. In this original in-house code, the transport equations (mass, momentum, energy balance equations) are solved based on a time-splitting numerical scheme including four stages sequentially referred to below as Phase A, B, C and D, as depicted in Figure 3. Phase A computes the effect of spray and combustion as source terms. In Phase B or Lagrangian phase, the cell is moved with the fluid and all the physical properties such as pressure, temperature and velocity except of the convection terms are calculated implicitly with the SIMPLE numerical scheme [START_REF] Patankar | Numerical heat transfer and fluid flow: Computational methods in mechanics and thermal science[END_REF], including a BICGSTAB and SOR preconditioners [START_REF] Bohbot | IFP-C3D: an Unstructured Parallel Solver for Reactive Compressible Gas Flow with Spray[END_REF]. Then, in Phase C (Eulerian stage), the grid cell boundaries are mapped back to their original position (in the absence of wall movement). The obtained solutions from Phase B are updated in the Phase C using a quasi-second order upwind (QSOU) explicit numerical scheme. The Minmod slope limiter is used for scalar fluxes, and Van Leer slope limiter is used for momentum fluxes (see [START_REF] Bohbot | IFP-C3D: an Unstructured Parallel Solver for Reactive Compressible Gas Flow with Spray[END_REF]). The final stage, phase D contains the stiff relaxation algorithms for the interfacial velocity (7-Equation systems only), pressure, temperature and chemical potentials (if with phase change) [START_REF] Habchi | a Gibbs Energy Relaxation (Germ) Model for Cavitation Simulation[END_REF]. Inside the original IFP-C3D code, the inherent EoS is the ideal gas EoS for gas phase and the SG-EoS for liquid phase. In the present work, all the subroutines in the code involving EoS have to be replaced by PR-EoS (see the stages marked with shadow, in Figure 3). Due to the fact that liquid and gas phases are both solved with multicomponent system, the involved dissolved gas in liquid phase makes the model totally distinct with previous research using SG-EoS [START_REF] Habchi | a Gibbs Energy Relaxation (Germ) Model for Cavitation Simulation[END_REF], [START_REF] Chiapolino | A simple and fast phase transition relaxation solver for compressible multicomponent two-phase flows[END_REF], [START_REF] Zein | Modeling phase transition for compressible two-phase flows applied to metastable liquids[END_REF]. Firstly, as for the initialization of the simulation and to obtain the phases composition at given temperature, pressure and component composition, a phase equilibrium computation or vapor-liquid equilibrium (VLE) calculation is carried out. The initial VLE computation has assumed the fluid is in saturated state at the beginning of the simulation. Next, if the simulated configuration includes inlets or outlets, a VLE computation is performed again to keep the boundary conditions consistent with the main field. Then, the transport equations are resolved as described above from Phase A to Phase C. Within the Lagrangian phase B, after the pressure and temperature are updated during the SIMPLE algorithm, a series of thermal properties must be updated correspondingly by solving the single-phase PR EoS. Following Phase C, with the known internal energy, specific volume and component composition, a new temperature, pressure and phase compositions must be calculated for a new time-step (or cycle). This is the role attributed to the UV flash model implemented in Phase D. In fact, the UV flash process with PR EoS has replaced the relaxation procedures for the pressure, the temperature and the Gibbs energy in the original 7-Equation model.

One noting point is that in the current coupled system, each phase is resolved with an independent PR EoS which ensures the validity of the speed of sound especially inside the vapor dome. The composite EoS can help avoid entering the negative sound speed region as discussed in previous work [START_REF] Yi | A multicomponent real-fluid fully compressible fourequation model for two-phase flow with phase change[END_REF] and Appendix A(2).

1D test cases simulation

Before performing the multidimensional ECN Spray A injector simulations, the 4EQ-PR-EQ and 4EQ-PR-Wo-EQ models are employed to investigate the spurious oscillations issue which is one of the main challenges in transcritical modelling [START_REF] Ma | An entropy-stable hybrid scheme for simulations of transcritical real-fluid flows[END_REF], [START_REF] Pantano | An oscillation free shock-capturing method for compressible van der Waals supercritical fluid flows[END_REF]. This study is firstly carried out using a 1D advection tube in transcritical conditions. Then, the 4EQ-PR-EQ model is further applied to simulate a 1D shock tube to verify the capability of current solver in dealing with strong shock relevant to the injector needle opening conditions.

1D transcritical advection tube

In the 1D advection tube, a 0.5 m long liquid n-dodecane is initially located in the center of a 1m long tube. The operating conditions are similar with the Spray A conditions. The initial liquid temperature is 363 K. The rest of the tube is full of gaseous nitrogen with a temperature of 900 K. A constant pressure of 60 bar is also specified. The advection velocity is set with 100 m/s with periodic boundary conditions. The grid size sensitivity for both 4EQ-PR-EQ and 4EQ-PR-Wo-EQ models is tested in this section, and the grid size ranges from 0.02 to 1.0 mm. The numerical profiles of pressure, temperature, velocity and density at an instant of 0.2 ms with the two models are presented in Figure 4 (a,b,c,d) and compared to the exact solution. As expected, the main spurious oscillations appear at the contact discontinuity locations, as shown in the profiles of pressure and velocity (see Figure 4 (a,b)). The amplitude of the oscillations reduces significantly as the grid is refined. It is well known that the spurious oscillation problem has posed a severe challenge to the transcritical injection modelling [START_REF] Ma | An entropy-stable hybrid scheme for simulations of transcritical real-fluid flows[END_REF]- [START_REF] Terashima | Investigations of trans / supercritical jet mixing using a high-resolution numerical method[END_REF]. The problem appears as the finite volume discretization strategy is applied to multicomponent system with high concentration gradient, even when ideal-gas EoS is used, and deteriorates with the non-linear real fluid EoS such as cubic series. Many researchers have contributed to unravelling this difficulty [START_REF] Ma | An entropy-stable hybrid scheme for simulations of transcritical real-fluid flows[END_REF], [START_REF] Pantano | An oscillation free shock-capturing method for compressible van der Waals supercritical fluid flows[END_REF], [START_REF] Abgrall | How to prevent oscillations in multicomponent flow calculations: A quasi conservative approach[END_REF]. Nevertheless, no extra efforts are paid to address this issue in current work since no stability problems are encountered during transporting such strong gradient and the spurious oscillations are controllable with enough refined mesh. The sharpness of the interface for these two models are well kept, as seen in the density profile (Figure 4 (d)). The other noting point is that the model with VLE has demonstrated much less oscillations compared with the frozen model as observed in the profiles of temperature and velocity (Figure 4 (b,c)). This trend can also be confirmed in the evolution profile of the averaged velocity in the tube, as depicted in Figure 5 (e). After one cycle period of the advected liquid (10 ms), the oscillations amplitude of the 4EQ-PR-EQ model has proved to be much smaller than that of the 4EQ-PR-Wo-EQ model. It implies that reaching the thermodynamics equilibrium in each cell at the end of each time step helps in reducing spurious oscillations. In addition, due to the fact that the internal energy instead of the total energy is transported in the energy equation Eq. (1.4), the overall energy loss in the whole computation has been checked in Figure 5 (f). The amount of total energy loss in one cycle for these two models is less than 1%. However, the energy loss error for the 4EQ-PR-EQ model (0.01%) is much less compared to the frozen model. Thereby, it is concluded that the error induced by the non-conservative formulation of the energy balance equation will not affect the overall computational accuracy in limited timescale and the model with phase change (4EQ-PR-EQ) is more stable than without phase change. (a,b,c,d) profiles present the evolution of pressure, temperature, velocity and density at 0.2 ms with the two varied mesh sizes of 1mm and 0.2 mm, respectively. 6

4 E Q -P R -E Q ( x = 1 m m ) 4 E Q -P R -E Q ( x = 0 .0 2 m m ) 4 E Q -P R -W o -E Q ( x = 1 m m ) 4 E Q -P R -W o -E Q ( x = 0 .0 2 m m ) E x a c t (

1D shock tube

To further verify the capability of current model in dealing with two-phase strong shocks, the 4EQ-PR-EQ model is used to solve the 1D shock tube problem with Spray A condition [START_REF]ECN-Engine Combustion Network (Spray A&B)[END_REF]. The initial discontinuity is set in the middle of the tube. In particular, the left side is filled with pure liquid n-dodecane with the initial conditions of 150 MPa, 363 K and the right side is full of high temperature nitrogen with the initial conditions of 6 MPa, 900 K. The simulation has been conducted with various mesh resolutions ranging from 0.02 mm to 1 mm. The evolution of pressure, density, speed of sound and mass fraction of vaporous dodecane are presented below in Figure 6 (a, b, c, d). Firstly, obvious spurious oscillations are detected in the pressure profile with the mesh size of 1mm. Whereas the oscillations are dumped significantly as the mesh is refined to 0.2 mm and further to 0.02 mm. Affected by the strong expansion wave, a prominent decrease of pressure is detected in the left part, which leads to the reduction of density and the speed of sound correspondingly, as shown in Figure 6 (b,c). One noting point is that the volume shift formula is not adopted to adjust the density and speed of sound for the concern of potential non-consistence in thermodynamics. Thus, it may explain the relatively low-density value (~693 kg/m 3 ) and extremely high speed of sound (~2768 m/s) for n-dodecane at the initial pressure, 150 MPa and temperature, 363 K. The contact discontinuity is also well shaped, as depicted in the zoom of the density profile Figure 6 (b). Since the model used here is with phase change, the evaporation wave appears between this contact discontinuity and the expansion wave, as proved by the generated vaporous dodecane in Figure 6 (d). In addition, due to evaporation, a local two-phase zone is formed which indicates the decrease of speed of sound from 612 m/s at the compressed gas state to a minimum value of 253 m/s following the Wood formula (Eq. (A.1) )) [START_REF]A Textbook of Sound[END_REF], as it could be observed in Figure 6 (c). At last, the effect of compression wave is displayed in the slightly elevated density and speed of sound behind the evaporation front. One noting point is that the increasing mesh resolution is not only helpful to reduce the spurious oscillations, but also favorable to improve the accuracy of results especially in the wave front Figure 6 (a,b). Hence it is

4 E Q -P R -E Q ( x = 1 m m ) 4 E Q -P R -W o -E Q ( x = 1 m m ) (e ) t [s ]
|u -u 0 |/u 0 0 0 .0 0 2 0 .0 0 4 0 .0 0 6 0 .0 0 8 0 .0 1

1 0 -7 1 0 -6 1 0 -5 1 0 -4 1 0 -3 1 0 -2 4 E Q -P R -E Q (  x = 1 m m ) 4 E Q -P R -W o -E Q (  x = 1 m m ) ( f )
concluded that keeping a high mesh resolution is crucial to achieve a reliable result, especially in real fluid 1 simulation as previously confirmed in Ref. [START_REF] Pantano | An oscillation free shock-capturing method for compressible van der Waals supercritical fluid flows[END_REF]. This recommendation is generally largely respected when 2 simulating in-nozzle and near-nozzle two-phase flow, as described in the next Section for Spray A 3 simulations. 4

ECN Spray A injection

In this section, the aforementioned two-phase flow fully compressible diffuse interface models (4EQ-PR-EQ, 4EQ-PR-Wo-EQ) are applied to simulate the real size diesel injector Spray A from Engine Combustion Network (ECN) [START_REF]ECN-Engine Combustion Network (Spray A&B)[END_REF] at HTHP evaporating conditions. The simulation results from both models are compared to experimental data, along with a detailed analysis of the numerical results in order to complement the ECN experimental database.

Numerical setup

The Spray A injector is a common-rail single-hole injector with a nominal diameter of 0.09 mm. Taking into account the effect of in-nozzle flow and expected pressure oscillations, the configuration has included the lower control chamber just above the needle, in addition to the sac and orifice, as shown in Figure 7 (a). Thus, the whole computation domain contains not only the classical spray inside the chamber, but also more challenging in-nozzle flow. In the simulation, the needle lift is stationary at its maximum value taken as 500 μm. The computational grid was generated by using the ANSYS ICEM Package with body-fitted multi-block hexahedral cells. In addition, the actual shape of the sac and orifice available at the ECN website [START_REF]ECN-Engine Combustion Network (Spray A&B)[END_REF] is used (stl* for the injector serial # 210675). To keep the computational costs tractable, a varying grid strategy is employed aimed at different zones of interest. Particularly, the in-nozzle part and main spray region is much more refined compared to other regions and the far downstream is distributed with relative coarse mesh (Figure 7(c)). The orifice is discretized with 24 cells with an average size of 3.75 𝜇𝑚, as depicted in the zoomed view in Figure 7 (b,d). The minimum mesh size (0.8 𝜇𝑚) is located in the near nozzle region. The coarsest mesh with the size of 8 mm exists in the far downstream region, which is not of interest in this work. (d) sac as recently observed using X-ray experiments [START_REF] Manin | Transient cavitation in transparent diesel injectors[END_REF]. Indeed, due to multiple successive injection, the orifice and sac may be full of fuel or a mixture of fuel and air bubbles. The study here has referred to the conclusion from Payri, et al. [START_REF] Payri | Assessment on Internal Nozzle Flow Initialization in Diesel Spray Simulations[END_REF], in which they have proposed to initialize the simulations with the orifice filled with gas especially in the situation with high back pressure. As stated in this study [START_REF] Payri | Assessment on Internal Nozzle Flow Initialization in Diesel Spray Simulations[END_REF], the presence of gas in the nozzle may significantly impact the initial transient phase of the jet. In this work, the initial pressure distribution is illustrated in Figure 8. For the convenience of computation, the initial gradient is set in the middle of the orifice with almost pure liquid n-C 12 H 26 (including an initial N 2 mass fraction equaling to 10 -5 ) in the left side and pure N 2 in the right. The working fluid is n-C 12 H 26 with the temperature of 363 K. The chamber is initialized with N 2 at the temperature and pressure of 900 K and 60 bar, respectively, which corresponds to the gas density of 22.07 kg/m 3 . In this study, a constant pressure (150 MPa) boundary condition is specified at the inlet canal of the lower control chamber, as depicted in Figure 7 (a). The inlet temperature and concentrations are assumed the same as the liquid fuel initial conditions. Other numerical settings can be found in Table 1. The simulations are conducted within the large eddy simulation (LES) framework. The involved subgrid scale model is the Smagorinsky model, as described in Section 2.1. Binary Interaction Parameter 0.19 [START_REF] García-Córdova | Vapor-liquid equilibrium data for the nitrogen + dodecane system at temperatures from (344 to 593) K and at pressures up to 60 MPa[END_REF] Needle lift 500 μm

Comparison with experimental results

In this section, the numerical results from the 4EQ-PR-EQ and 4EQ-PR-Wo-EQ models are compared to experimental data. Since no phase equilibrium computation is considered in the 4EQ-PR-Wo-EQ model, the computational efficiency is much higher than for the 4EQ-PR-EQ model. Therefore, the simulation with phase change were run only for 100 μs, due to limited computational resource. In contrast, the case without phase change has been run more than 230 µs.

Jet evolution comparison

During the simulations, no numerical instabilities are encountered for both models. The jet evolution snapshots from these two models at the early time (up to 89 µs) are illustrated in Figure 9. The figures with back illumination are experimental data snapped directly from the ECN website [START_REF]ECN Spray A DBI675 injector[END_REF]. The relevant theory is according to the Mie-scattering theory [START_REF] Pickett | Measurement Uncertainty of Liquid Penetration in Evaporating Diesel Sprays[END_REF]. Based on this theory, the recommended criterion to evaluate the liquid penetration length should be the liquid volume fraction (LVF) less than 0.0015 in the evaporation conditions [START_REF] Pickett | Measurement Uncertainty of Liquid Penetration in Evaporating Diesel Sprays[END_REF]. Thereby, this criterion (LVF (𝛼 𝑙 ) = 0.0015) is applied to the numerical results of the 4EQ-PR-EQ model. In contrast, the criterion based on the mixture mass fraction of n-dodecane (𝑌 𝑐12ℎ26 = 0.6) referred to the previous studies employing the no phase change model [START_REF] Lacaze | Analysis of high-pressure Diesel fuel injection processes using LES with real-fluid thermodynamics and transport[END_REF], [START_REF] Ma | Modeling and Simulation of Diesel Injection at Transcritical Conditions[END_REF], is adopted in Figure 9 (right row) for the sake of comparison. Therefore in this Figure 9, the liquid penetration length denoted with blue iso-surface is determined with a liquid volume fraction criterion (LVF = 0.0015) for the 4EQ-PR-EQ model results; and a criterion based on mixture mass fraction of n-C 12 H 26 ((𝑌 𝑐12ℎ26 = 0.6) for the 4EQ-PR-Wo-EQ model results. In fact, as shown in Eq.(1.1), the flow system can provide the information of liquid volume fraction (𝛼 𝑙 ) even without the equilibrium model. However, the liquid penetration length predicted with the criterion of LVF (𝛼 𝑙 ) = 0.0015 has presented significant difference compared to the prediction with the criterion of 𝑌 𝑐12ℎ26 = 0.6 for the model without VLE. The comparisons from both standards are displayed in Appendix A(3) in detail. In all the images, the background contour denotes the temperature variation in the range of 353 K-900 K.

At the early time sequences (10 μs-30 μs), the modelling results of both models have demonstrated slightly longer penetration than the experimental value. Since no vapor is generated by the 4EQ-PR-Wo-EQ model, the tip of the jet stays less spread than for the 4EQ-PR-EQ model. Thereby, evaporation seems to trigger earlier radial spreading that has led to a slightly larger angle in the tip of the jet. This trend becomes more evident from 43 μs , because the evaporation becomes severe after a short heating period, as may be seen in Figure 9 (middle row). At this moment, the vapor cloud formed in the tip of the jet with 4EQ-PR-EQ model has spread much wider than previous time instant. As time further elapses, both models can estimate the liquid penetration length (about 9 mm) very well compared to the experimental value. However, obvious differences are detected in the vapor penetration, as discussed further below based on Figure 10 and Figure 11. Moreover, the temperature predicted in the model without phase change is lower than that with phase change as observed in the background contour in Figure 9 which will be further discussed in the following thermodynamics analysis (see Section 4.3,below). Even if a good qualitative agreement has been achieved for the liquid penetration, the near-nozzle spray angle seems underestimated compared to the experimental results (Figure 9). Indeed, a much slim jet is observed for both simulation results. One important reason for this may be attributed to the initial setting of pressure gradient in the middle of the hole. As shown in Figure 8, the initial pressure specified in the sac (1500 bar) has proved to be not appropriate because it far exceeds the actual pressure which should be similar to the order of the chamber pressure (60 bar). This is also the reason why the injection velocity is soared to 650 m/s in less than 15 μs, also much higher than the average 600 m/s in the experiment. This justifies the longer liquid penetration in the earlier time sequences as noted above. Since the current Spray A injection condition is still in subcritical regime with weak but existent surface tension, the consideration of the primary atomization in the near-nozzle region may remedy the underestimation of spray angle. Other causes may be from the omission of mass diffusion in the numerical model as discussed in Section 2.1.

Further study is needed to clarify this problem.

Jet penetrations and flowrate comparison

As a further quantified validation, the rate of injection (ROI) and the liquid and vapor penetration from both models have been compared with experimental data as plotted in Figure 10 and Figure 11. As mentioned above, affected by the pressure initialization (Figure 8), the overall penetration length from both models have exceeded the experimental results in the early time (< 30 μs). The unappropriated initialization has also led to a much higher mass flow rate at the initial time compared to the ROI from CMT [START_REF]CMT-Motores Térmicos, Rate of Injection Model[END_REF], [START_REF] Payri | A new methodology for correcting the signal cumulative phenomenon on injection rate measurements[END_REF], as shown in Figure 10. Influenced by the initial strong shock inside the orifice, the mass flux has soared to the maximum value (2.48 g/s) in less than 15 μs. After this early injection time, the mass flow rate is approaching the CMT predictions (Figure 10), which denotes the pressure variation in the nozzle has resumed to the normal range. As for the penetration predictions, both models can predict liquid penetration very well since the liquid penetration is largely affected by the ambient condition instead of the injection pressure [START_REF] Siebers | Liquid-Phase Fuel Penetration in Diesel Sprays[END_REF]. However, as for the vapor penetration, the results from 4EQ-PR-EQ model can basically have a good agreement in the early 100 μs. It is worth noting that a 10 µs delay has been adopted in Figure 11 to adjust the initial penetration slopes and ensure fair comparisons. Another noting point is that the penetration length estimated by using the vaporous n-C 12 h 26 in 4EQ-PR-EQ model can have a very close result with the predictions by employing the mixture mass fraction criterion (𝑌 𝑐12ℎ26 = 0.01). Whereas, significant deviations are detected for 4EQ-PR-Wo-EQ after 80 μs by using the same criterion. Since no vapor is generated from the 4EQ-PR-Wo-EQ model because of lacking the phase change model, the empirical estimation by using mixture mass fraction appears to be not suitable with current model. Indeed, the liquid length predicted by LVF equaling 0.0015 in 4EQ-PR-Wo-EQ model has presented the same length as the vapor penetration predicted by using the criterion (𝑌 𝑐12ℎ26 = 0.01). This coincidence implies the so-called 'vapor' in the 4EQ-PR-Wo-EQ is not appropriate for the modelling of the actual vapor coming from evaporation process as generated in the phase change model. The overestimation of vapor penetration has also been reported in recent simulation of the Spray A modelling with phase equilibrium model [START_REF] Matheis | Multi-component vapor-liquid equilibrium model for LES of highpressure fuel injection and application to ECN Spray A[END_REF] in which they have attributed this deficiency to the shortcoming of PR EoS in predicting liquid density. Some researchers have attributed this deficiency to an underestimation of turbulent mixing in radial direction which induces an overestimation of axial momentum and a faster convection of the vapor downstream [START_REF] Knudsen | Compressible Eulerian needle-totarget large eddy simulations of a diesel fuel injector[END_REF]. In current study, the causes can be attributed mainly to the unappropriated initialization of the pressure distribution inside the nozzle, as noted above.

Figure 10 Illustration of the mass flux or rate of injection (ROI) from 4EQ-PR-EQ, 4EQ-PR-Wo-EQ models, as well as the predictive data by using the injection rate model [START_REF]CMT-Motores Térmicos, Rate of Injection Model[END_REF], [START_REF] Payri | A new methodology for correcting the signal cumulative phenomenon on injection rate measurements[END_REF].

C M T 4 E Q -P R -E Q Figure 11
Illustration of liquid and vapor penetration for 4EQ-PR-EQ and 4EQ-PR-Wo-EQ models compared to experimental data from ECN website (https://ecn.sandia.gov/bkldaal4liquid/), referred to [START_REF] Pickett | Measurement Uncertainty of Liquid Penetration in Evaporating Diesel Sprays[END_REF], [START_REF] Pickett | Relationship Between Diesel Fuel Spray Vapor Penetration/Dispersion and Local Fuel Mixture Fraction[END_REF], [START_REF] Manin | SP2-4 Evaluation of the liquid length via diffused back-illumination imaging in vaporizing diesel sprays(SP: Spray and Spray Combustion,General Session Papers)[END_REF]. The evaluation of vapor penetration is based on the mass fraction of n-C 12 H 26 (𝑌 𝑐12ℎ26 = {0.015,0.01,1e-5}) for 4EQ-PR-Wo-EQ model. The penetration length decided by vaporous n-C 12 H 26 is also evaluated for 4EQ-PR-EQ model. The evaluation of liquid penetration is based on the liquid volume fraction with a critical value of 0.0015 (LVF = 0.0015) for 4EQ-PR-EQ and the mass fraction of n-C 12 H 26 with a critical value of 0.6 (𝑌 𝑐12ℎ26 = 0.6) for 4EQ-PR-Wo-EQ model. The penetration length determined by the liquid volume fraction of 0.0015 (LVF = 0.0015) is also evaluated for 4EQ-PR-Wo-EQ model.

Mixture mass fraction and velocity comparison

More validations concerning mixture mass fraction of n-C 12 H 26 and actual velocity in radial direction are performed for 4EQ-PR-Wo-EQ model as shown in Figure 12. The instantaneous LES results of mass fraction are depicted in Figure 12 (left) at the time interval of 156-230 𝜇𝑠. A good agreement has been achieved in the amplitude of mass fraction even though there exist some deviations in the radial width which may be caused by the insufficiency of radial turbulent mixing. In addition, the velocity distribution at the time interval of 212-222 𝜇𝑠 is also shown in Figure 12 (right). It is noted that the velocity distribution has an excellent agreement with experimental results. These agreements imply the evaporation may not seriously affect the mass and momentum distribution in the radial direction.

Figure 12 Validations of mixture mass fraction (Left side) and velocity (Right side) with experiments at an axial distance of 19.84 mm and 17.5 mm from the nozzle exit respectively. The instantaneous LES results (green solid line) and the average value (red bold line) are from 4EQ-PR-Wo-EQ model. The experimental data at 217 𝜇𝑠 (blue dash line) are referred to [START_REF] Pickett | Relationship Between Diesel Fuel Spray Vapor Penetration/Dispersion and Local Fuel Mixture Fraction[END_REF], [START_REF] Meijer | Engine Combustion Network:" Spray A" basic measurements and advanced diagnostics[END_REF].

Thermodynamics analysis

The T-x diagram of n-C 12 H 26 and N 2 at the chamber pressure 60 bar is plotted in Figure 13. The adiabatic mixing temperature (AMT), usually called frozen temperature 𝑇 𝐹 and equilibrium temperature 𝑇 𝐸 are also calculated using PR EoS. The scattered points are obtained from current numerical simulations at the instant of 50 μs. The mixture temperature distribution from equilibrium 4EQ-PR-EQ model can follow very well the equilibrium temperature 𝑇 𝐸 in almost the whole N 2 concentrations range except a small temperature deviation that may be seen in the single-phase liquid region where the N 2 mole fraction is smaller than 5%. Indeed, since the basic 𝑇 𝐹 and 𝑇 𝐸 lines are calculated with the initial fuel temperature of 363K, the lower temperature detected in the simulations could prove the fuel has been through a cooling stage in the hole before entering the chamber, which will be confirmed in the following Section 4.4. Even if phase change is not considered, the trend of mixture temperature evolution from 4EQ-PR-Wo-EQ model can still correlate well with equilibrium temperature, particularly when the N 2 mole fraction is smaller than 0.6. Evident deviations start at N 2 mole fraction higher than 0.75, where the equilibrium temperature differs significantly from frozen temperature as shown in Figure 13. Indeed, a slightly lower temperature is detected by 4EQ-PR-Wo-EQ model near the two-phase vapor side. This temperature differences further enlarge with increasing N 2 mole fraction up to the two-phase vapor side limit. After transiting to single-phase gas region (with more than 90% N 2, see Figure 13), the temperature profile resumes to the equilibrium trend. One noting point is that even only the internal energy is conserved in the flow system (see Eq.(1.4)), no evident error is detected in the temperature prediction by these two models.

Radial distance [mm]

Effect of evaporation on the jet

As aforementioned in Section 2.2.2, one of the main advantages of the current model is to be able to predict the vapor fraction (𝜓 𝑣 including N 2 and C 12 H 26 ) in each cell after solving the transport equations. Indeed, this thermodynamic parameter represents the extent of fuel and gas mixing. The evolution of 𝜓 𝑣 for the model without VLE is mainly driven by the mixing of the fuel with N 2 as the jet develops. In this case, the evolution of 𝜓 𝑣 is primarily due to the gas entrainment . In contrast, the contribution of 𝜓 𝑣 will include the fuel evaporation and mixing process for the model considering phase change. However, even though the evaporation is very intense at high temperature conditions, the dominant gas in the vapor phase is still N 2, as shown in Figure 14(a). The amount of N 2 exceeds 90% in terms of the molar fraction in the vapor phase ( 𝑦 𝑛2,𝑣𝑎𝑝𝑜𝑟 ), as denoted with the isolines in Figure 14 (a). This implies the injection is dominated more by the mixing process instead of the evaporation as mentioned by previous researchers [START_REF] Dahms | Understanding high-pressure gas-liquid interface phenomena in Diesel engines[END_REF]. Nevertheless, the temperature and the jet shape still present significant differences between these two models in the downstream dilute zone. The overall temperature for the model considering VLE is slightly higher than without equilibrium as confirmed in above T-x diagram (Figure 12). Indeed, affected by the vapor flow rate, the tip of the jet is much wider and more spread than the model without phase change. It seems that in the model with phase change, the interaction due to the counter flows between gas entrainment and evaporation flow from the interface has triggered the instability of the jet slightly earlier which has led to stronger turbulence, as verified by the eddy-viscosity in Figure 14 (e,f) demonstrate the eddyviscosity distribution in the near nozzle region for these two models.

More flow information from VLE model

Except the computational time consumption issue, the VLE model can provide abundant and accurate information related with phase properties. As aforementioned, one important step of VLE computation is to determine the phase state. By solving TPD functions, the real phase number and state can be verified.

The study here has adopted three integers to represent phase state (TPD = 0: single-gas phase; TPD = 2: two-phase; TPD = 1: single-liquid phase) [START_REF] Yi | A multicomponent real-fluid fully compressible fourequation model for two-phase flow with phase change[END_REF]. To identify the real fluid state during the jet developing, an iso-surface of TPD value equaling 2 (two-phase region) is presented in Figure 15 (a). As shown in this Figure, the jet starts evolving from the exit of the orifice with single-liquid state corresponding to almost pure liquid n-dodecane and forms an intact liquid core. One noting point is the liquid boundary delimited with the criterion (LVF = 0.0015) can have a good agreement with the boundary of two-phase region (TPD = 2) except the near nozzle liquid core which confirms the previous predictions in the phase diagram. With phase equilibrium model, the density of each component at any phase is also available. In current study, the mixture density is computed with 𝜌 𝑚𝑖𝑥 = 𝛼 𝑔 𝜌 𝑔 + 𝛼 𝑙 𝜌 𝑙 . The maximum density (646 kg/m3) is in the liquid core where the liquid volume fraction approaches one (Figure 15 (b)). As more N 2 is mixed into the downstream jet, the mixture density of the jet decreases significantly. In Figure 15 (c), the mass fraction of each component in any phase (𝑌 𝑝,𝑘 ) is defined as the ratio of specific density with respect to mixture density 𝑌 𝑝,𝑘 = 𝛼 𝑝 𝜌 𝑝,𝑘 (∑ 𝛼 𝑝 𝜌 𝑝,𝑘 ) 𝑘 . Thereby, the mass fraction of vaporous dodecane in Figure 15 (c) represents the percentage of generated vaporous dodecane related with the total mixture. The generated vapor is accumulated in the jet front where the mass fraction of vaporous dodecane reaches 30%. The maximum specific density of vaporous dodecane is around 20 kg/m 3 , as shown in Figure 15 (d). Figure 15 ((e), (f)) present the variation of dissolved N 2 in the jet. Since the chamber pressure is very high (~60 bar), the dissolved N 2 part becomes non-negligible. The mass fraction of dissolved N 2 reaches 1% mostly located in the two-phase zone where the N 2 specific density arrives to 6 kg/m 3 .

In-nozzle flow analysis

The important influence of in-nozzle flow on the downstream spray development has been widely investigated in the community [START_REF] Xue | Eulerian CFD Modeling of Coupled Nozzle Flow and Spray with Validation Against X-Ray Radiography Data[END_REF], [START_REF] Battistoni | Large-Eddy Simulation (LES) of Spray Transients: Start and End of Injection Phenomena[END_REF]. Thereby, to know the instantaneous flow variation inside the sac and orifice, the evolution of pressure, temperature, compressibility factor and density stretching from the sac up to chamber (Figure 16) are illustrated in Figure 17 (a,b) based on the predictions of 4EQ-PR-EQ model. Firstly, an abrupt pressure drop is observed inside the orifice. The strong shock has brought in rapid increase of velocity from 0 m/s to 650 m/s in less than 15 𝜇𝑠. The sudden increase of velocity also induces cooling effect on the fuel (~10 K) which can be verified in the variation of temperature profile as shown in Figure 17 (b). Thus, it corresponds to the lower temperature in the single-phase liquid region of the T-x diagram, as illustrated in Figure 13. One noting point is the evolution of compressibility factor. Since it represents the repulse force between the molecules at high pressure condition, the compressibility factor also far exceeds 1 even in pure liquid condition as shown in Figure 17 (b). From the sac to the chamber, the fluid has been through the transition from the extreme dense liquid to less dense liquid which corresponds to the variation of the compressibility factor from 12 to 0.6. Accordingly, the strong expansion has brought in a reduction of density by around 40 kg/m 3 . The significant variation of compressibility factor and density imply the compressibility of liquid is definitely not negligible in the high pressure injection simulations.

Parametric study

A series of parametric studies including different ambient temperature and pressure are carried out in this section. Two test cases keeping the similar density (~22 kg/m 3 ) (Table 2) as the previous simulations are performed. Since the simulations show very high demand in CPU resource, both cases have not been run more than 80 𝜇𝑠. The early jet evolution for these two cases are displayed in Figure 18 and Figure 19.

Both simulations are performed with 4EQ-PR-EQ model. Both cases can achieve a good agreement with experimental results in the prediction of liquid penetration. With higher temperature and pressure in the chamber, the liquid penetration has shortened (case 1) compared to the case 2 with lower temperature and pressure. The cause is attributed to much less evaporation appearing in case 2 as proved in the contour of the mass fraction of vaporous dodecane (Figure 20 (d)) [START_REF] Naber | Effects of Gas Density and Vaporization on Penetration and Dispersion of Diesel Sprays[END_REF]. The region of the mass fraction of vaporous dodecane locating at the range of 0.3-0.5 in case 1 is obviously larger than in case 2 ( Figure 20 (b,d)). However, the two-phase zone in the low temperature case is much broader than in high temperature situation. In addition, the spray angle with high evaporation rate is also narrower than for the low temperature case, which implies the strong evaporation can bring cooling effect on the jet [START_REF] Pickett | Relationship Between Diesel Fuel Spray Vapor Penetration/Dispersion and Local Fuel Mixture Fraction[END_REF]. Since the initialization is the same as previous simulation, the vapor penetration length also slightly exceeds the experimental value as reported in previous case.

Summary, conclusions and future work

In this work, a fully compressible Eulerian-Eulerian based two-phase flow model combined with real fluid EoS is developed and applied to the simulation of ECN Spray A injector. To be more specific, the thermal solver can address the phase change problems with phase equilibrium theory and also be capable of solving the EoS without equilibrium assumption. Indeed, the use of the latter non-equilibrium thermodynamics model is justified as the time scale of relaxing chemical potential is much longer than relaxing pressure and temperature. Especially in extremely fast injection modelling, it is reasonable to assume the fluid has not reached the equilibrium state within each time step.

Firstly, the coupled flow and thermodynamics solvers have been used to investigate the spurious oscillations issue. The results have confirmed that the spurious oscillations problem exists in the current solver even though the flow solver is quasi-conservative in the sense of using the internal energy instead of total energy. However, the oscillations have been shown to be significantly dumped as the mesh is refined. Since no stability problem is encountered during all the simulations, the study here has not paid extra efforts to address this issue thoroughly. Then, the solver has been applied to model the strong shock in a 1D shock tube with the Spray A conditions. The involved expansion wave, contact discontinuity, compression wave and evaporation wave are shown all resolved well with the current model when the grid is sufficiently refined. Finally, the model is used to simulate a full 3D Spray A injector including the innozzle flow. Several conclusions are drawn from the modelling results:

1) High mesh resolution can enhance the accuracy of numerical results by effectively reducing the spurious oscillations. In addition, the model considering VLE has presented higher numerical stability than the model without VLE.

2) The initialization of in-nozzle flow plays a significant role in the early jet development, cone angle and therefore on vapor phase penetration.

3) A good agreement of the vapor, liquid penetration length and rate of injection predictions has been achieved for the model considering phase change.

4)

As for the frozen model, excellent agreements have been achieved in the prediction of radial mass and velocity distribution. However, some uncertainties and deviations are detected in the estimation of penetration length according to different measure criterions.

5) Both models have underestimated the spray angle which may be caused by the non-fully resolved shear stress layer and the omission of enthalpy diffusion and mass diffusion terms in transport equation. Since the current Spray A injection condition is still in subcritical regime with weak but existent surface tension, the consideration of the primary atomization in the near-nozzle region may remedy the underestimation of spray angle. Further study is needed to clarify this problem.

6) The error caused by the non-conservative energy equation has not adversely affected the temperature prediction during CFD modelling.

7) Even though the fuel vapor molar fraction contribution is lower than 10% with respect to the overall vapor fraction, the evaporation still presents significant effects on the temperature distribution, turbulence intensity and jet development, which may be important for combustion.

8) The extreme large variation of liquid compressibility factor and density in the nozzle implies the importance of considering the compressibility of liquid during injector simulations.

Since evaporation still plays a critical role in the accurate modelling of the spray development, it will be essential to turn to more efficient method such as tabulation method to replace the iterative phase equilibrium solver. Furthermore, the addition of subgrid models aimed at modelling the atomization and break-up phenomena will probably enhance the prediction of spray angle. More accurate EoS is also required. These are among our future work in order to be able to compute the ECN Spray C and G injectors. Finally, the new model proposed in this article paves the way for the investigation of other different industrial configurations like aeroengines, rockets and aero-space vehicles. phase composition (𝑦 𝑘 𝑜𝑟 𝑥 𝑘 ) at the constant temperature of 593.5 K. In this Figure, it is evident to find that the curve plotted with liquid phase composition (𝑥 𝑘 , blue triangles) has coincided totally with the curve plotted using the feed (𝑧 𝑘 ). The vapor curve (𝑦 𝑘 , red dots) also corresponds well to the right side of the feed (𝑧 𝑘 ) curve. The reconstruction process inside the vapor dome is therefore realized by defining mixture rules combining the liquid phase and vapor phase results. During computation, the selection of the liquid volume must exceed the minimum molecular volume 𝑏 (see Eq. (1.5)) based on the PR-EoS as depicted in Figure A.2 . Moreover, if three roots are returned when solving the phase cubic equation, the middle root must be discarded to ensure the stability of the flow system. In contrast, the maximum root is selected for the vapor phase. With these effective liquid and vapor molar volumes, the speed of sound is guaranteed to be real in each phase. Thereby, the mixture sound speed computed with Eq. (A.1) is also ensured to be real even in two-phase zone. However, the possibility of losing hyperbolicity exists in the Lagrangian Phase B stage during solving the transport equation. One strategy of solving this issue is to restart the current cycle of the simulation with reduced time step. In which, 𝐾 𝑖 is equilibrium constant. 𝑁 𝑐 denotes number of components.

3) Liquid penetration length predictions with different criterions

As discussed in Section 4.2.1, the liquid penetration length predicted with the empirical criterion 𝑌 𝑐12ℎ26 = 0.6 [START_REF] Lacaze | Analysis of high-pressure Diesel fuel injection processes using LES with real-fluid thermodynamics and transport[END_REF] can achieve a good agreement with experimental results. However, much longer penetration is detected as with the criterion of LVF (𝛼 𝑙 ) = 0.0015 as depicted in Figure A.3. Since the criterion of LVF (𝛼 𝑙 ) = 0.0015 is validated by the experimental theory [START_REF] Pickett | UNCERTAINTY QUANTIFICATION FOR LIQUID PENETRATION OF EVAPORATING SPRAYS AT DIESEL-LIKE CONDITIONS[END_REF] and the criterion of using mixture mass fraction is more prone to be empirical adjustment, the deviation of liquid penetration length predicted with the more reliable liquid volume fraction method seems imply the importance of the evaporation in the whole spray developing to some extent.

Figure 1

 1 Figure 1 Flow chart of 4EQ-PR-EQ model, the inputs denoted by (*) are the solution from the flow solver.

Figure 2

 2 Figure 2 Flow chart of 4EQ-PR-Wo-EQ model, the inputs denoted by (*) are the solution of the flow solver.

Figure 3

 3 Figure 3 Schematic of coupling of thermo-solver with flow solver in IFP-C3D

1 Figure 4

 14 Figure 4 1D advection tube operated with Spray A thermal conditions ( 𝑇 𝑓𝑢𝑒𝑙 = 363 𝐾, 𝑃 𝑎𝑚𝑏. = 6 𝑀𝑃𝑎, 𝑇 𝑎𝑚𝑏. = 900 𝐾).(a, b, c, d) profiles present the evolution of pressure, temperature, velocity and density at 0.2 ms with the two varied mesh sizes of 1mm and 0.2 mm, respectively.

Figure 5

 5 Figure51D advection tube operated with Spray A thermal conditions ( 𝑇 𝑓𝑢𝑒𝑙 = 363 𝐾, 𝑃 𝑎𝑚𝑏. = 6 𝑀𝑃𝑎, 𝑇 𝑎𝑚𝑏. = 900 𝐾). (e, f) profiles illustrate the variation of average velocity and error of total energy in one cycle (10 ms) with the mesh size of 1mm and a CFL value equaling 0.2, respectively.

Figure 6 8 x

 68 Figure 6 1D two-phase shock tube operated with Spray A conditions ( 𝑃 𝑖𝑛𝑗. = 150 𝑀𝑃𝑎, 𝑇 𝑓𝑢𝑒𝑙 = 5 363 𝐾, 𝑃 𝑎𝑚𝑏. = 6 𝑀𝑃𝑎, 𝑇 𝑎𝑚𝑏. = 900 𝐾) using the 4EQ-PR-EQ model. (a, b, c, d) denotes the evolution 6of pressure, density, speed of sound and mass fraction of vaporous dodecane at an instant of 0.2 ms with a 7 varied mesh resolution from 1 mm, 0.2 mm to 0.02 mm and a CFL value equaling 0.2. 8

Figure 7 (

 7 Figure 7 (a) Configuration of whole computational zone, note the lower control chamber below the fuel inlet canal. (b) Cut slice of the computational domain. (c) Mesh illustration in the cut plane, note that only the near nozzle region (20 mm axial length) is well refined. (d) Zoomed region in the near nozzle zone.

Figure 8

 8 Figure 8 Demonstration of the initial pressure distribution in the whole computational domain (needle lift = 500 µm).

Figure 9

 9 Figure 9 Illustration of jet evolution at different time instants. The left column is the experimental data snapped from (https://ecn.sandia.gov/dbi675/) website and the middle and right columns are the modelling results from 4EQ-PR-EQ model and 4EQ-PR-Wo-EQ model respectively. The liquid penetration length denoted with blue iso-surface is determined with a liquid volume fraction criterion (LVF = 0.0015) for the 4EQ-PR-EQ model; and a criterion based on the mixture mass fraction of n-C 12 H 26 (𝑌 𝑐12ℎ26 ) = 0.6 for the 4EQ-PR-Wo-EQ model. The background contour represents the temperature distribution ranging from 353 K to 900 K.

Figure 13 T

 13 Figure13T-X diagram of n-C 12 H 26 and N 2 system at a pressure of 60 bar. 𝑇 𝐸 denotes the equilibrium temperature between the mixture. 𝑇 𝐹 symbols the frozen temperature or adiabatic mixing temperature. The red scattered triangles are obtained from CFD modelling with 4EQ-PR-EQ model at t = 50 μs. The blue scattered squares are obtained from CFD modelling with 4EQ-PR-Wo-EQ model at t = 50 μs.

2 TFigure 14 Figures

 214 Figure 14 Figures (a, b) depicts the distribution of vapor molar fraction at 80 μs for the 4EQ-PR-EQ and 4EQ-PR-Wo-EQ models respectively. The two isolines denote the molar fraction of N 2 in the vapor phase (𝑦 𝑛2,𝑣𝑎𝑝𝑜𝑟 ) with the value of 0.95 (white line) and 0.99 (black line) respectively in Figure (a). Figures (c,

  d) demonstrate the temperature distribution at 80 μs for both models. Figures

Figure 15 (

 15 Figure 15 (a) depicts two-phase region in the jet by using an iso-surface of TPD = 2 at t = 62 μs. The dark isolines delimit the liquid penetration length based on the criterion (LVF = 0.0015). (b, d, f) depict the density of mixture fluid (𝜌 𝑚𝑖𝑥 ), vaporous dodecane (𝜌 𝑐12ℎ26,𝑣𝑎𝑝𝑜𝑟 ) and N 2 in liquid phase(𝜌 𝑛2,𝑙𝑖𝑞𝑢𝑖𝑑 ) in the two-phase zone of the jet at the time of 62 μs, respectively. Figures(c, e) demonstrate the mass fraction of vaporous n-C 12 h 26 (𝑌 𝑐12ℎ26,𝑣𝑎𝑝𝑜𝑟 ) and dissolved N 2 in the liquid phase(𝑌 𝑛2,𝑙𝑖𝑞𝑢𝑖𝑑 ) in the twophase region at the time of 62 μs, respectively. The background contour in all pictures represents the temperature variation from 353 K to 900 K.

Figure 16

 16 Figure 16 Demonstration of the axial direction from the sac to the chamber.

Figure 17 (

 17 Figure 17 (a) In-nozzle evolution of pressure and temperature from the sac to the chamber. (b) In-nozzle evolution of compressibility factor and density from the sac to the chamber.

Figure 18

 18 Figure 18 Illustration of early jet evolution for case 1 (𝑇 a𝑚𝑏. = 1200 K, 𝑃 a𝑚𝑏. = 8 MPa). The liquid penetration length marked with blue iso-contour is presented with the criterion of LVF = {0.0015}. The left column is the experimental images snapped directly from ECN website (https://ecn.sandia.gov/dbi675/). The right column corresponds to the results of 4EQ-PR-EQ model. The background contour represents the temperature field at the range of 353 K-1200 K.

Figure 19 Figure 20

 1920 Figure 19 Illustration of early jet evolution for case 2 (𝑇 a𝑚𝑏. = 700 K, 𝑃 a𝑚𝑏. = 4.6 MPa). The liquid penetration length marked with blue iso-contour is presented with the criterion of LVF = {0.0015}. The left column is the experimental images snapped directly from ECN website (https://ecn.sandia.gov/dbi675/). The right column corresponds to the results of 4EQ-PR-EQ model. The background contour represents the temperature field at the range of 353 K-700 K.

Figure A. 2

 2 Figure A.2 Schematic of Pv diagram based on the molar fraction feed (𝑧 𝑘 , black square), liquid molar fraction (𝑥 𝑘 , blue triangle). And vapor phase molar fraction (𝑦 𝑘 , red point) at a constant temperature of 593.5 K. The black arrow mimics the phase split process. 1 ρ 𝑚𝑖𝑥 C s,mix 2

Figure A. 3

 3 Figure A.3 Illustration of jet evolution from 4EQ-PR-Wo-EQ model at different time instants. The liquid penetration length denoted with blue iso-surface is determined with the mixture mass fraction of n-C 12 H 26 criterion (𝑌 𝑐12ℎ26 = 0.6) in the left column and the liquid volume fraction criterion (LVF = 0.0015) in the right column. The background contour represents the temperature distribution ranging from 353 K to 900 K.

Table 1

 1 Operating conditions for the ECN Spray A modelling

		4EQ-PR-EQ (With Equilibrium solver)
	Models
		4EQ-PR-Wo-EQ (Without Equilibrium solver)
	Compressibility	Fully compressible two-phase flow
		Large Eddy Simulation;
	Turbulence model
		Sub-grid scale Smagorinsky model
	Grid type	Hexahedral
		Total cells number: 2505255
	Mesh resolution	Minimum cell size: 0.8 𝜇𝑚
		Maximum cell size: 8 mm
	Time integration precision	First order
	Spatial discretization	Second order
	Time step	2E-10, CFL: 0.12
	Inlet Boundary Condition	Pressure, 150 MPa
	Outlet Boundary Condition	Pressure, 6 MPa

Table 2

 2 Operating conditions for ECN Spray A modelling

	Case No.	𝑇 a𝑚𝑏. [K]	𝑃 a𝑚𝑏.	𝜌 𝑎𝑚𝑏.
			[MPa]	[kg/m 3 ]
	1	1200	8	22.04
	2	700	4.6	22.79

Acknowledgement

This project has received funding from the European Union Horizon 2020 Research and Innovation program. Grant Agreement No 675528 for the IPPAD project. The author would like to thank Mr. Kumar Rajesh for his generous help in preparing the mesh and Dr. Gilles Bruneaux for his valuable guidance in the injector modelling and reviewing the paper.

𝑙,𝑘 * 𝑀𝑤 𝑙 * ∑ (𝛼 𝑙 * 𝜌 𝑙,𝑘 * ) 𝑘 * 𝑀𝑤 𝑘 * , 𝑦 𝑘 * = 𝛼 𝑔 * 𝜌 𝑔,𝑘 * 𝑀𝑤 𝑔 * ∑ (𝛼 𝑔 * 𝜌 𝑔,𝑘 * ) 𝑘 * 𝑀𝑤 𝑘 * , 𝜓 𝑣 * = (∑ 𝛼 𝑔 * 𝜌 𝑔,𝑘 * ) * 𝑀𝑤 𝑚𝑖𝑥 * 𝑘 (∑ 𝛼 𝑝 * 𝜌 𝑝,𝑘 * ) 𝑝,𝑘 * 𝑀𝑤 𝑔 * (1.12) 𝑌 𝑘 * = (∑ 𝛼 𝑝 * 𝜌 𝑝,𝑘 * ) 𝑝 (∑ 𝛼 𝑝 * 𝜌 𝑝,𝑘 * ) 𝑝,𝑘 , 𝑧 𝑘 * = (∑ 𝛼 𝑝 * 𝜌 𝑝,𝑘 * ) 𝑝 * 𝑀𝑤 𝑚𝑖𝑥 * (∑ 𝛼 𝑝 * 𝜌 𝑝,𝑘 * ) 𝑝,𝑘 * 𝑀𝑤 𝑘 * (1.13)

Appendix A

1) Validation of phase equilibrium computation

To validate the vapor-liquid equilibrium model, isothermal-isobaric flash (TP flash) is computed for two binary hydrocarbon and nitrogen system which are widely encountered in industry. The thermal properties for each component are summarized in Table 3. 𝑇 𝑐 , 𝑃 𝑐 are the critical temperature and pressure respectively. ω is the acentric factor. 𝑀 𝑤 denotes the molar weight of the species. BIP stands for binary interaction parameter which is taken from the references [START_REF] García-Córdova | Vapor-liquid equilibrium data for the nitrogen + dodecane system at temperatures from (344 to 593) K and at pressures up to 60 MPa[END_REF], [START_REF] García-Sánchez | High-pressure (vapor + liquid) equilibria in the (nitrogen + n-heptane) system[END_REF] . The calculation results are compared with experimental data [START_REF] García-Córdova | Vapor-liquid equilibrium data for the nitrogen + dodecane system at temperatures from (344 to 593) K and at pressures up to 60 MPa[END_REF], [START_REF] García-Sánchez | High-pressure (vapor + liquid) equilibria in the (nitrogen + n-heptane) system[END_REF] as presented in Figure A.1. The phase behavior of this kind of mixture is rather difficult to model with a cubic equation of state due to the distinct molecular structures between N 2 and hydrocarbons. It can be noticed from the plots that the match with experimental data is good at most involved pressures but fails around the critical point. Since the only tunable parameter for the equation of state model is the binary interaction parameter (BIP), it is expected that using temperature dependent interaction parameters could improve the accuracy. It is noteworthy to mention that, for type III mixture [START_REF] Konynenburg | Critical Lines and Phase Equilibria in Binary Van Der Waals Mixtures[END_REF], the mixture critical pressure rises significantly as the concentration of nitrogen is increasing.

Whereas the mixture critical temperature reduces as more nitrogen exists in the mixture.

2) Recovery of sound speed inside the vapor dome

The non-convexity of cubic equation inside the vapor dome may lead to the loss of hyperbolicity in Euler system as discussed by many researchers [START_REF] Saurel | A general formulation for cavitating, boiling and evaporating flows[END_REF], [START_REF] Ma | An entropy-stable hybrid scheme for simulations of transcritical real-fluid flows[END_REF], [START_REF] Petitpas | Diffuse interface model for high speed cavitating underwater systems[END_REF], [START_REF] Menikoff | The Riemann problem for fluid flow of real materials[END_REF]. To cure this deficiency, one efficient strategy is to use the saturation line in two-phase state. Since an unphysical region exists inside the twophase region, the possibility of finding the saturation state between liquid and vapor phase always exists. In current study, the speed of sound is approximated using the Wood formula (Eq. (A.1) ) [START_REF]A Textbook of Sound[END_REF]. The gas volume fraction (α 𝑔 ) can be expressed as function of the vapor mole fraction (𝜓 𝑣 ) and specific volume (𝑣) in each phase, as written in Eq.(A.2) . Whereas the vapor mole fraction of a single component system can be easily obtained through the relation of volume, 𝜓 𝑣 = 𝑣-𝑣 𝑙 𝑣 𝑔 -𝑣 𝑙 (where 𝑣 𝑙 , 𝑣 𝑔 are the specific volume of liquid and vapor phase, respectively); for multicomponent systems, the specific volume of each phase is obtained by solving the PR EoS with each phase p composition (𝑥 𝑘 , 𝑦 𝑘 ). The calculation of phase composition and vapor mole fraction are through the VLE computation. The fundamental equation solved in the VLE computation is the Rachford-Rice equation [START_REF] Rachford | Procedure for use of electronic digital computers in calculating flash vaporization hydrocarbon equilibrium[END_REF], formulated as Eq.(A.3). With the relation

, the phase composition 𝑥 𝑘 , 𝑦 𝑘 are formulated as Eq. (A.4) derived from Eq. (A.3). Thus, the mixture system is reconstructed with each phase EoS which is also named as composite EoS [START_REF] Yi | A multicomponent real-fluid fully compressible fourequation model for two-phase flow with phase change[END_REF]. To better comprehend the reconstruction process, a VLE computation for n-C 12 H 26 -N 2 system is conducted with the initial temperature, pressure and molar fraction of N 2 equaling 593.5 K, 1 MPa and 0.1, respectively. A good agreement with experimental results is achieved for this condition as shown in Appendix A,