High pressure/high temperature multiphase simulations of dodecane injection to nitrogen: Application on ECN Spray-A

A R T I C L E I N F O


Introduction
High pressure/high temperature fuel injection and mixing is a topic of interest for many transportation applications, involving classical Internal Combustion Engines (ICE) in the automotive & transportation sector [1], but even extending to aerospace and rocket propulsion [2,3]. Fuel injection is a rather complex phenomenon on its own and challenging from a computational perspective [4], due to the vast disparity of temporal and spatial scales, ranging from the intact fuel jet just right at the injector exit, the formation of irregular liquid structures during primary atomization, to eventually the formation of finely atomised droplets at the end of the secondary atomization. Further complexities of studying fuel injection involve, interactions of internal flow phenomena with the jet (e.g. cavitation [5,6], combined with moving control needle [7]) and thermodynamic effects (phase transitions) that occur downstream during the mixing process [8]. Whereas over the recent years there have been investigations contributing to the understanding of the link between internal flow phenomena, primary atomization (see indicatively [5,6,[9][10][11]) and phase change at mild conditions (indicatively [12][13][14], at pressures and temperatures well below the critical point, hence classical atomization regime), thermodynamic effects at extreme conditions like those appearing in modern engines, exceeding the critical conditions of the fuel/oxidizer, are not thoroughly understood. This phenomenon applies mainly to Diesel fuel injection [15], which will be the focus of the present work, though it is also relevant to cryogenic fuel injection and combustion in rockets [2,16].
In terms of physical mechanisms, as temperature and pressure increase close to the critical point of a material, liquid and vapour phases are more difficult to discern, so there is smaller difference between saturation densities, the latent heat of vaporisation is reduced and surface tension effects diminish [17]. Among the first works to quantify and provide an explanation of mixing processes at high pressure/temperature conditions were by Dahms and Oefelein [18][19][20]. In these works, the authors justify the transition from the classical to the dense fluid regime due to the increase of intermolecular forces at the liquid and gaseous sides of the interface, leading to the thickening of the interface and the disappearance of the classical atomization pattern, thus instead of small droplets and ligaments, there is diffusive mixing governed by turbulence and molecular diffusion.
Naturally, the representation of materials at such conditions necessitates the use of a relevant thermodynamic model, termed as Equation of State (EoS), relating pressure, temperature and density. In order to capture such phase transitions, it is necessary to employ an EoS of at least third order (cubic polynomial); perhaps the earliest such example is the van-der-Waals EoS, which has been vastly superseded by more modern variants [21]. The situation becomes more complicated with the introduction of multiple materials and the formation of mixtures, with the inclusion of mixing rules and binary interaction parameters. An overview of the current state of the art on numerical methods and thermodynamic models applied on supercritical or transcritical injection is provided in the review work of Ma et al. [22]; recent computational works employing diffuse interface methods to study jet mixing at such conditions mainly employ thermodynamic closures based on the Peng-Robinson (PR), Redlich-Kwong (RK) or Soave-Redlich-Kwong (SRK) EoS [21]. Indicatively, Ihme et al. [23] studied the jet disintegration for the Spray-A case, identifying vapor and liquid penetration through mixture mass fraction, further expanding to autoignition characteristics under reacting conditions. Knudsen et al. [24] used an explicit, density based Finite Volume solver to describe the dense fluid mixing between dodecane and nitrogen, however without an explicit comparison of the liquid penetration of the jet; instead the authors quantified roughly a range of liquid penetration in terms of the mixture density variation. The most recent works, on the estimation of liquid penetration involved complex thermodynamic modelling and vapor-liquid-equilibrium (VLE) calculations under a purely Eulerian framework for ECN Spray A cases. Matheis et al. [25] performed a comprehensive analysis of the effect of different p/T conditions on spray evolution, although ignored in-nozzle flow. Yang et al. [26] performed similar parametric studies, while also including the injector, needle and upstream fuel system parts. All the previous methods employ the assumption of dense fluid approximation, which implies that all materials are tracked in an Eulerian manner, as a diffuse-interface, homogeneous mixture. For the sake of completeness, other researchers may employ spray modelling as Lagrangian particle tracking [27,28], Lagrangian-droplet-Eulerian-fluid with VLE calculations [29], Eulerian surface area density (Σ-Y) models [30][31][32], high resolution sharp interface methods, such as Volume of Fluid (VoF) [33], or hybrid techniques based on VoF for primary atomization and transitioning to Lagrangian particle tracking towards the secondary atomization [34]. However such approaches are not of interest in the present work, as they commonly employ classical assumptions of constant properties/ incompressible flow, droplet parcels and the relevant modelling methods (e.g. evaporation models).
Despite the extensive use of the PR and SRK EoS in modern numerical investigations, it is well known that these models suffer from deficiencies [35]. In particular, these cubic EoS are known to underpredict the density of the liquid phase and overpredict the speed of sound [36][37][38] (for a more comprehensive comparison between different thermodynamic models the interested reader is addressed to the work of Perez et al. [39]). Such discrepancies can lead to important errors, indicatively see Fig. 1 showing the relative error for the prediction of density between different models comparing to NIST RE-FPROP properties for dodecane [40]; note that the traditional cubic PR EoS may have an error of up to 10-15% near critical point and at high pressures. In fact, this deficiency has been recognised by other researchers, forcing them to use corrections to match the correct mass flow rate [25] and proposing alternative models for future work, such as the volume translated PR [41], or the generalized RK-PR. Nevertheless, volume translated PR still suffers from thermodynamic inconsistencies [42] whereas the generalized RK-PR, even if its predictions are improved comparing to the reference, still suffer from large deviations at very high p/T conditions (see Fig. 1).
In general, more advanced thermodynamic models do exist, though they have barely been explored in the frame of high p/T injection, mainly due to complexity and computational cost. In particular, a promising candidate is the Perturbed Chain -Statistical Associating Fluid Theory (PC-SAFT) model. This model provides a better agreement comparing to reference properties (see Fig. 1), while also being able to model the effect of multi-component mixtures [43] and realistic fuels [44][45][46][47][48] with minimum input of only three parameters per component. Extensions of the PC-SAFT model can provide information on transport properties, such as dynamic viscosity and thermal conductivity [49]. Additionally, it can provide information on the Liquid-Vapor Equilibrium and phase co-existence [37,50].
A particular drawback of advanced thermodynamic models is their complexity which limits the on-the-fly performance of the flow simulation algorithm. Hence many authors have resorted to thermodynamic evaluation of properties and the construction of a thermodynamic table at a precursor stage of the simulation over a range covering the conditions to be analysed. This approach is common for real-fluid materials with complex behaviour such as water, using cartesian Adaptive Mesh Refinement [51] tables, or even Arificial Neural Networks [52]. Similarly, for hydrocarbons unstructured tables have been used combined with bilinear interpolation [53]. Recently, thermodynamic interpolation methods have been extended to multi-species applications, using weighted inverse distance interpolation [54].
In the present paper it is aimed to employ such advanced thermodynamic models and apply them in a realistic simulation of the ECN Spray-A test case, including the internal injector geometry. The novel elements of the present paper involve: -tabulation of properties in a structured table of log 10 (p)-T, can be generalised to log 10 (p)-T-y form. This tabulation is accurate and very fast to search. Searching speed is independent of table resolution (see Appendix 1), hence the only limit is memory consumption of the property table. It is applicable to any thermodynamic model, be it high accuracy NIST REFPROP, PC-SAFT [57], generalised cubic equation of state [56] or traditional cubic equation of state.
-thermodynamic modelling; to the best of the author's knowledge, it is the first time that NIST properties, the most accurate to date, have been used for thermodynamic modelling using full p/T dependence. The major advantage of this approach is that there is no need for density/velocity corrections as is the case for Peng-Robinson or existing cubic EoS (see Appendix B of Ref. [25]). -the present investigation discusses the mechanisms of vorticity generation, which, to the authors' knowledge, has not been discussed in the past for compressible jets at such conditions/transitions.
The present paper is structured as follows: first, the mathematical model of the governing equations is presented, justifying the use of selected models. Then, a description of the thermodynamic property tabulation is provided, along with a discussion on the liquid identification. Further to that, a description of the case set-up, conditions and non-dimensional numbers are provided. Next, a description of results will be made, involving in-nozzle flow instances and liquid and vapour penetration. Finally, a short discussion follows, discussing mechanisms of vorticity generation, and finally a conclusion section.

Mathematical model
The simulations are carried out by assuming a diffuse interface approach, under a homogeneous mixture assumption under mechanical and thermal equilibrium; that is both materials involved share the same velocity, pressure and temperature fields [58]. Hence, the model consists of four Partial Differential Equations, plus the thermodynamic closure, which will be described in the next section.
The governing equations are provided below: -Mixture mass conservation: where ρ is the mixture density and u is the velocity vector field.
-Dodecane mass fraction, y C12 , transport equation: Any diffusion of mass fraction is ignored, since convective phenomena have been found much more dominant than diffusion processes [43]. Quadratic Upwind Interpolation For Convective Kinematics (QUICK [59]) has been used for the transport equation of dodecane.
-Mixture Momentum equation: where p stands for the pressure and τ corresponds to the stress tensor eff T , with μ eff the sum of laminar, μ, and turbulent, μ t , dynamic viscosity. The momentum equation is resolved using a Bounded Central Differencing [60] scheme, to minimise numerical diffusion, while also maintaining stability.
-The mixture energy equation: where E is the total energy, defined as internal energy plus the kinetic where h is the enthalpy, provided as a function of pressure and temperature (see section of Thermodynamic Properties). The total effective thermal conductivity, λ eff , is equal to the thermal conductivity, λ, being a function of thermodynamic conditions p, T and composition, y C12 , plus the turbulent contribution, as follows: where c p is the heat capacity of the mixture and Pr t is the turbulent Prandtl number, assumed equal to 0.85, based on the average value obtained from multiple experiments using different materials [61]. The energy equation is discretized using a Second Order Upwind scheme [59]. The necessary turbulent closure for the turbulent viscosity is provided using the Wall Adaptive Large Eddy (WALE) model, as it has combined good performance near walls [62]. Turbulent viscosity, μ t , is calculated as: where S ij is the rate of strain tensor and S ij d is the traceless symmetric part of the square of the strain of the velocity gradient tensor, i.e.: with = g ij u x i j and δ ij the Kronecker delta. The length scale, L s , is based on the filter size and the cell to wall distance, d wall , as follows: where the used model constants are: κ the von Karman constant, 0.41, and C w = 0.325. Radiation heating is omitted for the present cases. This can be justified on the basis of estimating the temperature increase of the liquid jet from the surrounding hot environment. By employing the Stefan-Boltzmann law, the radiation heat flux can be estimated as: where: -σ is the Stefan-Boltzmann constant, 5.6703 . 10 8 W/m 2 K 4 -T ∞ is the farfield temperature. As a worst case scenario, the maximum examined temperature of 1200 K will be used in the present estimation. -T jet is the temperature of the emerging jet, which is 363 K.
-A is the surface area of the jet; here it will be assumed that it is perfectly cylindrical, hence equal to π . D jet .
Based on these values, the heat flux towards the jet is~33 W per meter length of the jet. The equivalent temperature rise of dodecane can be estimated assuming the heating of the liquid cylindrical mass, as ρ C12 . π . R jet 2 . L~4.5 μg, assuming a density of 699 kg/m 3 , at the lowest examined chamber pressure. Given the heat capacity of dodecane is 2440 J/kg . K, the temperature rise due to radiation for the time scale of the injection (1.5 ms) is~4.5 K. Based on this temperature increase, and considering the highly turbulent mixing between the dodecane and nitrogen, convective heat transfer and mixing is assumed to be the dominant mechanism of heat transfer, hence radiation is ignored.
Time advancement is performed with a second order backward differencing scheme, using an implicit, pressure-based, fully coupled solver [63,64]. The implicit nature of the solver allows to use time steps much larger than those imposed by the restrictive acoustic Courant number; in the present study all simulations were carried out with a time step of 25 ns. The fully coupled nature of the solver provides better stability compared to segregated-type solvers (e.g. SIMPLE, PISO) [59]. The Fluent v19.1 software [65] has been used to carry out the present simulations, customized with external UDFs to incorporate the material properties and thermodynamic model using the functionality of User Defined Real Gas Model (UDRGM) [65]. Validation of the energy conservation of the mixing process of the jet and indicative shock tube cases are presented in Appendix 2.

Thermodynamic properties
Properties for dodecane are based on NIST REFPROP [40,66]. Thermodynamic properties are sampled along a range of 0-2500 bar and 280-2000 K. Note, that due to the limited applicability range of dodecane properties, extrapolation has been used to properties behaving in a reasonable manner (constant monotonicity). Other properties, such as thermal conductivity, which displayed erratic behaviour, especially at temperatures above 700 K, have been adjusted using thermodynamic models based on the PC-SAFT EoS.
Tables are stored in a structured data table of constant intervals as a function of the decimal logarithm (log 10 ) of pressure and temperature. As demonstrated in Fig. 2, using the decimal logarithm of pressure provides superior reconstruction of the properties (here shown for density), both at low and very high pressures; this is attributed to (i) the non-linear distribution of sampling points, refined towards low pressures where phase change is expected to happen and (ii) logarithmic functions are an accurate way of expressing liquid EoS, see Ref. [67]. The Average Absolute Deviation (AAD) at the tested range is 81.5%, 150% and 48% for the linear interpolation at temperatures 400, 600 and 800 K respectively, comparing to NIST REFPROP formulas. For the log 10 interpolation at the same range, AAD is 0.28%, 0.15%, 0.19% at temperatures 400, 600 and 800 K respectively, showing much superior accuracy. The actual table used had 75 pressure interpolation points, hence AAD against NIST REFPROP should be even lower. This procedure has the advantage that it can be done independently of the flow Fig. 2. Comparison between the linear (dashed lines) and log 10 (continuous lines) interpolation of dodecane density as a function of pressure, for different temperatures. Note that for this graph, only 50 interpolation points have been used along a pressure range of 100 Pa to 2500 bar. Temperatures span from subcritical (400 and 600 K) to supercritical (800 K). Note the excessive smoothing and large errors of the linear interpolation at subcritical pressures. P. Koukouvinis, et al. Fuel 275 (2020) 117871 solver execution, hence one can avoid the time-consuming calculation of thermodynamic properties on the fly. Furthermore, this procedure can be generalised to any equation of state or mixing rules. The ambient gas is modelled as ideal gas, with variable heat capacity, conductivity and viscosity based on NIST REFPROP [66], as the deviation from the real-gas, both in terms of density and speed of sound, is rather small, see Appendix 3 for indicative comparisons over the pressure and temperature ranges examined. Mixture properties are assumed to be mass or volume weighted averages of the properties of the individual components involved, assuming thermodynamic equilibrium, hence metastable states and loss of hyperbolicity is avoided [68]. Speed of sound for mixtures is calculated based on the Wallis or Wood speed of sound formula [68]. The PC-SAFT model has also been used to obtain tables of pressure-temperature-composition of the materials involved, which can help identify liquid volume fraction, as will be described in the next subsection. The interested reader is addressed to the Appendix 4 for an outline of the modelling used for PC-SAFT calculations. Indicative tables are provided as supplementary data, see Appendix 5 for the table format.

Liquid identification
The identification of the liquid length is not straightforward to derive from dense fluid simulations, at elevated pressure/temperature conditions [69]. In the literature, researchers have used as an indication the mixture fraction [15,23,70], however this choice of the liquid presence is somewhat arbitrary and often is associated with parametric investigations based on various mixture fraction values. It should be highlighted that the sensitivity of liquid length on the mixture fraction value is quite strong; indicatively, the liquid length can vary by 100% over a mixture fraction range of 0.4-0.95. On the other hand, relying on single-phase thermodynamic parameters in the form of a simple criterion e.g. comparing local temperature, T, to the critical temperature, T c (in the form if T < T c , then liquid), is a rather crude assumption and overestimates the liquid length, as will be shown later.
Mixture thermodynamic modelling can provide information through Vapor-Liquid-Equilibrium (VLE) calculations. Due to the lack of documented applicability of NIST REFPROP properties for mixtures with largely different volatilities, as is the case here, the capability of the PC-SAFT model will be exploited in identifying VLE instead. The process is as follows: based on the PC-SAFT p-T-y diagram (details on derivation are provided in Appendix 4), for a given pressure, p, the composition-temperature iso-line, here denoted as T L (y), corresponding to the desired liquid fraction (here 0.15%, this value is based on previous experimental work [71] and has been used in similar CFD studies in the past [25]) is identified. If the local temperature of a computational cell, T, at the cell's composition, y, is lower than T L (y), then it is assumed that the state of fluid is liquid/gas mixture, else it is vapour/ gas mixture. The procedure is displayed schematically in Fig. 3. Note that, alternative criteria based on the path-integrated, projected liquid volume fraction have also been examined, though differences between the methods are less than 0.5% (see Appendix 6).

Case set-up
The computational mesh consists of the Spray A injector geometry #210675 [73] and the downstream injection volume. It is a single hole injector, with nominal orifice diameter of 90 μm, orifice length~1 mm and k-factor of 1.5 (see Fig. 4a). The injection volume extends from the exit of the injector orifice, up to 40 mm downstream and a radius of 10 mm as to reduce any interference from boundaries. To achieve high resolution in areas of interest, while also maintaining computational efficiency, an unstructured, telescopic refinement mesh was used, with five distinct refinement zones and near wall refinement (a view of the refinement zones is in Fig. 4b). The resolution in these refinement zones is: . The maximum mesh size at the farfield is 0.52 mm. Near wall boundary layers are used inside the orifice, with a near wall resolution of 0.5 μm, resulting to a max. y +~10. The total cell count is 9.5million cells. A complete structure of the mesh is shown in Fig. 5, demonstrating the transitions between the refinement zones.
The mesh sizing described in Fig. 5 and the time step used result to a convective Courant number of~5 in the orifice and the refinement zone #1, dropping to less than 0.25 at the rest zones for the time step used. Acoustic Courant number is less than 10; this is a beneficial aspect of using the implicit pressure-based solver, as there is no stability constraint due to acoustic waves [74,75].
The conditions analysed correspond to those examined by ECN [76] and are summarized in the table below: Fig. 3. Procedure for determining the liquid volume fraction: starting from the PC-SAFT pressure-composition-temperature, p-T-y, diagram: (a) at a given pressure, the composition-temperature y-T diagram is found; here indicated as a dark red slice. (b) the isoline of desired liquid fraction threshold (here 0.15%, based on Refs. [25,72]) separates the composition-temperature diagram to two regions; anything below the thick, black isoline, belonging to region 1, is liquid dominant, hence will be identified as liquid. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) P. Koukouvinis, et al. Fuel 275 (2020) 117871 Boundary conditions at the inlet are based on the mass flow rate, at fixed temperature (see Table 1). The fluid at the inlet is prescribed as pure dodecane (no mixing with non-condensable gas). Mass flow profile is based on the online tool provided by CMT [77], using 1500 bar upstream pressure, 20.4-80 bar downstream pressure depending on the case, outlet diameter of 89.4 μm, density of 786.7 kg/m 3 (corresponding to 1500 bar and 363 K, based on the properties used), discharge coefficient 0.9 and duration 1.5 ms. The resulting profile is shown in the figure below ( Fig. 6), max. variation between the mass flow rates for different backpressures is~2%. Fixed pressure is imposed at the farfield and zero gradient for all other quantities (velocity, mass fraction). Nozzle walls are considered adiabatic, as the injector is kept   at fixed temperature same as the fuel. Initially, the sac volume and 95% of the orifice length were patched to pure liquid dodecane at farfield pressure and fuel temperature (363 K). Based on the conditions above, the liquid dodecane is injected at max. pressure of~500 bar and temperature of 363 K, at downstream pressures of 20.4 bar−80 bar and the jet velocity is in the order of 600 m/s. At the aforementioned operating pressure/temperature range, ρ C12 = 699.7-786.9 kg/m 3 and μ C12 = 0.578-1.885 mPa.s, hence the maximum Reynolds number is, considering the injection conditions and properties of dodecane: The Taylor length scale for this Reynolds number range is 1.1-2 μm (2 times the cells size at the near wall region of the hole, or 1/3rd of the cells located at the core of the hole region), which is comparable to the finest resolution level used inside the orifice and at the vicinity of the exit of the injector.
The Weber number of the jet requires information on the surface tension coefficient between dodecane and nitrogen. In general, this coefficient is strongly dependent on the liquid/gas combination, though here it will be assumed similar to that of diesel/nitrogen for which measurements exist [78]. For the conditions examined, the surface tension coefficient ranges from 0.0166 N/m (at 80 bar and 363 K) to 0.019 N/m (at 20.4 bar and 363 K); note these calculations are based on the assumption that dodecane has just exited the orifice, hence it is still cold, at the upstream temperature of 363 K. As dodecane heats up, surface tension will diminish strongly, hence the surface tension coefficient values provided are the maximum possible.
Hence the minimum Weber number of the jet is (based on the maximum surface tension coefficient): Given the fact that Weber number expresses the ratio of inertial to surface tension forces, large values indicate dominance of inertial effects [79], at least for the large jet scales. Hence, it becomes apparent that surface tension plays a rather marginal role when considering the whole jet, justifying its omission from the governing equations. When considering smaller structures, the role of surface tension is not well understood. It has been found that as pressure increases, evidence of surface tension influence becomes smaller; indicatively, for ambient pressures of 20 bar and temperatures of 900 K, droplets or ligaments were still detectable [80], though their number was much smaller than lower pressures. Further experiments at higher pressures showed a diminishing effect of surface tension and transition to diffusive mixing. After careful observations using long range microscopy of droplet lifetime, a map has been derived [81], indicating regimes of classical evaporation, where surface tension is still present, transitional mixing and diffusive mixing, where surface tension may exist only for a short duration of the droplet's lifetime. It is highlighted that all the previous discussion is relevant after the end of injection, at characteristic droplet sizes of 10-20 μm and velocities of 2-10 m/s (We~0.1-5), which is not the interest of the present study.

