Statistical Analysis of Turbulent Flame-Droplet Interaction: A Direct Numerical Simulation Study

Turbulent combustion of mono-disperse droplet-mist has been analysed based on three-dimensional Direct Numerical Simulations (DNS) in canonical configuration under decaying turbulence for a range of different values of droplet equivalence ratio (ϕd), droplet diameter (ad) and root-mean-square value of turbulent velocity (u′). The fuel is supplied in liquid phase and the evaporation of droplets gives rise to gaseous fuel for the flame propagation into the droplet-mist. It has been found that initial droplet diameter, turbulence intensity and droplet equivalence ratio can have significant influences on the volume-integrated burning rate, flame surface area and burning rate per unit area. The droplets are found to evaporate predominantly in the preheat zone, but some droplets penetrate the flame front, reaching the burned gas side where they evaporate and some of the resulting fuel vapour diffuses back towards the flame front. The combustion process in gaseous phase takes place predominantly in fuel-lean mode even for ϕd > 1. The probability of finding fuel-lean mixture increases with increasing initial droplet diameter because of slower evaporation of larger droplets and this predominantly fuel-lean mode of combustion exhibits the attributes of low Damköhler number combustion and gives rise to thickening of flame with increasing droplet diameter. The chemical reaction is found to take place under both premixed and non-premixed modes of combustion and the relative contribution of non-premixed combustion to overall heat release increases with increasing droplet size. The statistical behaviours of the flame propagation and mode of combustion have been analysed in detail and detailed physical explanations have been provided for the observed behaviour.

of the flame propagation and mode of combustion have been analysed in detail and detailed physical explanations have been provided for the observed behaviour.
Keywords Direct numerical simulation · Turbulent combustion · Droplet diameter · Droplet equivalence ratio · Turbulence intensity

