Large-Eddy Simulation (LES) of Spray Transients: Start and End of Injection Phenomena

— This work reports investigations on Diesel spray transients, accounting for internal nozzle ﬂ ow and needle motion, and demonstrates how seamless calculations of internal ﬂ ow and external jet can be accomplished in a Large-Eddy Simulation (LES) framework using an Eulerian mixture model. Sub-grid stresses are modeled with the Dynamic Structure (DS) model, a non-viscosity based one-equation LES model. Two problems are studied with high level of spatial and temporal resolution. The ﬁ rst one concerns an End-Of-Injection (EOI) case where gas ingestion, cavitation, and dribble formation are resolved. The second case is a Start-Of-Injection (SOI) simulation that aims at analyzing the effect of residual gas trapped inside the injector sac on spray penetration and rate of fuel injection. Simulation results are compared against experiments carried out at Argonne National Laboratory (ANL) using synchrotron X-ray. A mesh sensitivity analysis is conducted to assess the quality of the LES approach by evaluating the resolved turbulent kinetic energy budget and comparing the outcomes with a length-scale resolution index. LES of both EOI and SOI processes have been carried out on a single hole Diesel injector, providing insights in to the physics of the processes, with internal and external ﬂ ow details, and linking the phenomena at the end of an injection event to those at the start of a new injection. Concerning the EOI, the model predicts ligament formation and gas ingestion,


INTRODUCTION
Research on fuel injection technology, sprays and combustion is continuously motivated by the increase in the global energy demand, especially for the transportation sector. Improvements in combustion engine efficiencies, the introduction of alternative fuel blends and the increase of energy demand for commercial transportation continue to propel research in the area of both light and heavy duty vehicles.
Direct injection engines, in the quest for higher efficiency and pollutant reduction, are making use of higher and higher injection pressures, up to 3 000 bar for compression ignition engines. At such elevated pressures, the flow develops high velocities within the injector that can lead to cavitation (Arcoumanis et al., 2000). Fast transients due to short injection durations are also of utmost importance as both direct injected gasoline and Diesel engines often adopt multiple injections, and possibly injection rate shaping, to cope with emission limits and fuel efficiency targets. Under these conditions, the transient part of spray behavior is most of the time the dominant part of an injection event. In Diesel common rail systems, commanded pulse duration can easily be as short as 200 ls, for pilot or post injections, even targeting injected mass as low as 1 mg/shot. Furthermore, the minimum dwell time between two consecutive injections has been progressively reduced, approaching nowadays 200 ls for solenoid injectors. As near nozzle region is concerned, most of the research has been devoted to the steady state part of a fuel jet, focusing on the steady liquid break-up, air entrainment, and mixing. Little information is available on the transient behavior during both opening and closing phases.
The shape of injection ramps affects spray penetration and mixing process, especially with short multiple injection strategies . It has been shown that after an injection pulse fuel dribble can be produced, and a significant amount of ambient gas is also ingested in the nozzle sac depending on factors like ambient back pressure, hole sizes, etc. (Swantek et al., 2013(Swantek et al., , 2014Battistoni et al., 2014b;Eagle and Musculus, 2014). The start of a new injection event is consequently affected by the presence of residual gas in the sac. During the first 100 ls after SOI liquid penetration is shorter and evaporation rate is altered . Additionally, droplets dibbled after the EOI create locally fuel-rich pockets within an overall fuel-lean region near the injector, the combination of which may contribute to soot and/or unburned hydrocarbon emissions (Musculus et al., 2013).
In the area of nozzle flow and sprays, the large variability of scales makes the two-phase flow modeling particularly difficult and computationally demanding. Some numerical methods track the liquid-gas interface explicitly like Level-Set (LS) or some formulations of Volume Of Fluid (VOF) methods, other methods do not pursue interface surface reconstruction (Gorokhovski and Herrmann, 2008). The first group of models requires very fine grids and large computational resources since it aims at directly resolving the mechanics of the liquid backup and drop formation. These approaches generally assume incompressible phases, no mass transfer due to phase change, relatively low Reynolds numbers, and generally limited to only external jet domain (Arienti and Sussman, 2014;Bianchi et al., 2007;Bode et al., 2014). To solve all scales these methods rely on extreme fine grids as no turbulence models can be easily incorporated. In addition, downstream they can be coupled with Lagrangian particletracking models to further simulate the motion of the smallest drops in the dilute spray region. The inclusion of complex geometries or moving boundaries is generally quite challenging, therefore these simulations, based on first-principles with high-fidelity, are limited to the external nozzle region, and require an inflow boundary derived from other type of in-nozzle flow calculations.
The second group of methods uses a diffuse-interface modeling approach. Indeed the liquid breakup is not directly resolved but it is modeled as diffusion, still thorough the aid of a fine supporting grid. These approaches offer the advantages of easy inclusion of the energy equation, with complex equation of states for the fluids, mass transfer due to phase change and compressibility. They generally allow the modeling of complex and moving geometries at high Reynolds numbers. Turbulence is closed either with Reynolds Average Navier-Stokes (RANS) for relatively inexpensive computations or with LES. All these approaches are Eulerian in nature, and can be formulated as single-fluid (Oefelein et al., 2014;Schmidt et al., 2010;Demoulin et al., 2007;Xue et al., 2015) or multi-fluid (Alajbegovic, 1999;Wang et al., 2014) two-phase models. Eulerian models can be coupled with Lagrangian particle-tracking models in the dilute spray regions as well, resulting in the method termed the "Eulerian Lagrangian Spray and Atomization (ELSA)" model (Demoulin et al., 2007;Vallet et al., 2001). The modeling of phase change is also particularly relevant within these approaches. First principles can never be applied to resolve locally the dynamics of single bubbles, therefore simplified models based on Rayleigh-Plesset equation or on the homogeneous assumption (either equilibrium or non-equilibrium) are generally employed (Wang and Greif, 2006;Habchi, 2014). Significant improvements have been demonstrated if noncondensable gases are accounted for in the multi-phase model (Battistoni et al., 2014a).
In this paper, we present an Eulerian method to simulate continuously the internal nozzle flow and near exit jet region, accounting for the needle motion, liquid compressibility, phase change due to cavitation and presence of non-condensable gas. The high pressure fuel flows through the injector needle valve, then through the orifice and is injected in an initially quiescent ambient gas. In this method, both the liquid and gas phases are solved in the same Eulerian coordinate system using a single fluid representation. The interface between these two phases is not tracked explicitly, but it is resolved only using the grid. The turbulence effects are modeled performing LES with the Dynamic Structure (DS) sub-grid scale model (Pomraning and Rutland, 2002).
Due to the high Reynolds number involved in practical Diesel injections (Re is greater than 500 000 in the present cases) addressing these problems is extremely challenging even using High Performance Computing (HPC) resources. In addition, because of the large differences in the length and time scales involved in the internal flow compared to the outward spray, the concurrent simulation of the internal and external domains is extremely demanding. This article presents an analysis of injection transient phenomena by means of high resolution LES, and, recognizing these challenges, part of the work focuses on monitoring the quality of these scaleresolving simulations and on quantifying the amount of resolved scales for assessing the model applicability. Despite the modeling challenges, the rest of the article highlights the physics of the injection transients in detail, and simulation results are compared with experimental data.
The remainder of the article is organized as follows. First, SOI and EOI phenomena are briefly reviewed. Then the Eulerian multi-phase model and the LES framework are described. Next, EOI and SOI cases are presented as applications. Results are compared against available synchrotron X-ray data and mesh requirements are discussed. Finally, the last section concludes the paper. Yu et al. (1983) conducted a pioneering study showing that EOI rate shape, sac volume and nozzle geometry significantly affect the production of unburned hydrocarbons in Diesel engines. They identified local over-mixing, poor end of injection and fuel emptying from the sac as important sources of hydrocarbon production. Optical measurements of Diesel injection showed that the mixing between fuel vapor and air is enhanced during the decelerating period after the end of injection, leading to low local fuel concentration (Bruneaux, 2005). Musculus (2009) showed that fuel-lean regions that form during the ignition dwell after the end of injection are likely a significant source of unburned hydrocarbons. A region of increased entrainment was identified. He also suggested that the deceleration rate could be shaped to achieve a desired mixing level after the end of a fuel pulse. It has also been shown that the EOI behavior affects soot formation and oxidation both for conventional Diesel combustion and Low Temperature Combustion (LTC) conditions (Musculus et al., 2013).

End of Injection
To the best of our knowledge, few attempts have been made on simulating end of injection process within the injector and in the near nozzle region. Hu et al. (2010) conducted a LES of a transient air jet focusing on the EOI behavior in the external jet region. They suggested that the deceleration of the jet enhances flow instabilities, increasing small-scale turbulence and causes higher entrainment.
Only recently, in regard to the EOI behavior, the focus has shifted toward the near nozzle region and the internal nozzle flow. Research at ANL using different X-ray diagnostic techniques has revealed new features of the EOI behavior (Kastengren et al., 2008). Figure 1 shows examples of ingested bubbles and dribble formation after an injection event for: a) a seven hole nozzle and, b), c) two single hole nozzles. Images are obtained using synchrotron X-ray phase contrast technique on metal Diesel injectors for the internal flow and the outside region. These results show that at the end of an injection pulse, rapid deceleration of the liquid jet eventually leads to its break up, forming poorly atomized spray and ligaments injected after the needle closes (dribbles). Correspondingly, a significant amount of ambient gas can be entrained in the orifice and in the sac of the injector nozzle, as the needle valve closes. Gas bubbles have been observed inside metal nozzles, and the volume of residual void after the injection process has been quantified for some injectors.
However, the details of the physics of this process still need to be fully understood. The dynamics of liquid ligament formation, and gas ingestion are very fast processes. Numerical simulations can help identifying the key parameters of the process, providing information in terms of pressure and void fraction evolutions inside the nozzle. These data should help answering questions like the following. Are internal bubbles due to local cavitation or due to gas ingested from the chamber? Is there a link between the extent of bubbles and the amount of external dribbles?
In this paper, numerical simulation were performed to characterize the EOI process in single hole injectors. Primarily, the focus is on the last part of the injection, encompassing the needle closure and extending for a large amount of time after it.