Internal flow
Internal flow depends relatively weakly on the downstream conditions, for the range of pressures used here, hence results will be only shown for the 900 K and 60 bar case. The velocity distribution inside the orifice is shown in Fig. 7. After 0.1 ms the flow practically stabilises; the injected dodecane has an average velocity magnitude of~630 m/s and average root mean square of velocity fluctuation of~4 m/s, with maximum of 70 m/s near the walls. As shown in the figure, flow velocity asymmetry is observable; near the entrance of the orifice, velocity is slightly higher at the lower side, due to the sharper turn that the flow has to follow. This causes a small detachment region locally, resulting to a thicker boundary layer at that side of the orifice.
An indicative instance, showing the temperature distribution inside the injector is shown in Fig. 8; here the effects of viscous heating and liquid cooling due to depressurization are visible. The strong shear along the orifice walls causes heating, due to the viscous effects (the last term of Eq. (4)). Locally, near wall temperatures may reach~500 K. On the other hand, isentropic liquid expansion imposes cooling of the liquid; an isentropic expansion from the inlet pressure and temperature  to the orifice exit pressure would result to a temperature of~341 K (temperature drop of~22 K) which is similar to the minimum value found in the orifice exit as shown in Fig. 8. In the whole cross-section, the average outlet temperature is reduced by~13 K (average outlet temperature of~350 K at max. flow rate). The observed temperature reduction corresponds to a velocity coefficient, c v , of~0.92, similar to that found in previous studies [82,69,70], whereas the discharge coefficient is~0.91. Note that, due to the injector orifice eccentric placement, the one side of the orifice shows more heating, see Fig. 8.