Introduction
Flame propagation into turbulent droplet-laden mixtures has a number of important engineering applications including Internal Combustion (IC) engines (e.g. Direct Injection (DI) and Compression Ignition (CI) engines) [1,2], gas turbines [2,3], and hazard prediction and control [4]. Thus, fundamental understanding of the physical processes which underpin these applications would be beneficial for the purpose of both designing highly efficient engines and gas turbines and in order to identify reliable and accurate measures of hazard predictions and control. In spite of the importance of flame-droplet interaction in the aforementioned engineering applications, this topic has received relatively limited attention. The pioneering experimental work performed by Burgoyne and Cohen [5] indicated that both the diameter and concentration of liquid fuel droplets can have significant influences on the resulting flame propagation and burning rate. The experimental analysis by Ballal and Lefebvre [6] demonstrated that the flame speed arising due to droplets small enough to evaporate completely before reaching the flame front was similar to that of a purely gaseous mixture and that the flame speed decreases with increasing initial droplet diameter. It was furthermore experimentally observed [6] that there exists an optimal initial droplet diameter for a given fuel and overall equivalence ratio, φ ov = φ g + φ d (where φ g is the contribution arising due to gaseous fuel and φ d the contribution arising due to liquid fuel), for which the laminar flame propagation was enhanced compared to that of a purely gaseous fuel. This enhanced flame speed was observed experimentally for both lean (φ ov < 1.0) and rich (φ ov > 1.0) overall equivalence ratios [7]. For an overall rich flame, it was suggested that the flame speed enhancement was due to the incomplete evaporation of large droplets on the unburned gas side which led to a local gaseous equivalence ratio close to unity (i.e. close to stoichiometry) in the region of the flame front [8][9][10]. Faeth and his co-workers [11,12] indicated that the evaporation characteristics of the droplets can also contribute to the post-evaporation flame propagation behaviour based on a number of experimental investigations. Lawes and Saat [13] showed that the inertia of the droplets may lead to enrichment of the fuel vapour, which, under quiescent or laminar flow conditions, may augment (reduce) the burning rate for fuel-lean (-rich) droplet-laden mixtures, though this effect is diminished at high values of turbulence intensity [13]. Lawes and Saat [13] furthermore showed that the equivalence ratio of the gaseous fuel-air mixture which develops due to the evaporation process can also have significant effects on both early and late stages of flame propagation, but, once again, these effects are diminished as a result of increasing turbulent velocity fluctuations. Thus, the entire body of experimental evidence suggests that flame propagation in turbulent droplet-laden mixtures is a complex physical process depending on the simultaneous interaction of evaporative heat and mass transfer, fluid dynamics and combustion thermo-chemistry and a thorough understanding of all these phenomena will be required in order to develop accurate models for the design and development of reliable, energy-efficient engines and combustors. Direct Numerical Simulations (DNS) have contributed significantly to the physical understanding and modelling of the combustion of turbulent droplet-laden mixtures in the recent past [14][15][16][17][18][19][20][21][22][23][24][25][26], where either a single-step [14][15][16][17][18][19][20][21][22][23][24] or a detailed [25,26] chemical reaction mechanism is employed. In the aforementioned DNS studies, the gaseous phase is treated in a typical Eulerian fashion and the droplets are considered as sub-grid particles which are tracked in a Lagrangian manner. In this numerical framework, appropriate source terms involving suitable relaxation time-scales have been used in the gaseous phase for the mass, momentum and energy conservations in order to account for the contributions arising due to droplet evaporation. The same treatment of droplets as sub-grid particles will also be adopted here to analyse the statistical behaviour of flame-droplet interaction. Although the atomisation of the liquid fuel in realistic systems gives rise to poly-disperse sprays, fundamental understanding of the factors controlling flame propagation can be obtained by studying mists of monodisperse droplets. The effects of volatility, droplet diameter and droplet equivalence ratio on burning velocity in one-dimensional flames where fuel is supplied in the form of monodisperse droplets have recently been analysed by Neophytou and Mastorakos [27] based on detailed chemistry n-heptane and n-decane laminar flame simulations. The effects of equivalence ratio, droplet diameter and droplet group number on n-decane spray flame structure have also been previously investigated numerically by two-dimensional unsteady laminar counter-flow [28][29][30] and turbulent jet [23] flame simulations. The analysis of Neophytou and Mastorakos [27] indicated that the high values of burning velocity and chemical reaction rate are obtained for 'small' droplets with overall equivalence ratio φ ov either less than or equal to 1.0 (i.e. φ ov ≤ 1.0) and for 'large' droplets with overall equivalence ratio greater than one (i.e. φ ov > 1.0). Furthermore, it has been found by Neophytou and Mastorakos [27] that combustion can be sustained even for φ ov values which are greater than the rich-flammability limit for gaseous fuel-air mixture. This behaviour is due to the incomplete evaporation of the droplets which leads to a gaseous equivalence ratio much lower than the overall equivalence ratio. Furthermore, numerical results based on two-dimensional unsteady laminar counter-flow [28][29][30] and turbulent jet [23] simulations suggest that the resulting spray flame structure shows a considerable extent of premixed combustion in addition to an expected diffusion mode of burning [23,[28][29][30]. The extent of premixed mode of combustion has been found to increase with decreasing droplet diameter. Moreover, the burned gas temperature can be different from a gaseous diffusion flame due to the combined effects of latent heat of evaporation and premixed combustion [28][29][30]. It is further demonstrated by Nakamura et al. [28] that the group combustion number which is defined based on a homogenous field of unburned droplets and air mixture is not sufficient to characterize the combustion process in spray flames due to the premixed-like combustion. The simulation results by Fujita et al. [23] and Watanabe et al. [30] demonstrated that radiative heat transfer can play an important role in spray combustion and flame propagation under some conditions. The present analysis extends the analysis of Neophytou and Mastorakos [27] for turbulent flames by carrying out three-dimensional compressible DNS of turbulent flame propagation into droplet-laden mixtures. The introduction of a turbulent flow field enhances mixing and influences the local evaporation rate, which in turn affects the local equivalence ratio in such a manner that extremely fuel-lean and fuel-rich regions are less likely to occur. This, in turn, is likely to have significant influences on heat release and reaction rate in turbulent spray flames. It is yet to be assessed if the aforementioned observations of Neophytou and Mastorakos [27] based on laminar flame simulations also remain qualitatively valid under turbulent flow conditions. Furthermore, the influences of root-mean-square (rms) velocity fluctuation on flame structure and its propagation in droplet mists are yet to be analysed in detail. The aforementioned gaps in the existing literature have been addressed by this analysis, in which the DNS data has been used to analyse the effects of initial droplet diameter, a d , droplet equivalence ratio, φ d , and root-mean-square (rms) value of turbulent velocity, u , on the flame structure and the statistical behaviour of the flame propagation in droplet-laden mixtures. The main objectives of the present study are as follows: 1. To demonstrate the effects of a d , φ d and u on the flame structure when fuel is supplied in the form of droplets on the unburned gas side of the flame. 2. To understand the physical processes which significantly influence flame propagation and burning rate in turbulent flame-droplet interaction.
The remainder of the paper takes the following form. The necessary mathematical background related to the current study is provided in the next section. This will be followed by a brief discussion of the numerical implementation. Following this, the results are presented and subsequently discussed. Finally, the main findings are summarised and conclusions are drawn.

Mathematical Background
DNS of turbulent reacting flows should ideally account for both the inherent threedimensionality of turbulence and the detailed chemical structure of the flame. However, in general, the limitations of computer hardware restrict most DNS studies to either twodimensional domains with detailed chemistry or three-dimensional domains with simplified chemistry. Although three-dimensional DNS with detailed chemistry is now possible, it remains prohibitively expensive (e.g. of the order of millions of CPU hours [31]) to carry out an extensive parametric analysis which has been performed in the current investigation. Thus, for the purpose of carrying out the present extensive parametric analysis at a modest computational cost, a single-step irreversible Arrhenius-type chemical mechanism has been used: Fuel where s is the oxidiser-fuel ratio by mass (i.e. the mass of oxygen consumed per unit mass of fuel). The fuel reaction rateω F is given by the Arrhenius-type expression: where ρ is the gas density and Y F and Y O are the fuel and oxygen mass fractions respectively. T, the non-dimensional temperature, β, the Zel'dovich number, α, a heat release parameter, and B * , the normalised pre-exponential factor, are defined as: 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 τ = T ad(φ g =1) − T 0 /T 0 is a heat release parameter. Here, the pre-exponential factor B is taken to be B ≈ 2.86 × 10 6 S 2 b(φ g =1) /α T0 (where α T0 is the thermal diffusivity of the unburned gas) for the present thermo-chemistry, and the activation energy, E ac , and the heat of combustion are taken to be functions of the gaseous equivalence ratio, φ g , following the suggestion of Tarrazo et al. [32], which correctly predicts the equivalence ratio φ g dependence of the unstrained laminar burning velocity S b(φ g) in hydrocarbon-air flames, especially for fuel-rich mixtures. Accordingly, the Zel'dovich number, β, is expressed as β = 6f φ g where: Furthermore, the heat release per unit mass of fuel H φ g = T ad(φ g) −T 0 C P Y F0(φ g) is [32], where α H = 0.18 and Y F0(φ g) is the fuel mass fraction in the unburned gas for a premixed flame of equivalence ratio φ g . Interested readers are referred to Refs. [32] (see Fig. 2 of Ref. [32]), where the chemical mechanism is shown to compare favourably with both experiments and detailed chemistry simulations for typical hydrocarbon-air mixtures. It was also shown in Ref. [33] that this chemical mechanism captures the realistic variation of laminar burning velocity with gaseous equivalence ratio φ g in the context of the thermo-chemistry used in this analysis. The Lewis numbers of all species are taken to be equal to unity and all species in 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 gaseous specific heats at constant pressure and volume respectively) and Prandtl number (Pr =μC g P /λ= 0.7, where μ and λ are the dynamic viscosity and thermal conductivity of the gaseous phase respectively).
The droplet transport equations employed by Reveillon and Vervisch [15] have been used for the purpose of this analysis. Amongst the quantities transported for each droplet are the position, x d , velocity, u d , diameter, a d and temperature, T d . The 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/decay timescales for droplet velocity, diameter and temperature respectively, which are defined as: where ρ d is the droplet density, C L p is the specific heat for the liquid phase, C u is the corrected drag coefficient and is given by: In Eqs. 9-10, Re d is the droplet Reynolds number, Sc is 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 defined as [15]: where Y s F is the value of the fuel mass fraction Y F at the surface of the droplet. Equations 11-13 also invoke the unity Lewis number assumption. The Clausius-Clapeyron equation for the partial pressure of the fuel vapour at the droplet surface, p s F , is used to evaluate the Spalding number B d , which leads to the following expressions: where T s ref is the boiling point of the fuel at pressure p ref , R is the gas constant, T s d is assumed to be T d , and W O and W F are the molecular weights of oxidiser and fuel respectively.
In this analysis, the drag coefficient is corrected (see Eq. 10) for large values of the droplet Reynolds number 1 ). This correction in was used in several previous analyses [15][16][17][18][19][20][21][22][23][24]34]. The droplet Reynolds number Re d remains either of the order of unity or smaller than unity (i.e. Re d ≤ 1) in most (∼95 %) locations (especially in the vicinity of the flame) and thus the exact nature of the corrections to the drag coefficient is not expected to have a significant influence on the simulation results. For this reason a simple drag coefficient correction [34] has been considered in this analysis. A more involved empirical expression for C u accounting for blowing velocity at the droplet surface was used in some analyses [14,[28][29][30] in the past, but the exact nature of empiricism is likely to have higher order effects for small values of Re d when DNS is carried out for the carrier phase and the droplets are considered as sub-grid point sources.
Droplet evaporation leads to mixture inhomogeneities, which are often characterised in terms of the mixture fraction, which is defined as: where Y F∞ = 1.0 is the fuel mass fraction in the pure fuel stream and Y O∞ = 0.233 is the oxygen mass fraction in air. The hydrocarbon fuel used for 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. One can furthermore define a reaction progress variable, c, based on a species mass fraction and mixture fraction such that c rises monotonically from 0.0 in unburned reactants to 1.0 in fully burned products. In droplet combustion it is advantageous to employ an oxygen-based reaction progress variable, which takes the following form [17,19]: From the definition of the reaction progress variable (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 ξ : where the first term on the right-hand-side arises due to molecular diffusion, the second represents reaction rate, the third is the source/sink term arising due to droplet evaporation, which originates due to a non-zero source term in the mass conservation equation and the last is the cross-scalar dissipation term arising due to reactant inhomogeneity [33,35]. The cross-scalar dissipation termȦ c in Eq. 18 arises due to mixture inhomogeneity [33,35], which in this case is induced by droplet evaporation. The evaporation of droplets increases the local fuel mass fraction at the evaporation sites, which in turn affects the local oxidiser mass fraction and, hence, the local mixture fraction (since ξ = ξ(Y F , Y O )) and, finally, the reaction progress variable (since c = c(Y O , ξ)). Due to the definition of c (see Eq. 17), the definitions ofω c ,Ṡ c andȦ c depend on the value of ξ . The reaction rateω c of the reaction progress variable is given as:ω The expressions forṠ c andȦ c are as follows: ,ξ>ξ st (20) 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 molecular diffusion term in Eq. 18 may be further broken down into its normal and tangential components to give: where N = −∇c/| ∇c| is the flame normal vector, κ m = 0.5(∇· N) is the arithmetic mean of the two principal curvatures of a given isosurface c = c * . In Eq. 22 the first term on the right-hand-side gives the component of molecular diffusion normal to the flame front and the second expression gives the tangential molecular diffusion component. The statistical behaviour of the termsω c , N · ∇ ρD N · ∇c , (−2ρDκ m |∇c|),Ṡ c andȦ c are discussed in detail in Section 4 of this paper.

Numerical Implementation
A three-dimensional compressible DNS code SENGA [19,21] is used for the present numerical investigation. A 10th order central difference scheme is used for spatial discretisation for the internal grid points and the order of differentiation drops gradually to a one-sided 2nd order scheme at the non-periodic boundaries [19,21]. A low-storage third-order explicit Runge-Kutta scheme [36] is used for the explicit time advancement.
Here a rectangular computational domain of size 30δ th ×20δ th ×20δ th has been considered, 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 discretised using a Cartesian grid of size 384 × 256 × 256 which ensures about 10 grid points within δ th . This grid is also fine enough to resolve the Kolmogorov length-scale η. Partially non-reflecting boundary conditions are used for the mean direction of flame propagation (i.e. x-direction), while the y-and z-directions are considered to be periodic. The non-periodic boundary conditions are specified using the Navier Stokes Characteristic Boundary Conditions (NSCBC) technique [37]. The droplets are initially distributed uniformly in space throughout the y-and z-directions and in the region 0.0 ≤ x/δ th ≤ 10.0. The reacting flow field is initialised using the steady laminar solution for the required initial values of droplet diameter a d and droplet equivalence ratio φ d . The steady laminar solution is generated using a commercial software COSILAB [38], where the one-dimensional forms of the 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 mono-disperse droplets on the unburned gas side of the flame. Interested readers are referred to Ref. [27] for the mathematical framework for generating steady state one-dimensional solutions of the laminar flame-droplet interaction for mono-disperse sprays. An incompressible homogeneous isotropic velocity field generated using a standard pseudo-spectral method [39] is used for initialising turbulent velocity fluctuations. This initial velocity field is superimposed on the steady laminar spray flame solution generated using COSILAB [38]. For the present analysis the unburned gas temperature is taken to be T 0 = 300K, which yields a heat release parameter τ = T ad(φ g =1) −T 0 /T 0 = 6.54. In the full turbulent three-dimensional DNS analysis, as in the one-dimensional laminar case, the fuel is supplied purely in the form of mono-disperse droplets with non-dimensional 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 10δ th from the point in the laminar flame at whichT= 400K, which corresponds to a nondimensional temperature T ≈ 0.05. The values of the laminar burning velocity for each droplet case, S b,d are listed in Table 1 after normalising them by the unstrained laminar burning velocity of the stoichiometric mixture S b(φ g =1) . The number of droplets at t = 0 varies between 1.16 ≤ (ρ N ) 1/3 δ th ≤ 2.27, where ρ N is the droplet number density in the region 0.0 ≤ x/δ th ≤ 10.0.. The liquid volume fraction remains much less than 0.01 for all cases considered here. The maximum droplet diameter experimentally analysed by Burgoyne and Cohen [5] did not exceed approximately 10 % of the unstrained laminar thermal flame thickness of the stoichiometric mixture. The ratio of the droplet to computational cell Table 1 Ratio of droplet laminar burning velocity to the unstrained laminar burning velocity of the stoichiometric mixture volume remains about 0.11, 0.27 and 0.5 for a d /δ th ≈ 0.06, 0.08, 0.1 respectively in the DNS simulations presented here Hence, the sub-grid point source treatment of droplets is expected to be valid for the flame-droplet interaction analysed here. In all cases droplets are supplied at the left hand side boundary to maintain φ d ahead of the flame. The droplets evaporate as they approach the flame front. Due to the high volatility of n-heptane, evaporation commences on entry. The droplets which interact with the flame are in reality much smaller than the initial size of the droplets which have been injected at a distance 10δ th from the location corresponding to T ≈ 0.05 The droplet diameter decreases by at least 25 % when it reaches the most reactive region of the flame. As the droplets which interact with the flame front are in reality much smaller than the initial size of the droplets, the modelling of sub-grid evaporation is not expected to significantly affect the statistics of flame-droplet interaction, which is the main focus of the analysis. For this reason, it can be expected the qualitative nature of the findings are not sensitive to the evaporation modelling used in this analysis.
The simulations are carried out under both laminar flow conditions and for normalised root-mean-square (rms) turbulent velocities u /S b(φ g =1) = 4.0 and 7.5 with a non-dimensional longitudinal integral length-scale L 11 /δ th = 2.5. Thus, the domain size of this analysis is given by 12L 11 × 8L 11 × 8L 11 i.e. the span-wise direction contains 8 integral length-scales. This ensures that enough number of integral eddies are retained when the statistics are extracted. The number of integral eddies within the domain for this analysis is either comparable to or greater than that used in several previous DNS based analyses [24][25][26][27][40][41][42][43]. It would indeed be desirable to have a larger number of eddies within the domain than that used here, but that would require a larger domain size which would greatly increase the computational cost for the parametric analysis carried out in this paper. The ratio of droplet diameters to the Kolmogorov scale is a d /η ≈ 0.3, 0.4, 0.5 in the cases of a d /δ th ≈ 0.06, 0.08, 0.1 respectively for the highest u /S b(φ g =1) case (i.e. the initial value of u /S b(φ g =1) = 7.5) considered in this paper. The ratios of droplet diameter to the Kolmogorov length scale and droplet volume to cell volume remain comparable to several previous analyses [15,16,18,19,21]. The mean normalized inter-droplet distance s d /η ranges between 0.0220 and 0.0432 (i.e. 0.0220 <s d /η< 0.0432) for the highest u /S b(φ g =1) case. All simulations are allowed to evolve until t final ≥ 3t turb ≈ 4t chem , where t turb =L 11 /u is the initial turbulent eddy turnover time and t chem = D 0 /S 2 b(φ g =1) the chemical timescale, where D 0 is the unburned gas diffusivity. The simulation time used in the current analysis remains comparable to a number of recent DNS analyses [19,21,[24][25][26][40][41][42][43], which significantly contributed to the fundamental understanding of turbulent combustion.
The Reynolds/Favre averaged value of a general quantity Q (i.e.Q andQ) are evaluated by ensemble averaging the quantity Q over the y − z plane at a given x location. The statistical convergence of the Reynolds/Favre averaged values has been ensured by comparing the values obtained on full sample size with the corresponding values based on distinct half of the available sample size. For the sake of conciseness, the values based on full sample size will be reported in the next section of this paper.

Global features of flame-droplet interaction
The instantaneous distributions of non-dimensional temperature, T, fuel mass fraction, Y F , magnitude of normalised fuel reaction rate, |ω F | × δ th /S b φ g = 1 , and reaction progress variable, c, at the central x − y mid-plane at t ≈ 4t chem are shown in Fig. 1 for each droplet size with droplet equivalence ratio φ d = 1.0 and also, for the sake of comparison, for a premixed stoichiometric flame. The droplets, depicted as black dots, are those residing in the cells immediately above or below the plane shown in the figure. It is evident from Fig. 1 (top row) that the droplets shrink due to evaporation as they approach the flame, although complete evaporation may not take place before a droplet passes through the flame front, as can be seen from those droplets which have penetrated the flame front. It is furthermore evident from Fig. 1 (top row) that the burned gas temperature is noticeably lower in the droplet cases in comparison to the stoichiometric premixed flame. This can be attributed to two factors: firstly, the evaporating droplets absorb latent heat from the background gas (this occurs on both sides of the flame, but is more noticeable on the burned gas side), secondly, since, in many cases, the evaporation of droplets is not complete on their arrival at the flame front, the reaction takes place predominantly under fuel-lean conditions and, thus, the heat release due to combustion and the resultant burned gas temperature are smaller than in the corresponding stoichiometric premixed flame. The reduction of burned gas temperature due to absorption of the latent heat of evaporation is consistent with previous findings based on two-dimensional unsteady laminar [28][29][30] and turbulent [23] flame simulations. The extent by which the gaseous fuel mass fraction field in droplet cases is reduced in comparison to the premixed stoichiometric case can be seen from Fig. 1 (second row). The premixed case benefits from a uniformly stoichiometric mixture, whereas the gaseous fuel field in the droplet cases is both fuel-lean and non-uniform. The degree of non-uniformity of fuel distribution (i.e. isolated islands of high fuel concentration) on the unburned gas side of the flame appears to increase with droplet size. This is due to the gaseous fuel clouds surrounding individual droplets coalescing to form a single large cloud in the case of small droplets. This does not occur for larger droplets, which are fewer in number (for constant φ d ) and also evaporate more slowly. The droplets continue to evaporate in the burned gas region and some of the evaporated gaseous fuel diffuses back towards the flame front, thereby augmenting the fuel mass fraction in the region of the flame front. As mentioned above, the fuel reaction rate magnitude is greatly reduced for the droplet cases in comparison to the gaseous premixed flame. This can be seen by comparing the magnitude of the reaction rate fields ( Fig. 1 (third row)) for the premixed and droplet cases. High reaction rates are absent along most of the flame front in the droplet cases. The notable exceptions can be seen in the vicinity of regions where the gaseous fuel mass fraction is close to stoichiometry (as indicated in the third row of Fig. 1). These regions arise both due to the presence of a droplet in the immediate vicinity of the reaction front, in which case the locally high temperature leads to rapid evaporation and augmentation of the gaseous fuel mass fraction, and due to post-flame evaporation, in which case the local fuel mass fraction is augmented by the back-diffusion of the newly-evaporated gaseous fuel from the burned gas region. Droplets residing at the flame front can be seen, for example, in the lower half of the medium-sized droplet case (a d /δ th = 0.08), in which case the ξ = ξ st isosurface is restricted to the immediate proximity of the reaction front. However, in the upper half of the same plot the ξ = ξ st isosurface can be seen to extend into the burned gas region because there exists a source of fuel beyond the flame front from which the fuel is back-diffusing (i.e. a droplet undergoing post-flame evaporation). A similar behaviour of the ξ = ξ st isosurface extending into the burned gas region can be seen in the large droplet case (a d /δ th = 0.10). It is noteworthy that both these cases have occurred in regions where the nearby flame front is convex to the reactants. This flame geometry does not expose the droplets approaching the flame to the focusing of heat as in a case where the flame front is concave to the reactants, thus increasing the probability of finding droplets beyond the flame front. On the other hand, the region of burned gas into which these droplets have entered is visibly hotter than neighboring regions and this contributes to a more rapid evaporation. It is worth noting that the penetration of droplets through the flame and the subsequent diffusion of fuel vapour back towards the flame due to the post-flame evaporation is consistent with previous one-dimensional detailed chemistry laminar flame simulations by Neophytou and Mastorakos [27].
A comparison of the fuel mass fraction, Y F , and reaction progress variable, c, (Fig. 1  (bottom row)) fields indicates that, almost without exception, the inhomogeneous mixture arising from evaporation remains fuel-lean (i.e. Y F <Y Fst = 0.0621) on the unburned gas side (see Fig. 5). Finally, it can be seen from Fig. 1 (bottom row) that, whereas in the premixed stoichiometric case, the c-isosurfaces lie in close proximity to one another, in the droplet cases the distance between them appears to increase with increasing droplet size, indicating a wider flame front and reaction zone. The reasons for this distinction and its significance are discussed later in this paper (see Fig. 7).
The global burning rate can be characterised in terms of the turbulent flame speed, S T , which is defined as S T =(1/ρ 0 A p ) ω c dV, whereω c is the reaction rate of reaction progress variable (see Eq. 19), A p is the projected flame area in the direction of mean flame propagation and dV is an infinitesimal volume element. Figure 2 shows the temporal evolution of the turbulent flame speed normalised by its initial value, S T,0 , which is identical to the laminar burning velocity of the corresponding one-dimensional laminar spray flame case. The effects of turbulence intensity, droplet size and droplet equivalence ratio on S T /S T,0 are clearly visible from Fig. 2. Evaporation of droplets absorbs latent heat from the surrounding gas which acts to reduce the temperature of the gaseous phase and, thereby, reduces the rate of fuel consumption. This behaviour is particularly prominent for high values of φ d . The smaller droplets evaporate more quickly, leading to a sharper drop in both temperature and fuel consumption rate for high values of φ d . This drop may be compounded further by a build-up of gaseous fuel in the vicinity of the droplets, resulting in a locally non-flammable fuel-rich mixture. This situation may be partially or even wholly alleviated either by the effects of turbulent mixing or (at a much slower rate) by molecular diffusion. These features are apparent most clearly from the evolution of the turbulent flame speed in the case of small droplets (i.e. a d /δ th = 0.06), in which the behaviour of S T can be seen to depend most strongly on the value of φ d and on the turbulent flow conditions. The timescale required for the burning rate in the φ d = 1.70 case to recover from the initial dip decreases with increasing turbulence intensity. The same can be said to be true in the case of larger droplets, but the effects, such as the initial dip when φ d = 1.70, are not as prominent due to the slower evaporation of the droplets. However, it can be seen from Fig. 2 that combustion can be sustained for φ d = 1.70 in the case of flame propagation in droplet mist even though a gaseous equivalence ratio of 1.70 is beyond the rich flammability limit, which is consistent with previous laminar flame simulations [27]. By the end of the simulation time (t ≈ 4t chem ), in most cases the turbulent flame speed has begun to level off, indicating that droplet evaporation and gaseous fuel consumption approximately balance each other. However, it is clear that throughout the simulation time the burning rate of the droplet cases does not reach that of the premixed cases, and that, furthermore, the difference between them increases with increasing turbulence intensity.
The quantity A = ∫ |∇c|dV provides a measure of the flame area, and, when normalised by its initial value, A 0 , it provides a measure of the global wrinkling of the flame surface. The temporal evolution of normalised flame area A/A 0 is shown here in Fig. 3. All droplet cases show an initial increase in A/A 0 with time, before levelling off at later times. Figure 3 reveals that A/A 0 (i.e. total amount of global wrinkling) for all sizes of droplets is principally affected by u /S b(φ g =1) , rather than φ d , and that A/A 0 increases primarily with increasing turbulence intensity u /S b(φ g =1) . Once again the premixed stoichiometric flame attains higher values of normalised flame area than the comparable droplet cases and the difference between them increases with increasing turbulence intensity.
The quantity = ( ω c dV/ |∇c|dV) can be used to characterise the amount of burning per unit area of flame. The temporal evolution of /ρ 0 S b(φ g =1) is shown in Fig. 4 and indicates that this quantity follows the trend of S T /S T,0 and A/A 0 , eventually approaching a quasi-steady state for most cases considered here since by this time the decay of turbulence is roughly compensated by turbulence generation due to heat release and thus the turbulent kinetic energy evaluated over the domain does not vary rapidly with time. As a result of this, the flame wrinkling also does not change rapidly with time, and the generation and  Figure 4 also shows that, after a short transient period, the ratio /ρ 0 S b(φ g =1) increases with increasing φ d for both medium and large droplet cases. The same trend develops at later times for small droplets under high turbulence intensity (i.e. u /S b(φ g =1) =7.5), after a considerably longer transience due to the effects of more rapid evaporation which, as was previously noted, inhibit the burning rate for cases of large values of φ d . The reason for this trend is that an increase in φ d provides larger amount of fuel in the gaseous phase, producing a locally non-flammable mixture. Turbulent mixing eventually produces a flammable mixture. This process takes place more rapidly under high turbulence intensity than under low turbulence intensity. For this reason, the case of small droplets under low turbulence intensity has not yet developed the aforementioned trend, although Fig. 4 indicates that it would eventually develop. It is evident from Fig. 4 that /ρ 0 S b(φ g =1) remains smaller than unity for all droplet cases (including initial values obtained from steady one-dimensional laminar spray flame), indicating that the burning rate per unit area for the flames considered here is smaller than the corresponding value in the stoichiometric laminar premixed flame. It is also shown in Fig. 4 that /ρ 0 S b(φ g =1) remains close to unity for turbulent stoichiometric premixed flames.

Gaseous phase equivalence ratio and its implications for flames propagating into droplet-laden mixtures
It is useful to understand the composition of the fuel-air mixture in the gaseous phase in order to explain the reduced burning rate of the spray flames where fuel is supplied in the form of mono-disperse droplets in comparison to that of the stoichiometric laminar premixed flame. A comparison between the fuel mass fraction Y F and reaction progress variable c fields indicates that, almost without exception, the inhomogeneous mixture arising from evaporation remains fuel-lean (i.e. Y F < Y Fst ) on the unburned gas side. The predominance of fuel-lean conditions observed in Fig. 1 can further be substantiated from By comparing the PDFs of the two regions one can identify the source of certain features of the PDFs as arising due to events in either the burned or unburned gases. Small droplets under laminar flow conditions (Fig. 5(a1, a2)), which are able to evaporate before reaching the flame front, show an acute dependence on the droplet equivalence ratio: too few droplets (φ d = 1.00) lead to an almost entirely fuel-lean distribution and too many droplets (φ d = 1.70) lead to an almost entirely fuel-rich distribution. It is only for mid-values of droplet equivalence ratio (φ d = 1.25, 1.50) that the droplets are plentiful enough to encourage the development of a significant likelihood of encountering a stoichiometric mixture in laminar flow conditions, but not so plentiful that the mixture is fuel-rich. In these cases a stoichiometric mixture develops in regions where c > 0.90, as can be seen by comparing Fig. 5(a1) and (a2) for small droplets. The effects of chemical reaction and heat release are strong for φ g ≈ 1.0 mixtures and, thus, the evaporation rate is expected to be high close to these burned gas regions. This eventually leads to the aforementioned peak value of the PDF at φ g ≈ 1.0 in Fig. 5(a1). The same phenomena can be observed in the presence of turbulent flow ( Fig. 5(b1, b2) and (c1, c2)), although the extreme leanness or richness of the mixture is ameliorated by the turbulent mixing. Medium and large droplets exhibit somewhat similar PDFs to small droplets, but the slower evaporation of these droplets precludes the development of locally fuel-rich mixture to the extent found in the case of small droplets. This implies that all cases of medium and large droplets exhibit a high probability of finding fuel-lean mixture. Once again, when the extended range (0.10 ≤ c ≤ 0.99) is taken into account, a high probability of finding stoichiometric mixture (i.e. ξ = ξ st ) develops, which, although always non-negligible, is sometimes more likely to occur than the fuel-lean mixture and sometimes less likely, depending on the specific choice of droplet size, equivalence ratio and turbulence intensity. In general, turbulent mixing has the effect of spreading the distribution more uniformly, thereby reducing both the high likelihood of finding fuellean mixture and, where appropriate, stoichiometric or fuel-rich mixture. Finally, the PDFs of φ g in the progress variable range 0.01 ≤ c ≤ 0.99 at time t = 0 (not shown) indicate an even greater probability (up to an order of magnitude larger) for fuel-lean mode of combustion. The predominance of fuel-lean mode of combustion reduces the burning rate per unit area for the flames considered here in comparison to the corresponding value in the stoichiometric premixed laminar flame in spite of having a peak of PDF(φ g ) at φ g ≈ 1.0. The difference between overall equivalence ratio φ ov and the gaseous equivalence ratio φ g in spray flames and its dependence on droplet size are found to be consistent with previous detailed chemistry laminar flame simulations [27]. Figure 1 (bottom row) shows that the initial droplet diameter a d has important influences on the distribution of the reaction progress variable, c. It can be seen from Fig. 1 (bottom row) that c-isosurfaces representing the preheat zone (e.g. c < 0.5) are more wrinkled than the isosurfaces in the reaction zone (e.g. 0.7 < c < 0.9) for the turbulent droplet cases considered here. This behaviour is typical of the thickened flame regime of combustion (e.g. thin reaction zones regime), where turbulent eddies can penetrate into the flame, but decay before entering into the reaction zone, which allows for the continuation of chemical reactions without extinction. The regime of combustion can be characterised with the help of Karlovitz number Ka = (u /S b(φ g =1) ) 3/2 (L 11 S b(φ g =1) /α T ) −1/2 (where α T0 is the thermal diffusivity of the unburned gas), which provides a measure of the ratio of flame thickness to the Kolmogorov length scale [44]. For a stoichiometric mixture one obtains a Karlovitz number of 9.0 for the values of u /S b(φ g =1) = 7.5 and L 11 /δ th = 2.5 (and Ka = 3.5 for the values of u /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 takes place well within the thin reaction zones regime [44]. Furthermore, the Damköhler number Da = t turb /t chem = L 11 S 2 b(φ g ) /u α T0 scales as Da ∼Re 1/2 t /Ka [44] and, hence, low Damköhler number combustion is likely for the high Karlovitz number (i.e. Ka 1) cases considered here. This can further be substantiated from the variation of the segregation factor g = c 2 /c (1−c) withc shown in Fig. 6 for the turbulent cases considered here, wherec=ρc/ρ is the Favre-averaged reaction progress variable and c = c −c is the Favre fluctuation of c. According to Bray et al. [45] (1−c)) where PDF(c) shows a bi-modal distribution with delta functions at c = 0.0 and 1.0. Since the chemical reaction in the gaseous phase takes place predominantly in a fuel-lean mode, Da is expected to be small (due to small values of S b(φ g ) for fuel-lean mixtures), which leads to a reduced value of g (i.e. g < 1.0) in turbulent droplet cases. The variations of g withc for turbulent premixed stoichiometric flames are shown in Fig. 6, which shows that g remains smaller than 1.0 even for these flames as combustion takes place under small values of Da (i.e. Da < 1). It is evident from Fig. 6 that, as has already been observed with regard to the c-isosurface contours in Fig. 1, in most cases g assumes a lower value in the preheat zone (e.g. c < 0.5) than in the reaction zone (e.g. 0.7 < c < 0.9), indicating that the effects of reduced Damköhler number are more dominant in the preheat zone of these turbulent flames. It is also evident that the low Damköhler number effects (i.e. effects of slow chemical reaction) are strong for all cases where PDF(φ g ≈ 1.0) assumes a low value (see Fig. 5). This is due to the fact that the highest reaction rates are experienced in regions where the mixture is close to stoichiometry. Thus, in low u cases small droplets, which evaporate more easily, give rise to favourable fuel-air mixtures, attain higher g-values and result in a thinner flame, whereas larger droplets, which evaporate more slowly, give rise to highly fuel-lean conditions under the same flow conditions and, hence, much lower reaction rates and thicker flames. Increasing the turbulence intensity has the effect of reducing both the difference arising due to droplet size and that due to droplet equivalence ratio, as can be seen by comparing the left and right columns of Fig. 6. When these results are compared with those of turbulent stoichiometric premixed flames, one can see that, in general, the droplet cases assume smaller g-values than the premixed cases, with the exception of some of the a d /δ th = 0.06 cases (the ones which enjoy the most favourable conditions in terms of the availability of the gaseous fuel among all droplet cases), showing g-values which are comparable to the corresponding stoichiometric premixed case towards the burned gas side of the flame brush. Although the magnitude of the segregation factor g is much lower for most droplet cases, the profile of the variation of g withc is similar, especially for high turbulence intensity cases (e.g. u /S b (φ g = 1) = 7.5) and especially for high values ofc (i.e. towards the burned gas side of the flame brush). The difference in the behaviour of c 2 and g at low and high values ofc arises due to the fuel-lean mixture which is prevalent at low c values (see Fig. 5), and the mixtures with φ g ≈ 1.0 which are more likely to be found at high c values (see Fig. 5). This allows the droplet cases to share some qualitative similarities with the stoichiometric premixed flame in that region.
It can be seen from Fig. 1 (bottom row) that for the small droplets most c-isosurfaces remain close to each other, whereas, for the medium and large droplets the isosurfaces are increasingly separated from each other, indicating a broadening of the flame in these cases. The thickening of the reaction zone can be confirmed from Fig. 7 which shows the variation of the normalised mean values of the surface density function (SDF) (i.e. |∇c|) conditional on c. The inverse of the normalised SDF, |∇c| ×δ th , can be understood as a measure of the local flame thickness relative to that of the unstrained laminar flame thickness of a stoichiometric premixed flame (i.e. |∇c| ∼ 1/δ, where δ is the local flame thickness). Figure 7 shows that the peak value of the conditional mean value of normalised SDF (i.e. |∇c| ×δ th ) remains smaller than unity throughout the flame for all droplet cases, indicating a thicker flame in comparison to the stoichiometric premixed laminar flame in spite of a peak value of gaseous equivalence ratio PDF at φ g ≈1.0. As combustion takes place predominantly under fuel-lean mode (see Fig. 5), the flame thickness δ ∼α T0 /S b(φ g ) for droplet cases is expected to be greater than that in the stoichiometric premixed flame. As a result, the maximum value of the mean normalised SDF |∇c| ×δ th conditional on c remains smaller than unity for all droplet cases considered here. The results for the stoichiometric premixed flames are also shown in Fig. 7, which shows that the maximum value of the mean value of normalised SDF conditional on c indeed attains a value close to unity (i.e. |∇c| ×δ th ≈ 1.0). It is also noteworthy that the peak value of the conditional mean value of normalised SDF in the stoichiometric premixed cases are obtained at lower values of c than in the droplet cases, due to the readily available gaseous fuel in the stoichiometric premixed flames as opposed to the droplet cases which depend on droplet evaporation to supply gaseous fuel. The availability of gaseous fuel increases towards the burned gas side of the flame brush for the droplet cases (see Fig. 5), which acts to reduce the flame thickness δ ∼α T0 /S b(φ g ) towards the burned gas side. This is reflected in the shift of the peak value of the conditional mean value of normalised SDF towards a higher value of c in the droplet cases than in the corresponding stoichiometric premixed flames. The availability of more reactive fuel-air mixture and thinner flame front in the burned gas side than in the unburned gas side, along with the decay of turbulence intensity within the flame, are expected to give rise to an increase in Da∼L 11 S b(φ g ) /δu from the unburned to the burned gas side of the flame brush. This is reflected in the increase in g with increasingc in Fig. 6 for the turbulent droplet cases. Figure 7 further shows that, under each flow condition considered here, the highest peak value of the mean normalised SDF (i.e. thinnest flame) is attained by one of the flames involving the smallest droplets (i.e. a d δ th = 0.06). Although, in the absence of turbulent mixing, it is the lowest value of φ d that attains the peak value of the mean normalised SDF, whereas, under sufficiently turbulent flow conditions, it is the highest value of φ d that does so. Furthermore, for most cases of medium and large droplets, the normalised SDF increases with increasing φ d . This shows that flames possessing lower droplet equivalence ratios tend to have a lower reaction rate and to be thicker. The same is true for small droplets under sufficiently turbulent flow conditions, although their rapid evaporation rate leads, under laminar flow conditions, to a very low reaction rate and a very thick flame (low SDF value) for the highest value of φ d considered here (i.e. φ d = 1.70). The maximum value of the mean normalised SDF in the absence of turbulent mixing is attained by the intermediate values of φ d (i.e. φ d = 1.25, 1.50). This is due to the incomplete evaporation of the droplets, which means that for the φ d = 1.00 case the mixture remains mostly fuel-lean (since not all of the liquid fuel has evaporated), whereas the medium cases of φ d are more likely to produce stoichiometric mixture for a system of partially-evaporated droplets. Under these circumstances, the reaction rate for the medium φ d cases will exceed that of the low φ d case.
In conclusion, it may be said that the role of turbulent mixing has a lesser effect on the flame thickness of medium and larger droplets (although some change in thickness is noticeable in Fig. 7) than for small droplets, even though it is apparent from Fig. 6 that turbulent mixing is essential for medium and large droplets to attain high values of g and c 2 in droplet-laden flames.
The effects of φ d , a d and u on the thermal expansion due to heat release can be understood by examining the statistical behaviour of the dilatation rate ∂u i /∂x i . Figure 8 shows the PDFs of the normalised dilatation rate, ∂u i /∂x i ×δ th /S b (φ g = 1), for all droplet cases with φ d = 1.00 and φ d = 1.70 and compares the results with that of a premixed stoichiometric flame. It is evident from Fig. 8 that ∂u i /∂x i assumes predominantly positive values throughout the flame, but there is a finite probability of obtaining locally negative value of ∂u i /∂x i . The probability of obtaining high values of ∂u i /∂x i decreases on both the unburned and burned gas sides of the flame brush for all cases considered here. Under laminar flow conditions, the most likely dilatation rate for the premixed flame increases significantly as the reaction zone (i.e. 0.7 < c < 0.9) is approached from the unburned gas side and then drops beyond the flame (e.g. c = 0.9). In contrast, the most likely dilatation rate is much smaller in all of the droplet cases. For φ d = 1.00, the positive tail of the PDF for the laminar droplet cases extends further for smaller droplets, indicating a greater likelihood of thermal expansion in those cases. However, when the droplet equivalence ratio is increased to φ d = 1.70, this has an adverse effect on the thermal expansion for the laminar cases with small droplets, but a noticeably positive effect is observed for the laminar cases with the medium-size droplets. This reversal in behaviour arises once again due to the evaporation characteristics of the different size droplets: the small droplets, which evaporate easily, lead to a mixture too fuel-rich to sustain a high reaction rate and thus experience low thermal expansion, whereas the medium size droplets, which evaporate more slowly, benefit from the increase in the droplet equivalence ratio and thereby enjoy an increased reaction rate = .
= . and thermal expansion. Under turbulent flow conditions the extent of the positive tail of the PDF tends to increase with decreasing droplet size. It is further extended by increasing both u and φ d , such that the longest positive tail is attained for small droplets under high turbulence intensity with high droplet equivalence ratio: conditions which encourage rapid evaporation and facilitate good mixing of the gaseous fuel. The stoichiometric turbulent premixed flames show predominantly positive values of dilatation rate (i.e. the peak of the PDF is to be found when ∂u i ∂x i > 0) throughout the flame, although the exact magnitude of ∂u i ∂x i depends on the position within the flame. The location of the peak value for the turbulent droplet cases, when φ d = 1.00, is also positive, but small in magnitude. However, some droplet cases, when φ d = 1.70, show the peak of the PDF at higher values of ∂u i ∂x i . The availability of combustible fuel-air mixture increases with increasing φ d for all turbulent droplet cases, which leads to an increase in the probability of finding high positive ∂u i /∂x i values.

Combustion modes in flame-droplet interaction
The mode of combustion can be characterised with the help of a flame index: [46]. A positive flame index (i.e. FI > 0) indicates a premixed mode of combustion and a negative value (i.e. FI < 0) is indicative of a non-premixed mode. It is expected that droplet evaporation not only leads to inhomogeneities in the fuel mass fraction field, but also induces a non-premixed mode of combustion. The relative contributions of premixed and non-premixed modes of combustion change with u , a d and φ d . This can be seen from Fig. 9 which shows the percentages of heat release arising from premixed and non-premixed modes of combustion for different values of a d and φ d . The total heat release is given by H R total = Vω T dV whereω T = H φ |ω F | is the heat release = .
( ) = .   [46]. Figure 9 shows that medium and large droplets exhibit significant percentages of both modes of combustion such that the heat release for large droplets may be over 50 % due to non-premixed mode, whereas small droplets exhibit predominantly a premixed mode, such that in the laminar flow with φ d = 1.70 the percentage heat release due to premixed combustion is 100 %. The existence of premixed mode of combustion and its predominance for small droplet cases are consistent with previous findings based on two-dimensional simulations of unsteady laminar [28][29][30] and turbulent jet [23] spray flames. The evaporation rate decreases with increasing a d so the clouds of evaporated fuels for large droplets induce greater mixture inhomogeneity, thereby strengthening the relative contribution of non-premixed combustion to overall heat release. An increase in φ d strengthens the premixed mode of the small droplets and the non-premixed mode of the large droplets, but the effect of φ d for the medium sized droplets is non-monotonic in nature. An increase in φ d increases the number of evaporation sites, but a d affects the separation between the droplets for a given value of φ d . The level of mixture inhomogeneity and its distribution depend on both φ d and a d , and thus the relative strength of non-premixed combustion does not show any monotonic trend for the medium sized droplets. The evaporation takes place more rapidly and the released fuel vapour mixes more readily with surrounding air for small droplets than the medium and large droplets. Thus the level of mixture inhomogeneity decreases with increasing droplet equivalence ratio for small values of a d , leading to an increase in the contribution of premixed combustion with increasing φ d . Finally, it can be seen from Fig. 9 that the percentage of heat release from non-premixed mode of combustion increases with increasing u /S b (φ g = 1) for a given set of initial values of droplet sizes and droplet equivalence ratios. Turbulent advection acts to transport the evaporated fuel from the evaporation sites and mix it with pure air to produce a conducive environment to support non-premixed combustion. Furthermore, alongside augmented mixing, turbulent fluid motion also increases the rate of heat transfer to the droplets and promotes their evaporation, which in turn increases the probability of obtaining non-premixed combustion.

Reaction progress variable transport
By analysing the relative strengths of the terms in the transport equation for the reaction progress variable (see Eq. 18) it is possible to identify the terms which play important roles and in which regions of the flame this is the case. Figures 10 and 11 show that the variation of the mean values of the termsω c , N · ∇ ρD N · ∇c , (−2ρDκ m |∇c|),Ṡ c andȦ c conditional on the reaction progress variable c itself and on the mixture fraction ξ . The effects of turbulence intensity, droplet size and droplet equivalence ratio are exhibited in both figures and are compared with relevant results from a stoichiometric premixed flame under similar flow conditions. It can be seen from Fig. 10 that for the stoichiometric premixed flames only the mean values of reaction rateω c and the flame normal molecular diffusion rate N · ∇ ρD N · ∇c play significant roles, since the mean value of the tangential molecular diffusion rate is negligible for statistically planar flames, and the other terms are absent. In the unburned gas region (c < 0.5), it is only the normal molecular diffusion term N · ∇ ρD N · ∇c that is significant. The contribution of chemical reaction rateω c becomes important for c > 0.5 and the system is driven by the balance between the (positive) reaction rateω c and the (negative) flame normal molecular diffusion rate term N · ∇ ρD N · ∇c . This pattern is followed in the droplet cases as well, although the magnitudes of the mean values ofω c and N · ∇ ρD N · ∇c remain smaller than the corresponding terms in the stoichiometric premixed flame. The mean contribution of the tangential diffusion rate term (−2ρDκ m |∇c|) remains small for all the statistically planar cases considered here but locally this term can be significant and play an important role. In the preheat zone (i.e. c < 0.5 ), the flame normal molecular diffusion rate N · ∇ ρD N · ∇c plays a dominant role even for droplet cases.
In the region c > 0.5, the reaction progress variable transport is driven by the competition of the normal molecular diffusion rate N · ∇ ρD N · ∇c and the reaction rate termω c . In  the remainder of the figures in this paper, the region 0.0 ≤ c ≤ 0.5 is not shown, becauseω c assumes negligible magnitude in the preheat zone and thus the statistics in this region are determined by the molecular diffusion and droplet evaporation and not by chemical reaction. It can be seen from Fig. 10 that the magnitude of the terms in the reaction progress variable c transport equation increases with increasing φ d and with increasing u , but with decreasing droplet size. The only exception being the case of small droplets under laminar flow conditions for which increasing φ d is detrimental to the reaction rate term due to the development of non-combustible fuel-rich mixture in the vicinity of the droplets. In the region c > 0.5 only the mean values of reaction rateω c and normal molecular diffusion rate N ·∇ ρD N · ∇c are significant, and all other terms are negligible in that region. The mean contribution of the cross-scalar dissipation termȦ c assumes non-negligible values towards the unburned gas side of the flame for the cases with high u /S b φ g = 1 (e.g. see initial u S b φ g = 1 = 7.5 cases in Fig. 10) but this term remains negligible in comparison to the mean values ofω c and N · ∇ ρD N · ∇c for the rest of the flame front, which is consistent with previous observations by Malkeson and Chakraborty [33] in the context of gaseous stratified mixture combustion. The mean contribution of the evaporation termṠ c remains negligible in comparison to the mean values ofω c and N · ∇ ρD N · ∇c for the major part of the flame front. This is consistent with previous findings by Neophytou et al. [21] in the context of edge flame propagation for droplet-laden mixtures.  (−2ρDκ m |∇c|),Ṡ c andȦ c conditional on mixture fraction ξ values for the reaction progress variable range where the effects of chemical reaction and heat release are significant (i.e. 0.5 < c < 0.9 ). Two observations are apparent from all cases reported in Fig. 11: firstly, the mean values of both flame normal and tangential molecular diffusion rates (i.e. N · ∇ ρD N · ∇c and (−2ρDκ m |∇c|)) and the cross-scalar dissipation contri-butionȦ c assume peak values at ξ = ξ st , and secondly, the peak value of the mean reaction rateω c is obtained at a location close to ξ = ξ st but at ξ > ξ st where the mean values of N · ∇ ρD N · ∇c assumes negative values in accordance to the observations from Fig. 10.
The location corresponding to ξ = ξ st is associated with high temperature which encourages evaporation (and thus the high mean value ofȦ c close to ξ = ξ st due to evaporation induced mixture fraction gradient) leading to a considerable increase in the mean values of N · ∇ ρD N · ∇c and (−2ρDκ m |∇c|) due to the newly evaporated fuel vapour. The focusing of heat at the negative curvature locations encourage evaporation of droplets and gives rise to high magnitudes of positive (−2ρDκ m |∇c|), which contributes to the positive peak mean values of the tangential molecular diffusion rate in Fig. 11. The magnitudes of N · ∇ ρD N · ∇c and (−2ρDκ m |∇c|) decrease with increasing φ d . The larger number of cooler droplets at higher values of φ d leads to greater heat loss in the form of latent heat, thereby reducing the background gas temperature and the amount of evaporation that can take place. This leads to the reduction in the magnitudes ofȦ c , N · ∇ ρD N · ∇c and (−2ρDκ m |∇c|) for high values of φ d . The peak in the reaction rate term occurs at a location ξ ≈ 0.07, which corresponds to φ g ≈ 1.1, at which equivalence ratio the highest flame speed (and reaction rate) occurs in gaseous flames. The effect of increasing the turbulence intensity on the peak mean values of the molecular diffusion terms and reaction rate of reaction progress variable is most noticeable in the case of small droplets, less so for medium sized droplets and hardly at all for large droplets. Finally, the magnitude of the cross-scalar dissipation termȦ c increases with increasing droplet size: there is less mixture inhomogeneity for small droplets than for large droplets. This is because smaller droplets evaporate more easily and are more numerous than larger droplets, thus the fuel mass fraction field is more uniform for small droplets that for larger droplets, especially under the effect of turbulent mixing. The statistical behaviour ofω c , N · ∇ ρD N · ∇c , (−2ρDκ m |∇c|),Ȧ c andȦ c in response to c and ξ determines the statistical behaviour of flame propagation rate.

Flame propagation and density-weighted displacement speed statistics
It is possible to define a displacement speed based on a kinematic description of the reaction progress variable, c (i.e. the speed at which the flame surface moves normal to itself with respect to an initially coincident material surface, c = c * ) [44]: where S d , is the displacement speed, which is defined in the following manner [21,33,[47][48][49][50][51][52]: The effects of subdividing this range of c reveals that the probability of finding positive (negative) values of S * r (S * n ) increases from c = 0.5 to c = 0.9, which is consistent with the variations of the mean values ofω c and N ·∇ ρD N · ∇c shown in Fig. 10. These termsω c and N · ∇ ρD N · ∇c (S * r and S * n ) largely balance each other throughout the entire region 0.5 < c < 0.9 such that the statistical behaviour of the total density-weighted displacement speed S * d including its curvature and strain rate dependence remains qualitatively similar across this region.
The variations of the mean values of S * d and its components conditional on normalised values of curvature κ m and tangential strain rate a T = δ ij − N i N j ∂u i /∂x j acting on the reaction progress variable isosurface are shown in Figs. 13 and 14 respectively for the region corresponding to 0.50 ≤ c ≤ 0.90 for a d δ th = 0.06, 0.08 and 0.10 in the case of droplet equivalence ratios φ d = 1.00 and φ d = 1.70. The same qualitative behaviour has been observed for other values of φ d . The results from the droplet cases are compared with those of a stoichiometric premixed flame under same flow conditions. It is evident from Fig. 13 that S * d is negatively correlated with κ m for all cases and that the correlation coefficient remains greater than −1.0. Figure 13 reveals that the deterministic negative correlation between S * t = −2ρDκ m ρ 0 and κ m is responsible for the net negative correlation between S * d and κ m . It can be seen from Fig. 13 that S * r is negatively (positively) correlated with κ m for negative (positive) curvature values, whereas S * n remains only weakly correlated with κ m . The curvature κ m dependence of S * r and S * n originate principally due to κ m dependence on |∇ c|, which also exhibits both positive and negative correlating branches (not shown here). The complex curvature κ m dependences of the combined contribution of (S * r +S * n ) is consistent with several previous DNS findings in the context of turbulent premixed and edge flame propagation [21,33,41,42,[47][48][49][50][51][52] and the κ m dependences of the combined contribution of (S * r +S * n ) leads to a negative correlation between S * d and κ m , which has a negative correlation coefficient greater than -1.0. As S * z and S * s do not play important roles in the curvature κ m dependence of S * d because of their small magnitudes, / ( ) Fig. 14 Variation of mean values of normalised density-weighted displacement speed and its components, S * i /S b(φg=1) , conditional on normalised tangential strain rate a T ×δ th /S b(φg=1) values. Details are as described in Fig. 13 the curvature dependence of density-weighted displacement speed S * d in droplet cases is qualitatively similar to the turbulent premixed cases. Figure 14 shows that S * t and a T are strongly positively correlated, whereas both S * r and S * n are less strongly positively correlated with a T . The strength of all correlations decreases with increasing φ d . The observed strain rate dependence of (S * r +S * n ) originates due to the (positive) correlation between |∇c| and a T (not shown here), whereas a T and κ m are negatively correlated with each other (not shown here), which leads to a positive correlation between a T and S * t = −2ρDκ m ρ 0 . The positive correlation between S * t and a T overcomes the weak correlation between (S * r +S * n ) and a T giving rise to a net positive correlation between S * d and a T , which is consistent with several previous findings in the context of turbulent premixed and edge flame propagation [21,33,41,42,[47][48][49][50][51][52]. The additional S * d components in droplet combustion (i.e. S * z and S * s ) assume negligible values in the reaction zone and thus the strain rate dependence of density-weighted displacement speed S * d in droplet cases show qualitative similarities to the corresponding statistics in turbulent premixed flames. Figure 15 shows the mean values of S * d and its components conditional on |∇ξ | for the region corresponding to 0.50 ≤ c ≤ 0.90 for a d δ th = 0.06 , 0.08 and 0.10 in the case of φ d = 1.00 and φ d = 1.70. It is clear from Fig. 15 that there is a strong positive correlation between S * d and |∇ξ |, which originates from a positive correlation on S * r . For low φ d (i.e. φ d = 1.00) |∇ξ | also seems to have some correlation with S * t , but this disappears for sufficiently high values of φ d . Other components do not play a noticeable part in the correlation with |∇ξ |. It has been shown in Fig. 9 that a considerable amount of non-premixed  Fig. 15 Variation of mean values of normalised density-weighted displacement speed and its components, S * i /S b(φg=1) , conditional on normalised magnitudes of mixture fraction gradient, |∇ξ |×δ th . Details are as described in Fig. 13. Inset (top left) shows close-up of data combustion takes place in the droplet cases. A positive correlation between the reaction rate magnitude |ω F | with |∇ξ | is expected in non-premixed combustion [54] and this leads to a predominant positive correlation between S * r and |∇ξ | for the droplet cases considered here.
by a modified one-step mechanism [32] for the purpose of computational economy. Furthermore, a simple sub-grid evaporation model has been used to account for point source droplet evaporation. Thus, further analysis based on detailed chemistry and more elaborate droplet evaporation modelling will be necessary for deeper insight into the flame structure in flame-droplet interaction and quantitative predictions, which will form the basis of future investigations.