Start of Injection
Multiphase flow is generally expected in fuel injectors due to cavitation, but also due to ingestion of ambient gas in the nozzle sac after an injection event, as discussed in the previous section. The sac can be partially empty at the start of injection. The presence of hot ambient gas in the sac can be responsible for delayed liquid injection and higher liquid vaporization thereafter (Crua et al., 2010;Pickett et al., 2013). Crua et al. (2010) and Stetsyuk et al. (2014) observed Diesel jets in hot conditions using ultra-high speed video and found that the initial stage of fuel injection is more complex than the accepted linear temporal evolution. A gaseous pre-jet, with spheroidal cap shape, was observed for a wide range of conditions and issues sub-sonically into hot and still compressed air. It may consist of cavitation pockets expanded after previous injection; ingested in-cylinder gases after previous injection; heated and vaporized fuel trapped in the injector holes between injections. The high speed liquid fuel jet initially trails the gaseous cloud but later its tip proceeds further ahead. This process is observed to occur within the first 1-2 mm of jet penetration. Pickett et al. (2013) reports that in vaporizing conditions the liquid penetration of the Engine Combustion Network (ECN) spray A can be about 2 mm shorter in the first 200 ls than the steady state value. The behavior is attributed to the initial presence of residual gas or fuel vapor in the sac, which enhances vaporization rate of the first-injected fluid, as the gas trapped from the chamber in the sac is most likely warmer than the incoming liquid fuel. This also suggests that this issuing fuel can be ignited by the hot air mixed with remaining burnt products in the chamber.

Governing Equations
The compressible Navier-Stokes equations, which express the conservation of mass, momentum, energy and species, for a single-fluid mixture read as follows: where q is the mixture density (or pseudo-density), u is the velocity, p is the pressure, and r ij is the mixture shear stress tensor assumed to behave as Newtonian. F sf is the surface tension, formulated as a volume force, and evaluated using the hypothesis of continuum surface force (Brackbill et al., 1992). The momentum equation does not include any explicit term for sub-grid phase slip, therefore this contribution is neglected. However, as a consequence of the turbulent mass dispersion closure (cf. first right-hand side term in Eq. 4) a diffusion velocity difference between the two phases does exist (Demoulin et al., 2007). Y i denotes the mass fraction of species i. The energy Equation (3) is formulated in terms of specific internal energy e of the mixture. T is the mixture temperature, and h k is the specific enthalpy of species k. In addition, k is the thermal conductivity and D is the mass diffusion coefficient expressing the diffusion of a component in the mixture. D is assumed constant for each component. With this formulation, summation over i of Equation (4) yields global continuity (1). Continuity is solved to obtain the mixture density, and at the same time all the species transport equations are solved for the species mass fractions. A reconciliation procedure at the end of each loop ensures that any discrepancy due to round-off error is adjusted by renormalizing the species mass fractions to satisfy global continuity. Lastly, S ik are source terms related to species mass transfer due to phase change. Cavitation/condensation sources are calculated according to the Homogeneous Relaxation Model (HRM) of Bilicki and Kestin (1990) and Schmidt et al. (2010), described later in Section 2.2.
The multi-phase system is assumed to be a multi-component two-phase mixture comprising of: the liquid fuel (subscript 1), and multiple gaseous species, i.e.: fuel vapor (subscript 2), non-condensable gas mixed with the liquid (subscript 3), other gaseous species, like ambient chamber species, etc. (subscripts 4, etc.). The mixture has a variable density given by the weighted average of the single components densities. Gas phase densities are determined using the equation of state for an ideal gas: where R i is the gas constant and T is the temperature. Liquid is considered compressible. The following barotropic (q liq ¼ f p ð Þ) equation of state is used: where K is the liquid bulk modulus, and subscript 0 refers to an arbitrary reference condition.
The sum of all gas species will be referred to with the subscript g. The mixture density q is expressed through the following equation: where a i is the volume fraction of each component. Volume fractions a i and mass fractions Y i are related through: (without summation over repeated index). Using Equations (7) and (8), the cumulative void fraction a g takes the form: This cumulative void fraction includes voids generated by phase change (cavitation), expansion of non-condensable gases, and mixing with ambient gas.
In this implementation, the void fraction a g is not calculated with an explicit transport equation, as found in other VOF formulations, but the interface is resolved only using the grid. Individual species mass fraction Y i of each fluid is tracked throughout the domain, by first solving the related species Equation (4). Then the void fraction is calculated using Equation (9).
This approach results in a diffused-interface modeling. As an example, the liquid break-up is modeled as diffusion and through the aid of a sufficiently fine grid. This approach, compared to geometric interface reconstruction methods usually employed in VOF, offers the advantages of easy coupling with energy and turbulence equations, inclusion of phase change, and modeling of complex and moving geometries for high Reynolds number flows.