ECN Spray A − 900 K, 60 bar and comparison with 700 K, 46 bar
An indicative comparison of the spray evolution for downstream conditions of 60 bar and 900 K is shown in Fig. 9. The cyan iso-line represents the liquid region, with a liquid volume fraction greater than 0.15%, whereas the greyscale colouring represents the mixture fraction of dodecane to obtain a representation similar to the experiment. The results show a qualitatively good agreement with similar propagation speed for both liquid and vapour. The strong shear forces and density gradients lead to instabilities which break the jet forming mushroomtype structures, see also Fig. 10. The computational results are similar to those obtained by other numerical investigations, e.g. see Refs. [23,25]. Fig. 10 shows details of the start-up of the jet. As expected, during the beginning of the injection, vortex roll-up occurs, leading to the formation of a hemispherical cap at the tip of the jet, see the instance at Fig. 8. Temperature distribution inside the injector. Near the wall viscous heating can be observed, whereas in the core of the flow temperature drop due to liquid depressurization is observed. The black isolines indicates the region of liquid, cooled below the inlet temperature (< 363 K), according to the denoted value. Dimensions in mm.

Fig. 9.
Indicative results of the liquid and vapor penetration at the start up of the jet for dodecane injection to 900 K, 60 bar nitrogen (left), comparing to experimental data [76] (right). The thick cyan iso-line represents region where the liquid volume fraction is higher than 0.15%; note that the average liquid length is 10.8 mm for the injector #210675 [76]. Colouring is according to dodecane mass fraction. Dimensions in mm. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) 10 μs. Further, as the jet propagates, it is disintegrated due to turbulent mixing (Kelvin-Helmholtz and Raylegh-Taylor instabilities), leading to its widening from 90 μm at the orifice exit to 400 μm in a distance of 4 mm. Note also the slight asymmetry in the jet formation, clearly visible near the vicinity of the orifice exit, due to the orifice eccentric placement. Around the jet surface, mushroom-like interfacial instabilities develop, due to density gradients (Rayleigh-Taylor instabilities).
The asymmetry of the jet induced by the orifice geometry can be clearly shown in Fig. 11. The deformation of the jet surface is shown through the isosurface of dodecane mass fraction for a value of 10%. Also, slices along the jet cross-section are provided, to demonstrate the progressive deformation of the spray surface. It can be seen that, downstream the exit from the orifice, the initially circular jet grows circumferential prolongations, mainly at opposing sides, see Fig. 11; the main ones are indicated along the jet direction as dashed red paths and numbers from 1 to 4. These are the result of tangential, to the jet crosssection, velocity components emanating from surface roughness/irregularities of the injector orifice and it has been demonstrated that they are strongly dependent on surface representation [33]. Whereas in the present study the jet cannot be resolved at the same level of detail as previous studies with interface capturing schemes [33], due to lower resolution (both temporal and spatial) and modelling assumptions (dense fluid, omission of surface tension), there are similarities with previous investigations in respect of jet characteristics. Indicatively, the level of tangential velocity magnitude at the exit of the orifice, being in the order of 25 m/s, and the arrangement of striations in opposing groups, are similar to the study of Agarwal et al. [33], who used Volume Fig. 10. Details of the start-up of the jet, ECN Spray-A, injection to 60 bar, 900 K. Left: vortical structures around the jet, indicated using the second invariant of the velocity gradient (also known as q-criterion, value 5 · 10 10 s −2 ). Right: Dodecane mass fraction isosurface, y C12 = 10%. Note the formation of a vortex ring due to vortex roll-up at the beginning of the jet and the formation of mushroom-like structures (indicated with arrows), outcome of Rayleigh-Taylor instabilities.
of Fluid method, albeit omitting any thermodynamic effects.
A quantitative comparison of the spray penetration is shown in Fig. 12. Here, experimental results of the vapour penetration from the #210677 injector and liquid penetration from the #210675 injector (obtained from the relevant web pages of the ECN website [83,76], respectively) are compared against the numerical results showing a good agreement. The vapour penetration has been estimated using a dodecane mass fraction value of 0.1%, as it is the recommendation from ECN [84]. Whereas vapor penetration is somewhat insensitive to the mixture fraction value used [23], using higher values of e.g. 10%, can lead to relatively slower penetration estimation by~3-4%, especially at later stages, when the jet becomes diluted and in areas of coarse resolution, see Fig. 12. As observed in the graph, vapour and liquid propagate at the same velocity, until a time of~0.04 ms, at a distance of~8 mm. Also, as a comparison, the injection at the lower temperature of 700 K and 46 bar pressure is shown in Fig. 13; both vapor and liquid penetration at the two conditions are similar, however vapor penetration at 700 K and 46 bar is slightly faster, since the ambient gas density is slightly lower (by~1%). Liquid penetration at the lower temperature extends further downstream, due to the lower heating rate, thus inducing slower vaporisation.
Finally, Fig. 14 shows the comparison of the radial distribution of dodecane mass fraction, between simulation and experiment, for downstream conditions of 900 K and 60 bar pressure. Experimental results are shown for two different radial directions; note the slight asymmetry at each direction. The error bars of the experimental results are presented using 95% of confidence level ( ± 2σ), as in the Ref. [85], and are shown only for the one experimental dataset (the blue points), for clarity. The numerical results are obtained by time-averaging instantaneous mass fraction distribution, after the spray has established, starting at 0.16 ms ASOI and using samples over 50 μs duration (2000 samples) and then performing circumferential averaging to obtain the radial distribution.

