Large-Eddy Simulation of ECN Spray A: Sensitivity Study on Modeling Assumptions

: In this study, various mixing and evaporation modeling assumptions typically considered for large-eddy simulation (LES) of the well-established Engine Combustion Network (ECN) Spray A are explored. A coupling between LES and Lagrangian particle tracking (LPT) is employed to simulate liquid n -dodecane spray injection into hot inert gaseous environment, wherein Lagrangian droplets are introduced from a small cylindrical injection volume while larger length scales within the nozzle diameter are resolved. This LES/LPT approach involves various modeling assumptions concerning the unresolved near-nozzle region, droplet breakup, and LES subgrid scales (SGS) in which their impact on common spray metrics is usually left unexplored despite frequent utilization. Here, multi-parametric analysis is performed on the effects of (i) cylindrical injection volume dimensions, (ii) secondary breakup model, particularly Kelvin–Helmholtz Rayleigh–Taylor (KHRT) against a no-breakup model approach, and (iii) LES SGS models, particularly Smagorinsky and one-equation models against implicit LES. The analysis indicates the following ﬁndings: (i) global spray characteristics are sensitive to radial dimension of the cylindrical injection volume, (ii) the no-breakup model approach performs equally well, in terms of spray penetration and mixture formation, compared with KHRT, and (iii) the no-breakup model is generally insensitive to the chosen SGS model for the utilized grid resolution.