Phase Change
Mass transfer due to phase change can occur at the liquid/gas interface. In general, in a multi-component mixture, it can occur between the liquid fuel and its vapor, or it can potentially occur when gas is released from or absorbed into the liquid solution. Different species are involved in the two cases and different mechanisms are involved.

Liquid-Vapor Mass Transfer Model
The model for mass exchange between liquid and vapor (cavitation) is based on the HRM model as developed by Bilicki and Kestin (1990). The model assumes a simple first order rate equation for the evolution of the instantaneous vapor quality x towards its equilibrium value x over a given time-scale H. This synthetic approach has been proved to be able to represent not only cavitation but also flash boiling (Schmidt et al., 2010). These are phase change processes driven by a pressure drop. Usually, the process is referred to as 'cavitation' when a tiny amount of superheating is required to evaporate the liquid and the process is very fast, whereas 'flash boiling' is used when a higher degree of superheating is achieved, and the process takes a relatively long time. The HRM model expresses in a synthetic way the extremely complex process by which the two phases exchange heat and mass. The rate equation takes this form: The instantaneous non-equilibrium vapor quality x (or dryness) is defined as: The equilibrium vapor quality " x is a function of the thermodynamic properties at the local pressure, i.e., " x p; h ð Þ is expressed by the following equation with bounds at zero and unity: where h is the actual fluid enthalpy of liquid and vapor mixture, excluding air content, and h 1,sat and h v,sat refer to saturated liquid and vapor, respectively. The time-scale H is evaluated using an empirical fit proposed by Downar-Zapolski et al. (1996): H ¼ 3:84 Á 10 À7 a À0:54 / À1:76 ðunits of sÞ ð13Þ where a is the gas phase volume fraction, and subscripts sat and crit denote saturation and critical pressures, respectively. Correlation provided by combining Equations (13) and (14) has been determined in experiments of flashing flows of water in pipes with upstream pressures of in excess of 10 bar. In the absence of specific values for other liquids and conditions, it has been applied in the present form in this study for representing Diesel fuel. The time integration of the rate Equation (10) allows to compute the source terms S 12 and S 21 (= ÀS 12 ) in Equations (4). Additional details on this source term evaluation are given in Battistoni et al. (2015).