ECN Spray A − 1200 K, 80 bar
A similar comparison of the spray evolution for downstream conditions of 80 bar and 1200 K is shown in Fig. 15. Due to the higher ambient temperature, the liquid core disintegrates much faster, at around~5-7 mm, leaving unconnected liquid blobs that may reach occasionally 10 mm. Note that such structures may have a size of 100 μm or less, which is comparable to the pixel size of the imaging system [76]; most likely this could affect the detection in the experiment. The jet propagation is similar in both cases, though the vapour penetration is somewhat over predicted towards the end of the provided instances.
A quantitative comparison of the spray penetration is shown in Fig. 16. Note that for the liquid penetration, the maximum coordinate of the identified liquid isosurface has been used, hence smaller than the experimental imaging system resolution structures have been considered in the liquid length calculation. Nevertheless, the predicted Fig. 11. Demonstration of jet surface striations, generated at the exit of the orifice. Slices are placed every 0.1 mm, starting from the exit of the orifice. Note that the vectors show only direction, their magnitude is constant; only one every 25 vectors is shown for clarity. Instead, tangential velocity magnitude of dodecane is shown as the slice contour coloring, only for dodecane mass fractions of 10% or more. Fig. 12. Indicative results of the jet liquid and vapor penetration. CFD stands for the present study; L is the liquid penetration (volume fraction > 0.15%), V is the vapour penetration (for 10% and 0.1% mass fraction), T c is the propagation of temperature isosurface at the dodecane critical temperature. Ref-liq corresponds to published numerical data [25]. Vapour -Exp. correspond to data from the #210677 injector [83] and Liquid -Exp to diffused back-illumination (DBI) of the #210675 injector [76]. The dashed black line corresponds to the average liquid length of 10.8 mm. trend is similar close to the average liquid length reported in the ECN website [76] and published data from previous numerical investigations employing LVE calculations [25]. Concerning the time mismatch between the simulation and the experimental data, it seems that it is related to timing anomalies of the injection system as both the present simulation and the reference simulation data have a very similar time shift. Nevertheless, the average liquid length is correctly predicted.