Introduction
Multiphase flow mixing and evaporation phenomena play pivotal roles in various industrial and energy-related applications, such as fuel spray in engines of stationary energy and transportation sectors. Fuel spray is described by the injection of a high pressure and momentum liquid fuel into ambient gas. Mixing of the liquid fuel and ambient gas under engine-relevant conditions with the presence of turbulence is an intricate process, which requires specific considerations and simplifying assumptions in numerical simulations. In this study, different modeling assumptions with relevance to fuel spray mixing and evaporation in large-eddy simulation (LES) within the context of compression ignition (CI) engines are examined.
Mixing and evaporation process in fuel sprays has been the subject of studies since several decades ago, for instance c.f. [1][2][3][4]. These fundamental studies have described the details of multi-stage mixing, droplet transport, atomization, and evaporation of fuel sprays. In a typical fuel spray, the liquid core has similar properties to the potential core of a single-phase jet [5][6][7][8] and it is under a strong shear stress. This strong shear results in primary breakup or atomization, leading to the formation of irregular liquid ligaments along the surface of the liquid core. The rate of primary breakup tends to control the liquid core length, which is a function of the nozzle diameter (d n ), and the spray-to-ambient relative density [8]. In engine-relevant fuel sprays, such as diesel spray, an elevated ambient gas density or pressure results into a shorter liquid core and a rapid atomization process [9]. After the onset of primary atomization, the resulting droplets along the surface of the liquid core are intrinsically unstable and may be subject to secondary breakup due to aerodynamic forces. The breakup pattern is governed by droplet Weber number (We) defined as We = ρd p u 2 rel /σ, where ρ, d p , u rel , and σ denote the gas density, droplet diameter, droplet-to-gas slip velocity, and droplet surface tension, respectively. During the secondary breakup phase, droplets breakup until reaching a stable We below 12 [10,11], after which they start to evaporate.
In droplet-laden (multiphase) fuel sprays, droplet atomization process plays an important role in mixing and evaporation. However, from the numerical modeling perspective, resolving near-nozzle length scales is challenging and there is insufficient experimental data available. We note that there have been some attempts to numerically resolve the primary breakup of dense sprays using e.g., the level-set method (LSM) [12] or volume of fluid (VoF) [13][14][15][16] approach, in addition to further strategies reviewed in [17]. However, in computationally expensive simulations such as reactive fuel sprays or even evaporative sprays with relatively long penetration lengths, it is a common practice to model the dispersed phase in the Lagrangian framework without resolving the liquid core and primary breakup. In these simulations, the gaseous phase is resolved in the traditional Eulerian framework. In particular for LES, a coupling between the Lagrangian particle tracking (LPT) and LES is usually utilized to model the atomized Lagrangian fuel droplets which lie in the subgrid range and to take their interactions with the surrounding carrier phase into account [18,19]. In addition, to avoid extreme computational costs, one common approach is to implicitly consider the primary breakup by introducing Lagrangian droplets with particular injection parameters and size distributions into a small cylindrical injection volume which is larger than actual fuel nozzle dimensions [20,21].
According to literature, LES/LPT approach has been studied in non-reacting fuel sprays by several researchers, e.g., [22][23][24][25], wherein droplet secondary breakup model was utilized under various spray and turbulence conditions. Moreover, limitations of the LPT approach in specific spray applications were addressed in [26]. Additionally, Vuorinen et al. [27][28][29][30] performed systematic studies using implicit LES (ILES) to model subgrid scales (SGS) for various droplet parameters, wherein LPT computational parcels were introduced into the quiescent gas through an injection disk area (not a cylindrical volume). Particularly, they demonstrated significant effects of droplet dispersion in terms of Sauter mean diameter (SMD), volume fraction, and Stokes number (St) on the spray shape, mixing quality, and global flow characteristics, compared to single phase jets. Additionally, in [31], the influence of droplets breakup due to the resonance were taken into account. In high ambient density conditions, Kaario et al. [32,33] studied fuel spray behavior using LES/LPT with a cylindrical injection volume approach and one-equation eddy diffusivity SGS model (K sgs -model), concluding that SGS effects on particle dispersion are small. Furthermore, Kaario et al. [20] investigated mixing and evaporation of an engine-related selective catalytic reduction (SCR) system, modeled using LPT and a hybrid methodology using LES and Reynolds-averaged Navier-Stokes (RANS) models to account for spray-wall interactions. However, these numerical studies have investigated arbitrary or specific system configurations due to the lack of a unified benchmark data.
The Engine Combustion Network (ECN) collaboration provides well-documented baseline experimental/computational data for model validations in engine-related research. With relevance to this study, in ECN Spray A, a liquid n-dodecane is injected through a nozzle hole of diameter d n = 90 µm at high injection pressure (150 MPa) into hot quiescent ambient gas of temperature T = 900 K and density ρ = 22.8 Kg/m 3 . Under non-reacting conditions and by assuming an inert ambient, several numerical studies utilized the LES/LPT approach against the ECN Spray A data. For instance, Wehrfritz et al. [34] and Xue et al. [35] investigated different secondary breakup models, including the Kelvin-Helmholtz Rayleigh-Taylor (KHRT) and the enhanced Taylor analogy breakup (ETAB) models, and different SGS models, including Smagorinsky, K sgs -model, and ILES, respectively, for various grid resolutions. Tsang et al. [36] developed and validated a model that predicts the spray sub-grid dispersion velocity. Senecal et al. [37] additionally studied LES multi-realization uncertainties on global and local spray characteristics. Moreover, uncertainties in fuel properties, ambient conditions, and injection parameters on spray-related objective functions were investigated via the global sensitivity analysis method by Pei et al. [38]. Further studies with relevance to spray modeling parameter sensitivity can be revised in [39,40]. More recently, Kaario et al. [21] investigated fuel properties effects on the liquid spray formation by replacing n-dodecane in Spray A with other fuels. Seeking a more accurate modeling approach, Koukouvinis et al. [41] employed advanced thermophysical closure models to account for capturing complex details of the flow. Following a similar pathway, Matheis and Hickel [42] and Yang et al. [43] both accounted for real gas thermodynamics to investigate the co-existence of transcritical effects in the spray, while showing that surface tension effects still persist in the sub-critical regime. Under reacting conditions, various works employed different combustion models along with the LES/LPT approach to investigate ignition-related phenomena such as flame structure and stabilization [44][45][46][47][48][49][50], radiative heat transfer [51], as well as ignition modes and dynamics [52]. Additionally in [53][54][55][56][57], the cylindrical injection volume filtering concept was utilized together with the KHRT and ILES models, as discussed in [34].
Besides the aforementioned numerical studies, several experiments with relevance to the ECN Spray A conditions were conducted. For instance, the works by Pickett et al. [58], Payri et al. [59], and Kastengren et al. [60] delivered quantifications of various global spray metrics in the near-nozzle region using optical visualization techniques. Focusing on the internal nozzle flow effects, Daly et al. [61] and Payri et al. [62] considered, respectively, the nozzle temperature and geometry as design parameters. Also, Agarwal and Trujillo [63] demonstrated the effects of nozzle imperfection on generating stronger radial velocity components with shorter liquid core and faster primary atomization. With relevance to the high temperature and pressure environment in Spray A, Falgout et al. [64] reported evidences on the existence of supercritical effects in the mixing layer, thereby supporting the hypothesis which was previously investigated in [65][66][67] and theoretically established in [68,69]. On the other hand, recent observations by Crua at al. [70] demonstrated that droplets surface tension is yet clearly present, which is consistent with numerical observations in [42,43] as discussed above, hence the existence of spray classical atomization and evaporation processes. Additionally, angular dispersion effects, in terms of spreading angle variation, on mixing and evaporation processes were investigated by Jung et al. [71]. Spray shape and flame structure under reacting conditions were also discussed by several authors such as [72][73][74].
In the present numerical study utilizing LES, and following the baseline operating conditions of the ECN Spray A, we shift our attention towards some of the modeling assumptions that have been extensively utilized in recent numerical simulations while their details have never been elaborated. In particular, motivations of the present research are as follows. First, the volume filtering approach has been frequently utilized in literature to implicitly mimic the primary breakup in the near-nozzle region, but topological details of the cylindrical injection volume and the spray sensitivity to variations of the volume dimensions is missing. Second, although utilization of the no-breakup model approach for droplets secondary breakup has been reported in literature, a comparison between a classical breakup model (such as KHRT) and the no-breakup model on the spray liquid, vapor penetration and droplet evaporation has not been performed. Third, a study on the subgrid scale modeling effects when the no-breakup model approach is employed is absent from literature. Accordingly, the overarching goal of this study is to perform a systematic analysis to understand the sensitivity of simulations to such modeling parameter variations for the non-reactive baseline ECN Spray A. Therefore, the main objectives of the present study are to: (i) investigate details of the near-nozzle modeling for evaporative ECN Spray A, using the injection volume concept, and report the effect of cylindrical injection volume dimension variations on the spray characteristics, (ii) assess the importance of droplet secondary breakup model, by comparing the no-breakup model against the well-established KHRT model, and (iii) study the effect of SGS modeling (Smagorinsky, K sgs -model, ILES) on the mixture formation, when the droplet no-breakup model is utilized.

