Large Eddy Simulation of cavitation effects on reacting spray flames using FGM and a new dispersion model with multiple realizations

In the present study Large Eddy Simulations (LES) are performed to investigate the ignition behavior of cavitating and non-cavitating n -dodecane flames. Inspired by the work of Tsang et al. (2019) [1] an LES specific dynamic dispersion model, that evaluates the dispersion velocity locally, is proposed. Although the model discards tuning of global model constants, excellent mixing predictions are obtained for all cases. The resulting model is extensively validated using inert Spray A conditions as defined by the Engine Combustion Network (ECN). Subsequently, it is applied to the larger orifices of Spray C and Spray D. A Flamelet Generated Manifolds (FGM) approach that takes the effect of scalar dissipation into account is adopted for combustion modeling. The coupling of the turbulence modeling approach and FGM shows excellent predictions of ignition characteristics on Spray C and Spray D, suggesting a minor effect of cavitation on ignition development. For the sake of understanding the injection-to-injection variations of LES, multiple realizations are performed. Based on the analysis of structure similarity index (SSI), it is found that a single realization is sufficient for global parameters such as ignition delay time (IDT) and lift-off length (LOL). However, different number of realizations are needed for different scalar fields. It is suggested that the temperature-based IDT is preferred for a single realization while a radially integrated intensity is needed for an OH-based IDT or LOL. © 2021 The Author(s). Published by Elsevier Inc. on behalf of The Combustion Institute. This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ )


Introduction
The complex turbulent spray processes taking place in an internal combustion engine are driven by the mixing of fuel and air, leading to the chemical reactions that release a fuels energy. With the demand for clean and efficient internal combustion engines, studying these topics still needs detailed and considerable effort s to achieve a deep knowledge of the underlying physics.
The spray development, i.e., the mixing of fuel vapor and air, directly influences the flame stabilization, i.e., the lift-off length (LOL). A longer ignition delay time (IDT) prolongs, and thus enhances the mixing process, reducing the formation of harmful emissions [2] . The impact of spray development on soot emission has been quantitatively investigated in several studies [3,4] . It is pointed out that the equivalence ratio distribution at the lift-off length is a key indicator for soot formation. More importantly, soot formation was shown to be avertable when the centerline equivalence is kept below a value of two, regardless of the experimental setup [5][6][7] . This phenomenon can be attributed to the variations tions. Although dispersion (i.e., the turbulence fluctuation induced term in drag force here) has been successfully modeled and widely applied in Reynolds-averaged Navier-Stokes (RANS) implementations, it is only recently adapted to Large-Eddy Simulation (LES) using adjustable global model constants [1] . In reference [1] , the importance of adapting the RANS-type dispersion modeling for LES has been addressed, while the effect of dispersion modeling on spray development was extensively studied. In this paper, attention is thus paid to the dispersion modeling specifically for LES.
The Engine Combustion Network (ECN) [11] aims to collaboratively address and reduce uncertainties of experimental setups, providing a solid and complete characterization of the combustion process such that improvements on predictive computational fluid dynamics (CFD) models can be achieved. ECN Spray A is a well-defined single-orifice spray into a meticulously designed ambient regarding temperature, pressure, and oxygen concentration closely resembling those in diesel engines. An extensive database has been built using advanced techniques as validation with respect to both inert (liquid and vapor penetration, and spatial fuel distribution [4] ) and reacting (IDT and LOL [12] , major species distribution [13,14] and soot production [15] ) characteristics. More recently, Spray C and Spray D were introduced as larger orifice injectors to study the effects of cavitation, and particularly its influence on soot emission.
Cavitation initiates where the local liquid pressure inside an large-nozzle-size injector is lower than its vapor pressure, leading to a phase change and thus a "choked" flow inside the nozzle. Injected mass is directly affected by the cavitation, generally followed by an increased spreading angle and enhanced fuel atomization, i.e. the fuel-air mixing [16,17] . It is thus expected that cavitation can be one of the driving forces determining both the ignition behavior and emissions. Indeed, this is confirmed and studied extensively in both experimental investigations and numerical analyses [4,[18][19][20][21] . Cavitation is directly a consequence of injector geometry, flow separation and depressed pressure gradient are enhanced especially when passing a sharp edge [21,22] . Spray C is thus designed as a straight-hole nozzle with sharp edges to promote cavitation, while Spray D adopts converging nozzle as well as a hydro-eroded inlet to avoid cavitation. The different configurations lead to significantly different flow development inside the nozzle, and consequently have an impact on the morphological spray characteristics. Subsequently, in this investigation, both the inert and reacting cases of Spray C and Spray D are simulated and compared by directly implementing spray parameters obtained in experiments, such as the discharge coefficient, spreading angle, etc.
The flamelet concept adopted in this work is a combustion modeling approach in the turbulent flame locally described as laminar flames (flamelets) [23] . This approach is widely applied in turbulent spray combustion simulation due to its significant advantage in reducing the number of transport equations that need to be solved [24][25][26][27][28][29] . The Flamelet Generated Manifolds (FGM) [30] applied in this article can be regarded as a typical tabulated flamelet approach that retrieves information from a prepared database. The concept uses several independent controlling variables and has been successfully applied to several LES studies of ECN cases [31,32] . Turbulence-chemistry interaction (TCI) accounts for the influence of subgrid scale fluctuations on chemical reactions. Dedicated investigations suggest that neglecting TCI leads to the failure of capturing the important physics during combustion processes, particularly when higher temperatures or fuel-rich conditions are encountered [33][34][35][36] . In this sense, here a presumed probability density function (PDF) using the top-hat filter is adopted to describe the local flame structures.
An LES is a statistical model that provides a more precise energy spectrum compared to RANS by filtering only the small scales. However, more uncertainties arise as the filtering operation highly depends on the mesh resolution. In spray simulations, the uncertainties increase since both breakup and turbulence-spray interaction are influenced by a random number that depends on its seed. The initialization of droplets during breakup is realized by allocating velocities whose direction and magnitude are randomly chosen. A statistical average is thus preferred to analyze the global results. However, due to the high computational cost of LES, it is common practice to analyze a single realization. Only a few of investigations performed and addressed the difference of ensembleaveraged LES. This motivates to perform multiple simulations to understand the potential limitations and uncertainties of a single LES realization, especially when using FGM together with the turbulence modeling approach presented in this investigation.
Following the work by Tsang et al. [1] , the main contribution of the present paper is an improved dispersion model such that dispersion velocity can be automatically determined locally. Five main objectives of this investigation are addressed. (1) To validate the new dispersion model. (2) To evaluate the use of the turbulence modeling approach for non-cavitating and cavitating sprays with a different nozzle diameter. (3) Analysis of the required number of realizations for specific combustion parameters. (4) A comparison of the characteristics of different phases during ignition between Spray C and Spray D. (5) To bring insights to the soot formation for reacting Spray C and Spray D by providing the ensemble-averaged spatial vapor fuel distribution of the "quasi-steady" state.
The outline of this paper is as follows. The modeling approach and general equations for both liquid and gaseous flow are described in Section 2 , where the improved dispersion model is introduced. In Section 3 a case description and computational configuration is given. The non-reacting spray simulations are presented in Section 4 . The modeling approach is first extensively validated via ECN Spray A, and further applied to Spray C and Spray D. Section 5 introduces the coupled FGM with top-hat filtering, and presents the prediction of reacting conditions for Spray C and Spray D. In this section, a detailed analysis of ignition behavior is carried on. Subsequently a single LES realization is compared to the ensemble-average to eliminate the injection-to-injection variation of LES. The findings are finally concluded and summarized in Section 6

Basic equations and numerical methods
The present work applies a Lagrangian-Eulerian approach to simulate the spray process. A customized solver in OpenFOAM 2.4.x [37] is created. The finite-volume method is adopted for the gaseous phase simulation using a PISO (Pressure-Implicit with Splitting of Operators) corrector to satisfy mass conservation. An implicit second-order backward differencing time scheme is used for temporal integration, where the time-step is locally adjustable during the simulation according to the Courant-Friedrichs-Lewy (CFL) number (set as 0.3 in this study, leading to in-situ simulation time steps of approximately 6 ×10 −3 and 7 ×10 −3 second for Spray C and Spray D, respectively). All the diffusive terms are discretized by a second-order linear central differencing scheme, while the convective terms are treated by a blend of the first-order upwind scheme and a central differencing scheme.

Governing equations
In the present study, LES equations govern the gas-phase physics, including mass and momentum conservation with unclosed terms that need to be modeled: The source terms, ¯˙ S ρ and ¯˙ S u,i , arise from either evaporation or displacement between phases and are listed in Table 1 . The Favre filtered values are obtained by applying the density-weighted spatial-filtering operation: φ = ρφ/ ρ. Any variable is thus decomposed into a resolved part, φ, and a subgrid one, φ , that is unresolved. ρ represents the fluid density, u i the velocity, p the pressure, and τ i j the viscous stress tensor. i j is the turbulent stress tensor that needs to be modeled.
It is worth mentioning that the FGM treats the species differently compared to detailed chemistry. In an FGM, species information is directly retrieved from the database rather than applying a transportation equation. This in conjunction with the treatment of thermophysical properties that is described in detail later.

Subgrid modeling
i j in Eq. (2) is treated as a non-viscosity shear stress that is closed via a dynamic structure model (DS) [38] : ν nozzle is introduced as an artificial kinematic viscosity that is induced by the high-speed spray jet to allow for a relatively coarser mesh [39] . For the k sgs , the subgrid kinetic energy, an additional transportation is needed which is given by: where ν sgs is the subgrid scale kinematic viscosity, P and D represent the production and dissipation effect respectively. The detailed description of these two terms are described in reference [38] and [1] . Following [1] , in this study the model parameters for P and D are set to 0.1 and 0.5, respectively. The source term ¯˙ S k,sgs is found to be prominent in improving the prediction of penetration. For more details the reader may refer to Tsang et al. [1] .

Liquid phase modeling
The liquid phase is described using standard Discrete Droplets Model (DDM) [40] . Lagrangian parcels are tracked instead of individual droplets (for convenience in the following part "droplet" refers to a parcel unless a specific explanation is added), and are governed by a momentum and a position equation: u d,i denotes the i th velocity component of a droplet. In OpenFOAM, the "face-to-face" method [41] is applied for particle tracking, the cell index to which a particle belongs is stored as particle information and continuously updated. During simulation, the program iterates particles rather than cells. The drag force on a droplet, F i,d , appears in the modeling of both the liquid ( Eq. (5) ) and the gaseous phase (the source terms of Eqs. (2) and (4) ). The so-called "two-way coupling" between the gaseous and the liquid phase is thus realized. Drag force can be considered as a friction that is induced by the velocity difference between phases. A spherical drag force model [42] is applied in this study: u i represents the dispersion velocity. ρ c is the carrier-phase density while r d is the droplet radius, and C D is the drag coefficient. The carrier-phase velocity u c,i is an interpolated velocity at the droplet position considering the velocities at its surrounding cells. In this study, a first-order interpolation is adopted. The mass and heat transfer of Lagrangian particles is modeled by the Ranz and Marshall correlation [43] . The injected fuel is distributed near the nozzle outlet according to the Rosin-Rammler distribution. The droplet breakup is modeled by the Kelvin-Helmholtz and Rayleigh-Taylor (KH-RT) instabilities including a restriction of RT instabilities in the near-nozzle zone [44] . Table 1 lists all the sources and sinks that appear in the gaseous phase equations that are induced by the Lagrangian parcels. The effect of parcels on the fluid are considered by summing over the number of parcels, np, in each cell. Each parcel contains n d droplets sharing same thermophysical properties. The volume of a computational finite-volume cell is marked by V f .

The dispersion model
It is evident that turbulence has a direct influence on momentum, mass, and heat transfer between phases in the scope of dense sprays [45] . Dispersion modeling is to include the effect of turbulent fluctuations, i.e. the subgrid motion in LES, on spray development. The Approximate Deconvolution Model (ADM) approximates the subgrid velocity via the truncated Van Cittert series expression [46] , which reads: Ideally, ADM reconstructs 95% of the sub-grid fluctuations [47] . However, the information of those that exceed the cut-off length (which is highly dependent on the mesh quality) cannot be rebuilt [48] , and the reconstructed ratio can be much lower in many applications. Hence, there is a residual part to fully cover the sub-grid motion. In analogy with the classic stochastic dispersion model adopted in RANS, the residual part can be created randomly by a Gaussian distribution as it is a reasonable assumption to consider the residual turbulent motion as locally homogeneous and isotropic [49] . The subgrid dispersion velocity is thus decomposed into a deterministic term and a residual term: where the residual part is Here, σ is the variance of Gaussian distribution.
The subgrid kinetic energy k sgs can be explicitly expressed by 1 / 2 | u 2 i | . By applying Eq. (9) , the subgrid kinetic energy then equals: | u sgs,ADM + u sgs,RES | 2 = k sgs,ADM + k sgs,RES + u sgs,ADM · u sgs,RES , k sgs,ADM is the only explicit term on the right-hand side of Eq. (11) . The effect of residual induced terms can be either an augmentation or an abatement of k sgs,ADM . Tsang et al. [1] proposed a straightforward solution to determine the u sgs,RES by applying two case-dependent model constants.
In the framework of RANS, the variance of a Gaussian distribution is controlled by 3( 1 2 σ 2 ) = k . In this way, the dispersion is assumed to be isotropic and the variance of each component of the dispersion velocity accounts for 1/3 of the kinetic energy. Both the magnitude and the direction of a velocity component is created randomly. Tsang et al. [1] interpreted the role of the variance from controlling the velocity component on the kinetic energy level ( k ) to the level of subgrid kinetic energy residuals ( k sgs,RES ). As k sgs,RES is an unknown term, the first parameter is introduced to identify u 2 sgs,RES , by assuming its variance follows 3( 1 2 σ 2 ) = C sig k sgs . As the reconstructed deterministic term k sgs,ADM can be either greater or less than the sub-grid kinetic energy k sgs , the second parameter . C sgs and C sig are defined globally such that the sum of the two velocities, u sgs,ADM and u sgs,RES , hold the magnitude of k sgs on a correct level. Although the model constants are globally defined, an extensive study on different combinations of the two factors by Tsang et al. [1] justifies the importance of correctly modeling the sub-grid dispersion velocity.
In this paper, a new model is proposed to dynamically define u sgs,RES following the concept of the typical stochastic model within RANS. As u sgs,RES can either lead k sgs,RES being an augmentation or an abatement of k sgs,ADM , its variance can be determined by considering the two extreme conditions. In the case of an augmentation, the maximal condition is that the direction of the residual sub-grid velocity is exactly the same as the reconstructed u sgs,ADM given by Eq. (9) , saying that u sgs,RES = γ · u sgs,ADM . As the sum of the two velocities ( u sgs,ADM + γ · u sgs,ADM ) explicitly determines the subgrid kinetic energy k sgs , the constant γ can be calculated and thus the k sgs,RES . The Gaussian distribution randomly creates the velocity components via Eq. (10) by keeping its variance to obey 3( 1 2 σ 2 ) = k sgs,RES . After derivations, by applying σ 2 = 2 3 k 2 sgs,RES = 1 3 | γ · u sgs,ADM | 2 , yields: The new dispersion model is validated on the inert cases of Spray A, Spray C, and Spray D, and the results are presented in Section 4 . Subsequently, the model is applied in reacting Spray C and Spray D cases in Section 5 .

Case description
The benchmark ECN case, Spray A, is used in this study to validate the turbulence modeling approach. The approach is then applied to the simulation of both inert and reacting cases of Spray C and Spray D, respectively. The spray development until 3 ms is simulated. The injection starts from the beginning and lasts until the end of the simulation. All cases feature the single-orifice injection of n -dodecane at 150 MPa and share the same ambient conditions. The injectors are designed such that they can be used to study the cavitation effect. The details are summarized in Table 2 .

Computational setup
In the present work, a cylindrical domain with a diameter of 47 mm is used to simulate the part of the combustion chamber in which the flame is present. The injector is located on the base. The computational domain is discretized into a structured mesh with hexahedral elements. The mesh is refined in a quadruple prism zone near the nozzle outlet to improve the accuracy of prediction in this region. Inside the prism, the cells are equally distributed along the radial direction and only stretched alongside the spray downstream. The rest of the mesh applies expansion ratios of 1.015 and 1.01 for radial and axial directions, respectively. Two meshes are compared. One applies the commonly adopted 62.5m μm as the minimum grid spacing following the recommendation in literature [50] , while the other one doubles the smallest cell size to 125 μm. This leads to a total number of 3.6 and 2 million cells, respectively. This study is carried out using the Yao chemistry mechanism [51] with 54 species and 269 reactions.

Simulation of non-reacting spray
Following the described modeling approach, this section presents the prediction of the inert Spray A, Spray C, and Spray D cases.

Governing equations for the non-reacting spray
Transport equations of species as well as the sensible enthalpy transportation are added, which are used in conjunction with species information to calculate the temperature. The added equations read: Y i indicates the mass fraction of species, ˜ h s is the sensible enthalpy. λ and C p represent the heat capacity and thermal conductivity of gas mixture, respectively. Sc, P r denote Schmidt and Prandtl number, respectively. The subscript sgs refers to the subgrid scale properties. In this study, Schmidt and Prandtl number for both the resolved and unresolved parts are set to 1. ¯˙ S Y i and ¯˙ S h s are sources due to phase change. In both equations, no reaction induced source exists as this approach is only used to study the inert cases. For the reacting spray we apply the FGM approach, which does not require transportation of species and enthalpy.

Spray A validation
In this section, the standard ECN Spray A case is validated as it has the most extensive database including liquid penetration, vapor penetration, and spatial fuel distribution. Two dispersion models are tested. The first one follows Tsang et al.'s approach [1] with elaborately tuned model parameters as is reported in Bao et al. [74] . For brevity this dispersion modeling approach is indicated by u global . The other one adopts the approach described in Section 2.3 and is marked by u local . They are tested on both the refined (62.5 μm) mesh and the coarser mesh (125 μm). It is found that the new dispersion model, by locally defining the model parameters, allows a coarser mesh. For brevity, in the following parts, only results using the refined mesh for u global and for the approach with u local using the coarser mesh are shown.
In Fig. 1 the performance on the prediction of penetrations is shown. The vapor and liquid penetration are defined as the farthest axial position from the injector outlet where gaseous fuel mass fraction has decreased to 0.001, and where 90% of the injected liquid-phase fuel are contained, respectively. The u local model gives even a better prediction on a coarser mesh than the u global model on the finer mesh. Moreover, no tuning of C sig and C sgs is needed for u local while multiple simulations were performed for different cases for the global one to identify the best-tuned model parameters. A more detailed quantitative comparison can be found in Fig. 2 . In Fig. 2 a spatial and time-averaged quantities of the steady-state mixture fraction along the spray axis are shown. Figure 2 b shows the radial distribution of the time-averaged mixture fraction and temperature fields. Both are averaged from 1.2 ms until 3 ms. The value at each radial location is achieved by computing an azimuthal average of the time-averaged field around the spray axis. Again, it is seen that the new dispersion model allows for a coarser mesh. According to the standard deviation of mixture fraction, u local gives even better performance. In most of the LES studies, under-prediction of the downstream mixture fraction profile is a common difficulty [54,55] . The new dispersion model gives good agreement when compared to the experiments for all available experimental profiles. The radial distribution of mixture fraction even suggests that u local gives a better agreement in the upstream region (0.018 and 0.035 m) while the downstream prediction of both u local and u global are in good agreement with the experiment. This illustrates that the tuned global model parameters are not the best for all regions. Figure 3 shows the instantaneous temperature distribution in mixture fraction space at 3 ms after the start of injection for the u local model by considering all the cells in the computational domain. The symbols are colored by the local liquid volume fraction. The adiabatic mixing line is shown by the dashed line. The adiabatic mixing line in this case is obtained by setting the fuel temperature to 155 K to implicitly include the effect of latent heat. This procedure is also described in Akkurt [56] . The solution of the adiabatic mixing line is obtained from the simulation of a one-dimensional counterflow diffusion flame configuration using CHEM1D with the chemistry artificially turned off. It is seen that temperature mostly distributes around the adiabatic mixing line, while the liquid parcels occupy only a limited volume of a cell (seldom more than 0.1 of cell volume, indicated by bright blue). This can be partially attributed to the use of a coarser mesh. Hence, it is expected that especially near the nozzle outlet the extreme changes in physical properties can be avoided since a coarser mesh can be applied. More features of the explicit LES-type dispersion model can be found in Appendix A .

Performance on Spray C and Spray D
In this section, the modeling approach of Spray A is further applied to Spray C and Spray D. The inlet configurations of Spray C and Spray D are chosen following Westlye et al. [19] . Following the suggestions given, the spreading angles are set as 21.5 and 19 de-   grees for Spray C and Spray D, respectively. The liquid and vapor penetration can be found in Fig. 4 . Both vapor and liquid penetration are well captured in all cases. Following ECN guidelines, a more reliable definition for the liquid penetration that is consistent with the experimental procedure is also adopted as a reference. The liquid penetration can be evaluated via the projected liquid volume (PLV), which is defined as: where LV F refers to the liquid volume fraction. The distance from the nozzle outlet where a threshold P LV F = 0 . 2 × 10 −3 mm 3 liquid/mm 2 is found, is defined as the liquid penetration [57] . Figure 5 shows the liquid penetration defined by the time-averaged projected liquid volume fraction (from 1.2 ms to 3 ms). The predicted values are 21.1 and 24.3 mm for Spray C and Spray D, respectively. This is very close to the experimentally calculated val-ues of 20.7 and 24.2 mm. Furthermore, this illustrates that a commonly used definition based on the furthest axial distance in the computational domain that contains 90% of the liquid provides similar values to the projected liquid volume fraction method recommended by the ECN. It must be noted that for the u global model, not only the global model constants in the dispersion model but also the model parameter of breakup model (KH-RT in this case) were tuned depending on each case as well. For the u local model, a single set of uniform breakup model parameters were adopted for all three cases. Besides, the dispersion model discards the fitting of model parameters ( C sgs and C sig ).

Simulation of reacting spray
This section presents the reacting Spray C and Spray D results based on the validations using non-reacting spray simulations in Section 4.3 . An extended FGM that takes the effect of scalar dissipation into account is adopted to model combustion [74] .

Governing equations for the reacting spray
FGM is a tabulated chemistry reduction approach that retrieves the thermophysical properties from the prepared flamelet database. Each flamelet contains information of a representative laminar flame (counterflow diffusion flame in this case) that is solved in physical space using the software CHEM1D [58] . The details can be found in Bao et al. [74] . The FGM considers the combustion as a function that is parameterized by controlling variables. Generally, the mixture fraction Z and a reaction progress variable Y c are chosen. Z is defined following Bilger's approximation [59] that indicates the fuel-air ratio locally, while Y c is typically a linear combination of several representative species whose value has a monotonically increasing (or decreasing) dependency on the reaction progress. The FGM applied in this study includes the effect of the scalar dissipation rate χ. χ is directly related to the strain rate a in the flamelet configuration and is used as an extra parameter to represent the state of individual flamelets as described in Bao et al. [74] . The FGM applied in this article automatically selects the strain rate during the simulation according to the local flow characteristics. For more details please refer to [74] .
In turbulent combustion modeling, the pre-computed FGM library based on laminar flames needs to be extended to include the TCI. This can be realized through a presumed PDF as described in Section 5.2 . The present work applies a top-hat PDF [60] to the mixture fraction and a δ-PDF on the progress variable and scalar dissipation. For each controlling variables a conservation equation needs to be solved for LES: and an equation for the variance of ˜ Z , ˜ Z v to determine the shape of the PDF: In this work, Y c is the unscaled progress variable which is defined as: The chosen definition of Y c was validated in our previous study and was targeted at not only capturing the ignition delay, but also C 2 H 2 development. The details are described in Bao et al. [74] . In Eq. (18) , the filtered scalar dissipation rate, χ , indicates to which extent the variance of mixture fluctuations decay by mixing. It can be decomposed into two parts: The resolved part is directly calculated through χ res = ν/Sc|∇ Z | 2 , while the subgrid part is modeled following the local equilibrium hypothesis [61] : χ sgs = ν sgs /Sc sgs / 2 ˜ Z v . Note that in the case of reacting spray simulation in the current study, temperature is directly retrieved from the database, so that enthalpy transportation is not needed.

Turbulence-chemistry interaction
In this investigation, the FGM is extended to take the effect of turbulence fluctuations on local flame structures into account using a top-hat filtering PDF [60] , considering its great advantage in memory reduction. A filtered scalar can be expressed as a function of independent controlling variables: where C χ is the normalized scalar dissipation, and P is the PDF. In practice, it is common to assume that the controlling variables are statistically independent of each other. By applying a δ-PDF for C Y c and C χ , Eq. (21) reduces to: A top-hat filtering approach assumes the distribution probability of scalars in a LES cell to be equal within the bounds of the independent variables (i.e. the mixture fraction, Z, and in this case, the lower and upper bounds are respectively Z l and Z u ). Integration in Eq. (22) can thus be expressed as: where Z is the width ( Z u − Z l ) of the Z interval that is directly determined by its variance: C is a model constant. The lower and upper bounds by definition are thus Z l = Z − Z/ 2 and Z u = Z + Z/ 2 , respectively. However, depending on the variance, Z l and/or Z u can exceed the allowable mixture fraction range. Therefore, predicted filtered scalars cannot be calculated by a top-hat function alone. Three categories are identified: 1. Both Z l and Z d are within the mixture fraction limits.
Z l ≥ 0 Z u ≤ 1 ( Fig. 6 a). The top-hat filter is applied to the whole bounding range, leading to: 2. Both Z l and Z u are outside of the mixture fraction limits. Z l < 0 Z u > 1 ( Fig. 6 b). The top-hat filter is applied only within the mixture fraction limits, while a Dirac peak is adopted at Z = 0 and Z = 1 : 3. Either Z l or Z u is out of the mixture fraction limits. a Z l < 0 Z u ≤ 1 ( Fig. 6  b Z l ≥ 0 Z u > 1 ( Fig. 6 The top-hat filtered FGM in this study only stores the mean of a thermochemical property as function of mixture fraction and progress variable and its integrated value over that range in mixture fraction. Thus, the top-hat filter only increases the database size by a factor of 2. A beta-pdf formulation with e.g. 10 levels for the variance of mixture fraction involves 10 fold increase (one table for each variance level). The performance of the top-hat filters is evaluated via the unsteady progress of a laminar flame at a = 500 s −1 . In Fig. 7 the filtered source of the progress variable is presented as it is the most intensively changing variable among all those included in the database. A constant C = 1 / 6 for the top-hat PDF is adopted in this case. It is seen that the top-hat gives very similar results to a β filter at both low and high C Y c (close to the chemical equilibrium), especially at a lower variance. Meanwhile, a higher variance leads to a smoother scalar distribution on of the progress variable source in mixture fraction space for both PDFs, as expected.

Multiple realizations
Regarding the ECN cases, different numbers of realizations are recommended in various investigations. Hu et al. [62] suggested 14 realizations to achieve a converged flow, while Senecal et al. [54] declared that 9 realizations are sufficient to represent the flow characteristics. Note that for an LES, this can be dependent on mesh resolution. Pei et al. [63] concluded that by refining the mesh, the minimum number of realizations can be reduced to 5 and 6 for mixture fraction and temperature, respectively, while the  number is increased to 14 for a scalar with more fluctuations such as OH. It is thus expected that the number of realizations required must be evaluated based on different quantities for each case.
In this work, multiple realizations are taken for both Spray C and Spray D by perturbing the random seed. The convergence of representative scalars is quantitatively analyzed with the help of the so-called relevance index [64] . 8 realizations are carried out for both the Spray C and the Spray D flame. Following Liu et al. [64] , the relevance index, i.e. the structure similarity index (SSI), is defined by projecting one field onto another: where u and v are two vectors of an evaluated quantity such as temperature, representing the ensemble-average of a certain number of realizations and the ensemble-average of all realizations performed as a target (8 in this study), respectively. (u, v ) represents the inner product over the whole domain and || u || = (u, u ) 1 / 2 . It is typical to evaluate the relevance index at a random time instance since variation in time does not change the trend of the relevance index [63] . In Fig. 8 the relevance index of representative scalars is plotted for Spray C and Spray D, respectively, at 1.5 ms and 2 ms after the start of injection. It is seen that relevance index increases with the number of realization as expected. Besides, the tendency of the relevance index performs differently for different scalars, while it is hardly influenced by the chosen time instance. It is noted that for a certain scalar of interest, Spray C and Spray D have similar relevance index tendencies. In both Spray C and Spray D, the relevance index of temperature barely changes with ensemble size, even a single realization leads to a relevance index of 0.99 which can be considered as a "converged" result. Contrarily, 8 realizations seem insufficient to provide a reliable ensemble-averaged OH field, as the gra-dient of the relevance index is still non-zero and the value of SSI is much lower than 0.99 even when the ensemble size approaches the total number of simulations. The relevance index of mixture fraction and formaldehyde, however, increases sharply within the first 3 and 4 realizations and subsequently levels off. To be precise, 4 realizations are needed for mixture fraction while an average of over 6 realizations is sufficient for formaldehyde. The result suggests that temperature should be preferred when analyzing ignition behavior in consideration of computational cost.

Ignition delay and flame lift-off
In this section, the ignition delay time (IDT) and flame lift-off length (LOL) for Spray C and Spray D are quantitatively compared to the measurements. 6 LES realizations for both Spray C and Spray D are performed to investigate the sensitivity of IDT to shot-to-shot variation.
In the present work, the IDT is based on temperature. It is defined as the time between the start of injection (SOI) and the time the greatest rate of maximum temperature rise in the domain ( dT max (t) dt | max ) is observed [65,66] . The definition used here depends on temperature rather than the mass fraction of OH. This choice is a direct consequence of the observation presented in Section 5.3.1 , and leads to IDT of 0.4425 and 0.4457 ms respectively for Spray C and Spray D (both are the average of 6 realizations). For both Spray C and Spray D, a very limited variation of IDT across realizations is observed (the variability reaches approximately 2%). The standard deviations of IDT for Spray C and Spray D are 5 . 52 × 10 −3 and 7 . 23 × 10 −3 ms, respectively. This illustrates that a single realization is qualified to be representative. As a global parameter, the variation of IDT is in accordance with the observation by Xue et al. for non-reacting spray, that global parameters can be reflected by a single LES realization [67] . Furthermore, the results show that IDT for Spray C and Spray D are very similar. This corresponds to what has been observed experimentally [68] . According to the ECN, the experimental IDT for Spray C and Spray D are 0.561 ms and 0.563 ms, respectively [11] . The details regarding measurements are discussed by Maes et al. [68] , where the differences between the chemiluminescence-and pressure-derived IDT are addressed. Note that depending on the experimental facility, IDT values that are 10-20% lower are also observed.
The LOL indicates the stabilization location of a flame, and is generally analyzed using OH (e.g. the closest axial distance from nozzle orifice where 14% of the maximum OH mass fraction during simulation is found). In literature, different value of OH thresholds (from 2% to 50%) for the LOL can be found [28,63,69] . Pei et al. concluded that applying various thresholds for the LOL only leads to minor changes [63] . In this study, the LOL is evaluated as function of time following our previous work [74] . For each time instance, the averaged OH mass fraction of transverse planes (with regards to the spray axis, noted by OH a v e,sur f i ) are considered. By evaluating the average over a plane, the high variations of local OH mass fractions generated by the LES can be avoided. In this investigation, we take the closest axial distance from the nozzle outlet where 2% of the maximum OH a v e,sur f i is observed as the LOL. From the instantaneous LOL, a time-average value is calculated, which can be compared to the experimental results. While it was shown above that the OH field is more susceptible to variations, 6 realizations for each case are considered to discard the shot-to-shot variations when comparing LOL values. The results are shown in Fig. 9 .
The instantaneous LOL varies significantly during the quasisteady stage, and the variations are quite different between different realizations. However, a time-averaged LOL from a single realization shows negligible difference compared to the time-averaged LOL obtained from the ensemble-average LES. The fluctuation of the instantaneous LOL is significantly reduced by performing more realizations. For Spray C the standard deviations of a single LES realization and the ensemble-average of 6 realizations are 1 . 2 × 10 −3 and 5 . 5 × 10 −4 m, respectively, while for Spray D these values are 8 . 1 × 10 −4 and 3 . 3 × 10 −4 m. Given these standard deviations, the predicted difference between LOL (marked by * ) of Spray C and Spray D is significant and compare well to the experimental values and trend.

Temporal development of flame
In this section, the temporal behavior of Spray C and Spray D during ignition are analyzed. For both cases, a single LES realization is compared to the ensemble-average of 8 realizations. Representative time instances, covering the cool flame stage until the stage in which the lifted flame can be considered as quasi-steady are sampled and discussed. For brevity, the IDT is denoted by τ * .
Five representative time instances for Spray C and Spray D are selected: (1) 0.1 τ * that is prior to the first stage ignition, (2) 0.5 τ * characterizes the low temperature combustion stage, (3) 0.9 τ * that is approaching to the IDT when both low and high temperature zones exist, (4) 1.01 τ * , slightly after the IDT when ignition kernels can be observed (this is the nearest output time instance in simulation) (5) and finally 3 τ * as a representative of the quasi-steady flame. In Figs. 10 and 11 , the spatial ignition characteristics are analyzed. The scatter plots in Fig. 12 show the distribution of temperature in mixture fraction space, where the distance from the injector orifice is denoted by color from blue to red. To enhance visibility, only mixture fraction values greater than 0.001 are shown in the scatter plots to discard the regions that are not influenced by the spray. Figures 10 and 11 show that the ensemble-average reduces the fluctuations as expected, resulting in much smoother local flame structures. It is clear that the flame depicted by an ensembleaverage is characterized by an uninterrupted stoichiometric mixture fraction contour (black solid line) as found in RANS studies. Nevertheless, the variations originating from the LES still exist. Although in Section 5.3.1 it is proven that 8 realizations are more than enough to achieve a "converged" temperature prediction according to the relevance index, local fluctuations can still be observed. At 0.1 τ * , the spray is prior to the first-stage ignition. Temperature change in this period depends on mixing rather than combustion. In simulations of diesel-like ECN sprays, it is common practice to differentiate between the first-stage (cool flame) and second-stage ignition (high temperature combustion) as the time when a temperature rise of 400 K is achieved [36,63,65] . 0.5 τ * and 0.9 τ * are thus respectively attributed to first-and secondstage ignition according to their temperature rise. The flame develops quickly during the second-stage ignition. A significant change in spatial temperature distribution is observed in a relatively short period from 0.9 τ * to 1.01 τ * , while the area indicated by the stoichiometric mixture fraction only shows slight variation. This illustrates that in contrary to the cool flame stage, chemical reactions dominate the temperature evolution rather than mixing. Following Kahila et al. [66] , the ignition kernels are defined as regions with a temperature higher than 1550 K inside. According to the contour plots, multiple ignition kernels appear at 1.01 τ * . It is noted that in both Spray C and Spray D, the ignition kernels are found at the periphery of the jet, upstream from the head-vortex of the spray like observed by Desantes et al. [35] , Ong et al. [65] . This is different from what is found for the standard Spray A flames where ignition occurs at the head of the jet [69] . The ensemble-averaged contour plots show that ignition kernels appear at around 30 mm and 33 mm away from the nozzle outlet for Spray C and Spray D, respectively. This is roughly 5 mm farther downstream than their corresponding LOLs. After IDT, the ignition kernel distribution is wider compared to 1.01 τ * . Although the flame continues to develop farther downstream, the behavior within the ignition kernel region remains stable in terms of temperature so that the flame can be considered as "quasi-steady". As can be seen in Figs. 10 and 11 , stoichiometric mixture fraction becomes quasi-steady very quickly. In each case, the stoichiometric mixture fraction contour appears at the same distance from the injector tip at selected time instances. The location seems to be similar to quasi-steady LOL or CH 2 O inception point. This is similar to observations of the LOL and CH 2 O inception in experiment [70] . An interesting phenomenon is that in a single LES realization, the highest temperature is located at a very thin layer along the stoichiometric mixture fraction. However, for the ensemble-averaged LES, temperature is more widely distributed. This is due to the averaging of the different spatial distributions induced by different realizations. For the same reason, the peak temperature decreases by a value of the order of 150 K in an ensemble-averaged LES.
To substantiate previous findings, and to gain additional insight into the ignition and flame stabilization processes of Spray C and Spray D, Fig. 12 shows the early temperature evolution in mixture fraction space. At 0.1 τ * , the temperature distribution follows the mixing line for both cases. However, a significant difference between Spray C and Spray D is seen. In both the single LES realization and the ensemble-averaged LES of Spray D a leaner mixture is established. The difference in temperature distribution seen in contour plots is definitely due to the spray development. Increasing temperatures in the following two panels (before τ * ) are first predicted at fuel-rich mixtures. After τ * , the highest temperature shifts back towards oxidizer side until it approaches the stoichio-   metric mixture fraction. This is also observed in our previous study without the inclusion of TCI [74] .
A detailed analysis regarding the flame location can be visualized via Fig. 13 , where the temperature distribution is depicted against radial distance from the spray axis, and color coded by mixture fraction. Before the quasi-steady stage ( t < 3 τ * ), higher temperatures always appear away from the spray core region, showing that ignition in Spray C and Spray D occurs at the periphery of spray. Also denoted here is that the relatively higher temperatures typically appear where the mixture fraction is lower than 0.2. At the same time, a relatively low temperature is observed in the near nozzle region (i.e. this area is primarily represented by a magenta color at 1.01 τ * in Fig. 12 ).
Following Maes et al. [13] , the transient development of the flames is evaluated through I xt -plots. These plots show the radially integrated intensity of a certain species along the spray axis x against time t: where i represents a species. In this work, the intensity corresponds to the mass fractions of the species. In Fig. 14 the I xt -plots of CH 2 O, an important intermediate species formed during firststage ignition, and OH, an important marker of high-temperature combustion for Spray C and Spray D (single realization LES) are depicted. In both cases, CH 2 O formation is slightly earlier than the temperature-based IDT as described in Section 5.3.2 (marked by the vertical dashed line), while OH appears at about the same time of the IDT. An OH-based IDT is introduced as a comparison, which is defined as the time 2% of maximum I xt,OH is achieved. The ignited location is marked by the magenta circle. It is seen that the OH-based IDT defined by I xt is consistent with the temperaturebased IDT.
The intensity of CH 2 O experiences a decrease after the IDT due to the consumption of formaldehyde in the pre-mixed burn phase, observed experimentally by Sim et al. [71] , Maes et al. [72] . Meanwhile, OH is always more luminous in the right-upper corner due to the developing spray flame in the downstream region. Regions of high-intensity CH 2 O and OH hardly coincide due to their different phase in the combustion process, in line with what has been observed experimentally for spray flames [13] . The high intensity OH in Spray D is clearly farther downstream than that of Spray C, corresponding to a longer LOL. It is noted that the high-intensity CH 2 O region expands not only downstream, but also closer to the nozzle outlet around the IDT. This illustrates that first ignition starts farther downstream from the location where it will eventually stabilize. Subsequently, during the ignition process, the cool-  flame wave described by Dahms et al. [36] , reaches upstream with a velocity that is higher than that of the spray. The phenomenon can also be observed through the ignition kernels (marked by the magenta dashed line) depicted in Figs. 10 and 11 . According to the I xt -plots, in both Spray C and Spray D, combustion starts around their respective LOL locations (indicated by CH 2 O) and finally stabilize at the LOL (indicated by OH).
In general, the analysis of temporal evolution reveals very similar behavior for Spray C and Spray D especially after the secondstage ignition. This is consistent with the observation that cavitation does not affect ignition behavior considerably in experiments of Spray C and Spray D [18,68] .

Quasi-steady flow
It is evident that soot emission is highly dependent on local conditions. As is pointed out by Pickett and Siebers [73] , soot formation is favorable in turbulent flames when the equivalence ratio at the lift-off length (LOL) is higher than 2. Our turbulence modeling approach shows very good agreement for mixture fraction distribution for Spray A. Thus, in this section the mixing prediction of reacting Spray C and Spray D are analyzed to quantitatively understand the soot emissions found in experiments. Figure 15 depicts the equivalence ratio distribution of reacting Spray C and Spray D conditions. The ensemble-averaged LES consists of 6 realizations and is time-averaged between 2 ms to 3 ms. The radial equivalence ratio distribution at the lift-off length locations are shown in Fig. 16 . As is shown, in both cases, the equivalence ratio at the LOL is above 2, implying that soot will most certainly be formed. Although the equivalence ratio of Spray D initiates farther downstream from the injector exit compared to that of Spray C, it is seen that the distribution at the LOL is similar to that of Spray C ( Fig. 16 ). Experiments at various ambient temperature levels concluded that Spray C generally produces more soot than Spray D [68] . But when gradually reducing the ambient temperature, a turning point is revealed roughly around the ambient temperature of 900 K (the condition investigated in this work), where Spray C starts to form slightly less amount of soot. Indeed, our results imply similar equivalence ratio distributions for Spray C and Spray D while clear differences are seen upstream. Further investigations on soot formation concerning the mixing differences at the lift-off length due to oxidizer temperature are needed to prove that the trend of reducing soot formation is captured appropriately. This will be part of a future study.

Conclusion
In the present work LES have been carried out within the FGM framework to understand the cavitation effect on ECN target conditions. A new dispersion model specifically designed for LES is implemented, which discards the need for tuning of model parameters. A non-viscosity LES applying this new dispersion modeling shows promising mixing prediction by a detailed validation of the nonreacting ECN Spray A target condition. Consequently, it was applied to Spray C and Spray D, achieving excellent agreement on liquid and vapor penetrations with the experiments. The spray modeling approach is further applied in reacting conditions. A top-hat filtering PDF is introduced in conjunction with an FGM that takes the effect of scalar dissipation into account.
Multiple realizations are performed to discard LES variations. The minimum number of LES realizations needed for several representative fields are explicitly investigated by using the structure similarity index. Furthermore, a detailed analysis of both global and local ignition behavior is carried out by comparing single LES realizations to the ensemble-averaged LES.
The main conclusions of this work can be summarized as follows: • Predictions for the baseline non-reacting ECN spray A target condition show excellent agreement to the experiments regarding both liquid and vapor penetrations, and spatial distribution of mixture fraction. Considering the significant impact of the global dispersion model parameters on turbulence [1] , the effect of our locally defined model constants is justified. The new dispersion model avoids elaborate tuning of the dispersion model parameters.
• The same breakup model constants are directly applied to Spray C and Spray D and show good agreement regarding both liquid and vapor penetration. By allowing a coarser mesh, the modeling approach shows the potential of reducing the commonly observed excessive local liquid void fraction which may induce irregular fluctuations during ignition. • For the reacting conditions the simulation perform well for Spray C and Spray D conditions. The potential of discarding injection-to-injection uncertainties are considered via multiple realizations. Analysis of global ignition behavior, including LOL and temperature-based IDT, suggests that a single realization is sufficient to target the experiments. Predictions provide similar IDTs between Spray C and Spray D which is consistent with experiments. The trend that the LOL of Spray C is approximately 3 mm shorter than Spray D is captured as well. Regarding local ignition behavior, it is noted that different number of LES realizations are needed for different quantities according to their similarity indexes. Temperature and mixture fraction need less than 4 realizations while 8 realizations are not sufficient for quantities with a large spatial variation such as OH. It is thus suggested that an radially integrated intensity is preferable when OH-based IDT and LOL are needed. • It is found that although the low-temperature combustion stage for Spray C and Spray D shows significant differences, their flame development during the second-stage ignition are very similar. In both cases, ignition kernels appear on the wings of the flame rather than at the spray tip such as for Spray A. This is also reported in previous simulations by other researches [35,65] . • By analyzing multiple realizations, it is found that an ensembleaveraged LES results in a more widely distributed reacting zone which is similar to a RANS, and its peak temperature is approximately 150 K lower than that a single realization. • In the experiments, the ambient temperature of 900 K investigated in this study seems to be the turning point where Spray C starts to form comparable amount of soot to Spray D. This correlates well with the observed similarity in the equivalence ratio distributions for Spray C and Spray D at their corresponding LOLs.
To conclude, the new dispersion model captures the mixing of a high-pressure spray very well, without any fitting of model parameters. Besides, relatively 605 larger grid cells are allowed. The coupling of FGM with the top-hat filtering predicts all combustion phenomena with high fidelity in both cavitating and noncavitating flow in ECN configurations. Multiple LES realizations are preferred regarding the detailed analysis of local flame behavior.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
In this section, the role of the LES-type dispersion model is discussed by comparing three cases. The first one does not include dispersion, the second one adopts the u local model, and the third one applies the dispersion model that is directly adapted from RANS. All three cases are based on a 0.125-mm mesh. Figure 17 depicts the predicted results for liquid and vapor penetration, while   t SOI Z x =20 mm (t ) dt / (t − t SOI ) (Spray A [52] ). Fig. 18 shows the spatial distribution of mixture fraction and temperature on the spray axis. The time-averaged liquid penetration for no u , u local and RANS-type u are 13.33, 14.97 and 16.05 mm, respectively. It is observed that both liquid and vapor penetration are over-predicted when dispersion is ignored. Accordingly, the mixture fraction distribution in the upstream region is significantly different from that given by the u local model. As is frequently done in literature, a solution can be fitting the KH-RT model parameter. However, a shorter liquid penetration, illustrating a confined breakup, corresponds to a lowered value of mixture fraction on the spray axis. Note that the LOL of standard Spray A is around 17 mm [11] . This will lead to an even more severe underprediction of mixture fraction at the LOL as is already seen with the current model constant applied. It is expected the stabilization of the flame is different from the prediction by using the u local model, therefore directly impacting soot formation. An alternative solution can be refining the mesh. In this way, the effect of dispersion as a consequence of the subgrid fluctuations is reduced. However, a well-known problem is that Lagrangian-Eulerian simulation hampers the mesh refinement, especially where parcels cluster together. Compared to not including the dispersion, the difference between a RANS-type dispersion and LES-type dispersion is minor in terms of Figs. 17 and 18 . Although the u local model gives the best vapor penetration, the prediction of a RANS-type dispersion is still acceptable. According to Fig. 18 , very similar fuel distributions are seen in the period between 1.2 to 3 ms. The slight difference in liquid penetration supports the observation in the paper by Tsang et al. [1] that the dispersion modeling hardly influences vapor behavior in the downstream area. However, an evaluation of the quasi-steady state may not be sufficient. The work by Tsang et al. proves that the LES-type dispersion provides reasonable stochastic droplet trajectories compared to a RANS-type model. They declared that the former is able to accurately capture the anisotropic feature seen in both mean and variance of subgrid dispersion velocity. Figure 19 focuses on a probing location located at 20 mm downstream of the nozzle outlet (on the spray axis). In this case, the mixture fraction is averaged in time up to the ac-tual time ( Z x,a v e (t) = t t SOI Z x =20 mm (t ) dt / (t − t SOI ) ) and is depicted in the figure. For visualization, the development of mixture fraction around its inception time is magnified. The results show a clear difference between the two types of dispersion model. The inception time at 20 mm from the injector tip is approximately 0.02 ms, which is on the order of 1/10 of the IDT in standard Spray A. This may result in different IDTs. As soot is very sensitive to its inception location, the LES-type dispersion model is expected to promote soot formation due to the accelerated mixture fraction development.