Comparison between 900 K, 60 bar and 900 K, 20.4 bar -effect of ambient density
The purpose of this comparison is to investigate the spray penetration under the same temperature conditions, but different ambient gas density, as imposed through the ambient pressure. To be specific, the cases compared here involve a density of 22.5 kg/m 3 and 7.6 kg/m 3 , see also Table I, hence it is expected that the lower ambient density should present a reduced spray drag and faster spray penetration. Indeed, as shown in Fig. 17, the reduced ambient density of 7.6 kg/m 3 shows a much faster jet penetration, both in terms of liquid and vapour. The same observation can be made by comparing the spray propagation in respect to time, as shown in Fig. 18; indeed, from 40 μs and onwards, the vapour spray has a faster propagation by~40-50% at the low density/low pressure condition.

Discussion
As shown in the previous figures, the break-up process of the jet shows quite intricate and complex patterns of mixing. The compressed liquid jet disintegrates as it becomes heated, forming the characteristic spray cone. A general observation from the present results is that the simulated liquid core thickness appears smaller than the one observed in the experiments; this is a common pattern in previous similar numerical investigations, see e.g. [25] or [15]. It is likely that this observation is due to simplifications of the surface representation of the orifice geometry; in fact, it has been shown that differences in the resolution of Spray-A geometry can lead to dramatic changes of the breakup pattern of the spray, see Ref. [33]. A further reason, relevant during the start-up, is the associated uncertainties of initial conditions inside the injector and rate of injection profile. In the present study, the injector (sac and 95% of the orifice length) was assumed to be initially  Fig. 15. Indicative results of the liquid and vapour penetration at the start up of the jet for dodecane injection to 1200 K, 80 bar nitrogen, comparing to experimental data [76]. The thick cyan iso-line represents region where the liquid volume fraction is higher than 0.15%; note that the average liquid length is 7.4 mm. Colouring is according to the dodecane mass fraction. Dimensions in mm. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) P. Koukouvinis, et al. Fuel 275 (2020) 117871 full of liquid, however this does not necessarily reflect the actual conditions in an injector before injection. In particular, there are experimental results demonstrating that sac and orifice may be partially occupied by gas bubbles, sucked in from the previous injection cycle, or by movement of the needle valve on the same cycle [86]. The exact amount of gas in the sac volume can lead to deviations of the spray penetration timing, whereas the spatial distribution of such gaseous pockets can lead to shot-to-shot variations and increased spray dispersion/penetration [87]. Another potential reason for the mismatch is the omission of turbulent fluctuations, originating upstream from the injector. It is reminded here, that the injector is a complex electro-hydraulic device, that includes a flow control needle valve, just upstream of the simulated section. This valve, as well as the upstream injector components, induce flow disturbances, identified in recent a study [88] as Schlieren-like structures, due to the density-gradient changes observed through high-speed imaging. Such flow disturbances will undoubtedly have an effect on the secondary flows in the sac and, consequently, the spray. Finally, even though the mass profile used aims to replicate the needle opening, in fact the actual presence of the needle induces viscous heating, as it limits the discharge coefficient at low lift operation; indicatively, a temperature increase of 60-80 K is expected at lifts below 100 μm, depending on injection pressure [89]. Hence, it is reasonable to expect that such heating can contribute to mismatch during early transients effects, especially at the beginning of injection, explaining the slight lag observed in Figs. 12, 13 and 16. Nevertheless, the aim of the present study is to demonstrate the capability of the current model in explaining fundamental features and spray patterns. Despite the simplifications concerning upstream effects (such as those mentioned in the previous paragraph), it is of interest to Fig. 16. Indicative results of the jet liquid and vapour penetration. CFD stands for the present study; L is the liquid penetration (volume fraction > 0.15%), V is the vapour penetration (for 10% and 0.1% mass fraction), T c is the propagation of temperature isosurface at the dodecane critical temperature. Ref-liq. corresponds to published numerical data [25]. Liquid -Exp. corresponds to diffused back-illumination (DBI) measurements of the 210,675 injector [76]. The dashed black line corresponds to the average liquid length of 7.4 mm reported in the experiments.  P. Koukouvinis, et al. Fuel 275 (2020) 117871 provide an explanation on the mechanisms of jet break-up during interaction with ambient gas. This can be done by examining sources of vorticity generation. In particular, the vorticity evolution equation, in the absence of external body forces (e.g. gravity, which is omitted here), can be expressed as follows [90]: The first term in the right hand side represents vortex stretching/ tilting, the second term represents vortex dilation, the third term represents the baroclinic torque and the last term expresses vortex diffusion due to viscous stresses, as defined in Eq. (3). Vortex stretching or tilting is due to the effect of velocity gradient on vorticity; it is a crucial mechanism in the generation of complex vortical structures and is considered responsible for the kinetic energy cascade process in turbulence [91]. Vortex dilation is due to the volumetric expansion/contraction (the velocity divergence), which describes how fluid compressibility affects the vorticity. The baroclinic torque occurs due to the different alignment of density and pressure gradients and is responsible for the formation of Rayleigh-Taylor instabilities, like those observed in Fig. 10.
In Fig. 19a comparison of the strength of different factors affecting vorticity evolution is shown, on a slice along the length of the jet, for downstream conditions of 900 K and 60 bar; results are indicative for the rest cases and time instances, as well. As observed from the magnitude of the terms involved, viscous stresses have the lowest contribution in vorticity; high values are mainly located just outside the orifice exit, where velocity gradients are very strong. Further downstream, the contribution of viscous stresses diminishes, as the jet decelerates. The strength of baroclinic torque and vortex dilation terms is comparable, though the baroclinic torque is slightly stronger; both have a magnitude in the order of 10 11 1/s 2 , i.e. an order of magnitude stronger than viscous effects. Finally, the strongest contribution in vorticity comes from the vortex stretching/tilting term, which leads to the formation of a cascade of smaller vortical structures. These observations denote the relatively weak influence of viscosity, which however can contribute at the start-up of the jet (t < 5 μs), in the initial vortex roll-up or the formation of Kelvin-Helmholtz rolling instabilities. At the maximum flow rate, the break-up process is dominated by vortex stretching/tilting, interfacial instabilities, such as Rayleigh-Taylor instabilities, and compressibility effects (i.e. vortex dilation), creating hairpin-like structures. Similar observations on the nature of break-up mechanisms at comparable density ratios (ρ min / ρ max~0 .05) and moderate Reynolds/high Weber numbers have been found in previous numerical investigations on idealized, incompressible jets [92].