Gas Phase Governing Equations
Turbulent gaseous flow is described using the compressible Navier-Stokes equations in the Eulerian framework. The Favre filtered LES formulations for conservation of mass, species mass fractions, momentum, and energy are given as, ∂ρ ∂t where ρ, u i , p, Y k , h s , andτ ij denote density, velocity component in x i direction, pressure, chemical specie k mass fraction, sensible enthalpy, and conventional viscous stress tensor assuming Stokes' hypothesis for isotropic Newtonian fluid, respectively. The overbar and tilde (∼) are filtering operators denoting unweighted and density-weighted ensemble averages, respectively. Unity Schmidt number is assumed for all species, hence the mass diffusivity D = µ/ρ, wherein µ denotes the viscous diffusion rate of the gas mixture. The total enthalpy h t is defined as the summation of h s and specific kinetic energy. The filtered source terms S ρ , S u i , S Y k , and S h correspond to coupling between liquid and gaseous phases with respect to mass, momentum, species and energy, respectively. Equations (1)-(4) together with the thermodynamic and caloric state equations for ideal gas, are numerically solved using the finite volume method (FVM) within the OpenFOAM-6 framework [75]. Temporal advancing is based on a transient, second order implicit scheme. The coupling between pressure and momentum is resolved via the reacting PISO algorithm [76]. Pressure and diffusion terms are discretized by a second-order central scheme, whereas the interpolation of the convective flux terms are relevant to the subgrid terms implementation in LES, i.e., turbulence SGS modeling, which is further discussed in the following section. Further details on the source terms are provided in Spray Modeling section (Section 2.3).

Subgrid Scale Modeling
The present study employs three different closure models, namely the Smagorinsky model, one-equation eddy diffusivity model (K sgs -model) and ILES model, for the unresolved turbulent subgrid scales. The residual SGS term for a quantity φ ∈ {u j , Y k , h s } is defined as τ sgs ij = u i φ −ũ iφ , which is incorporated for each mentioned variable in Equations (2)-(4). A brief discussion of the implemented SGS models is given in the following.

Smagorinsky Model
The Smagorinsky model [77] is one of the most popular SGS models in turbulent flows. The residual SGS stress tensor in the momentum conservation equation (Equation (2)) is decomposed into isotropic and deviatoric parts, with D sgs denoting the SGS eddy diffusion coefficient, ) is the resolved strain rate tensor. The approximation in Equation (5) is based on the eddy diffusivity assumption, postulating a linear relationship between the SGS shear stress (τ sgs ij − 1 3 τ sgs kk δ ij ) and dev(˜ ) ij . The subgrid scale kinetic energy, k sgs , is defined as: In OpenFOAM, the subgrid scale diffusivity is computed as where C k = 0.094 is a default value for the model constant, and ∆ is the filter width, computed as the cubic root of the cell volume, which is a suitable choice considering the isotropy of the utilized hexahedral mesh. In Smagorinsky model, the SGS kinetic energy k sgs is computed by assuming a balance between the SGS energy production and dissipation, i.e., : where (:) denotes the double dot product, and C = 1.048 by default. The convective fluxes of the SGS terms as well as the resolved terms are herein discretized by a second-order central scheme.

One-Equation Eddy Diffusivity Model (K sgs -Model)
As the name suggests, the present model utilizes the eddy diffusivity approximation, similar to Smagorinsky model, and the SGS stress tensor τ sgs ij is evaluated from Equation (5). The difference between Smagorinsky model and the one-equation eddy diffusivity model (K sgs -model) is that while the Smagorinsky model assumes a local equilibrium between the SGS energy production and dissipation to compute k sgs , the K sgs -model solves a transport equation for k sgs from which it accounts for the SGS energy production, dissipation, and diffusion [78]. The transport equation for k sgs is given as wherein the terms on the left-hand side denote temporal and spatial (convection) variations, and those on the right-hand side denote diffusion, production, and dissipation, respectively (from left to right).
We note that as in the Smagorinsky model, a second-order central difference scheme is utilized herein for numerical discretization of the convective flux terms, including the resolved and SGS terms, in addition to the k sgs fluxes emerging from the present model.

