Statistical Analysis of the Reaction Progress Variable and Mixture Fraction Gradients in Flames Propagating into Droplet Mist: A Direct Numerical Simulation Analysis

ABSTRACT The statistics of reaction progress variable, , and mixture fraction, , and their gradients (i.e., and ) in flames propagating in droplet mist, where the fuel was supplied in the form of monodisperse droplets, have been analyzed for different values of turbulent velocity fluctuations (), droplet equivalence ratios (, and droplet diameters ( based on three-dimensional direct numerical simulations (DNS) in a canonical configuration under decaying turbulence. The combustion process in the gaseous phase has been found to take place predominantly in fuel-lean mode, even for . The probability of finding fuel-lean mixture increases with increasing initial droplet diameter due to slower evaporation of larger droplets. It has been shown that the joint probability density function (i.e., joint PDF) of and (i.e., ), cannot be approximated in terms of discrete delta functions throughout the flame brush for the cases considered here. Furthermore, the magnitude of cannot be adequately approximated by the product of marginal PDFs of , and variable, (i.e., ). The statistical properties of the Favre probability density functions (Favre-PDFs) of the mixture fraction, , and oxidizer-based reaction progress variable, , have been analyzed at several locations across the flame brush and a -function distribution has been found to capture the Favre-PDFs of and obtained from the DNS data. Furthermore, a log-normal distribution has been shown to capture the qualitative behaviors of the PDFs of the gradient of the mixture fraction and the gradient of the reaction progress variable, and , respectively, but discrepancies between the log-normal distribution and the DNS data were observed at the tails of PDFs. In addition, the interrelation between and was examined in terms of the PDFs of the cosine of the angle between them (i.e., and it was observed that most droplet cases exhibited much greater likelihood of positive values of than negative values. Finally, the joint PDF of and , , has been compared with that of P( (i.e., assuming statistical independence of and ) and a good level of agreement has been obtained. The bivariate log-normal distribution has been considered both assuming correlation between and and assuming no correlation for the purpose of modeling , and the variant with no correlation has been found to be more successful in capturing qualitative behavior of although quantitative discrepancies have been observed due to inaccuracies involved in parameterizing P( and P(by log-normal distributions.

the gaseous phase has been found to take place predominantly in fuellean mode, even for ϕ d >1. The probability of finding fuel-lean mixture increases with increasing initial droplet diameter due to slower evaporation of larger droplets. It has been shown that the joint probability density function (i.e., joint PDF) of and c (i.e., P ; c ð Þ), cannot be approximated in terms of discrete delta functions throughout the flame brush for the cases considered here. Furthermore, the magnitude of P ; c ð Þ cannot be adequately approximated by the product of marginal PDFs of , and variable, c (i.e., P ð Þ Á P c ð Þ). The statistical properties of the Favre probability density functions (Favre-PDFs) of the mixture fraction, , and oxidizer-based reaction progress variable, c, have been analyzed at several locations across the flame brush and a β-function distribution has been found to capture the Favre-PDFs of and c obtained from the DNS data. Furthermore, a log-normal distribution has been shown to capture the qualitative behaviors of the PDFs of the gradient of the mixture fraction and the gradient of the reaction progress variable, Ñ j j and Ñc j j, respectively, but discrepancies between the log-normal distribution and the DNS data were observed at the tails of PDFs. In addition, the interrelation between Ñ and Ñc was examined in terms of the PDFs of the cosine of the angle between them (i.e., cos θ ð ÞÞ and it was observed that most droplet cases exhibited much greater likelihood of positive values of cos θ ð Þ than negative values. Finally, the joint PDF of Ñ j j and Ñc j j, P Ñ j j; Ñc j j ð Þ , has been compared with that of P( Ñ j jÞÁ P Ñc j j ð Þ (i.e., assuming statistical independence of Ñ j j and Ñc j j) and a good level of agreement has been obtained. The bivariate log-normal distribution has been considered both assuming correlation between Ñ j j and Ñc j j and assuming no correlation for the purpose of modeling P Ñ j j; Ñc j j ð Þ , and the variant with no correlation has been found to be more successful in capturing qualitative behavior of P Ñ j j; Ñc j j ð Þ although quantitative discrepancies have been observed due to inaccuracies involved in parameterizing P( Ñ j jÞ and P ( Ñc j jÞby log-normal distributions.

Introduction
Flame propagation into turbulent droplet-laden mixtures is of pivotal relevance to several engineering applications ranging from Internal Combustion (IC) engines to aero-gas turbines to the prediction and control of hazards (Aggarwal, 1998;Heywood, 1998;Lefebvre, 1998). A number of experimental (Aggarwal and Sirignano, 1985;Ballal and Lefebvre, 1981;Burgoyne and Cohen, 1954;Faeth, 1987;Hayashi et al., 1976;Lawes and Saat, 2011;Nomura et al., 2000;Szekely and Faeth, 1983), analytical (Greenberg et al., 1998;Greenberg et al., 2009;Silverman et al., 1993), and computational (Fujita et al., 2013;Wacks et al., 2015;Watanbe et al., 2007Watanbe et al., , 2008 analyses concentrated on flame propagation in turbulent droplet-laden mixtures and indicated a complex physical interaction of evaporative heat and mass transfer, fluid dynamics, and combustion thermo-chemistry is responsible for flame droplet interaction. The experimental evidence suggested that there could be significant differences between the overall equivalence ratio (considering fuel in both liquid and gaseous phases) and gaseous equivalence ratio due to incomplete evaporation. Furthermore, droplet inertia along with the difference between overall and gaseous equivalence ratios could give rise to augmentation/ reduction of burning rate in quiescent and low turbulence conditions, but that these effects disappear for sufficiently high turbulence intensity (Aggarwal and Sirignano, 1985;Lawes and Saat, 2011;Nomura et al., 2000). Recently, direct numerical simulations (DNS) have made significant contributions to both the physical understanding and modeling of the combustion of turbulent droplet-laden mixtures (Luo et al., 2011;Miller and Bellan, 1999;Neophytou et al., 2010Neophytou et al., , 2011Neophytou et al., , 2012Reveillon and Demoulin, 2007;Reveillon and Vervisch, 2000;Sreedhara and Huh, 2007;Wandel, 2013Wandel, , 2014Wandel et al., 2009;Wang and Rutland, 2005;Xia and Luo, 2010). Neophytou and Mastorakos (2009) recently analyzed the effects of volatility, droplet diameter, and droplet equivalence ratio on burning velocity in one-dimensional (1D) laminar flames where fuel is supplied in the form of monodisperse droplets. Wacks et al. (2015) recently extended the analysis of Neophytou and Mastorakos (2009) for turbulent flames by carrying out 3D compressible DNS of freely propagating turbulent flame propagation into droplet-laden mixtures. The statistical behaviors of the gradients of mixture fraction and reaction progress variable c have been analyzed in several previous analyses on partially premixed and stratified flames (Bray et al., 2005;Malkeson et al., 2013;Ruan et al., 2012), where mixture fraction and reaction progress variable c are taken to be representative passive and active scalars. Although mixture fraction is not strictly a passive scalar in the case of spray combustion, the mixture inhomogeneity in the gaseous phase in spray flames is often characterized in terms of mixture fraction (Ge and Gutheil, 2006;Sadiki et al. 2012). Furthermore, the progress of chemical reaction in inhomogeneous mixtures in droplet combustion can be characterized by reaction progress variable c, which can be defined in terms of a suitable species mass fraction so that it assumes a value equal to zero in the unburned gas and monotonically increases to assume a value equal to unity in the fully burned gas. The interdependence of mixture fraction and reaction progress variable c plays a key role in the modeling of spray combustion using flamelet generated manifold (FGM) (Ma et al., 2014;Sadiki et al., 2012). Furthermore, the characterization of probability density functions (PDFs) of mixture fraction and reaction progress variable c and their respective gradients (i.e., Ñ and Ñc) in terms of presumed functions play pivotal roles in flamelet and conditional moment closures (CMC; Borghesi et al., 2011;Ge and Gutheil, 2006;Sadiki et al., 2012;Tyliszczak et al., 2014). Furthermore, the statistics of Ñand Ñc play key roles in determining the statistical behaviors of scalar dissipation rates N ¼ DÑ Á Ñ and N c ¼ DÑc Á Ñc and cross-scalar dissipation rate N c ¼ DÑc Á Ñ, which are necessary for the purpose of modeling turbulent combustion in partially premixed mixtures (Bray et al., 2005;Malkeson et al., 2013;Ribert et al., 2005;Robin et al., 2006;Ruan et al., 2012). These modeling issues will also be valid for droplet combustion because the evaporation of droplets will lead to partially premixed combustion in gaseous phase. Although different presumed PDF approaches have been used in the past in the context of spray combustion (Ge and Gutheil, 2006;Ma et al., 2014;Sadiki et al., 2012;Tyliszczak et al., 2014), the statistics of the interdependence of and reaction progress variable c and their gradients (i.e., Ñand Ñc) have received limited attention (Luo et al., 2011;Wandel, 2013Wandel, , 2014Xia and Luo, 2010) in existing literature. This deficit has been addressed here by analyzing the statistical behavior of and c, and their respective gradients (i.e., Ñand Ñc) based on 3D DNS data of freely propagating statistically planar turbulent flames propagating into droplet mist for different values of root-mean-square velocity fluctuation u 0 ; droplet diameter a d , and droplet equivalence ratio Þ st in which m F;d is the mass of fuel in droplet form, m O is the mass of oxidizer, and m F =m O ð Þ st is the ratio of the mass of fuel to the mass of oxidizer in the stoichiometric mixture. In this respect, the main objectives of the present analysis are: (1) To analyze the distributions of and c and their joint PDF within the flame front for turbulent spray flames for different values of u 0 ; a d , and ϕ d .
(2) To analyze the influences of u 0 ; a d , and ϕ d on the statistical behaviors of Ñ and Ñc and their interdependence. (3) To provide physical explanations for the aforementioned statistical behaviors and indicate modeling implications.
The remainder of the article takes the following form. The necessary mathematical background and information related to the numerical implementation are provided in the next two sections. Following this, the results are presented and subsequently discussed, and finally, the main findings are summarized and conclusions are drawn.

Mathematical background
A single-step irreversible Arrhenius-type chemical mechanism has been used for the purpose of carrying out the present extensive parametric analysis without an exorbitant computational cost: where s is the oxidizer-fuel ratio by mass (i.e., the mass of oxygen consumed per unit mass of fuel). The fuel reaction rate _ ω F is expressed as: where Y F and Y O are the fuel and oxygen mass fractions, respectively, and ρ is the gas density. In Eq.
(2), T, the nondimensional temperature; β, the Zeldovich number; α, a heat release parameter; and B Ã , the normalized pre-exponential factor, are given by the following expressions: whereT is the instantaneous dimensional temperature, T 0 is the unburned gas temperature, T ad ϕ g ¼1 ð Þ is the adiabatic flame temperature for the stoichiometric mixture, E ac is the activation energy, R 0 is the universal gas constant, B is the pre-exponential factor, and τ ¼ Here, the modified single step chemical mechanism proposed by Tarrazo et al. (2006) has been considered and in this framework the activation energy, E ac , and the heat of combustion are taken to be functions of the gaseous equivalence ratio, ϕ g , which results in more accurate predictions of the equivalence ratio ϕ g dependence of the unstrained laminar burning velocity S b ϕ g ð Þ in hydrocarbon-air flames, especially for fuel-rich mixtures. According to Tarrazo et al. (2006), the Zel'dovich number, β, is expressed as β ¼ 6f ϕ g where: f ϕ g ¼ 1:0 þ 8:250 ϕ g À 1:00 2 ; ϕ g 0:64 1:0 ; 0:64<ϕ g <1:07 1:0 þ 1:443 ϕ g À 1:07 2 ; ϕ g ! 1:07 According to Tarrazo et al. (2006), the heat release per unit mass of Þ is given by H ϕ g =H ϕ g ¼1 ¼ 1 for ϕ g 1 and H ϕ g =H ϕ g ¼1 ¼ 1 À α H ϕ g À 1 for ϕ g >1, where α H ¼ 0:18, and Y F0 ϕ g ð Þ and Y Fb ϕ g ð Þ are the fuel mass fraction in the unburned gas and fully burned gas, respectively, for a premixed flame of equivalence ratio ϕ g . For the present investigation, the Lewis numbers of all species are taken to be equal to unity and all species in the gaseous phase are considered to be perfect gases. Standard values have been taken for the ratio of specific heats (γ ¼ C g P =C g V ¼ 1:4, where C g P and C g V are the specific heats at constant pressure and volume for the gaseous mixture, respectively) and Prandtl number (Pr ¼ μC g P =λ ¼ 0:7, whereμ and λ are the dynamic viscosity and thermal conductivity, respectively).
The droplet transport equations used in several previous analyses (Luo et al., 2011;Neophytou et al., 2010Neophytou et al., , 2011Neophytou et al., , 2012Reveillon and Demoulin, 2007;Reveillon and Vervisch, 2000;Sreedhara and Huh, 2007;Wandel, 2013Wandel, , 2014Wandel et al., 2009;Wang and Rutland, 2005;Xia and Luo, 2010) have been adopted for this analysis. The position,x d , velocity,ũ d , diameter, a d , and temperature, T d , for individual droplets are tracked in Lagrangian manner, and the relevant transport equations for these quantities are given by: where L v is the latent heat of vaporization, and τ p d , τ u d , and τ T d are relaxation timescales associated with droplet velocity, diameter, and temperature, respectively, which are given by: where ρ d is the droplet density, C L p is the specific heat for the fuel in liquid phase, C u is the correction for drag coefficient and is expressed as (Clift et al., 1978): In Eq. (10), Re d is the droplet Reynolds number, Scis the Schmidt number, B d is the Spalding mass transfer number, Sh c is the corrected Sherwood number, and Nu c is the corrected Nusselt number, which are expressed as (Clift et al., 1978;Luo et al., 2011;Neophytou et al., 2010Neophytou et al., , 2011Neophytou et al., , 2012Reveillon and Vervisch, 2000;Wandel, 2013Wandel, , 2014Wandel et al., 2009;Xia and Luo, 2010): where Y s F is the value of fuel mass fraction Y F on the droplet surface. It is worth noting that the unity Lewis number assumption has been implicitly invoked in Eqs. (11)-(13). Furthermore, the Clausius-Clapeyron relation for the partial pressure of the fuel vapor at the droplet surface, p s F , is utilized to calculate the Spalding number B d : where T s ref is the boiling point of the fuel at reference pressure p ref , the droplet surface temperature T s d is taken to be T d , and W O and W F are the molecular weights of oxidizer and fuel, respectively.
The conservation equations in the gaseous phase can be generically written as (Neophytou et al., 2010(Neophytou et al., , 2011(Neophytou et al., , 2012Reveillon and Vervisch, 2000;Wandel, 2013Wandel, , 2014Wandel et al., 2009): for the transport equations of mass, momentum, energy, and mass fractions, respectively; and ψ ¼ e, respectively, with u j and e ¼ òT being the velocity in the jth direction and the specific stagnation internal energy respectively. The quantity σ ψ is an appropriate Schmidt number corresponding to ψ. The _ ω ψ term in Eq. (15a) arises due to chemical reaction rate, and _ S g and _ S ψ are the appropriate source terms in the gaseous phase and are due to droplet evaporation, respectively. The droplet source term _ S ψ is tri-linearly interpolated from the droplet's sub-grid position,x d , to the eight surrounding nodes. The droplet source term for any variable ψ may be expressed as (Neophytou et al., 2010(Neophytou et al., , 2011(Neophytou et al., , 2012Wandel, 2013Wandel, , 2014Wandel et al., 2009): where V is the cell volume, m d ¼ ρ d 1 6 πa 3 d is the droplet mass and the summation is carried out over all droplets in the vicinity of each node. As indicated in Eq. (15a), the variable ψ is identified as ψ ¼ 1; u j ; e; Y F ; Y O È É , however, since within the droplets Y F ¼ 1:0, the source term for both the continuity equation and the fuel mass fraction equation are identical.
Evaporation of droplets leads to mixture inhomogeneities, which is often quantified in terms of the mixture fraction, which is defined as: where Y F1 ¼ 1:0 is the fuel mass fraction in the pure fuel stream and Y O1 ¼ 0:233 is the oxidizer mass fraction in air. The hydrocarbon fuel used in this DNS analysis is n-heptane, C 7 H 16 , for which s ¼ 3:52 and the stoichiometric fuel mass fraction and mixture fraction values are given by Y Fst ¼ st ¼ 0:0621. Furthermore, it is possible to define a reaction progress variable,c, based on a species mass fraction and the mixture fraction such that c rises monotonically from zero in the unburned reactants to one in the fully burned products. Here, the progress variable, c, has been defined in terms of oxidizer mass fraction in the following manner following several previous analyses (Neophytou et al., 2010(Neophytou et al., , 2011(Neophytou et al., , 2012Wandel, 2014;Wandel et al., 2009): From Eq. (17) it is possible to derive a transport equation of c based on the transport equations for the oxidizer mass fraction Y O and the mixture fraction (Bray et al., 2005;Wacks et al., 2015): where D is the mass diffusivity and the first term on the right-hand side of Eq. (18) arises due to molecular diffusion, the second term represents reaction rate, the third is the source/sink term arising due to droplet evaporation, and the last is the cross-scalar dissipation term arising due to reactant inhomogeneity (Bray et al., 2005;Malkeson and Chakraborty, 2010;Wacks et al., 2015). The cross-scalar dissipation term _ A c in Eq. (18) arises due to mixture inhomogeneity (Bray et al., 2005;Malkeson and Chakraborty, 2010;Wacks et al., 2015), which, in this case, is induced by droplet evaporation. According to the definition of c (see Eq. (17)), the definitions of _ ω c , _ S c , and _ A c depend on the local value of . The expressions for _ S c and _ A c are given by the following expressions (Wacks et al., 2015): Þis the droplet source/sink term in the mixture fraction transport equation and _ are the droplet source/sink terms in the fuel and oxidizer transport equations, respectively, and Γ m is the source term in the mass conservation equation due to evaporation.
The reaction rate _ ω c of the reaction progress variable may be expressed as (Bray et al., 2005;Malkeson et al., 2013;Ruan et al., 2012): where N ζ 1 ζ 2 ¼ DÑζ 1 Á Ñζ 2 is the (cross-)scalar dissipation rate, for scalars ζ 1 and ζ 2 (here c or ). It is evident from Eq. (21) that the evaluation of (cross-)scalar dissipation rates is essential to accurately model _ ω c . In the context of Reynolds averaged Navier-Stokes (RANS) simulations, the Favre-averaged (cross-)scalar dissipation rates, ζ 1 ζ 2 ¼ ρDÑζ 00 1 Á Ñζ 00 2 = p (where the Favre mean and Favre fluctuation of a scalar q are given by q ¼ ρq= ρ and q 00 ¼ q Àq, respectively, with the overbar indicating a suitable Reynolds averaging operation) must be modeled. For RANS, the components of scalar dissipation and cross-scalar dissipation rates arising from the mean gradients are negligible in comparison to the contributions from gradients of fluctuations (i.e., ρDÑc:Ñc ) ρDÑc Ñc, ρDÑ:Ñ ) ρDÑ Á Ñ, and ρDÑc Á Ñ ) ρDÑc Á Ñ), which gives rise to the following expressions: cc % ρDÑc 00 Á Ñc 00 = ρ ¼Ñ cc % ρDÑ 00 Á Ñ 00 = ρ ¼Ñ and c % ρDÑc 00 Á Ñ 00 = ρ ¼Ñ c The scalar dissipation ratesÑ cc ;Ñ , andÑ c can alternatively be expressed as: where P q 1 ; q 2 ð Þ is the joint PDF of quantities q 1 and q 2 . Equations (23a)-(23c) indicate that the statistical behaviors of c, , Ñc j j, and Ñ j j are of fundamental interest for the purpose of modeling of cc , , and c (Malkeson et al., 2013;Ruan et al., 2012). The statistical behaviors of c, , Ñc j j, and Ñ j jwill be discussed in the Results and discussion section of this article.

Numerical implementation
A widely-used 3D compressible DNS code SENGA (Jenkins and Cant, 1999;Neophytou et al., 2010Neophytou et al., , 2011Neophytou et al., , 2012Wandel, 2013Wandel, , 2014Wandel et al., 2009;Wacks et al., 2015), which solves the standard transport equations of mass, momentum, energy, and species in nondimensional form has been employed for this analysis. In this framework, the spatial discretization for the internal grid points has been carried out using a 10th-order central difference scheme, where the order of differentiation drops gradually to a one-sided 2nd-order scheme at the nonperiodic boundaries (Jenkins and Cant, 1999). A low-storage third-order explicit Runge-Kutta scheme has been used for time advancement (Wray, 1990). A rectangular computational domain of size 63: Þ has been considered for the current investigation, where D 0 is the unburned gas diffusivity. For the present thermo-chemistry the Zel'dovich flame thickness D 0 =S b ϕ g ¼1 ð Þ is equal to about 0:625δ th where δ th ¼ T ad ϕ g ¼1 ð Þ À T 0 = max ÑT À Á L is the unstrained thermal laminar flame thickness of the stoichiometric laminar flame, and the subscript L refers to the values in an unstrained laminar premixed flame for the stoichiometric mixture. The simulation domain for the present analysis is discretized using a Cartesian grid of size 384 Â 256 Â 256, which resolves both the flame thickness, δ th , and the Kolmogorov length-scale, η. The boundaries in the mean direction of flame propagation (i.e., x-direction) are taken to be partially nonreflecting, whereas the transverse (i.e., y-and z-) directions are assumed to be periodic. The Navier-Stokes characteristic boundary conditions (NSCBC) (Poinsot and Lele, 1992) technique has been used for specifying the nonperiodic boundary conditions. The droplets are distributed uniformly in space throughout the y-and z-directions and in the region 0:0 ð Þ =D 0 16:53 ahead of the flame. The initial conditions for the reacting flow field have been generated based on the steady laminar solution obtained for the desired initial values of droplet diameter, a d , and droplet equivalence ratio, ϕ d . This initial condition has been generated using COSILAB (Neophytou and Mastorakos, 2009), where the 1D governing equations for the gas and liquid phases are solved in a coupled manner for spray flames where fuel is supplied in the form of monodisperse droplets on the unburned gas side of the flame. The turbulent velocity fluctuations have been initialized using a standard pseudo-spectral method (Rogallo, 1981), and this field is superimposed on the steady laminar spray flame solution generated using COSILAB. For the present analysis the unburned gas temperature is For all simulations, the fuel is supplied purely in the form of monodisperse droplets with nondimensional diameters a d =δ th ¼ 0:06; 0:08; 0:10 for different values of droplet equivalence ratio: ϕ d ¼ 1:0; 1:25; 1:5; 1:7 at a distance 16:53D 0 =S b ϕ g ¼1 ð Þ from the point in the laminar flame at whichT ¼ 400K, which corresponds to a nondimensional temperature T % 0:05. The droplet number density ρ N at t ¼ 0 varies between 1:16 ρ N À Á 1=3 δ th 2:27 in the region 0:0 xS b ϕ g ¼1 ð Þ =D 0 16:53. In all cases the liquid volume fraction remains much smaller than 0:01. In all cases droplets are supplied at the left-hand-side boundary to maintain a constant ϕ d ahead of the flame. The droplets evaporate as they approach the flame front and the droplet diameter decreases by at least 25% by the time it reaches the most reactive region of the flame, such that the volume of even the largest droplets is now less than half that of the cell volume, which validates the sub-grid point source treatment of droplets adopted for flame-droplet interactions analyzed here since this study is concerned primarily with regions where reaction rate is non-negligible. This droplet to cell size for the present analysis is consistent with several previous analyses (Neophytou et al., 2010(Neophytou et al., , 2011(Neophytou et al., , 2012Reveillon and Vervisch, 2000;Sreedhara and Huh, 2007;Wandel, 2013Wandel, , 2014Wandel et al., 2009;Wang and Rutland, 2005).
The simulations are carried out for normalized root-mean-square (rms) turbulent velocities u 0 =S b ϕ g ¼1 ð Þ ¼ 4:0 and 7:5 with a nondimensional longitudinal integral length-scale L 11 =δ th ¼ 2:5. The ratio of droplet diameter to the Kolmogorov scale is a d =η % 0:3; 0:4; 0:5 for a d =δ th % 0:06; 0:08; 0:1, respectively, for initial u 0 =S b ϕ g ¼1 ð Þ ¼ 7:5, and these ratios are larger in magnitude for the cases with u 0 =S b ϕ g ¼1 ð Þ ¼ 4:0. The ratio of a d =η remains comparable to several previous analyses (Neophytou et al., 2010(Neophytou et al., , 2011(Neophytou et al., , 2012Reveillon and Vervisch, 2000;Sreedhara and Huh, 2007;Wandel, 2013Wandel, , 2014Wandel et al., 2009;Wang and Rutland, 2005). The mean normalized inter-droplet distance s d =η ranges between 0.0220 and 0.0432 (i.e., 0:0220 < s d =η < 0:0432) for initial u 0 =S b ϕ g ¼1 ð Þ ¼ 7:5 cases. All simulations have been carried out until t final ¼ max 3t turb ; 4t chem ð Þ , where t turb ¼ L 11 =u 0 is the initial turbulent eddy turnover time and t chem ¼ D is the chemical timescale. This simulation time is either comparable to or greater than the simulation duration used in a number of recent DNS analyses (Neophytou et al., 2010(Neophytou et al., , 2011(Neophytou et al., , 2012Reveillon and Vervisch, 2000;Sreedhara and Huh, 2007;Wandel, 2013Wandel, , 2014Wandel et al., 2009;Wang and Rutland, 2005), which significantly contributed to the fundamental understanding of turbulent combustion. It was shown by Wacks et al. (2015) that the volume-integrated reaction rate, flame surface area, and burning rate per unit area were not changing rapidly when the statistics have been extracted. This information is not repeated here for the sake of conciseness. The Reynolds/Favre averaged value of a general quantity Q (i.e., Q andQ) is evaluated by ensemble averaging the relevant quantity Q over the y-z plane at a given x location. The statistical convergence of the Reynolds/Favre averaged values has been assessed by comparing the values obtained on full sample size with the corresponding values based on distinct half of the available sample size and a satisfactory level of agreement has been obtained. The values based on full sample size will be reported in the next section for the sake of conciseness.

Flame turbulence interaction
A selection of instanteneous distributions of nondimensional temperature T, mixture fraction , and reaction progress variable c in the central x À z midplane at t % 4t chem are shown in Figure 1 for each droplet size with the droplet equivalence ratio ϕ d ¼ 1:0 and initial turbulent intensity u 0 =S b ϕ g ¼1 ð Þ ¼ 7:5. The droplets, which are depicted by black dots, are those residing in the cells immediately above or below the plane shown in the figure. It is evident from Figure 1 that the droplets shrink due to evaporation as they approach the flame. It is also noteworthy that several droplets of each size can be seen to have penetrated the flame front. The droplets continue to evaporate in the burned gas region and some of the evaporated gaseous fuel eventually diffuses back towards the flame front. Neophytou and Mastorakos (2009) reported pyrolysis of the droplets in the burned gas side in the absence of sufficient oxygen for laminar 1D calculations. However, Kuo and Acharya (2012) indicated that the flames with small value of group number G usually exhibit high temperatures, which promote pyrolysis at the fuel-rich core, whereas the temperature values in the external group combustion are usually not high enough to give rise to a significant amount of pyrolysis. The DNS methodology adopted here is for large values of group number (Reveillon and Vervisch, 2000) and Figure 1 shows that the temperature of the burned gas in the droplet cases shown here is considerably lower than that of the corresponding stoichiometric turbulent premixed flame (T % 1). The reduction in temperature originates predominantly due to the fuel-lean mode of combustion and due to extraction of latent heat by the droplets. This is consistent with previous findings based on unsteady laminar and 2D turbulent simulations (Fujita et al., 2013;Nakamura et al., 2005;Watanabe et al., 2007Watanabe et al., , 2008. Due to high values of group number and low values of burned gas temperature, the effects of pyrolysis are kept beyond the scope of the current analysis, which employs only a modified single-step Arrhenius-type chemical mechanism (due to the exorbitant computational costs involved in more detailed chemical mechanisms), which is not sufficiently detailed to model the pyrolysis process. Furthermore, pyrolysis does not directly affect the statistics of scalar gradients, which is the main focus of the current analysis. A similar assumption was made in several previous analyses (Luo et al., 2011;Neophytou et al., 2011Neophytou et al., , 2012Reveillon and Demoulin, 2007;Reveillon and Vervisch, 2000;Sreedhara and Huh, 2007;Wandel, 2013Wandel, , 2014Wandel et al., 2009;Wang and Rutland, 2005;Xia and Luo, 2010).
A comparison between the mixture fraction field and reaction progress variable c field indicates that, almost without exception, the inhomogeneous mixture arising from evaporation remains fuel-lean (i.e., < st % 0:0621) on the unburned gas side. The reaction progress variable c fields in Figure 1 show important differences that arise due to the change of initial droplet diameter a d , although turbulent intensity and ϕ d remain constant. It can be seen from Figure 1 that for the smallest droplets most c isosurfaces lie very close together, whereas, the isosurfaces are increasingly separated from each other for the medium and large droplets, indicating a broadening of the flame in these cases.
The simulation methodology adopted here is restricted to situations where the droplet size is smaller than the Kolmogorov length scale. Thus, the simulation methodology adopted here is valid only for the external group combustion and external sheath group combustion according to the regime diagram by Chiu et al. (1982), which was discussed in detail by Reveillon and Vervisch (2000). Chiu and Liu (1977) defined a group number, LeN 2=3 a d =s d ð Þ (where Le and Sc are the Lewis and Schmidt numbers, respectively; N is the number of droplets in a specified volume; and s d is the mean inter-droplet distance) in order to distinguish between individually burning droplets (G ( 1:0) and external sheath combustion (G ) 1:0). All droplet cases considered here come under the category of external group combustion (i.e., have values of G much greater than unity). While Figure 1 gives an impression that individually burning droplets are present in the burned gas region, it should be noted that Figure 1 shows only one plane and that other droplets may reside in adjacent planes. The regime of combustion can be characterized with the help of the Karlovitz number , which provides a measure of the ratio of flame thickness to the Kolmogorov length scale (Peters, 2000). For a stoichiometric mixture with values of u 0 =S b ϕ g ¼1 ð Þ ¼ 7:5 and L 11 =δ th ¼ 2:5 one obtains a Karlovitz number of 9.0, and Ka ¼ 3:5 for the values of u 0 =S b ϕ g ¼1 ð Þ ¼ 4:0 and L 11 =δ th ¼ 2:5. Thus, the value of the Karlovitz number is likely to be large under fuel-lean conditions (since S bðϕ g <1Þ <S b ϕ g ¼1 ð Þ ), which suggests that combustion in all cases considered here takes place well within the thin reaction zones regime (Peters, 2000). Furthermore, the Damköhler number Da ¼ t turb =t chem ¼ L 11 S 2 b ϕ g ð Þ =u 0 D 0 scales as Da , Re 1=2 t =Ka (Peters, 2000) and, hence, low Damköhler number combustion is likely for the high Karlovitz number cases considered here. This can further be substantiated from the variation of f c 00 2 withc shown in Figure 2.
According to Bray et al. (1985) one obtains: f c 00 2 ¼c 1 Àc ð Þþ O 1=Da ð Þwhere the last term is negligible for high values of Damköhler number. It is possible to obtain maximum possible values of the variance f c 00 2 ¼c 1 Àc ð ) where the PDF of c (i.e., P c ð Þ) shows a bimodal distribution with delta functions at c ¼ 0:0 and 1:0, which takes the following form: where α c ; β c , and γ c are the weights associated with the PDF contributions and f b c ð Þ is the burning mode PDF, which originates from the interior of the flame. The third term on the right-hand side of Eq. (24) is of the order of 1=Da ð Þ (Bray et al., 1985) and thus this contribution can be ignored for Da ) 1. However, since the chemical reaction in the gaseous phase takes place predominantly in a fuel-lean mode for the droplet cases considered here, Da is expected to be small (due to small values of S b ϕ g ð Þ for fuel-lean mixtures), which leads to a significant departure of f c 00 2 fromc 1 Àc ð Þin turbulent droplet cases (i.e., it is expected that P c ð Þ departs significantly from the aforementioned bi-modal distribution). It is evident from Figure 2 that the low Damköhler number effects (i.e., effects of slow chemical reaction) are stronger for the cases with large and medium droplets and weaker for the cases with small droplets. This is once again due to the fact that smaller droplets evaporate more quickly, allowing more fuel vapor to be released than in the case of larger droplets. Similarly, the effects of increasing ϕ d are less noticeable for medium and large droplets, although this does lead to slight increases in f c 00 2 for these droplets. Whereas the effects for small droplets are much more noticeable, in most cases an increase in ϕ d leads in turn to an increase in f c 00 2 (and also faster chemical reaction) due to the extra gaseous fuel that is released. The exception being ϕ d ¼ 1:70 (for small droplets) where the abundance of gaseous fuel leads to fuel-rich mixture, which reduces the rate of chemical reaction.

Distributions of c and and their modeling
It is often necessary to evaluate f c 00 2 , f 00 2 , and g c 00 00 for the purpose of modeling2 cc ,2 , and c in turbulent partially premixed combustion according to the following expressions (Malkeson et al., 2013;Ruan et al., 2012): wherek ¼ ρu 00 i u i 00 =2 ρ and ¼ μð@u 00 i =@x j Þð@u 00 i =@x j Þ= ρ are the turbulent kinetic energy and its dissipation rate, respectively, and C cc ; C , and C c are the model parameters. Irrespective of the performance and validity of the expressions given by Eq. (25) for spray combustion, it is worth noting that f c 00 2 , f 00 2 , and g c 00 00 are unclosed terms and they can be expressed as: whereP q ð Þ ¼ ρ ð Þ q P q ð Þ= ρ is the Favre-PDF of a general quantity; ρ ð Þ q is the mean gas density conditional on q with P q ð Þ being the marginal PDF of q. Equations (23), (25), and (26) indicate that the knowledge of P c; ð Þ,P ð Þ, andP c ð Þ is necessary in addition to the PDFs of Ñc j j and Ñ j j for the purpose of modeling cc , , and c . Figure 3a shows the joint PDF of c and forc ¼ 0:1; 0:3; 0:5; 0:7; 0:9, which is denoted as P c; ð Þ, for the droplet cases a d =δ th ¼ 0:06; 0:08; 0:10 with ϕ d ¼ 1:00 and under turbulent flow conditions u 0 =S b ϕ g ¼1 ð Þ ¼ 7:5. Although P c; ð Þ shows considerable peaks near c % 0:0 onc ¼ 0:1 and near c % 1:0 onc ¼ 0:5; 0:7; 0:9, it is clear from Figure 3a that the distribution of c departs significantly from the presumed PDF of c given by Eq. (24). This is especially true for lower values ofc, which reside predominantly in regions of fuel-lean mixture, as has been noted with regard to Figure 1, resulting in a much slower combustion process with a low Damköhler number. This phenomenon is apparent for all droplet sizes considered here. The reason for this is that even the smallest droplets considered here are not sufficiently small to evaporate quickly enough to produce significant quantities of stoichiometric mixture on low values ofc: Figure 3a shows that the joint PDF P c; ð Þ cannot be accurately approximated in terms of discrete delta functions throughout the flame brush as was used previously by Ribert et al. (2005), Robin et al. (2006), and Mura et al. (2007) for approximating the joint PDF between Y F and (i.e., P Y F ; ð Þ). The joint PDFs between Y F and (i.e., P Y F ; ð Þ) have also been found to be inadequately represented by discrete delta functions (not shown here), which is qualitatively similar to P c; ð Þ. However, the profiles of P c; ð Þ begin to show some attributes of delta functions asc increases. This process occurs most rapidly for the small droplets, which evaporate more readily and, hence, are more likely to produce regions of stoichiometric mixture than larger droplets, such that onc ¼ 0:9 the profile of P c; ð Þ for the small droplets is confined to c % 1:0 and a narrow band of values bounding ¼ st . The larger droplets show the same trend, but due to their slower evaporation, do not resemble delta functions as much as the small droplets.
It can further be seen from Figure 3a that the peak value of the joint PDF P c; ð Þ is obtained at % st ¼ 0:0621 and c % 1:0 forc ¼ 0:5; 0:7; 0:9. The peak of P c; ð Þ at = st % 1:0 arises principally from the region c > 0:90 due to the high rate of droplet evaporation as a result of the high temperature in the burned gas region. Furthermore, the peak value of P c; ð Þ is obtained at % st ¼ 0:0621 and c % 1:0 due to the diffusion-type flame which develops as a result of droplet evaporation in the burned gas region and the subsequent back-diffusion of the evaporated fuel (Wacks et al., 2015).
The joint PDFs P c; ð Þ for the variation of ϕ d and u 0 =S b ϕ g ¼1 ð Þ are not shown here for the sake of brevity. It can be reported, however, that as ϕ d increases, for the same initial turbulent flow conditions u 0 =S b ϕ g ¼1 ð Þ ¼ 7:5, there is a shift in the position of the peak which is located at c % 1:0 (onc ¼ 0:5; 0:7; 0:9) towards higher values of , as may be anticipated due to the extra availability of fuel in liquid form, which in turn leads to more gaseous fuel. In addition, the peak at c % 1:0 is less well defined with the values being spread over a wider range of . As the turbulent intensity is decreased to u 0 =S b ϕ g ¼1 ð Þ ¼ 4:0, the joint PDFs P c; ð Þ show higher values of across c due to the lack of sufficient turbulent mixing. This effect is compounded by increasing the value of ϕ d and is most noticeable for small droplets, which evaporate more easily than the medium or large droplets. The relatively high likelihood of attaining > st signifies the development of regions in which the reaction proceeds slowly due to the fuel-rich mixture.
As was mentioned above, the evaluation of g c 00 00 can be expressed as a function of P c; ð Þ (see Eq. (26c)). In the cases of statistical independence of c and , the joint PDF P c; ð Þ can be approximated by P c ð Þ Á P ð Þ, which would then simplify the evaluation of g c 00 00 . Figures 3b and 3c compare the contours of P c; ð Þ and P c ð Þ Á P ð Þ onc ¼ 0:1 and c ¼ 0:3 for a representative droplet case. The contours on higher values ofc are not shown as they do not provide useful information due to the sharp peak of the joint PDFs. It can be seen from Figures 3b and 3c that P c; ð Þ and P c ð Þ Á P ð Þ are somewhat similar in appearance, but different in magnitude on these values ofc. However, it has already been discussed that high values of c (i.e., c % 1:0) are more likely for % st ¼ 0:0621, thus it cannot be expected that c and are statistically independent of each other. Accordingly, it may not be appropriate to model P c; ð Þ by P c ð Þ Á P ð Þ for the purpose of simplifying the integral given by Eq. (26c).
The Favre-PDFs of the mixture fraction (i.e.,P ð Þ) are shown in Figure 4 for all droplet sizes (a d =δ th ¼ 0:06; 0:08; 0:10) and for both turbulent intensities (u 0 =S b ϕ g ¼1 ð Þ ¼ 4:0; 7:5) for the lowest and highest values of droplet equivalence ratio (ϕ d ¼ 1:00; 1:70) at several locations across the flame brush (c ¼ 0:1; 0:3; 0:5; 0:7; 0:9). Figure 4 shows that the nature of the Favre-PDFs varies significantly with position in the flame and that in all cases shown here the most likely mixture shifts from fuel-lean on the unburned gas side (c ¼ 0:1) to stoichiometric ( st ¼ 0:0621Þ/fuel-rich on the burned gas side (c ¼ 0:9). In the mid-region (c ¼ 0:3; 0:5; 0:7), the progression from fuel-lean to stoichiometric/fuel-rich is achieved most swiftly by the smaller droplets due to a higher rate of evaporation, whereas larger droplet cases remain fuel-lean for a greater range ofc. Furthermore, the Favre-PDFs of for cases with lower turbulence intensity (u 0 =S b ϕ g ¼1 ð Þ ¼ 4:0), although very often following closely the Favre-PDFs of the cases with higher turbulence intensity (u 0 =S b ϕ g ¼1 ð Þ ¼ 7:5) both in shape and magnitude, show consistently a higher probability of a more fuel-lean mixture at most locations in the flame. The noteable exception to this is for small droplets at the higher droplet equivalence ratio (ϕ d ¼ 1:70). Under these conditions, the abundance of small droplets evaporate quickly and are not mixed sufficiently well by the flow of lower turbulent intensity, giving rise to a greater likelihood of fuel-rich regions than for the higher turbulent intensity, which is strong enough to mix the evaporated fuel. In general, the higher droplet equivalence ratio (ϕ d ¼ 1:70) gives rise to a greater spread of values of than the lower droplet equivalence ratio (ϕ d ¼ 1:00) and, in particular, to non-negligible probabilities of attaining fuel-rich mixtures. This phenomenon arises naturally from the higher probability of droplets being in closer proximity (due to higher droplet number density) for higher droplet equivalence ratios, thereby giving rise to higher concentrations of gaseous fuel, despite the effects of turbulent mixing. In order to model the Favre-PDF of , a presumed β-function PDF approach is often considered in the context of single gaseous phase combustion (Poinsot and Veynante, 2001). The β-function PDF of a scalar ζ is defined as: Figure 4 shows that the assumptionP ð Þ ¼ f ζ; a; b ð Þreasonably captures the general shape and magnitude of the Favre-PDFs obtained from the DNS data.
Indeed, the agreement of the β-function PDF withP ð Þ is consistent with the modeling assumptions made by Sadiki et al. (2012) and Ma et al. (2014) in the context of FGMbased simulations of turbulent spray combustion. However, it must be emphasized that the use of the β-function PDF to modelP ð Þ may not be suitable in all cases of droplet combustion, as has already been shown by Ge and Gutheil (2006) and Luo et al. (2011). Ge and Gutheil (2006) proposed a modification to the standard β-function PDF: f ζ; a; b ð Þ¼1=B a; b ð Þ ζ max À ζ min ð Þ 1ÀaÀb ζ À ζ min ð Þ aÀ1 ζ max À ζ ð Þ bÀ1 , where ζ max and ζ min refer to maximum and minimum values of ζ, which has recently been used by Tyliszczak et al. (2014) for LES of spray combustion. It is also worth noting that Borghesi et al. (2011) parameterizedP ð Þ in terms of the β-function PDF in the context of CMC modeling and assumed thatP ð Þ can be approximated as a β-function PDF when the variance of mixture fraction f 00 2 is greater than a threshold limit given by: ðY Fs ÀÞ 2:5 ffiffiffiffiffiffi ffi f 00 2 q whenÞY Fs and max is considered to be Y Fs (see Eq. (13)), whereas min can be considered to be zero (Tyliszczak et al., 2014). In order words, the modified β-function PDF proposed by Ge and Gutheil (2006) (i.e., f ζ; a; b ð Þ¼1=B a; b ð Þ ζ max À ζ min ð Þ 1ÀaÀb ζ À ζ min ð Þ aÀ1 ζ max À ζ ð Þ bÀ1 ) reduces to Eq. (27) for small magnitudes of 1 À max ð Þ¼ ð1 À Y Fs Þ and min À 0 ð Þ . In the cases considered here the magnitude of ð1 À Y Fs Þ well below unity (i.e., ð1 À Y Fs Þ ( 1:0) within the flame brush for n-heptane combustion and, thus, the standard β-function PDF predictsP ð Þ satisfactorily for the cases considered here. However, Eq. (27) may not be sufficient for a different fuel and flow conditions where Y Fs is much smaller than unity. Figure 5 shows the Favre-PDFs of c (i.e.,P c ð Þ) for all droplet sizes (a d =δ th ¼ 0:06; 0:08; 0:10) and for both turbulent intensities (u 0 =S b ϕ g ¼1 ð Þ ¼ 4:0; 7:5) for the lowest and highest values of droplet equivalence ratio (ϕ d ¼ 1:00; 1:70) at several locations across the flame brush (c ¼ 0:1; 0:3; 0:5; 0:7; 0:9). All locations in the flame shown here exhibit an almost monomodal distribution. Onc ¼ 0:1 it is most likely to achieve a value c ¼ 0. This is reflected to some extent onc ¼ 0:3, where the most likely outcome is slightly greater than c ¼ 0. In contrast, onc ¼ 0:5; 0:7; 0:9 it is most likely to achieve a value c ¼ 1, with the likelihood increasing asc increases fromc ¼ 0:5 toc ¼ 0:9. The likelihood of achieving a value c ¼ 0 onc ¼ 0:1 is greater for small droplets than for large droplets. This feature is reversed onc ¼ 0:9, where the likelihood of achieving c ¼ 1 is greater for large droplets than for small droplets. This is true at both values of droplet equivalence ratio shown here and for both values of turbulent intensity. It is clear from Figure 5 that the Favre-PDF of c cannot be approximated by Eq. (24), whereas the β-function PDF (see Eq. (27)) parameterized in terms ofc and f c 00 2 (i.e.,P c ð Þ ¼ f c; a; b ð Þ with a ¼ e c e c 1 À e c ð Þ= f c 00 2 À 1 h i and b ¼ a=e c À a) has been found to reasonably capture the general shape and magnitude of the reaction progress variable Favre-PDFs obtained from DNS data. This is consistent with the modeling assumption made by Ma et al. (2014) in the context of FGM based modeling of turbulent spray combustion. It is worth noting a number of analyses on single-phase homogeneous (Ahmed and Swaminathan, 2013;Kolla and Swaminathan, 2010) and inhomogeneous (Chen et al., 2015;Ruan et al., 2014) mixture combustion used β-function PDF for parameterisingP c ð Þ in the past.
Statistical behavior and the modeling of the PDFs of Ñ j j and Ñc j j Figures 6 and 7 show the PDFs of Ñ j j and Ñc j j, respectively. Often a presumed lognormal distribution is used to model these PDFs. The log-normal distribution of a scalar, ζ, is given by: where μ 0 and σ 0 are mean and standard deviation of the natural logarithm of ζ. Figures 6  and 7 show that both qualitatively and quantitatively the log-normal distribution succeeds in capturing the PDFs of the DNS data for all values of ϕ d ; a d , and u 0 =S b ϕ g ¼1 ð Þ . The only exception takes place at high values of Ñ j j and Ñc j j where there is insufficient data in the DNS dataset to accurately predict the value of the PDF and the regions with high probability of finding Ñc j j % 0, which log-normal distribution fails to predict. The discrepancy between the log-normal distribution and the tail of the scalar gradient PDFs is consistent with several experimental (Geyer et al., 2005;Karpetis and Barlow, 2002;Markides and Mastorakos, 2006;Sreenivasan and Antonia, 1997;Su and Clemens, 2003) and computational (Hawkes et al., 2007;Knaus et al., 2012;Malkeson et al., 2013;Pantano, 2004;Ruan et al., 2012) findings in the context of passive scalar mixing, non-premixed, and partially-premixed combustion. Despite some inaccuracies, the log-normal distribution appears to do a reasonable job for the purpose of parameterising Ñ j j and Ñc j j, which in turn could be used for modeling the scalar dissipation rates and cc .

Interrelation between Ñ andÑc
It is important to understand the relative alignment between Ñ and Ñc so that the scalar inner product of Ñ and Ñc (i.e., Ñ Á Ñc) can be parameterized in terms of Ñ j j and Ñc j j in Eq. (23c) for the purpose of evaluatingÑ c (and c ). The relative alignment between Ñ j j and Ñc j j can be quantified in terms of the angle betweenÑ and Ñc: A distinction may be drawn between flames in partially-premixed mixtures where the flame-front is propagating towards a leaner mixture, also known as back-supported flames, and where the flame-front is propagating towards a richer mixture, also known as front-supported flames. Hence, it follows that front-supported flames arise when the signs of Ñ and Ñc are opposite to each other and back-supported flames arise when Ñ and Ñc possess the same sign. According to Eq. (29), cos θ ð Þ<0 implies the presence of front-supported flames and cos θ ð Þ>0 implies the presence of back-supported flame-fronts. Figure 8 shows the PDFs of cos θ ð Þ for all droplet sizes (a d =δ th ¼ 0:06; 0:08; 0:10) and for  Figure 6. Comparison of PDFs of normalized mixture fraction gradient Ñ j j Â δ th across the flame brush between DNS data (solid lines) and lognormal distribution (dashed lines) atc ¼ 0:1; 0:3; 0:5; 0:7; 0:9 (red-green-blue-magenta-cyan) for droplet sizes (left to right) a d =δ th ¼ 0:06; 0:08; 0:10 at ϕ d ¼ 1:00; 1:70 with (a) both turbulent intensities (u 0 =S b ϕ g ¼1 ð Þ ¼ 4:0; 7:5) for the lowest and highest values of droplet equivalence ratio (ϕ d ¼ 1:00; 1:70) on several c-isosurfaces across the flame front (c ¼ 0:1; 0:3; 0:5; 0:7; 0:9). The PDFs of cos θ ð Þ for 0:01<c<0:99 are also shown. It is apparent from the PDFs of the entire range of cos θ ð Þ (i.e., subplot (b)) that in most droplet cases, although there is a non-negligible probability of finding front-supported flames, back-supported flames are far more likely to be found. This is consistent with the observation that can be made from Figure 1, which shows, based on the comparison of c Figure 7. Comparison of PDFs of normalized reaction progress variable gradient Ñc j jÂ δ th across the flame brush between DNS data (solid lines) and lognormal distribution (dashed lines) atc ¼ 0:1; 0:3; 0:5; 0:7; 0:9 (red-green-blue-magenta-cyan) for droplet sizes (left to right) a d =δ th ¼ 0:06; 0:08; 0:10 at ϕ d ¼ 1:00; 1:70 with (a) and fields, that flame predominantly propagates into fuel-lean mixture and thus positive values of cos θ ð Þ principally contribute to the evaluation ofÑ c (and c ) according to Eq. (23c). The likelihood of back-supported flames increases significantly with increasing droplet size and with increasing turbulent intensity, whereas the likelihood of frontsupported flames is largely unaffected by increases in either droplet size or turbulent intensity. The only exception to this trend can be seen in cases of small droplets with high droplet equivalence ratio and under low turbulent intensity: conditions, which give rise to the greatest inhomogeneities, due to the rapid evaporation of a large number of small droplets combined with a lack of turbulent mixing. Also included in Figure 8 are the subplots that focus on the peak values cos θ ð Þ % À1:0 (i.e., subplot (a)) and cos θ ð Þ % 1:0 (i.e., subplot (c)). These subplots show that, for small droplets, it is the c-isosurfaces towards the burned gas side (e.g., c ¼ 0:7; 0:9 under initial u 0 =S b ϕ g ¼1 ð Þ ¼ 4:0 and c ¼ 0:9 under initial u 0 =S b ϕ g ¼1 ð Þ ¼ 7:5) that are most likely to produce front-supported flames for ϕ d ¼ 1:70 due to relatively large amounts of evaporated fuel vapor, whereas for medium and large droplets, it is the unburned gas side (e.g., c ¼ 0:1 isosurface), which exhibits the high probability of finding front-supported flames because of localized islands of fuel-rich zones arising from evaporation on the unburned gas side of the flame front. Conversely, the c-isosurfaces on the burned gas side (e.g., c ¼ 0:7; 0:9) are most likely to produce back-supported flames in most droplet cases due to the predominantly fuel-lean nature of the unburned gas of the flame front, except for small droplets under low initial turbulence intensity, for which the unburned gas side (e.g., c ¼ 0:1) is most likely to exhibit backsupported flames because of greater probability of finding richer mixture on burned gas side in these cases due to conducive conditions for higher evaporation rate.
Finally, Figure 9 shows joint PDFs of Ñ j j and Ñc j j, P Ñc j j; Ñ j j ð Þ , and P Ñc j j ð ÞÁP Ñ j j ð Þ across the flame brush for several droplet cases (a d =δ th ¼ 0:06, ϕ d ¼ 1:00; 1:70, u 0 =S b ϕ g ¼1 ð Þ ¼ 4:0; 7:5). Other cases display the same qualitative behavior and thus are not shown here in the interest of conciseness. No clear correlation is apparent between Ñ j j and Ñc j j, which is consistent with a previous analysis (Malkeson et al., 2013;Ruan et al., 2012). Although the absence of correlation does not necessarily mean statistical independence, the similarity between P Ñc j j; Ñ j j ð Þ and P Ñc j j ð ÞÁP Ñ j j ð Þ suggests that statistical independence of Ñ j j and Ñc j j might be a valid assumption for these flames, and that the joint PDF can be modeled as P Ñc j j ð ÞÁP Ñ j j ð Þ. It should be noted that the discrepancies in the magnitude of the joint PDFs follow naturally from the definition of joint PDFs (i.e., that the volume under the joint PDF equals unity), with the result that minor differences in spread of the joint PDFs can give rise to noticeable differences in their magnitudes.
A comparison between Figure 9a and Figures 9c-9d reveals slightly better agreement between the joint PDF P Ñc j j; Ñ j j ð Þobtained from DNS data and the lognormal bivariate distribution without correlation (i.e., Eq. (30d)) than the predictions of Eq. (30a), which accounts for the correlation between Ñ j j and Ñc j j. The better agreement is apparent only at intermediate values ofc (i.e.,c ¼ 0:5; 0:7) for which the spread of the lognormal bivariate distribution without correlation is smaller than the corresponding spread with correlation, which more closely matches the spread of the joint PDF. The degree of similarity between Figures 9c and 9d depends on the specific value of the correlation coefficient, p, in each case, such that cases where the correlation coefficient is small, Figures 9c and 9d appear similar, since under such circumstances Eq. (30a) reduces to Eq. (30d). Thus, it may be said that the bivariate lognormal distribution with correlation better approximates the joint PDF P Ñc j j; Ñ j j ð Þwhen the correlation coefficient is small. In other words, the joint PDF P Ñc j j; Ñ j j ð Þmay be better approximated by the bivariate lognormal distribution without correlation. This is consistent with the findings from Figures 9a and 9b, which indicate that P Ñc j j; Ñ j j ð Þcan be approximated by P Ñc j j ð ÞÁP Ñ j j ð Þ because the approximation P Ñc j j; Ñ j j ð Þ%P Ñc j j ð ÞÁP Ñ j j ð Þ implicitly assumes statistical independence and thus zero correlation between Ñ j jand Ñc j j. The statistical independence between Ñ j jand Ñc j jis also implicitly assumed in the derivation of Eq. (30d). As log-normal distribution reasonably models P Ñc j j ð Þ and P Ñ j j ð Þ (see Figures 4 and 5), the prediction of Eq. (30d) also behaves qualitatively similarly to that of P Ñc j j ð ÞÁP Ñ j j ð Þextracted from DNS data. Figures 4 and 5 show that log-normal distributions do not adequately capture the tails of the PDFs of P Ñc j j ð Þand P Ñ j j ð Þ, and thus there are quantitative differences between the predictions of Eq. (30d) and P Ñc j j ð ÞÁP Ñ j j ð Þ. A similar behavior was reported earlier in the context of single phase gaseous partially-premixed combustion (Malkeson et al., 2013;Ruan et al., 2012). It is worth noting that the correlation coefficient provides a measure of linear dependence between the variables in question and thus the prediction of Eq. (30a), which explicitly accounts for the correlation coefficient between Ñ j j and Ñc j j, does not necessarily capture the significantly nonlinear relation between these quantities. Thus, based on the evidence presented in Figure 9 the bi-variate log-normal distribution without correlation (i.e., Eq. (30d)) seems to be slightly better suited for modeling P Ñc j j; Ñ j j ð Þthan the bi-variate lognormal distribution with correlation (i.e., Eq. (30a)).

Conclusions
The effects of droplet size, droplet equivalence ratio and turbulence intensity on the statistical properties of the scalars and c and their gradients, Ñ j j and Ñc j j, in statistically planar flames propagating into a liquid fuel droplet mist under decaying turbulence have been analyzed using modified simple chemistry (Tarrazo et al., 2006) 3D DNS data. The combustion process in the gaseous phase has been found to take place predominantly in fuel-lean mode, even for ϕ d >1, nevertheless droplets which were able to penetrate the flame front completed their process of evaporation in the burned gas region, releasing gaseous fuel which, by means of diffusion, increased the gaseous equivalence ratio in the region of the flame front. The probability of finding fuel-lean mixture decreases with decreasing initial droplet diameter due to faster evaporation of smaller droplets. An increase in droplet equivalence ratio ϕ d also reduces the probability of having fuel-lean combustion, whereas an increase in u 0 encourages the mixing of evaporated fuel with surrounding air and produces a more flammable mixture, which plays a crucial role for small droplets at high values of ϕ d . It was shown that the joint PDF of and c cannot, in many cases, be accurately modeled in terms of discrete delta functions, and its magnitudes cannot be predicted by P ð Þ Á Pðc) due to statistical dependence between and c at some locations within the flame. It has further been shown that a β-function distribution reasonably predicts the shape and magnitude of the Favre-PDFs of and c. A presumed log-normal distribution captures both qualitative and quantitative behaviors of the PDFs of Ñ j j and Ñc j j obtained from the DNS data across the flame brush, but there exist discrepancies between the log-normal distribution and DNS data at the tails of the PDFs. The PDFs of cosine of the angle between Ñ and Ñc (i.e., cos θ ð Þ) have been found to exhibit much greater likelihood of positive values than negative values for most droplet cases. This phenomenon arises due to the predominantly fuel-lean nature of the unburned gas region of the flame. A greater probability of negative values of cos θ ð Þ has been observed in cases of small droplet than in cases of larger droplets and detailed explanations for this behavior have been provided. It has been found that the joint PDF of Ñ j j and Ñc j j (i.e., P Ñc j j; Ñ j j ð Þ ) can be approximated by P( Ñ j jÞ Á P Ñc j j ð Þ, which implicitly assumes statistical independence of Ñ j j and Ñc j j. Two variants of bivariate log-normal distribution have been considered for the purpose of modeling P Ñ j j; Ñc j j ð Þ , assuming both correlation and no correlation between Ñ j j and Ñc j j. Significant quantitative discrepancies have been found between P Ñ j j; Ñc j j ð Þobtained from DNS data and the predictions of the parameterization based on bivariate log-normal distribution, which originate principally due to inaccuracies incurred during parameterising P( Ñ j jÞ and P ( Ñc j jÞby log-normal distributions. The variant with no correlation has been found to be more successful in capturing qualitative behavior of P Ñ j j; Ñc j j ð Þthan the variant which accounts for correlation between Ñ j j and Ñc j j. It is worth noting that the present analysis has been carried out for modified single step chemistry (Tarrazo et al., 2006) DNS data with moderate values of turbulent Reynolds number and thus experimental and detailed chemistry DNS data at high values of turbulent Reynolds numbers will be necessary for more comprehensive analysis and validation of the parameterizations, which have been found to exhibit the most promising behavior when compared to DNS data. Furthermore, the aforementioned parameterizations need to be implemented in actual RANS and LES simulations for the purpose of a-posteriori assessment.