Conclusion
The simulations presented so far, demonstrate the capacity of the discussed methodology, employing complex thermodynamics and accurate property libraries to capture details of the flow. The liquid length was found to be similar between simulations and experiments for the examined cases. The numerical results provide very fine details of the jet formation, multiphase phenomena and flow instabilities. Also, the analysis performed, provides indications on the driving mechanisms behind the jet break-up and association to Rayleigh-Taylor instabilities, considering the relative weak importance of viscosity comparing to the influence of baroclinic torque and compressibility.
In terms of thermodynamic model implementation, the tabulation technique employed decouples the evaluation of thermodynamic properties during execution of the flow solver, simplifying portability of models and greatly speeding-up execution of the simulation, which does not need to be done concurrently with the thermodynamic evaluation. This is especially important with complex thermodynamic models, as their computational cost per cell evaluation can be considerable and prohibitive for practical applications. NIST REFPROP formulations are the most accurate available so far, however these are limited only to a limited selection of materials and material combinations. Moreover, as demonstrated in Appendix 1, their applicability range is limited and extrapolation is not always straightforward or could be totally unreliable. On the other hand, the PC-SAFT model is capable of providing accurate information on the liquid penetration, through the pressure-temperature-composition diagram. Moreover, it can be directly applied to different components and multi-component mixtures, as databases for the needed coefficients already exist [93], without the tedious calibration needed for Helmholtz based EoS, like those in NIST REFPROP. It is highlighted that the formulations used by NIST, while unprecedented in terms of accuracy, may have more than 40 coefficients, depending on the material, specifically optimized to match the experimentally measured properties. Aim of further investigations is to exploit this capability of PC-SAFT to examine the effect of different components and component/gas combinations.
The focus of the present work was mainly placed at high p/T conditions, as the liquid/vapor separation occurs early, in a short distance from the orifice exit, hence relatively small simulation time is needed. Despite the great detail such simulations can provide, their computational cost is rather expensive, considering that to obtain the current results more than 50000CPU-hours were needed for a single simulation (60 CPU cores for each case), thus rendering parametric investigations and optimisation difficult on the industrial level. Nevertheless, the indications obtained from the present study, concerning the nature of instabilities that govern break-up, can lead to the development of RANS models tailored towards the description of relevant mechanisms.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix 1. Dodecane properties and tabulation
Fuel properties are described using structured log 10 (p)-T tables, based on REFPROP NIST (for the time being for dodecane, there is the possibility to use other sources or PC-SAFT models for other fuel components/mixtures). The constant resolution, structured format of the table enables very fast search irrespectively of the table resolution, as the table index over which to interpolate can be found through a simple division, as: where x represents one of the independent variables to be used for the table search and dx the constant discretization used for that specific independent variable. Linear interpolation is then performed over the [indx, indx + 1] interval. The choice behind using linear interpolation, instead of higher order alternatives, is the preserved monotonicity of the interpolation, which may be important in any noise produced due to evaluated Fig. 20. Tabulated properties according to NIST REFPROP; indicatively here density (left) and speed of sound (right). All properties (e.g. density, enthalpy, speed of sound, viscosity, etc) are expressed in a similar way, as a structured grid of the decimal logarithm of pressure, log 10 (p), and temperature, T. thermodynamic properties. Resolution can be increased for improved accuracy in sharp transitions without any impact on the execution time. The log 10 (p) format and interpolation provides good accuracy both at low and high pressures. An indicative representation of such a table is provided in Fig. 20.
Whereas most of thermodynamic properties behave reasonably over the pressure/temperature range examined (temperature range 280-2000 K, pressure range 0-2500 bar), thermal conductivity behaved non-monotonically at high temperatures; this is most likely related to the limited applicability range provided by NIST (up to~700 K). Since the conditions examined involve temperatures beyond 700 K, it was necessary to amend the thermal conductivity formulation. This was done partly by linearly extrapolating the thermal conductivity and by weighted average contribution of the PC-SAFT predictions at higher temperatures; NIST is used at temperatures below 700 K, linear blending between NIST and PC-SAFT between 700 and 900 K, then switching entirely to PC-SAFT at temperatures higher than 900 K. In this way the continuity and smoothness of variation of thermal conductivity is maintained in all the temperature range. An indicative comparison of the initial and modified thermal conductivity is shown below Fig. 21.