Implicit Model
As previously indicated, both the Smagorinsky and K sgs -models introduce an additional dissipation term to the system via explicit SGS diffusivity D sgs . However, in the present implicit model, such a dissipative contribution is implicitly introduced and it is controlled locally by choosing a dissipative numerical scheme for discretization of convection terms. Following the previous works on spray A [34,[53][54][55][56][57], ILES is implemented through a non-linear flux limiter termed as Gamma scheme [79].
The control parameter of the flux limiter is set to 0.3 for discretization of the momentum convection, whereas a parameter value of unity is considered for scalars convection to ensure a bounded total variation diminishing (TVD) solution.

Spray Modeling
Diesel engine sprays are characterized by a short liquid core and rapid atomization [9]. Hence, it is more efficient to treat the dispersed liquid phase in the Lagrangian framework, without resolving the length scales below the nozzle diameter in the near-nozzle region. Moreover, the dilute suspension feature of liquid droplets allows for a two-way coupling between the two phases without inter-particle collisions [80]. From a modeling perspective, the computational LPT parcels (single phase n-dodecane) are injected from a cylindrical volume geometry, as depicted in Figure 1, with an initial droplet size distribution (IDSD) to implicitly consider the droplets after primary atomization. The injection cylindrical volume is detached few millimeters from the upstream wall of the ECN vessel, hence allowing an upstream hot ambient entrainment. In addition, the flow in the under-resolved near-nozzle region (maximum cells/injection cylinder diameter ≈ 14) will become 3-dimensional downstream of the nozzle, with emerging vortical structures that allow further air entrainment into the spray. By considering such an injection volume approach, the near nozzle singularity is avoided with a computational gain. The injection conditions, in terms of pressure and mass rate profile, for ECN Spray A are assumed to hold in the injection volume. The Eulerian source terms S ρ and S Y k from Equations (1) and (3) are computed from the droplet evaporation model [81] utilizing the ideal gas assumption. The momentum source term S u (Equation (2)) is evaluated from sphere drag forces with a coefficient based on [82]. The energy source term S h (Equation (4)) is based on the Ranz-Marshall correlation for heat transfer [83].
It is worth noting that the employed injection volume approach models a dilute spray and it is considered as the injection source. Moreover, it assumes that the nozzle is located within a fixed distance after the dense spray region, which is rather short under Spray A conditions. Therefore, it is expected that altering the location of the injection volume does not affect the global spray behavior, although it could have a significant upstream impact depending on the detachment distance from walls.

Secondary Breakup Modeling
The present study employs two different models for droplet secondary breakup, namely the KHRT model and the no-breakup model. The KHRT model is a hybrid approach comprising the Kelvin-Helmholtz (KH) wave model and the Rayleigh-Taylor (RT) accelerative instability normal to the droplet surface. The KH mechanism considers droplets as liquid jets that suffer stripping after being injected into a gaseous environment, whereas the RT mechanism is driven by the density variation along normal directions of the liquid-gas interface. Considering the KH model, on the basis of the perturbation growth rate Ω KH and the wavelength Λ KH corresponding to the fastest Ω KH , a characteristic breakup time τ KH can be defined as where r denotes the radius of initial droplets. The stable diameter is then computed as d c = 2B 0 Λ KH . The KH model constants B 1 and B 0 control the breakup time and wavelength, respectively, and the choice of such parameter values is addressed later. The size of injected droplets under ECN Spray A conditions is typically within few micro-meters, thereby, the corresponding KH wavelengths for the droplets surface waves are separated from grid scales, which are much larger (e.g., 62.5 µm in this study). Considering the RT model, assuming a linear disturbance growth, a growth rate Ω RT and wavelength Λ RT can be determined. The RT breakup time is then obtained by the multiplicative inverse of the growth rate, i.e., τ RT = C τ 1 Ω RT , where C τ is a correction factor to delay breakup under certain conditions. Numerically, the initial droplet size is modeled within the size range of d p,max = d n /5 = 18 µm and d p,min = d n /90 = 1 µm. Such a droplet size range is chosen large enough to avoid immediate evaporation of the droplets. The IDSD is modeled according to the stochastic Rosin-Rammler distribution function [84]. Another approach to model the near nozzle region is the no-breakup model, which assumes that the droplets have already undergone primary and secondary breakup, and they are injected from V cyl with a constant initial droplet size. In fact, utilization of the no-breakup model idea originates from the observation in [34] for Spray A conditions, wherein droplets were noted to undergo breakup only up till ≈ 20d n distance from the nozzle exit, after which We becomes too small for droplet breakup to take place. Moreover, numerical evidence reported that SMD of the droplets reduces rapidly till below 1 µm [34,85], and a similar observation was experimentally noted in [86,87] as well. It is worth noting that these observations could be clarified in the light of experimental [64][65][66][67]70] and numerical [42,43] evidences showing the co-existence of transcritical and supercritical states in Spray A conditions in which they affect the spray atomization process. Following the previous works by Kaario et al. [20,21], the droplets are herein injected from V cyl with a constant initial size for all droplets which corresponds to a stable Weber number (We < 12), hence no further breakup, and only evaporation is considered.
In both models (KHRT and no-breakup), the implementation of the LPT method in OpenFOAM incorporates the parcel approach, which groups physically similar droplets into a parcel. Therefore, such an algorithm reduces the computational cost, efficiently. The total number of injected parcels (about one million for an injection duration of 1.5 ms) is chosen to correspond to an average parcel mass of 3.6 × 10 −6 mg and a total injected mass of ≈ 3.7 mg which follows injection rate recommendations by ECN [88], while a minimum parcel mass of 6 × 10 −7 mg is assumed. The parcels are randomly distributed within an injection cylinder of dimensions D cyl and L cyl and a half spray opening angle of α = 5 • . Such a chosen value for α is based on numerical testing which revealed, as demonstrated in Appendix A, that global spray metrics are rather insensitive to small opening angle variations.