Liquid-Air Mass Transfer Assumption
The model also account for the presence of non-condensable gas in the mixture (third component), however, the related species mass transfer source term is neglected, i.e., S 13 = 0. Within the time scale of the injector flow, events are assumed to be very rapid for significant mass transfer of contaminant gas to occur between a bubble and the surrounding liquid. The mass transfer of air from the liquid solution to the gaseous form (gaseous cavitation), and vice versa, is neglected. Indeed, mass transfer due to diffusion or dissolution is extremely slow compared to the time scales involved in the main flow (Brennen, 1995;Battistoni et al., 2015), especially during SOI or EOI transients.
In the current implementation, the amount of air added to the liquid fuel is introduced directly in gaseous phase, homogeneously mixed with other components. In other words, we have not introduced a specific model for air adsorption and release from the solution, and in Equation (4), the related source term S 13 is set equal to zero. As a result, as pressure in the fluid drops, air promptly expands, because it is already assumed to be in gaseous form. Conversely, if a cloud of expanded air exists and pressure rises back, gas will be compressed, but not reabsorbed in the liquid medium.
An effective constant amount of non-condensable air in gaseous form is assumed mixed with the liquid. In the current work, the value is set to Y 3 = 2910 -5 by mass, according to previous works (Battistoni et al., 2014a, and the assumption is based on an analysis of the characteristic time-scale of air mass diffusion in the liquid solution.

Sub-Grid Models for Filtered Equations
LES have been performed covering the inner nozzle region and the outside near exit region. The code uses a modified cut-cell Cartesian method for generating the grid at run time (Senecal et al., 2007). In this framework, the filter size D and shape are fixed by the computational cell size and shape, therefore the filter type is a box.
The set of Equations (1-4) is recast according to the LES decomposition into resolved fields and sub-grid fields, i.e., for any flow variable / ¼ " / þ / 0 . For compressible flows Favre averaging is conveniently applied, yielding / ¼ q/=" q. The procedure results in additional Sub-Grid Scale (SGS) terms that need to be modeled.
As far as the momentum equation is concerned, following Pomraning and Rutland (2002), the DS model exploits the significant similarity between the SGS and the Leonard stresses that has been observed in real as well as numerical experiments. The sub-grid stress tensor, where the Leonard stress tensor L ij can be calculated from the resolved velocity field. The hat symbol c ðÁÞ designates the test level filter, which is larger than the grid level filter, and in this implementation it is created from the center cell and the neighboring cells.
The DS model also involves an additional transport equation for the sub-grid turbulent kinetic energy, k sgs ¼ g u i u i Àũ iũi ð Þ =2, which takes the form: The sub-grid dissipation rate is given by: where the model constant C e is set equal to 1. As a consequence of the non-eddy viscosity nature of the DS model, the SGS energy diffusion due to turbulent viscosity is not present in the transport Equation (16): i.e., Pr sgs = 1. However, other sub-grid terms resulting in the energy and species conservation equations are modeled using an eddy viscosity concept. As an example, a turbulent diffusion concept is used and the turbulent diffusion is given by: Here, C k is a model constant set to 0.5 and Sc is the turbulent Schmidt number set to 0.71. Further details about this implementation are provided in Senecal et al. (2013Senecal et al. ( , 2014, Som et al., 2012b andXue et al. (2013).

Numerical Method
The overall set of equations is solved using the finite volume method. All variable values are co-located at the cell center and the Rhie-Chow algorithm (Rhie and Chow, 1983) is used to prevent pressure field decoupling (checker-boarding). First, the conservation equations of mass and momentum are solved using a pressure-velocity coupling iteration method, specifically the Pressure Implicit with Splitting of Operators (PISO) method of Issa (1986). Other transport equations, such as the energy and the species, are solved after the momentum predictor and the first corrector have been completed. Previous studies on nozzle flow simulations using this algorithm for a mixture model can be found in the work by Zhao et al. (2014) and Xue et al. (2014).
Spatial discretization is second-order accurate using central differencing schemes for all quantities wherever is possible. A minmod type of flux limiter is used to capture sharp gradients as well as to maintain high accuracy in the smoother regions . Grid refinements are managed through time-dependent fixed embedding located in the key regions. Time integration is first order accurate running fully implicit Euler scheme, with fine time-stepping controlled by a velocity-based Courant-Friedrichs-Lewy (CFL) below 0.1 and a speed-of-sound-based CFL below 2.0. This multiphase flow model has been recently implemented in the Computation Fluid Dynamics (CFD) code Converge.

ASSESSMENT OF LES RESULTS QUALITY
It has to be appreciated that the practice of building a grid and applying LES using the grid itself as a turbulence resolution length scale makes the model incomplete (in the sense of Pope, 2004), i.e., the specification is subjective and the grid becomes a parameter of the LES equations. Several efforts have been devoted to implement solution-adaptive gridding methodology (Pope, 2004;Celik et al., 2005;Brusiani and Bianchi, 2008;Davidson, 2009;Piscaglia et al., 2014;Xiao et al., 2014), however, here we only undertook an assessment of the results quality after the grid was prescribed.
Many metrics have been proposed in the literature to evaluate the quality of LES results. Among them, two main types are reviewed here. One is based on the resolved fraction of turbulent kinetic energy and the other is based on the length-scale resolution. The first type of criteria, following Pope (2000Pope ( , 2004, requires that the modeled turbulent kinetic energy be less than 20% of the total turbulent energy. It is expressed as: It is worth noting that the evaluation of the mean resolved turbulent kinetic energy, k res , requires temporal averaging or ensemble averaging over multiple realizations of the experiment (the mean operator is denoted here by angle brackets Á h i). In addition, to some extent the value of M is non-local, as it is affected by the grid resolution at other locations and at earlier times (Pope, 2004).
Another type of metric is based on the turbulence-resolution length scale. Following again Pope (2000Pope ( , 2004 and Brusiani and Bianchi (2008), the criterion requires that the resolved scales be larger than the dissipation range scales. The demarcation length-scale, between the inertial and dissipative ranges, l DI , is often estimated, according to Pope (2000), as l DI ¼ 60g, where g is the Kolmogorov length scale. The latter, g ¼ m 3 =e ð Þ 1=4 , can be evaluated using the sub-grid energy dissipation rate, e, which in this work is given by Equation (17). In synthesis, a Length Scale Resolution (LSR) parameter can be defined as: It is worth observing that expanding out its definition, LSR depends on k sgs , therefore on the turbulence closure model.
In addition, the evaluation of the demarcation length-scale l DI is based on scaling arguments strictly valid for equilibrium, isotropic turbulence.
With that in mind, if LSR is below 1, all the scales in the inertial sub-range and above are resolved, while only the dissipative range is modeled. In practical applications, as suggested by Brusiani and Bianchi (2008), a value of LSR < 3-5 can still be accepted for a good LES calculation at a reasonable cost. One advantage of this criterion, compared to the previous one, is that it does not require statistical averaging of the resolved quantities.
Lastly, one should always keep in mind that regardless of the metrics being used (either Eq. 19 or 20, for example), the numerical discretization errors can play a significant role in the assessment. In fact, the numerical viscosity generates an additional contribution, which should be evaluated explicitly (Celik et al., 2005). Therefore, the modeled SGS turbulent kinetic energy, which appears in both parameters (19) and (20), should be increased by the numerical viscosity counterpart. In both cases, the indexes would give a worse assessment. However, the estimation of the numerical error is beyond the scope of the present work and will not be included in the analysis.

APPLICATIONS TO TRANSIENT DIESEL JETS
A validation study of the proposed multiphase mixture model within the LES framework has been shown by , for the case of Diesel sprays producing shock waves at relatively high injection velocities and low temperature ambient conditions.
In this work, the model is applied to study the transient behavior of jets issuing from Diesel injectors, with focus on the initial and final part of a typical injection event. The in-nozzle flow is part of the computation, since the physics mostly depends on what occurs inside the injector. Recently, these topics have taken on particular importance as short and multiple injections are often used in advanced engine combustion strategies as discussed in the Introduction section. However, a combination of complex physics (multi-phase, compressibility, phase change) and uncertainties in boundary conditions (accuracy of the internal geometry, moving boundaries, etc.) make the prediction of these flows a particularly challenging task.

Problem Description
First, the multi-phase model is applied to simulate a singlehole mini-sac Diesel injector nozzle, for which synchrotron X-ray imaging data, collected at the Advanced Photon Source (APS) at ANL, are available (Swantek et al., 2013(Swantek et al., , 2014. In the experimental campaign, different injection pressure (from 500 bar to 1 500 bar), back-pressure (from 1 bar to 20 bar) and orifice diameter (from 110 lm to 180 lm) have been investigated. Here, we present the simulation of the case with maximum expected ingested gas among the above mentioned conditions, i.e., 1 500 bar injection pressure, 1 bar chamber pressure, and 180 lm orifice diameter. A summary of conditions and fuel properties, along with the main injector features is provided in Table 1. Also, the needle lift profile is shown in Figure 3.
Focusing on the simulation of the EOI process, we skipped the initial and central part of the injection event, and started the calculation at 1 900 ls where the lift is still high for this injector, being 150 lm. Results (shown later) demonstrate that starting from this lift the flow rate reaches a steady condition before the effective ramp-down starts.
The objective is to reduce the computational time, yet provide a good estimate of the initial flow for the simulation of the EOI behavior. The needle closure occurs in between 2 150 ls and 2 175 ls, as visible in Figure 3 (experimental data have a temporal resolution of 25 ls). To deal with very low needle lifts, we assumed the needle valve closed below 10 lm. This occurs at 2 140 ls, therefore the flow is simulated as a single domain up to this point. After the closure we switch to two disconnected domains, and the simulation continues up to 2 500 ls, to ensure that all the relevant dynamics is covered. The minimum cell size is 5 lm, and it has been maintained throughout the simulation.
The main focus for the gridding strategy has been on the sac, the orifice, and the near nozzle regions. A view of the computational grid is shown in Figure 4. Grid is temporally refined and changed for better adapting to the relevant dynamics. The region around the orifice is discretized with the finest cell size throughout all the ramp-down phase and it is retained after the needle closure. To save computational time, during the start-up phase of the calculation (before 2 000 ls, at high lift) the finest cell size has been set to 10 lm.
Associated with the transient needle movement profiles, pressure and temperature are specified at the inlet boundary. The turbulent intensity is set to 0.1% at the inlet, in terms of k sgs value. The Werner and Wengle (1991) wall model is employed at the wall boundaries. The outlet chamber is smaller than the real one, but its size is large enough to accommodate all the jet development in the time interval of interest. The length is 20 mm and the diameter is 16 mm. The bottom surface of the chamber is closed and the jet will eventually imping on it. The side walls are open and the static pressure is applied as boundary condition. The near nozzle region, below 10 mm, is unaffected by the  Figure 3 Fuel properties
choice of the chamber size, yet the savings in terms of total number of cells are substantial. Overall, the mesh size for a fine grid (Tab. 2, Grid 2) ranges from 10 to 16 million cells over the simulated period. Sections 4.1.2 and 4.1.4 provide more details on the grid suitability and simulation cost.

Discussion of Results: Dynamics of the Flow
The main global parameters predicted by the LES are shown in Figure 5. The top chart displays the average velocity at the nozzle exit. The middle chart displays the average and minimum pressure in the sac and orifice regions. The bottom chart shows the global (integrated) void fractions in the same regions. Observing the behavior of these quantities, we have a synthetic view of the dynamics occurring during the EOI as predicted by the model.
The first part of the simulation provides the initialization of the fuel jet while the needle lift is still high. The main flow features are basically steady for t < 2 100 ls, where the lift is above 50 lm, signifying that the flow losses are located in the orifice and not across the needle. After this point, the restriction imposed by the needle becomes relevant and ultimately causes the strong flow deceleration, which takes place in less than 40 ls, i.e., from 2 100 ls to 2 140 ls (instant of closure). At 2 170 ls, the average orifice pressure reaches a minimum, dropping by four orders of magnitude in about 25-30 ls. The average sac pressure shows a slightly smoother drop, due to the averaging over a larger region. The minimum pressure drops very fast in both regions. Spots of extremely low local pressures are found in the sac and in the orifice, immediately after the needle closure. The minimum pressure level detected after the closure reaches 1-10 Pa (absolute), therefore locally a strong cavitation is initiated, since the fluid has a vapor pressure equal to 1 000 Pa at the reference temperature of 25°C. Vapor phase is formed and at the same time the dissolved gas in the fuel expands. The global result of this gases being formed or expanded is the global void fraction a g shown in Figure 5c (integrated over each region). The effect of the abrupt deceleration of the liquid column at the needle closure is the formation of approximately 20% and 40% of void fraction in the sac and in the orifice, respectively. Later, from 2 250 l to 2 350 ls, a back flow takes place, as displayed by negative velocities at the hole exit in Figure 5a. During this time interval, gas from the chamber is ingested, which eventually fills entirely the orifice and partly the sac (Fig. 5c after  2 400 ls). Residual gases are therefore predicted inside the nozzle after the EOI.
In order to have a clear view of the results, contour and iso-surface plots from the LES are shown in Figure 6, side by side with X-ray phase contrast imaging data. Contours are shown in a plane slicing the domain through the hole  Table 2. axis, and refer to liquid volume fraction, LVF (¼ 1 À a g ). The iso-surfaces are built at LVF = 0.1 and are colored by pressure. Different timings are presented, ranging from the high lift condition to the end of the process. The quasisteady injection is shown at 2 000 ls, the ramp-down of the mass flow at 2 100 ls. In the model, the needle closure is at 2 140 ls, while the available experimental image is at 2 150 ls. Here, three calculated frames are shown, i.e., at EOI and after 5 ls and 10 ls. Shortly after the EOI, other comparisons are provided at 2 160, 2 200 and 2 230 ls. Finally, results at 2 300 and 2 500 ls refer to a later stage.
In addition, Figure 7 shows turbulent flow structures visualized in terms of Q-criterion at 2 100 ls.
Observing Figure 6, a satisfactory qualitative agreement is found between the calculated LVF and the experimental imaging. At high lift, an intact liquid core exits from the orifice (cf. results at 2 000 ls). During the ramp down phase (2 100 ls) liquid decelerates and the liquid-gas interface shows some wrinkling. When the needle touches the seat, the fuel is subject to the highest deceleration. The positive tensile stresses on the liquid ultimately cause its breaks-up into ligaments (results from 2 140 to 2 200 ls). At 2 150 ls two cavitation spots are observed in the sac near the valve seat. Vapor in the sac then recondenses, as the local pressure rises back. Large areas of void fraction are also visible in the orifice as a result of the liquid jet rupture. Soon after, a pressure wave travels from the chamber towards the sac and brings back a liquid pocket. Behind the residual liquid, also some gas from the chamber flows toward the sac (cf. results from 2 200 ls to 2 300 ls). As a result of this phase, ambient gas will partially fill the sac.
The formation of ligaments in the near nozzle region is predicted to occur from the time of closure up to approximately 200 ls later. The shape and structure of dribbles are very similar to those observed in the experiments, however, the temporal persistence of these structures is not as long as in the experiments. This is likely due to the diffuse interface method used, which does not track the interface geometry, but rather models the mass dispersion as a diffusive mixing process. Mass content and distribution are therefore correctly reproduced, but the atomization is not resolved in details.
At the end of the process, a fraction of the sac is empty. An interesting comparison is given in Figure 8, where the last available simulation time-step is compared to experimental images taken before the start of an injection pulse. It is found that the predicted void fraction left inside the injector matches quite well the experimental evidence of a residual gas bubble lying at the tip of needle. The experimental image is taken long time after the EOI, which cannot be practically simulated, therefore the last available time-step is used for the comparison. The visual comparison with the bubble detected in the experiments is certainly encouraging. The location and the amount of void fraction predicted by the model are also in good agreement with the X-ray image.
Overall, the case analyzed here clearly shows that a large amount of ambient gas can be pulled into the sac after the EOI, as the liquid jets decelerates and eventually breaks up forming dribbles. The mechanism of liquid deceleration, breakup and gas ingestion is well captured. The model needs to be further improved with respect to the preservation of the liquid-gas interface, and this is highlighted by the shorter   EOI behavior at different timings: X-ray imaging on the left, simulation results on the right. Each row shows X-ray phase contrast imaging of the external near nozzle jet with exposure time of 5 ls (Swantek et al., 2013), computed LVF on a cut plane, and iso-surfaces at LVF = 0.1 colored by pressure. Results shown refer to Grid 2 (fine) in Table 2. duration of the ligament phase compared to the experiments. An High Resolution Interface Capturing (HRIC) scheme has been used to prevent interface smearing (Zhao et al., 2014), however, this test case appears particularly challenging in respect to this, due to the very long physical time that needs to be simulated.

Grid Sensitivity and Evaluation of LES Quality
We investigate here the quality of the LES results, after the flow features have been already presented in the previous section. Two grids with different level of refinement have been built (Tab. 2). The evaluation of the simulation accuracy is done comparing the energy budget resolution parameter (Eq. 19) and the LSR parameter (Eq. 20). Results in terms of M and LSR indices are shown in Figure 9 for the fine grid (Grid 2) at high needle lift. The evaluation of M requires the calculation of turbulent statistics. Due to the computational cost of a single realization, ensemble averaging is not practical for this kind of analysis on injector closing transients. Therefore, the only viable option is temporal averaging: this has been accomplished by averaging the resolved field for 20 ls, from 1 980 ls to 2 000 ls, where the flow rate is substantially steady, as demonstrated in Figure 5. Temporal averaging is done at run time from the specified start-timing by tracking non-transport scalars.
Values for M are meaningful in the sac, orifice and near exit regions, while they may not be accurate close to the needle moving surface and at the spray tip for the intrinsic unsteadiness of the problem.
Focusing on the LES quality assessment, the contours of the turbulent energy resolution parameter, M, show that the grid is fine enough in the sac, in the core of the orifice and downstream up to 3-4 mm from the exit, where M is well below 0.2. On the contrary, observing the results in terms of LSR parameter, globally a better picture is perceived. The needle gap region seems well resolved (the lift here is above 100 lm), as LSR is well below 1. In addition, the mesh resolution at the orifice inlet edge is satisfactory based on the LSR metrics. We can comment on the comparison between the two quality criteria, by noting that the temporal averaging required for M makes this parameter not local (Pope, 2004)  Vortices visualized as Q-criterion iso-surface, marking a threshold that localizes the most coherent turbulent structures, at t = 2 100 ls. Top view shows details in the nozzle region; bottom view shows the whole field. Residual gases ingested at the end of the process. X-ray phase contrast image on the left, simulated LVF on the right (simulation time t = 2 500 ls; grid 2 -finein Table 2). Measures of LES quality for Grid 2 (fine) at high needle lift (t = 2 000 ls). Top plot displays Pope's criterion, M, (Eq. 19)threshold is 0.2; bottom plot shows contours of LSR criterion (Eq. 20),threshold is 1.
the flow through the gap region is highly dependent on the grid resolution. Also, the moving boundary may affect the results in terms of M, which appears to be very high, i.e., close to 1. On the other hand, the LSR value in the needleseat gap is very small. Based on this opposed answers, we might conclude that the use of these parameters in this narrow area needs to be further checked and tested; at the hole inlet edge, there is some recirculation and low velocities and at the same time high SGS turbulence is generated. This leads to very high M values. Additional cell refinement would likely be needed in this area; streaks of high M values are observed extending from the inlet edge down to the hole exit: this is probably due to the advection effect intrinsic in the definition of M. Conversely, the LSR parameter, which does not suffer from this non-locality problem, does show good grid quality in the orifice. After this assessment on the comparison of the energy budget and the length-scale resolution parameters, the grid suitability at different timings during the closing transient can be conveniently appreciated by using the LSR index, which is deemed a proper tool for dynamic resolution evaluations, especially in the presence of boundary movement and temporally changing flows. Figure 10 displays the contours at 2 100 ls, 2 140 ls and 2 170 ls. For each timing, a close-up and a global view are shown. Far from the nozzle (i.e., at distances greater than 5 mm) the LSR index becomes greater than 1, but these spots are located outside the region of interest for the current analysis, where the mesh is made coarse to reduce the total cost of the simulation. The focus is therefore inside the injector and close to the nozzle. The first contour plot (2 100 ls) refers to the jet deceleration phase, where velocities are still quite high (around 400 m/s). Here, the LSR results suggest again that near the hole inlet edge, the grid can be further refined. Elsewhere, the grid spacing is adequate. The second timing refers to the instant of needle closure. Here, the grid is satisfactory in the whole region of interest. At most, some refinements could be applied in the tail of the detaching jet and at needle seat exit. The last image refers to a timing when ligaments are being formed. Here, the turbulent flow is well resolved due to the low liquid speeds (few tens of m/s) and the large flow structures.
Computations with different grid resolutions have also been considered. A 10 lm grid is compared to a 5 lm grid (min cell size in the orifice and near nozzle region, Tab. 2). Results in terms of average pressure in the sac are shown in Figure 11. The prediction using the fine resolution shows that the pressure rises back earlier, therefore the backflow after the EOI occurs in advance. In addition, Figure 12 shows contours of LVF at three timings during ligament formation (2 160, 2 180 and 2 220 ls). The fine grid predicts an appreciably longer and sharper liquid dribble compared to the lower resolution mesh. This is probably due to a lack of small turbulence structures in the jet using the 10 lm grid, which are dissipated faster. As a matter of fact, the Eulerian mixture model is based on a diffuse interface assumption. Alternatively, a front tracking method could provide a better reconstruction of the phase interface for capturing ligament formation, however generally at the cost of increased complexity due to the absence of SGS turbulence or with additional limitations like phase incompressibility. Overall, this mesh resolution analysis suggests that the prediction of ligaments being formed after the needle closure requires a very fine grid in the nozzle region, and an even smaller cell size could be adopted. Instead, the global parameters describing the dynamics of the process do not change drastically with the resolution, therefore, trends are essentially well predicted within the accuracy of the model assumptions.

RANS-LES Comparison
The assessment of LES results can be complemented by the comparison with RANS results on the same grid and with similar numerical framework. Figure 13 shows the results of a RANS calculation with the standard k-e turbulence model on the low-resolution mesh (grid 1). The LVF contours can be directly compared to LES results shown in Figure 12 left column. All the time-steps shown have a common feature: the RANS turbulence model makes the flow extremely axisymmetric, as expected, but the consequence is that the liquid column rupture occurs with a very different mechanism compared to the LES prediction. The RANS model tries to capture the ensemble average behavior, which might be fully correct, but it cannot show the instantaneous features of a single realization as obtained with LES. Focusing inside the nozzle and the hole, the RANS model predicts a dynamics that has a one-dimensional character, therefore mass transfer, pressure, or other quantities propagate through quasi-planar wave fronts. In the RANS predictions, the residual liquid back flow occurs without fragmentation of the bulk mass, while in LES liquid breaks up generating complex three dimensional structures.
Another feature of interest is observed slightly after the needle closure (2 160 ls). In the RANS approach, an annular void region stems from the needle seat, and then re-condenses as the pressure rises back. This is not observed in LES, where spots of low pressure regions are produced even in areas detached from the needle seat.
The near exit region is also different. No ligaments are being formed using RANS, but rather a diffusive process takes place, which disperses the liquid in the chamber. With the same level of grid refinement, LES captures some low speed jet waving and eventually liquid islands. However, in the multi-phase mixture model the geometric reconstruction of the liquid-gas interface is not performed, therefore, even using LES neat pockets are not detected.  Effect of mesh refinements (Tab. 2) on a global quantity: average sac pressure versus time. Effect of mesh refinements (Tab. 2). Contours of LVF on a cut-plane.

LES
Lastly, Figure 14 shows the comparison of RANS and LES in terms of global quantities like pressure and void fraction averaged/integrated over the sac region. The pressure recovery in the sac takes more time with the RANS model. This is likely due to the "one-dimensional" character of the flow that has been observed in the RANS results. This point certainly requires further investigation, but it is our opinion that a fully three-dimensional pattern as observed with LES might be more realistic for a single realization. It is very likely that if several LES realizations had been run and averaged, results would be close to those obtained with RANS.

Computational Cost
Simulations were run using a time step controlled by a velocity-based CFL set to 0.1. With reference to the fine grid with 15.6 million cells (grid 2 in Tab. 2), the resulting timestep varied from 0.8910 -9 to 3.0910 -9 s during the simulated time. The total number of time steps was in the order of 300 000. The LES required about 33 days of wall clock time on 256 cores. Using the coarse grid the wall clock time was about 21 days on the same number of cores.

Problem Description
In this section, we present an analysis of the start of injection of the "spray A" single-hole mini-sac Diesel injector, which is widely characterized within the ECN. The multi-phase model is applied to simulate the needle opening phase starting from a very small needle lift (2.5 lm). The simulation aims at capturing the initial ROI transient, the internal nozzle flow features occurring at very low needle lift, the sac filling process with liquid fuel, and also the initial spray features. In particular, one of the objectives of this demonstration study is to assess the initialization of the internal nozzle regions and its effect on the initial jet behavior.
As discussed in the Introduction and shown in the previous EOI test case, there are several evidences that confirm presence of gas in the sac at the start of an injection event.
Here, as initial conditions we have assumed that the residual gas completely fills the sac and the hole volume at the beginning of the simulation, while liquid fuel only fills the region upstream of the needle-seat contact area.
A summary of nozzle geometry, operating conditions, and fuel properties is provided in Table 3, Figure 15, and Figure 16 (Tab. 4). Taking advantage of current efforts by the ECN to deeply characterize automotive injectors, the spray A internal geometry is thoroughly investigated. A high-resolution geometry reconstructed based on X-ray tomography is available, with surface features in the order of 1 lm (Fig. 16). The nozzle is tapered and has an equivalent outlet diameter of 89.4 lm. The hole axis has a significant offset from the nominal centerline. Needle motion, with three-dimensional trajectory, is available from high speed X-ray imaging through the metal injector operated in "cold" conditions, as tested at ANL . The lift curve is shown in Figure 15, and focuses on the initial phase. The needle opening starts at 150 ls after the electrical command. The first part of the rising ramp has a small velocity and after 150 ls the needle has only reached 9 lm lift. After this point, the main fast rise begins. In this study, we considered only the vertical lift motion and neglected the off axis components. In addition, the needle lift prescribed in the model is slightly different from the experimental data and it is shown by the linearized ramp in Figure 15. As a matter of fact, due to the uncertainty on the lift measurements (estimated to be around 2 lm) and to the noisy shape of the initial ramp, a smoothed linear rise from 150 ls to 310 ls is preferred. In addition, since a zero gap cannot realistically be modeled, the simulation starts at t = 200 ls, thus at point in time where the lift is about 2.5 lm. The simulated period extends up to t = 350 ls. i.e., encompassing only the initial injection stage. The fuel is n-dodecane introduced at 338 K. The injector body is at 303 K. Ambient conditions are non-evaporating (different from standard spray A specifications), but the chamber density matches the nominal density value of 22.8 kg/m 3 (pressure and temperature are correspondingly reduced). These conditions are those used for X-ray imaging and radiography at Argonne.
The high resolution geometry (Fig. 16) encompasses the orifice and the bottom part of the sac, while the needle-seat contact is not covered. Therefore, this part has been reconstructed in the CAD using two coaxial conical surfaces, with the same 30 degree nominal inclination with respect to the injector axis.
The computational grid (Fig. 16) has a base cell size of 40 lm and uses 4 levels of embedding to obtain a minimum cell size of 2.5 lm. The finest resolution is employed in the needle-seat gap and for 3 layers on the hole walls. Refinements are permanent and Adaptive Mesh Refinement (AMR) is not used. The outlet chamber is 16 mm 9 8 mm (length 9 diameter) and allows an adequate simulation of the jet up to 350 ls, as the tip penetration is around 10 mm at that time. Overall, the mesh size is 15.9 million cells. An additional "coarse" grid has been tested, which has 80 lm base size and 5 lm minimum cell size. Details are provided in Table 4.
The available experimental pressure trace has been applied at the inlet surface as total pressure. The chamber is large enough to be modeled as a closed vessel.  Spray A needle lift profile versus time, from X-ray imaging experiments . The main chart focuses on the starting phase. The insert shows the whole injection duration.

Discussion of Results
The representation of the moving needle and the use of the inlet pressure boundary condition enable a realistic calculation of the mass flux variation in time through the orifice. To show this, the global features of the transient simulation are shown in Figure 17. The top chart displays the liquid and gas mass flow rates exiting from the orifice. The bottom chart reports the liquid tip penetration, compared against available experimental data. The liquid outlet flow rate is negligible until the SOI is approached. Instead, the gas phase flow rate is basically always positive from the beginning of the simulation (200 ls), denoting the ejection of the residual gas from the sac. Evidently, the gas phase mass flux at the orifice exit is quite small because of the density difference compared to the liquid. Based on the liquid mass flow rate at the hole exit, the predicted fuel injection starts around t = 300 ls, when the curve has a steep increase. Experimental ROI curve is not available with sufficient time resolution and accuracy at the SOI, as needed for this validation.
A complementary assessment method is based on the analysis of the tip penetration, which is presented in the bottom chart of Figure 17. Experiments based on X-ray radiography (Argonne data, square symbols, by Kastengren et al., 2014) indicate that the SOI is approximately at 305 ls. Sandia penetration data, based on diffuse backillumination imaging, are also considered, but the available time information is only relative to the apparent SOI, which we assumed to be the same as in the X-ray dataset. Numerical penetration is calculated based on a radiographic projection of density distribution (Battistoni et al., 2014a), in analogy with the experimental methodology. The numerical predictions and the experimental data match quite well in terms of both start timing and slope of the liquid penetration curve.
The main result of these analyses is that if the sac is initially filled with gas, the liquid SOI, which occurs several tens of ls after the start of needle movement, is in good agreement with the experimental evidence. This delay is in the order of 100 ls, and it is compatible with the duration of the first slow-rising part of the needle movement. It is evident that any other types of sac initialization, e.g. partially filled or fully filled with liquid, would cause a premature SOI, which would be severely off compared to the experimental data.
We now switch the focus to the analysis of the flow features. The evolution of the sac filling process is displayed in Figure 18, with cut planes at different timings showing contours of Liquid Volume Fraction (LVF) and of residual gas fraction (calculated as transported passive scalar). The dynamics depicted in the sequence of images shows that at t = 210-220 ls, that is 10-20 ls after the start of the simulation, the liquid coming from the upstream high pressure   region begins to fill the sac. The liquid exiting from the gap area diffuses and occupies the space available. The predicted outlet velocities from the needle-seat gap are in the order of 150 m/s (in this initial phase). At the nozzle outlet, the residual gas jet emerges in a sort of 'mushroom' shape. Moving forward in time, at t = 250 ls a substantial amount of residual gas is still present in the sac. A small fraction of void is created at the needle-seat exit, after the expansion. Because of the off-axis location of the orifice, the liquid front does not proceed symmetrically. At this point in time, only gas phase is detected at the orifice outlet. Fine turbulent structures are observed in the gaseous jet.
At t = 270-290 ls, the residual gas concentration in the outlet chamber decreases, signifying that some high speed liquid has approached the hole exit and is passing through the residuals cloud. The near exit region, where fuel arrives, shows less turbulent structures because of the high density of the liquid phase.
At the time of t = 310 ls the liquid injection has newly started, and the jet is in an early developing stage. The features of an highly under-expanded jet are visible at the orifice outlet, with the formation of the Mach disk.
In addition, at this time liquid tip reaches the residual gases tip. Later in time, at t = 330 ls the fuel jet penetrates well beyond the gas plume and residuals are pushed away and diffused into the ambient air.
From a physical perspective, it is evident that the fuel mixes with the residual ejected gases first and with the ambient air later. It is expected that this mixing process can change the vaporization rate of the first part of the injection. In particular, since residual gases might be hotter than the ambient gases in a real engine combustion chamber, the initial liquid tip is expected to vaporize faster and produce a shorter liquid length. Further investigations are needed to support this interpretation, even if literature data do show a reduced liquid length for the first 100 ls after SOI for spray A in hot chamber conditions . Even other injectors reported in the literature do behave in the same way under hot conditions: e.g., spray H experiments with n-heptane show that initial liquid length is shorter than the steady value (Som et al., 2012a). The simulations shown here are restricted to a cold environment, therefore the observations based on the model can only refer to the mixing process and not to the vaporization rate. Simulations of vaporizing jets are planned for the near future, to further explore these physical processes. Figure 19 displays the evolution of turbulent structures both inside and outside the nozzle. Two grids are compared, the main features have already been shown in Table 4. Iso-surfaces mark a suitable Q-criterion threshold that identifies vortices and maximizes their rendering. Iso-surfaces are also colored by liquid mass fraction, Y liq , that gives an immediate perception of where the fuel is actually located at each time.

Grid Sensitivity
Looking at the outlet region, we notice that generally the two grids give a similar representation of the vortical structures, but the fine grid is able to provide more details (the residual gas jet tip at t = 250 ls, and the very near exit region at t = 340 ls).
At the bottom of Figure 19 a close-up view of the orifice is also provided. This comparison shows quite clearly that the coarse grid, using 5-10 lm resolution, is not adequate for capturing the smallest turbulent structures in the holes. On the contrary, the fine grid, that uses 2.5 lm cells near the wall and 5 lm cells in the core of the hole, is able to capture small turbulent structures in the orifice itself.
An interesting additional feature that the simulations capture is the formation of a shock wave at the tip of the liquid jet during the early stage of jet penetration, as the fuel penetrates in to the chamber. A snapshot is shown in Figure 20 for both grids. Density contours are plotted in the restricted range 20-30 kg/m 3 for better visualization. The prediction of  Table 4. shock wave occurrence matches quite well with the experimental X-ray radiographic projections discussed by Kastengren et al. (2014), as they found shock waves at the early stage of the jet penetration, around 20 ls after the SOI, as it is seen in the simulations. LES predictions of shock location and occurrence timing are therefore in very good agreement with experiments. The formation of shock waves at the jet tip are due to a relatively cold environment and to the high liquid tip velocity which slightly exceeds the speed of sound of the ambient gas. Concerning the grid effect, we observe a slightly more intense shock wave predicted by the coarse grid compared to the fine grid at t = 320 ls. At t = 325 ls the situation appears to be reversed. Some minor differences in the liquid tip velocities with the two grids cause this behavior. It is also interesting to observe that the fine grid allows capturing several secondary waves located behind the main shock wave.
The mesh resolution analysis shown here is not complete, as grid converged results are not clearly achieved yet (cf. turbulent vortices inside the orifice in Fig. 19). However, from these results we believe it is possible to get some useful indications for future directions in terms of mesh requirements. In a seamless simulation of the internal nozzle flow and outward jet, the main challenges remain the strong variation in the relevant length scales involved in the various regions. If 5 lm grid resolution seems adequate to provide a reasonable description of the turbulent flow in the free jet region, the wall bounded flow that takes place inside the orifice certainly requires a much finer grid spacing, as a 2.5 lm grid size may still not be sufficient for a proper flow representation. In addition, the flow through the needleseat passage at very low lifts (below 5 lm) is also critical and constitutes an area for future investigation. It is also acknowledged that wall-models are essential for enabling LES of realistic problems at high Reynolds numbers, as in this work. An effective wall model should be able to give reasonably accurate predictions of the wall shear stress irrespective of the grid resolution near the walls. The Warner and Wengle model here used has been proven to work quite well in various applications (Ameen et al., 2015), as with different near wall resolutions the predicted wall shear stress is substantially the same. Nevertheless, inclusion of more SGS physics, like usage of non-equilibrium wall models, will be considered for future work on this subject. The total computational cost of the simulation covering the SOI with the 15.9 million cells grid (cf. fine grid in Tab. 4) is 12 days of wall clock time on 256 cores. Simulations were run using a time step controlled by a velocitybased CFL set equal to 0.2. The resulting time-step is in the order of 1.0910 -9 s to 2.0910 -9 s. The total number of time steps required to cover 150 ls of physical time was 80 000. Turbulent structures inside and outside of the nozzle. Iso-surfaces mark a suitable Q-criterion threshold that maximizes vortex visualization. Colors refer to liquid mass fraction, Y liq . Coarse grid results are displayed on the left column, fine grid on the right column.

CONCLUSION
The objective of this study is to evaluate the feasibility and the possible benefits of multiphase LES applied to the simultaneous modeling of in-nozzle and near exit flow in Diesel injectors. Transient part of injections are simulated in single-hole injectors using a two-phase three-components Eulerian mixture model. High-resolution (2.5 lm) LES are devised and carried out to capture phenomena at small needle valve lifts and on small orifices. The multiphase model accounts for compressibility of both liquid and gas phase, and it also includes phase change due to cavitation. Transient jet features at the SOI and EOI are simulated and results compare well with synchrotron X-ray imaging and radiography data. Simulation results demonstrate the potential to comprehensively understand the physics of SOI and EOI transients complementing experimental data. Part of the article focuses on the description of the fast dynamics of these processes, which is resolved in great detail. Concerning the EOI behavior, liquid jet breakup, gas ingestion, in-nozzle bubble formation and dribbles are predicted. For the SOI behavior it is reported the initial ejection of residual gases trapped in the sac. Prescribing the needle motion from a very low lift, the timing of the first injected fuel is also captured. In addition, some potential effects on the initial evaporation rate and liquid length of the spray are highlighted.
Part of the article also focuses on the mesh assessment and suitability for LES of nozzle-flow and jet region. Two criteria are considered, one based on the percentage of resolved turbulent kinetic energy, and one based on the ratio between the SGS turbulent length-scales and grid spacing. Due to the high computational cost of the simulations, the collection of statistical turbulence quantities is performed in quasi-steady conditions via temporal averaging, and not as ensemble averaging over multiple realizations. This poses some specific challenges that will be subject of new work. Instead, the length-scale resolution index allows a local dynamic evaluation of the mesh resolution that can be used either for a posteriori grid assessment or for adaptive LES in the future. The grid resolutions employed, of about 2.5-5 lm results adequate for the resolution of the external spray, but may still not suffice for the extremely challenging requirements of the internal wall-bounded injector flow.
Based on these findings, future work will be devoted to the prediction of these phenomena at higher temperatures, more relevant to engine conditions, by developing a more comprehensive phase change model which is able to account for fuel evaporation. Improved wall models, particularly in wall-regions where cavitation/roughness is present and heat transfer is relevant, will be considered as well. Lastly, an equation of state for fuel density with simultaneous pressure and temperature dependency will be implemented.

ACKNOWLEDGMENTS
The submitted manuscript has been created by UChicago Argonne, LLC, Operator of Argonne National Laboratory (Argonne). Argonne, a US Department of Energy Office of Science laboratory, is operated under Contract No. DE-AC02-06CH11357. The US Government retains for itself, and others acting on its behalf, a paid-up nonexclusive, irrevocable worldwide license in said article to reproduce, prepare derivative works, distribute copies to the public, and perform publicly and display publicly, by or on behalf of the Government. This research was funded by DOE's Office of Vehicle Technologies, Office of Energy Efficiency and

Figure 20
Shock waves at the liquid jet tip, during spray A SOI. Coarse grid, left. Fine grid, right.