Appendix 2. Validation for planar cases -shock tube and jet
A simplified case is presented here to serve as a comparison in terms of the adiabatic mixing. The case resembles the ECN Spray-A case, but is simulated as a planar jet, since in this way much less computational resources are needed, hence the simulation is tractable on a simple desktop computer. Comparatively, the cell count is 1 M cells for planar 2D, instead of~10 M cells for the actual 3D case presented in the main investigation.
The planar jet has a width of 100 μm, injected in a computational domain of length 12 mm and width 6 mm. The planar jet is dodecane at 363 K, injected at 600 m/s in a nitrogen atmosphere at 110 bar and 973 K. The jet velocity is slightly less than the speed of sound of the ambient gas (M~0.94 in respect to the stagnant gas), hence the flow is transonic. In the indicative images shown below, shock waves can be seen emerging at the start-up of the jet, represented as shading according to the pressure gradient magnitude.
Indicative instances of the jet evolution are shown below, in Fig. 22. As shown below, at the start-up of the jet a shock wave is formed due to the sudden introduction of dodecane in the computational domain. The shock wave expands forming a bow shock, note also the temperature increase in the vicinity of the bow shock due to rapid gas compression; temperature locally can reach 1100 K (ΔT > 127 K). Behind the shock, the dodecane jet propagates and starts to disintegrate, due to the immense shear instabilities. These lead to the formation of mushroom-like structures, that cause further mixing and emission of secondary shockwaves, due to the complex reflections of the primary bow shock on the side boundaries and jet structures.
Concerning the jet heating, it is important to note that the path followed during dodecane/nitrogen mixing is close to the isobaric, adiabatic line, see Fig. 23. In the figure, the calculated states from the flow solver are shown as red points, scattered across the y C12 -T plane; scatter is expected to exist, as the pressure is not uniform at the whole computational domain, instead pressure wave transients can cause localised heating/cooling (see also the bow shock heating mentioned in the previous paragraph). Nevertheless, the general trend of the isobaric, adiabatic line is followed; note that isobaric, adiabatic lines are calculated using real-fluid NIST REFPROP or PC-SAFT properties, showing that the ambient gas simplification used in the present study does not affect much the resulting heating of the jet. Also, as a reference, the isobaric, isochoric mixing line [94] is also shown -this is the process followed by quasi-conservative methodologies [22], prone to energy conservation errors and mixing temperature over-prediction. Isobaric, adiabatic and isobaric, isochoric curves have been calculated using NIST REFPROP and PC-SAFT real-fluid mixing tables provided in Appendix 5. A schematic illustration of the procedure is shown in Fig. 24, using NIST REFPROP properties; the isobaric, adiabatic mixing curve passes from the intersection of mass fraction iso-lines with respective isenthalpic lines, calculated as h = (1 − y C12 )h N2 + y C12 h C12 . A similar calculation can be done for isobaric, isochoric lines; the process is identical, with the exception of using the intersection of mass fraction iso-lines with respective isochoric lines as v = (1 − y C12 )v N2 + y C12 v C12 , where v represents the specific volume. Further to the planar jet case, two shock tube cases are presented, at the conditions of the jet above, to demonstrate that the heating observed at the jet front is not a numerical artefact, but rather a physical effect of the supersonic gas compression. The shock tube is a one-dimensional case, involving the evolution of the flow in the presence of an initial discontinuity, examined here for the following cases: (1) with a large pressure discontinuity across the material interface: The shock tube spans from −2 to 2 m and the initial discontinuity is located at x = 0 m. It is resolved with 1000 cells and will be examined at 1 ms. The solution for these cases can be found utilising the exact Riemann problem solver for arbitrary EoS [53], combined with multi-material handling [95]. However, deriving the exact solution is rather cumbersome, given that it will involve the solution of two exact Riemann problems for tabular EoS, coupled across the material interface. Instead, an approximation of the exact Riemann problem solution will be provided as a reference, based on the characteristics-based approximate Riemann solver [96], to give an indication of the validity of the CFD predictions. For case (1), the predicted approximate solution is: u*~100 m/s, p*~134 bar, ρ* L~7 41.2 kg/m 3 , ρ* R~4 2.6 kg/m 3 , resulting to temperatures of T* L~3 18 K (dodecane side -cooling due to expansion), T* R~1 060 K (nitrogen side -heating due to compression). As shown in Fig. 25, the CFD prediction results to similar patterns as the approximate Riemann problem solution, with a clear formation of a shockwave at the gas side and rarefaction wave at the liquid side. The only exception is a slight overshoot of fluid velocity near the material interface, but in general the wave pattern is properly represented without serious oscillations.
For case (2), the predicted approximate solution is: u*~588 m/s, p*~251 bar, ρ* L~7 13.4 kg/m 3 , ρ* R~7 0.3 kg/m 3 , implying temperatures of T* L = 374 K (dodecane side), T* R = 1202 K (nitrogen side), both resulting as heating from compression. As shown in Fig. 26, the CFD prediction results to similar patterns as the approximate Riemann problem solution, with a clear formation of a shockwave towards both the gas and the liquid sides, due to the formation of a compressed state in between the two materials. Similar to before, a slight overshoot near the material interface is found both for pressure and velocity fields, but in general the wave pattern is properly represented without serious oscillations. In Fig. 26 (1). The dashed line indicates the material interface between C 12 (located at the left) and N 2 (located at the right), advected at the speed of the contact discontinuity wave. P. Koukouvinis, et al. Fuel 275 (2020)

Appendix 3. Gas properties
Below, a comparison between the real fluid properties of nitrogen, based on NIST REFPROP [66], and the ideal gas assumption are presented below. As shown, over the range of 50-125 bar and 300-1200 K the deviation is not that significant, with maximum error of~6% for density and 9% for speed of sound, mainly at low temperatures, with Average Absolute Deviation (AAD) of 2.6% for both. The AAD is defined as:   Fig. 26. Shock-tube results for case (2). The dashed line indicates the material interface between C 12 (located at the left) and N 2 (located at the right), advected at the speed of the contact discontinuity wave.