Numerical Setup
In the present work, simulations of liquid n-dodecane spray injection into an inert quiescent gaseous environment under the ECN Spray A operating conditions as noted in Table 1 are carried out. The computational domain, depicted in Figure 2, represents a circular cylinder of volume equivalent to the combustion vessel utilized for the ECN Spray A standard experiments by Sandia [9]. The domain is discretized into a background hexahedral mesh with a cell size ∆ = 1000 µm (R1). Two additional spatially refined regions are considered, the coarse region with a grid size of 250 µm (R2), and the dense region with a 62.5 µm resolution (R3).  Three main studies are considered herein. The first study uses KHRT secondary breakup model together with ILES to investigate the sensitivity of the near-nozzle modeling using the injection cylinder concept. In particular, five test cases are considered to vary the axial and radial dimensions of the injection cylinder. In the first three cases, the cylinder diameter (D cyl ) is varied while holding a constant volume (V cyl = 25πd 3 n ). In the remaining two cases, the cylinder length L cyl is altered with a constant diameter (D cyl = 5d n ). In the second study, the no-breakup model is compared to the KHRT secondary breakup model, with a fixed V cyl and using ILES. Moreover, two cases investigate the variation of KHRT characteristic breakup time (τ KH ), with the model constants presented in (Table 2) to show the possible tunability of such parameters for the spray validation. The third study utilizes the no-breakup approach together with a fixed V cyl , but varying the SGS model between ILES, Smagorinsky, and K sgs -model. A summary of the considered cases in these three studies are listed in Table 3.
It is worth mentioning two important notes. First, the KH model constant B 1 in Table 2, which controls the stripping rate, does not hold a recommended value by default and several works have tuned this parameter according to the problem setup [89,90]. However, in diesel sprays, the value B 1 = 40 has been extensively used [34,[91][92][93][94], hence the KHRT model utilizing such parameter value is herein referred to as the default KHRT model. Second, parameter values of V cyl = 25πd 3 n in Table 3, which correspond to cylinder dimensions D cyl = 5d n and L cyl = 4d n , are selected based on the authors' preference. However, these selected values are relevant to the a priori knowledge of the spatial information on primary breakup and droplets formation. The ratio N cell /D cyl in the table represents the number of grid cells within D cyl . All the simulations are carried out for evaporative sprays under non-reacting conditions for 2.0 ms simulation time. In the present numerical implementation, liquid penetration length is defined as the maximal axial distance to the location of 0.1 % of the liquid volume fraction averaged over a volume of 1 mm characteristic length. Similarly, vapor penetration length denotes the axial distance from the nozzle to where fuel mass fraction is 0.1 %. Validations for liquid and vapor penetration lengths are realized by comparing against the available experimental data for the ECN Spray A by Sandia [95][96][97], utilizing high-speed imaging techniques such as Mie scattering and diffused back-illumination (DBI) for the liquid length measurement, and the schlieren imaging for vapor length measurements. Experimental data utilized for fuel mixing validations are reported in [98,99] wherein Rayleigh scattering technique is employed.

Study-I-Cylindrical Injection Volume Dimensions
The first study investigated the effect of varying dimensions of the cylindrical injection volume, denoted by D cyl and L cyl , on the spray mixing and penetration lengths. As shown in Table 3, the first three cases A1-3 consider a constant V cyl with D cyl varying between 2.5d n and 10d n . Figure 3 depicts cutplanes of the mixture fraction (Z) field for each case at t = 0.25, 0.5, and 1.0 ms (top to bottom) as well as the mean Z field averaged over the time interval of 1.5-2.0 ms (the bottom row). Two general observations can be noted: (1) the jet width downstream the nozzle became larger with larger D cyl , and (2) the jet penetration length became shorter with larger D cyl as can be seen from the first two time samples. To further quantify mixing and evaporation in the mentioned sprays, liquid and vapor penetration lengths were plotted in time, and compared with the respective ECN Spray A data, as shown in Figure 4. It is observed that the introduced V cyl dimensions in cases A1 and A2 resulted in over-predicted liquid and vapor penetration lengths throughout the entire injection duration. Moreover, a slightly longer jet was noted for case A1, with the smallest D cyl , in comparison with case A2, while the liquid length was observed to be equivalent between both cases. Finally, when D cyl was enlarged to 10d n in case A3, a significant reduction in both liquid and vapor penetration lengths was observed, with a better convergence to the ECN Spray A data. Therefore, increasing D cyl was observed to reduce the jet penetration and it led to a faster liquid evaporation process.  The previous observation, noted from Figures 3 and 4, suggests the following. First, under the same injection parameters (injection pressure, V cyl , and total injected mass), the induced mass and momentum source terms, i.e., S ρ and S u i in Equations (1) and (2), are also constant. Second, by radially extending the injection cylinder with a fixed volume, the droplets momentum distribution varies accordingly. In particular, for a larger D cyl , the axial magnitude of the momentum becomes weaker when both V cyl and S u i are fixed. As a result, the jet axial penetration becomes shorter. The same holds for the liquid phase, where the collective mass flow rate of the droplets in the axial direction for larger D cyl and shorter L cyl results in a smaller St of the cloud droplets in the axial direction, hence a shorter liquid length and a faster evaporation.
To understand the relationship between D cyl and the axial magnitude of the momentum with more details, the axial streamwise component of the mean velocity, averaged over t = 1.5-2 ms, is plotted in Figure 5 along the axial direction, normalized by d n , for cases A1-A3. In particular, a Favre averaging procedure is performed along the circumferential data points within a radial distance of ∆ = 62.5 µm. The spatial averaging procedure is done in both the radial and the circumferential directions, for each axial location of a grid size interval away from the nozzle. Here, as expected, a monotonic drop in the peak streamwise mean velocity is observed in the near-nozzle region (x/d n ≈ 30-40) from case A1, with the smallest D cyl and the highest peak, to case A3, with the largest D cyl and the smallest peak. The present plot is a clear evidence that variation of D cyl has a significant impact on the axial magnitude of the spray momentum, and thereby on the liquid and vapor penetration lengths. The previous discussion elaborates on the significant influence of the radial change of the injection cylinder dimensions with fixed V cyl on the liquid and jet penetration lengths. We complemented this study by further investigating the spray metrics when V cyl is axially varied. In particular, two additional test cases are herein considered, A4 and A5 from Table 3, where L cyl was varied between 2d n and 8d n while D cyl = 5d n was fixed, hence V cyl was axially prolonged. Figure 6 depicts temporal evolution of liquid and vapor penetration lengths for cases A4 and A5, as compared to case A2 with the same D cyl and an intermediate value of L cyl . The figure shows a reduction in the jet length when V cyl was axially prolonged, while a mild reduction in the liquid length was attained. This observation can be referred to the droplet cloud distribution which became more relaxed with increased L cyl , hence a weaker momentum intensity since S u i became distributed over a larger volume. According to the presented discussion, it is concluded that a radial change of V cyl had a stronger influence on both liquid and vapor penetration lengths when compared to the axial variation of V cyl . This behavior corresponded to the direction of D cyl change being aligned with the direction of the largest gradients in the cross stream direction, hence enhancing the mixing process, while smaller changes in scalar gradients were noted in the axial direction. Overall, varying D cyl with a fixed V cyl impacted the axial magnitude of S u i , which played a significant role on the spray metrics.
According to these results, the injection volume characteristics, particularly the associated spatial parameters, play a significant role on the observed liquid length and vapor penetration. The relatively low amount of cells in the radial direction, c.f. cases A1-A3 in Table 3, leads to a resolution sensitivity in the near-nozzle region. Here, a reasonable compromise was observed with D cyl = 5d n , having ≈ 7 cells along the D cyl while avoiding the under-prediction of the vapor penetration and the relatively large jet width in the near-nozzle region with D cyl = 10d n . Therefore, for the remainder of this work, the case D cyl = 5d n with a constant V cyl = 25πd 3 n is selected.

Study-II-Breakup Model Comparison
The second study investigated spray characteristics when the classical KHRT breakup model was utilized with two different model parameter sets (see Table 2). In addition, the KHRT model was compared to an approach without a breakup model, as described in Section 3. Accordingly, the considered test cases correspond to a slower atomization process (B1), a faster atomization process (B2), and a complete atomization (B3) of the spray.
Before progressing with analyses, a comparison of liquid droplet size distributions in the considered cases was depicted through the global SMD temporal evolution within the injection duration, Figure 7a, in addition to the probability density function (PDF) of droplets diameters at time 0.5 ms after the start of injection (SOI), Figure 7b. First, the non-fluctuating nature of the temporal droplet size profile suggested that the corresponding PDFs did not alter their mean distributions during the entire injection duration. Second, it is observed that the average droplet sizes doubled between the three cases, particularly when moving from a complete atomization assumption in case B3 (0.5 µm with constant droplet size) to the droplets breakup with a faster stripping rate assumption in case B2 (smaller B 1 value in Equation (10)) until the slower stripping rates in case B1. In addition, from the reported experimental [87,100] and numerical [101] works, the measured SMD (≈ 1 µm sufficiently far from the nozzle) was found to be in close agreement with cases B2 and B3, compared with case B1. As a first analysis of the breakup model-or rather the droplet size-comparison, Figure 8 depicted liquid and vapor penetrations considering the respective three test cases B1, B2 and B3 in Table 3. The over-predicted spray penetration lengths could be adjusted by controlling the KHRT breakup time. In particular, a faster stripping rate implied a lower droplet St and faster evaporation, with a better convergence to the ECN Spray A experimental data. More importantly, when the no-breakup model was used in case B3, the liquid penetration was observed to be similar to the modified KHRT model (case B2) and both cases presented a very good agreement with the experimental data. To investigate the variation of mixture formation off the axial centerline, the radial distribution of Z at 1.5 ms after SOI was plotted at a representative location of 17.8 mm away from the nozzle. The field data points were interpolated and circumferentially averaged for the mentioned location. In addition, such an axial location was rather close to the liquid penetration in case B1 and thereby significant effects on the vapor phase were expected. As depicted in Figure 9, for case B1, there was a high concentration of Z near the axial centerline which then diffused away from the nozzle, i.e., an over-prediction compared with the ECN Spray A data. Such an over-prediction near the axial centerline was significantly improved in cases B2 and B3, i.e., with tuning B 1 parameter to a lower value or with the no-breakup model, respectively. Again, this observation can be justified from the droplet St perspective. In particular, in case B1 with a large B 1 value, there was a slower KH stripping rate and delayed breakup, hence larger droplets as demonstrated earlier. As a result of the large droplets formation, the corresponding St of the droplet cloud became larger and thereby a stronger influence on the turbulence with longer penetrations was achieved which, in turn, affected axial vapor spray concentrations. When a faster stripping rate was assumed (B 1 = 7) in case B2, the smaller obtained St led to a stronger impact on the droplets cloud by turbulence, hence more fuel was diffused off the centerline with an improved mixing with ambient. The complete atomization of droplets (no-breakup model in case B3) resulted in a radial profile which was rather close to that in case B2, thereby closer mixing characteristics between both cases were recognized within this spatial region.  Table 3. Data points are interpolated and circumferentially averaged for the considered axial location. Solid black lines and filled regions denote the ECN experimental counterparts of Z and its uncertainty, respectively.
After discussing the breakup model effects on spray penetrations and on radial distributions of Z, a further analysis is performed with relevance to the droplet heating and vaporization. As depicted in Figure 10, scattered data points of the temperature field are plotted in Z space for the spray region (conditioned by Z > 10 −4 ) at time 0.5 ms after SOI, and overlayed with conditional means (red curves) and standard deviations (filled regions). Two main observations can be herein noted. The first observation is the occurence of low temperature data points in case B1 within lean regions of the spray. This can be explained by the large droplet sizes, c.f. Figure 7, for case B1 with relatively long liquid lengths which implies stronger cooling effects away from the nozzle. Accordingly, having smaller droplet sizes in cases B2 and B3 would enhance mixing with the hot ambient and thereby the tendency of those lean regions to rapidly heat up, i.e., low temperature data points disappear, and approach the ambient temperature value. Interestingly, this breakup model-induced cooling effects could play an important role in reacting simulations, when most reactive zones and subsequently high temperature ignition kernels are realized within the spray periphery to identify the ignition delay. Following previous notes, the second observation is the larger Z magnitude in the case B3 with no-breakup model (smallest droplet sizes) which implies richer mixtures due to fast evaporation, hence the tendency to have a narrower Z distribution with larger droplet sizes, as noted in the case B1. It is important to mention that these observations were consistently noted afterwards throughout the injection duration, and only the data at 0.5 ms after SOI is herein presented for brevity. Figure 10. Scattered data points of temperature fields in Z space for spray region (Z > 10 −4 ), overlayed with conditional means (red curves) and standard deviations (filled regions) for cases B1-B3 in Table 3.
It is worth acknowledging that ECN Spray A conditions may support a supercritical state [42,43,[64][65][66][67]70], which is not directly accounted for in the present study. In fact, the rationale behind choosing small droplet sizes for this study is two-fold: (i) a modeling perspective, to match the ECN data on global spray metrics, and (ii) a physical perspective, particularly the transcritical/supercritical effects which indicate rapid atomization and phase change. The analysis in Study-II concludes that the no-breakup approach performs equally well in terms of predicting global Spray A metrics, compared with the KHRT model with tuned parameter. Moreover, from a numerical perspective, a no-breakup model approach offers parameter-free framework which only considers initial, constant, droplet size assumption, hence it avoids further uncertainties compared with the KHRT which require tuning various model constants. Based on that, we resume the analysis without a droplet breakup model in the next section for Study-III.

Study-III-Subgrid Scale Model Comparison
The third study investigates the spray characteristics when the SGS model is varied between ILES, Smagorinsky, and K sgs -model, while utilizing the no-breakup model and a constant V cyl . Figure 11 depicts spray liquid and vapor penetrations for cases C1, C2 and C3 in Table 3. Overall, it is shown that all the mentioned cases with different SGS models perform equally well on predicting the liquid length, while a slight discrepancy is observed for the jet length, in comparison with the experimental data. We note that the similarity of results with different SGS approaches is related to the considered grid resolution of 62.5 µm with which most of the important length scales are already resolved. In order to assess the spray mixing characteristics in the mentioned cases, first, the PDF of the mean mixture fraction was plotted in Figure 12. The sampled data considered only the spray region, which was identified by grid cell centers of value Z > 10 −5 , and the PDF bins were weighted by normalized cell volumes. Here, the PDF analysis shows that the considered test cases with varied SGS models demonstrated equivalent distribution of the mean Z. This observation suggests that SGS effects were insignificant on the spray mixture characteristics since the measured quantity essentially depended on the resolved scales, and it was a function of the scalar mass transfer corresponding to droplet evaporation. To further analyze the spray mixing with direct relevance to SGS, statistics of the scalar dissipation rate, defined by χ = 2D eff ∇Z 2 where D eff denotes the effective diffusivity, was examined by PDF calculations at t = 1.5 ms in dense and coarse mesh resolutions of the spray, as shown in Figure 13. The dense and coarse mesh regions of the spray were conditioned by grid cells of value Z > 10 −5 intersecting with those holding a cell volume below or above a threshold of 1.1 × 62.5 3 µm 3 , respectively. At this point, before discussing the plots, it is important to recall the following. First, the ILES model does not supplement system mass diffusivity D with any closure terms for turbulent subgrid scales (i.e., D sgs = 0). Second, χ represents a molecular mixing phenomenon which occurs mainly at the finest turbulent scales. Therefore, in contrast to kinetic energy dissipation, most of the scalar dissipation is usually modeled and not resolved [102]. Considering previous notes, it is necessary to choose a proper D eff for computing χ such that the results of implicit and explicit SGS models are comparable. Accordingly, we first considered the effective diffusivity as the resolved one, i.e., D eff = D, for evaluation of χ in explicit SGS cases (C2 and C3). The resulting PDFs utilizing the resolved D were plotted (dotted lines) and compared against the results from ILES (the yellow continuous line) in both spray mesh regions. Then, the SGS diffusivity contribution was considered in explicit models by utilizing the effective mass diffusivity as D eff = D + D sgs for evaluation of χ and the corresponding PDFs are plotted in both figures (continuous lines). In the dense mesh region of the spray, the corresponding PDF resulted in Figure 13a show that, first, without SGS contribution (D eff = D), all the cases exhibit equivalent distribution of χ. In particular, since D held a similar value in the three cases and by considering the same mixture composition and thermal properties, the resulting PDFs demonstratde statistics of the magnitude of the resolved scalar gradients, ∇Z , which were observed to be insensitive between the cases. Second, when accounting for SGS contributions (D eff = D + D sgs ), it is observed that there was not a large difference between ILES and explicit models in most of the larger χ values, whereas notable variations were observed within a small range of low dissipation rates. Except for this narrow window of small χ values, this plot indicated that SGS effects on the prediction of χ were rather small in this dense resolution.
On the other hand, for the PDFs in the coarse mesh region of the spray, shown in Figure 13b, rather larger differences were observed when accounting for SGS contributions (D eff = D + D sgs ). In particular, explicit SGS model cases showed emerged dissipation rates over a wide range of values. As stated earlier, scalar dissipation rate was a molecular phenomenon which occurs at the finest length scales and it was usually modeled. Additionally, SGS could be generally viewed as the sum of numerical diffusion and explicit model diffusion. Considering the existing coarse resolution, less resolved length scales impacted the value of the resolved D. Furthermore, in ILES (D sgs = 0), the numerically diffusing schemes herein considered mainly affected the convection of the flow variables, without altering D. Thereby, in the present implementation, it was expected for scalar dissipation to be stronger in the explicit SGS models than the ones from ILES, for the specified grid resolution.
Overall, the combined observations from Figure 13a,b can be summarized as follows. First, without considering SGS effects, the resolved scalar gradients demonstrated similar statistics in all the cases. Second, in the dense resolution near the nozzle, SGS effects on the scalar dissipation rate were observed to be small. Third, in the coarse resolution further away from the nozzle, SGS contribution was noted to be larger with respect to ILES. However, such effects were located far away from the nozzle, and the corresponding variation in χ seemed to have an insignificant impact on the spray metrics such as vapor penetration, while the liquid length was already resolved within the dense resolution region.
In addition to the previous analysis on spray mixing characteristics, according to Khani et al. [103,104], the mixing character and efficiency seemed to be quite sensitive to the turbulent Prandtl number. In the present analysis, in explicit LES, it is noted that the turbulent Prandtl number was in the second decimal order, although the value of the resolved counterpart approached unity. Motivated by the mentioned works, a thorough analysis might be useful to conduct in the future for a further insight within the spray modeling framework.
After demonstrating SGS effects on molecular mixing and scalar dissipation rates, their effect on the dissipation of kinetic energy (E k ) is another criteria to be explored in assessing the LES overall quality. Here, E k spectra are demonstrated through a power spectral density (PSD) plot, depicted in Figure 14. The velocity signal, recorded for 2 ms simulation time with the sampling frequency 10 6 Hz, was probed at a location 6 mm (≈ 67d n ) away from the nozzle and 1 mm (≈ 11d n ) off the spray centerline. The PSD was computed based on the Welch estimate [105] with modified Bartlett-Hann window [106] along multiple overlapping segments of the signal. Power spectral density of the kinetic energy. Velocity probe is located at distance 67d n away from the nozzle and 11d n off the spray centerline. Signal is recorded for 2 ms simulation time with sampling frequency of 10 6 Hz. PSD is computed based on Welch estimate.
It is worth noting that according to Taylor's frozen turbulence hypothesis [107], wherein a linear relation between the signal's frequency and streamwise wavenumber (κ) is assumed, the dissipation pattern of E k , herein depicted, in the frequency domain can be compared with the Kolmogorov spectrum [108]. Accordingly, the inertial range of scales, which separate the large energy-containing scales from the small dissipative ones and is governed by the relation E k ∝ ε 2/3 κ −5/3 with ε denoting the E k dissipation rate, was demonstrated in the figure for the three cases with varied SGS models. The imperfect coincidence of the slope of a signal's inertial range with respect to the −5/3 power law could be attributed to the accuracy of the employed numerical schemes, which would contribute to noise contaminations in the probed signal. Furthermore, in the higher frequency region corresponding to smaller dissipative range of scales, a coincidence between the three curves implied that all the SGS models were numerically stable without accumulations in E k and enstrophy, hence reliable quality LES was achieved. However, in the lower frequency range, significant discrepancies between the three cases were noted. In particular, the Smagorinsky SGS model was observed to damp E k at a wide range of low frequencies, which corresponded to larger length scales by Taylor's hypothesis, compared with the remaining SGS models. In fact, this observation could rather elaborate on the relatively longer vapor length in the Smagorinsky case, c.f. Figure 11. Particularly, the dissipation in the corresponding larger length scales would imply less flow turbulence and subsequently higher E k concentration in the mean flow, hence the impact on fuel vapor mixing and penetration.

Conclusions
The Lagrangian particle tracking (LPT) approach coupled with large-eddy simulation (LES) was utilized to investigate the sensitivity of various modeling assumptions on the global characteristics of an evaporative fuel spray under the ECN Spray A baseline conditions. A volume filtering approach was employed to model the near-nozzle region using an injection cylinder of volume V cyl and dimensions larger than those of the nozzle. The droplet secondary breakup was modeled using the KHRT model with two different model parameter sets, and the no-breakup model with constant droplet sizes. The subgrid scales (SGS) were modeled using the implicit LES (ILES), Smagorinsky, and one-equation model (K sgs -model). The sensitivity towards (i) varying the dimensions of V cyl while using ILES and KHRT, (ii) droplet no-breakup model, in comparison with KHRT model for droplet secondary breakup, and (iii) the SGS effects from different models when utilizing a droplet no-breakup approach, were all comprehensively investigated and compared to the ECN Spray A data. A summary of the main findings is highlighted as follows: 1.
Spray metrics are sensitive to the injection cylinder dimensions. The cross-stream dimension of V cyl has a stronger influence on the spray characteristics, including spray penetration lengths, momentum intensity, and droplets evaporation.

2.
The no-breakup model approach performs equally well compared with the tuned KHRT model.

3.
With respect to global spray metrics, the no-breakup model approach is shown to be insensitive to the chosen SGS closure model.

4.
According to statistical analyses on mixing indicators, explicit SGS models predict a wider range of scalar dissipation rates in coarse resolution regions compared with the ILES.

5.
According to the PSD analysis, dampened energy spectra in the lower frequency range, which correspond to larger length scales by Taylor's hypothesis, are produced by the Smagorinsky SGS model compared with the K sgs -model and the ILES. Acknowledgments: The authors acknowledge CSC (Finnish IT Center for Science) for providing the computational resources.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The  Figure A1 depicts the temporal evolution of liquid and vapor penetration lengths of the spray when the injection half cone angle (α) is varied between 2-10 • . Spray sub-models for SGS and droplet breakup follow those in cases B3 or C1 from Table 3. The plots show that liquid length is insensitive to α within the considered interval, whereas a mild variation is noted for vapor length, which is expected to be rather attributed to LES realization effects. Therefore, the present analysis justifies the choice of α = 5 • employed in this manuscript, since minimal effects on global spray metrics are noted for small variations in α. Results are compared with ECN data, wherein the filled region denotes experimental uncertainty. Spray sub-models follow those in cases B3/C1.