Appendix 4. Perturbed Chain -Statistical Associating fluid Theory (PC-SAFT) properties and Vapor-Liquid-Equilibrium (VLE) calculation
The PC-SAFT EoS [57] is a theoretically derived model, based on perturbation theory [97][98][99][100], that splits the intermolecular potential energy of the fluid into a reference term accounting for repulsive interactions and a perturbation term accounting for attractive interactions. The reference fluid is composed of spherical segments comprising a hard sphere fluid that then forms molecular chains to create the hard-chain fluid. The attractive interactions, perturbations to the reference system, are accounted for with the dispersion term. Intermolecular interaction terms accounting for segment self-or cross-associations are ignored. Hence, each component is characterized by three pure component parameters, which are a temperature-independent segment diameter, σ, a segment interaction energy, ε, and a number of segments per molecule, m. For dodecane and nitrogen these parameters are: where R is the universal gas constant, T is the temperature. Once the different contributions to the residual molar Helmholtz free energy have been defined, every other thermodynamic property can be calculated by its derivatives, as the Helmholtz free energy is a thermodynamic potential. In the case of thermophysical properties, i.e. thermal conductivity and viscosity, correlations need to be used as shown below.

Calculation of thermal conductivity
In this study, the residual entropy scaling technique is used for the calculation of the thermal conductivity [ In their original work, the authors applied the entropy scaling technique only for pure components, therefore an empirical mixing rule have been applied in this study to the coefficients and the reference thermal conductivity by the following formula: where X is any of the coefficients [A λ , B λ , C λ , D λ ] or λ ref . In case the conditions for the mixture falls within the Vapor-Liquid Equilibrium regime, the mixing rule is a simple weighted average using the vapor volume fraction, α g , and the thermal conductivities for the liquid and gas phases, i.e. λ l and λ g respectively. The present method is found to perform better than Chung's method [49] with an overall lower deviation in the prediction of thermal conductivity when comparing against NIST REFPROP at a range of 280-700 K; indicatively the present method has an average deviation from NIST of 10.7%, 8.4% and 4.8% at 1, 10, 100 bar respectively, whereas Chung's method 15.8%, 16.5% and 27% respectively.

Calculation of viscosity
The entropy scaling model was also used by Lötgering-Lin and Gross [102] for the calculation of dynamic viscosity [103]. For a pure component, the reduced dynamic viscosity, is obtained with the following expression: (4.9) The process for the calculation of the mixture value is explained in Vidal et al. [44]. In case the mixture is in VLE state, the following mixing rule is used, from Beattie and Whalley [104]

VLE calculation
For the calculation of the Vapor Liquid Equilibrium, the minimum of the molar Helmholtz Free energy A is calculated, defined in terms of density ρ, temperature Τ and composition z as: where R g is the universal gas constant and f i is the fugacity of the component i. This optimization problem is solved via a combination of the successive substitution iteration (SSI) and the Newton minimization method with a two-step line-search procedure, and the positive definiteness of the Hessian is guaranteed by a modified Cholesky decomposition [105]. The algorithm consists of two stages: first, the mixture is assumed to be in a single-phase state and its stability is assessed via the minimization of the Tangent Plane Distance (TPD), see also [106]. The stability is tested by purposely dividing the homogeneous mixture in two phases, one of them in an infinitesimal amount and called 'trial phase'. For any feasible twophase mixture, if a decrease in the Helmholtz free energy is not achieved, then the mixture is stable. The so-called tangent plane distance (TPD) as function of the density times the composition of the trial phase ρ'x i ' is where the tildes over the variables mean those calculated at the trial conditions and the asterisk those calculated at the feed conditions. In case the minimum of the TPD is found to be negative, the mixture is considered unstable and a second stage of flash, i.e. phase splitting, takes place consisting on the search for the global minimum of the Helmholtz Free Energy. As a result, the pressure of the fluid and the compositions of both the liquid and vapor phases are calculated. P. Koukouvinis, et al. Fuel 275 (2020) 117871 Appendix 5. Properties used in the current study The following property tables are provided: -the NIST derived log 10 (p)- T table for dodecane. This table includes the correction for thermal conductivity at high temperatures from the PC-SAFT  model. The format of the table is: decimal logarithm of pressure, temperature, density, enthalpy, entropy, heat capacity at constant pressure, speed of sound, thermal conductivity, viscosity, partial derivative of density in respect to temperature, partial derivative of density in respect to pressure, partial derivative of enthalpy in respect to temperature, partial derivative of enthalpy in respect to pressure. -the NIST derived log 10 (p)-T-y table for dodecane/nitrogen mixture. As the previous table, it includes the correction of thermal conductivity. It is calculated based on the weighted average of properties between real-fluid dodecane and nitrogen (no non-ideal mixing rules). The format of the table is: decimal logarithm of pressure, temperature, density, enthalpy, entropy, heat capacity at constant pressure, speed of sound, thermal conductivity, viscosity, partial derivative of density in respect to temperature, partial derivative of density in respect to pressure, partial derivative of enthalpy in respect to temperature, partial derivative of enthalpy in respect to pressure, nitrogen mass fraction. -the PC-SAFT table log 10 (p)-T-y used for liquid identification. It is calculated based on the weighted average of properties between real-fluid dodecane and nitrogen, including all non-ideal effects modelled using PC-SAFT. The format of the table is: decimal logarithm of pressure, temperature, mass fraction of dodecane, density, enthalpy, entropy, heat capacity at constant pressure, speed of sound, thermal conductivity, viscosity, partial derivative of density in respect to temperature, partial derivative of density in respect to pressure, partial derivative of enthalpy in respect to temperature, partial derivative of enthalpy in respect to pressure.

Appendix 6. Liquid volume fraction projection calculation and comparison with the present data
An alternative threshold for liquid penetration has been also identified through recent ECN meetings, based on the projected liquid volume fraction, see [107]. According to this criterion, the liquid volume fraction has to be integrated along the perpendicular, to the jet axis, direction, following the beam path, to obtain the path integral as follows: Liquid identification can be correlated to an integral value greater than two proposed threshold values, as established using experimental uncertainties; these threshold values are 0.2 · 10 −3 mm 3 /mm 2 and 2 · 10 −3 mm 3 /mm 2 .
An indicative instance is selected to demonstrate that the liquid penetration results presented here, calculated using the isosurface of 0.15% liquid volume fraction, are equivalent to the specified threshold values, see Fig. 27. In this figure (Fig. 27a), an indicative instance of the injection to 900 K, 60 bar is shown, at a slice of crossing through the midplane of the injector. The slice shows the liquid volume fraction, along with the process of calculating its path integral along the axial direction. Results of the path integral along the axial direction for the given slice are shown in Fig. 27b. Liquid penetration can be predicted as follows, with the different criteria discussed so far: -Threshold value of 0.2 · 10 −3 mm 3 /mm 2 : 10.262 mm -Threshold value of 2 · 10 −3 mm 3 /mm 2 : 10.217 mm -Isosurface of LVF 0.15%: 10.257 mm As shown, the deviation between min/max penetration is less than 0.5% for the different estimation methods, with the LVF = 0.15% isosurface being in between the 0.2 · 10 −3 and 2 · 10 −3 mm 3 /mm 2 threshold values. The projected liquid fraction, as calculated by using the path integration, along the axial direction (thick black line). Also, the two different threshold criteria are shown (blue dashed line and red dotted line, at 2 · 10 −3 and 0.2 · 0 −3 mm 3 /mm 2 ). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)