Statistical Behaviour and Modelling of Fuel Mass Fraction Dissipation Rate Transport in Turbulent Flame-Droplet Interaction: A Direct Numerical Simulation study

Three-dimensional Direct Numerical Simulation (DNS) data of statistically planar turbulent spray flames propagating into mono-disperse droplets for different values of droplet diameter a d and droplet equivalence ratio ϕ d has been used to analyse the statistical behaviour of the fuel mass fraction dissipation rate f ε Y and its transport in the context of Reynolds Averaged Navier-Stokes (RANS) simulations. Closures previously derived for high Damköhler number turbulent stratified mixture combustion have been shown not to capture the statistical behaviour of f ε Y for turbulent spray flames, because the underlying assumptions behind the original modelling are invalid for the cases considered in this analysis. The modelling of the unclosed terms of the fuel mass fraction dissipation rate f ε Y transport equation (i.e. the turbulent transport term T 1 , the density variation term T 2 , the scalar turbulence interaction term T 3 , the reaction rate term T 4 , the evaporation contribution terms T 5 and T 6 , and the dissipation rate term − D 2 ) has been analysed in the context of RANS simulations. The models previously proposed in the context of turbulent gaseous stratified flames have been considered here to assess their suitability for turbulent spray flames. Based on a-priori DNS analysis, suitable model expressions have been identified for T 1 , T 2 , T 31 , T 32 , T 33 , [ T 4 − D 2 + f ( D )] and [ T 5 + T 6 ], which have been shown to perform generally satisfactorily for all cases considered here.


Introduction
The closure of the mean reaction rate in the context of Reynolds Averaged Navier Stokes (RANS) simulations in turbulent combustion often requires the knowledge of the scalar dissipation rate of the fuel mass fraction Y F (i.e. e ε Y ¼ ρD∇Y ′′ F Á ∇Y ′′ F =ρ [1][2][3], where q, q e ¼ ρq=ρ and q ′′ ¼ q−q e are Reynolds average, Favre mean and Favre fluctuation of a general quantity q, D is the mass diffusivity and ρ is the gas density). Algebraic and transport equation based closures of e ε Y have previously been considered in the context of purely gaseous phase combustion where variations in equivalence ratio exist [2][3][4][5][6]. Analyses of simulations of turbulent droplet-laden combustion have also been carried out where the scalar dissipation rate transport equation has been considered [7]. Moreover, previous studies on droplet combustion analysed the statistical behaviour of the terms of the scalar dissipation rate transport equation [8], but the statistical analysis of e ε Y is yet to be addressed in detail. Furthermore, the validity of existing closures of e ε Y and the unclosed terms of its transport equation, which were originally proposed for purely gaseous phase combustion, is yet to be assessed for turbulent spray flames. These gaps in the existing literature have been addressed here by analysing the statistical behaviours of e ε Y and the terms of its transport equation using a three-dimensional compressible Direct Numerical Simulations (DNS) database [9][10][11][12] of statistically planar turbulent flames propagating into droplet-laden mixtures where the fuel is supplied in the form of mono-disperse droplets ahead of the flame. The current study considers selected cases from a large database [9][10][11][12] such that the effects of droplet diameter a d and droplet equivalence ratio ϕ d (i.e. fuel in liquid droplets to air ratio by mass, normalised by fuel to air ratio by mass under stoichiometric condition) on the statistical behaviours of e ε Y and its transport can be analysed in detail. The main objectives of this study are: (a) To analyse the statistical behaviours of e ε Y and the various unclosed terms of its transport equation for turbulent spray flames in the context of RANS. (b) To assess the validity of the existing models for the unclosed terms of e ε Y transport equation for turbulent droplet combustion.
The rest of the paper will be organised as follows. The mathematical background and numerical implementation pertinent to this analysis are presented in the next section. This will be followed by the presentation of results and their subsequent discussion. Finally, the main findings will be summarised and conclusions will be drawn.

Mathematical Background & Numerical Implementation
A modified single-step irreversible chemical mechanism [13] was used to perform the present analysis: Fuel + s · Oxidiser → (1 + s) · Products, where s is the oxidiser to fuel ratio by mass under stoichiometric condition. The activation energy and heat of combustion are taken to be functions of the gaseous equivalence ratio, ϕ g , so that a realistic ϕ g -dependence of unstrained laminar burning velocity S b ϕ g ð Þ can be obtained [13]. It has been shown by Tarrazo et al. [13] that the mechanism compares favourably with both experiments and detailed chemistry simulations for all hydrocarbon-air flames. It has been demonstrated by Swaminathan and Bray [14] based on experimental data that the normalised laminar burning velocity S b ϕ g ð Þ = S b ϕ g ð Þ n o max dependence of equivalence ratio ϕ g is not sensitive to the choice of fuel for hydrocarbon-air mixtures.
Moreover, it has also been found that there is no significant difference between ϕ g dependences of normalised laminar burning velocity S b ϕ g ð Þ = S b ϕ g ð Þ n o max and non-dimensional adiabatic flame temperature θ ϕ g ð Þ ¼ T ad ϕ g ð Þ −T 0 = T ad ϕ g¼1 À Á −T 0 À Á (with T ad ϕ g ð Þ and T 0 are the adiabatic flame temperature at gaseous equivalence ratio of ϕ g and unburned gas temperature respectively) obtained from the modified single-step [13] and multi-step detailed [15] chemical mechanisms.
Furthermore, the variations of S b ϕ g ð Þ = S b ϕ g ð Þ n o max and θ ϕ g ð Þ with ϕ g have been found to be in good agreement with experimental data [16].
In this analysis, all species are taken to have unity Lewis number and are assumed to be perfect gases. Standard values have been taken for the ratio of specific heats (γ = 1.4) and Prandtl number (Pr = 0.7) for the gaseous phase. The individual droplets are tracked in the Lagrangian sense and the quantities transported for each droplet are the position, x d , velocity, u d , diameter, a d and temperature, T d . The transport equations of x d , u d , a d and T d are given as [9][10][11][12][17][18][19][20][21][22][23]: where T is the instantaneous dimensional temperature, 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 [17][18][19][20][21][22]: where ρ d is the droplet density, C L p is the specific heat for the liquid phase, C g p is the specific heat at constant pressure for the gaseous phase, C u is the corrected drag coefficient and is given by [18][19][20][21][22][23]: Furthermore, 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 [9][10][11][12][17][18][19][20][21][22][23]: where Y s F is the value of Y F at the surface of the droplet. Equation (1iv) implicitly invokes the unity Lewis number assumption. The Clausius-Clapeyron relation 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: where T s ref is the boiling point of the fuel at pressure p ref , R 0 is the gas constant, T s d is assumed to be T d , and W air and W F are the molecular weights of air and fuel, respectively.
Droplet evaporation leads to mixture inhomogeneities, which are characterized by the mixture fraction: ξ = (Y F − Y O /s + Y O∞ /s)/(Y F∞ + Y O∞ /s), where Y F∞ = 1.0 (Y O∞ = 0.233) is the fuel (oxidiser) mass fraction in the pure fuel (air) stream. The fuel used here is n-heptane, C 7 H 16 , for which s = 3.52 and the stoichiometric fuel mass/mixture fraction is: Y Fst = ξ st = 0.0621. Using ξ, a reaction progress variable c can be defined in the following manner according to several previous analyses on droplet combustion [9][10][11][12][21][22][23][24][25] so that c increases monotonically from 0 in unburned reactants to 1.0 in fully burned products. It should be noted that the mixture fraction is not a passive scalar in the strict sense as there is an extra term in the transport equation of ξ due to evaporation. However, the Burke-Schumann relations where subscripts u and b refer to values in unburned reactants and fully burned products, respectively) remain reasonably valid for the oxygen mass fraction Y O because the evaporation contributions in the mixture fraction transport equation remain mostly negligible in both fully unburned and fully burned gases. Accordingly, using Eq. (1vi), one obtains the following transport equation for e ε Y [6]: Flow, Turbulence and Combustion ð2ivÞ where Γ is the source term in the mass conservation equation due to evaporation, The term T 1 is the turbulent transport term, T 2 arises due to density variation, T 3 originates due to the alignment of the scalar gradient with the fluiddynamic strain rates and essentially signifies generation/destruction of the scalar gradient by the velocity gradients, T 4 arises due to the chemical reaction rate, T 5 and T 6 arise due to droplet evaporation, whereas the term D 2 arises due to molecular dissipation. In Eq. (2), the terms of T 1 , T 2 , T 31 , T 32 , T 33 , T 4 , T 5 , T 6 and (−D 2 ) are unclosed terms in the context of second-moment closure and their modelling will be discussed in Section 3.
The present study uses a three-dimensional compressible DNS code SENGA [9][10][11][12][21][22][23]. High-order finite-difference (i.e. 10th central difference scheme for the internal grid points and the order of differentiation gradually reduces to a 2nd order one-sided scheme at the non-periodic boundaries) and explicit 3rd order low storage Runge-Kutta schemes are used for spatial differentiation and time advancement, respectively. A rectangular domain of size 63:3 Þ has been considered, where D 0 and S b ϕ g ¼1 ð Þ are the unburned gas diffusivity and the unstrained laminar burning velocity of the stoichiometric mixture, respectively. For the present thermo-chemistry is the unstrained thermal laminar flame thickness of the stoichiometric laminar premixed flame, where subscript 'L' refers to the values in the unstrained stoichiometric laminar premixed flame. The simulation domain is discretised using a Cartesian grid of size 384 × 256 × 256, ensuring that both the flame thickness, δ th , and the Kolmogorov length-scale, η, are adequately resolved. The boundaries in the mean direction of flame propagation (i.e. x-direction) are considered to be partially non-reflecting, whereas the other boundaries are taken to be periodic. The boundary conditions are specified using the well-known Navier Stokes Characteristic Boundary Conditions (NSCBC) technique [26]. The droplets are distributed uniformly in space throughout the y-and z-directions and in the region 0:0≤xS b ϕ g ¼1 ð Þ =D 0 ≤16:53 ahead of the flame. The reacting flow field is initialised based on the steady laminar solution generated using COSILAB [27] for desired values of a d and ϕ d , as done previously by Neophytou and Mastorakos [28] for one-dimensional laminar spray flame simulations. Initial turbulent velocity fluctuations, generated using a standard pseudo-spectral method [29] following Batchelor-Townsend spectrum [30], have been superimposed on top of the steady laminar spray flame solution. For the present analysis, the unburned gas temperature is taken as T 0 = 300K, which leads to τ ¼ T ad ϕ g ¼1 ð Þ −T 0 =T 0 ¼ 6:54 where T ad ϕ g ¼1 ð Þ is the adiabatic flame temperature of the stoichiometric mixture. The fuel is supplied purely in the form of mono-disperse droplets with a d /δ th = 0.06, 0.08, 0.10 for different values of ϕ d = 1.0, 1.7 at a distance 16:53D 0 =S b ϕ g ¼1 ð Þ from the point in the laminar flame at which T ¼ 400K. It should be noted that the droplet diameter a d has been non-dimensionalised using the welldefined thermal flame thickness of the stoichiometric mixture δ th , which is consistent with Buckingam's Pi Theorem. The initial droplet number density ρ N varies between 1.16 ≤ (ρ N ) 1/3 δ th ≤ 2.27 in the region 0:0≤xS b ϕ g ¼1 ð Þ =D 0 ≤16:53 and the liquid volume fraction remains well below 0.01. Droplets are supplied at the left-hand-side boundary to maintain a constant ϕ d ahead of the flame. Due to the high volatility of n-heptane, evaporation commences on entry and the droplet diameter decreases by at least 40%, 30% and 25% by the time it reaches the most reactive region of the flame for the initial a d /δ th = 0.06, 0.08, 0.10 cases respectively, such that the volume of even the largest droplets remains smaller than half that of the cell volume, which validates the sub-grid point source treatment of droplets adopted for flame-droplet interactions analysed here. The droplet diameter to grid size used in the current analysis remains comparable to several previous DNS analyses [18][19][20][21][22][23]31].
The cases considered here have initial values of normalised root-mean-square (rms) turbulent velocities u 0 =S b ϕ g ¼1 ð Þ ¼ 7:5 and non-dimensional 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. All simulations have been carried out until t final = max(3t turb , 4t chem ), where t turb = L 11 /u ′ and t chem ¼ are the initial eddy turnover time and chemical time, respectively. The simulation time remains either greater than or comparable to several previous analyses [18][19][20][21][22][23][32][33][34][35]. The volumeintegrated reaction rate, flame surface area and burning rate per unit area were not changing rapidly when the statistics were extracted [9]. The Reynolds/Favre averaged values of a quantity q (i.e. q and q e) are evaluated by ensemble-averaging q over the y-z plane at a given x-location. It should be noted that as the flames are statistically planar, the Favre averaged reaction progress variable c e and the scalar dissipation rate e ε Y and the terms of its transport equation are unique functions of the x 1direction, which is aligned with the mean direction of flame propagation. The spatial distribution of the terms of the scalar dissipation rate transport equation depending on the thickness of the flame brush, which changes from one case to another. Thus, to generalise the results, all the terms and model predictions are shown as a function of c e ¼ f x 1 ð Þ.

Figures 1a, b, and c present the instantaneous distributions of normalised fuel mass fraction
fields in the central x − z plane for a d /δ th = 0.08 and ϕ d = 1.0 at t = t chem , where the black dots indicate the droplets, which reside immediately adjacent to the plane shown. As they approach the flame, the droplets shrink due to evaporation. However, the droplets do not necessarily complete their evaporation until after passing through the flame. The evaporating droplets absorb latent heat from the background gas. This occurs on both sides of the flame, but is most noticeable on the burned gas side of the flame. For all cases considered in the current analysis, the reaction takes place predominantly under fuel-lean conditions and, therefore, the heat release due to combustion and the resultant burned gas temperature are lower than the adiabatic flame temperature of the stoichiometric mixture (i.e. T < 1.0) [9,10]. The predominantly fuel-lean combustion in these cases suggests a slow combustion process and low values of Damköhler number It should be noted that the fraction of total number of droplets that pass through the flame remains extremely small and furthermore, the evaporated vapour around the droplet on these occasions remain extremely fuel-rich due to rapid evaporation as a result of the combination of small droplet diameter and high temperature. Thus, the ignition delay of the evaporated vapour is likely to be much larger than the chemical time scale. This can be substantiated from the fact that the Damköhler number for the stoichiometric mixture Da about 0.67 for the cases considered here and, thus, the characteristic Damköhler number for fuel-rich mixture is expected to be even smaller. This suggests that the chemical time scale remains much greater than the mixing time scale, which scales with the eddy turn over time. Therefore, turbulent mixing enables mixing of the evaporated fuel vapour and leads to premixed and non-premixed mode of combustion before localised ignition takes place.
Furthermore, Chiu and Liu [36] proposed a group number, where Le and Sc are the Lewis and Schmidt numbers, respectively, s d is the distance between droplets, 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 in the current study come under the category of external sheath combustion (i.e. have values of G much greater than unity). The flames with small value of G usually exhibits high temperature, which promotes pyrolysis at the fuel-rich core, whereas the temperature values in the external sheath combustion are usually not high enough to give rise to significant amount of pyrolysis [37]. For the present configuration, the burned gas temperature remains mostly smaller than the adiabatic flame temperature of the stoichiometric mixture, and thus 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 sufficient to mimic the pyrolysis process. Furthermore, pyrolysis does not directly affect the statistics of scalar gradients and scalar dissipation rate which is the main focus of the current analysis. A recent DNS analysis [31] with comparable simulation parameters as that of the current analysis demonstrated with the help of Kerstein-Law parameter that either group combustion or individual burning are unlikely under these conditions. It should also be noted that detailed discussion of the flame-droplet interaction under both laminar and turbulent flow conditions for the current DNS database can be found elsewhere [9][10][11][12] and thus is not repeated here. The aforementioned references [9][10][11][12] also compared the global flame behaviour of the turbulent spray flames to the corresponding purely premixed gaseous case and interested readers are directed to this study for comprehensive discussions on these issues.
The variations of the fuel mass fraction dissipation rate e ε Y with c e across the flame-brush are shown in Figs. 2a-f. The scalar dissipation rate of fuel mass fraction, e ε Y , assumes non-zero values close to c e ¼ 0, which decrease before rising and attaining a maximum roughly at the middle of the flame-brush for all cases. The magnitude of e ε Y then decreases as c e approaches unity. A discussion of the variation of e ε Y across the flame brush has been provided in the following sections.
The fuel mass fraction dissipation rate e ε Y can be modelled using an algebraic expression following the linear relaxation approach as e where C YLR is a model parameter, k e is the turbulent kinetic energy and ε e is the dissipation rate of turbulent kinetic energy. Alternatively, Mura et al. [3] proposed algebraic expressions for e ε Y for turbulent flames with variations in equivalence ratio based on a presumed probability density function (pdf) approach with the Favre joint PDF between Y F and ξ (i.e. P e Y F ; ξ ð Þ¼ρP Y F ; ξ ð Þ=ρ) for turbulent stratified gaseous mixture combustion as: where P e ξjY F ð Þis the Favre PDF of ξ conditional on Y F and the quantities Y max (ξ) = ξ and Y min (ξ) = A(ξ)(ξ − ξ st ) are maximum and minimum values of Y F according to the Burke-Schumann relations [38] where A(ξ) = H(ξ − ξ st )/(1 − ξ st ) with H(ξ − ξ st ) being a Heaviside function. For Da ≫ 1, the last term on the right hand side of Eq. (3) disappears and λ w is unlikely to depend on ξ, which yields P e ξjY max ð Þ¼P e ξjY min ð Þ [3] such that for e ε Y in the following manner: Mura et al. [3] further proposed the following model for e ε Y : where C Y is a model parameter and S is a segregation factor that is defined as [3]: Flow, Turbulence and Combustion The performances of the models given by e ε Y ¼ C YLR ε e=k e À Á Y 00  interaction [39]) are shown in Figs. 2a-f in comparison to the corresponding quantity obtained from the DNS data for all cases. It is evident from Figs. 2a-f that none of the models considered here perform satisfactorily and, therefore, a modelled transport equation may need to be considered for e ε Y in the turbulent combustion of droplet-laden mixtures. It should be noted, however, that none of the models considered here were proposed in the context of low Damköhler number spray combustion and, therefore, shortcomings in the models are to be expected and there are no criticisms levied against them. However, as a starting point, it is important to assess their performance before analysing the transport equation based closure of e ε Y .

Statistical Behaviours of the Unclosed Terms of the e ε Y Transport Equation
The variations of the unclosed terms T 1 , T 2 , T 3 , T 4 , T 5 (here multiplied by 0.05 for visualisation purposes), T 6 , −D 2 (here multiplied by 0.05 for visualisation purposes) and f(D) with c e across the flame-brush for all cases considered are shown in Figs. 3a-f. The turbulent transport term T 1 is shown to exhibit large negative values towards the unburned gas side of the flame-brush (i.e. c e≈0), but remains small in comparison to other leading order terms for the majority of the flame-brush, obtaining both negative and positive values. The density variation term T 2 remains positive throughout the flame-brush with a peak value that moves from the unburned gas side to the burned gas side of the flame-brush for increasing initial droplet size a d /δ th . Moreover, the peak value for T 2 increases with increasing ϕ d and increasing a d /δ th . The scalar turbulence interaction term T 3 is shown to exhibit large positive values towards the unburned gas side of the flame-brush (i.e. c e≈0) before reducing in magnitude, but remaining mainly positive across the flame-brush. The relative importance of T 3 can be seen to increase with increasing ϕ d for medium and large droplets, and with increasing a d /δ th for the cases considered here. The statistical behaviour of the scalar-turbulence interaction term T 3 is determined due to the competition between the turbulent straining and the chemical heat release induced strain rate. Through scaling analysis, it has been postulated by Swaminathan and Bray [40] that the contribution of T 32 is likely to dominate other components of T 3 for large and moderate values of Damköhler number which can be observed from Figs. 4a-f where the contributions to the T 3 term (i.e. T 31 , T 32 and T 33 ) across the flame-brush for all cases considered in the current study are shown. Figures 4a-f show that T 31 exhibits large positive values and T 33 large negative values towards the unburned gas side of the flame-brush (i.e. c e≈0). However, in most cases considered here both T 31 and T 33 assume negligible values across the entire flame brush (i.e. c e > 0), with the exception of ϕ d = 1.7 for which T 33 assumes noticeably negative values for small droplets. This effect decreases with increasing droplet size. It has previously been shown that the scalar turbulence interaction term T 3 and its components, particularly T 31 and T 32 , are influenced directly by the inner products between the velocity and the scalar gradients. Therefore, the behaviour of T 3 is principally determined by the alignment between the scalar gradient ∇Y F and the local principal strain rates. For the currently considered cases, Wacks and Chakraborty [11] have shown that ∇Y F is predominantly aligned with the most compressive principal strain rate across the flame brush, where straining due to turbulence dominates the straining induced by heat release. These observations are consistent with the behaviour of T 31 , T 32 and T 33 in turbulent stratified flames under globally fuel-lean conditions [6].
The reaction rate contribution term T 4 is shown to assume positive values across the flamebrush but, given that the combustion is predominantly fuel-lean, remains small in comparison to the other leading order terms, which is consistent with earlier analysis of globally fuel-lean flames [6]. The droplet evaporation term T 5 has been found to be a leading order contributor in all cases exhibiting large positive values towards the unburned gas side of the flame-brush (i.e. c e≈0) before reducing in magnitude, but remaining positive across the flame-brush. The second peak value in the flame-brush is shown to shift towards the burned gas side of the flame-brush with increasing a d /δ th . The additional evaporation term T 6 has been found to be positive across the flame-brush, but small in magnitude in comparison to T 5 with a peak value in the flamebrush shown to shift towards the burned gas side of the flame-brush increasing a d /δ th . Furthermore, the magnitudes of T 5 and T 6 have been shown to increase with increasing ϕ d (a d / δ th ) for a given value of a d /δ th (ϕ d ). The behaviour of both evaporation terms T 5 and T 6 , and more noticeably of T 5 , close to c e ¼ 0 can be attributed to the aforementioned high volatility of n-heptane, which allows for evaporation to commence immediately upon entry of the droplets into the flame brush. The rate of evaporation is high at the entry point (c e ¼ 0) of the droplets into the flame brush due to the paucity of gaseous fuel in this region (see Fig. 1a). The rate of evaporation subsequently slows down (i.e. the magnitudes of T 5 and T 6 drop) as the droplets traverse the wide region of unburnt gas in which the value of c e remains negligible (see Fig.  1b). This effect is due to a drop in the Spalding mass transfer number (B d ) because its magnitude depends on the difference between the local value of Y F and the value of Y F on the droplet surface, Y s F (see Eq. (1iv)), which diminishes as evaporation progresses. Smaller droplets evaporate at a faster rate than larger ones and thus the local value of Y F approaches that of Y s F more quickly for smaller droplets. Consequently, the gradient of the decrease in T 5 and T 6 for small c e is noticeably greater for smaller droplets than for larger ones. The molecular dissipation term (−D 2 ) has been found to be a leading order contributor in all cases. It exhibits large negative values towards the unburned gas side of the flame-brush (i.e. c e≈0) before reducing in magnitude, but remaining negative across the flame-brush. The second peak value in the flame-brush is shown to shift towards the burned gas side of the flame-brush with increasing a d /δ th . The magnitude of (−D 2 ) has been shown to increase with increasing ϕ d (a d / δ th ) for a given value of a d /δ th (ϕ d ). The additional term f(D) is shown to be small in all cases relative to the other terms. The increasing trend of the magnitudes of T 5 , T 6 and (−D 2 ) with increasing a d /δ th (ϕ d ) for a given value of ϕ d (a d /δ th ) is consistent with increasing values of e ε Y under these conditions. The physical reasons discussed earlier for explaining the increasing trend of e ε Y with increasing a d /δ th (ϕ d ) for a given value of ϕ d (a d /δ th ) are also applicable for the increases in magnitudes of T 5 , T 6 and (−D 2 ).
It should be noted that the scalar dissipation rate transport of fuel mass fraction e ε Y in the purely gaseous premixed flame and its modelling are both qualitatively and quantitatively similar to the results shown Malkeson and Chakraborty [6] for turbulent premixed and stratified flames. These results are not repeated here for the sake of brevity but interested readers are directed to Ref. [6] for detailed discussion of the modelling of e ε Y transport in the purely gaseous and stratified flames.

Modelling of T 1
According to the scaling arguments of Swaminathan and Bray [40], the modelling of T 1 translates to the modelling of T 11 for large turbulent Reynolds numbers Re t . Therefore, in the context of statistically planar flames, the turbulent transport term becomes equal to The quantity ρu ′′ 1 ε Y is often modelled using the gradient hypothesis 1 where μ t ¼ 0:09ρk e 2 =ε e is the eddy viscosity and σ is an appropriate turbulent Schmidt number. An alternative model for ρu ′′ 1 ε Y , capable of addressing both gradient and counter-gradient behaviour was proposed for flames with varying model parameter ψ m in the following manner [6]: The model given by Eq. (7) was originally proposed in the context of turbulent stratified flames [6], allowing for both gradient and counter-gradient behaviour which had previously been reported by Veynante et al. [41] and Chakraborty and Cant [42]. The form of the model provided is consistent with previously proposed models for the turbulent flux of the reaction progress variable in the context of turbulent premixed flames [42][43][44]. A detailed explanation of the formation of the model given by Eq. (7) has already been provided by Malkeson and Chakraborty [6] and will not be repeated here for the sake of brevity. The performance of − μ t =σ ð Þ∂ e ε Y =∂x 1 (where σ = 1) and Eq. (7) (where ψ m = 0.25 × Y Fst ) compared to ρu ′′ 1 ε Y obtained from DNS is shown in Figs. 5a-f. It is evident from Fig. 5a- neither the qualitative nor the quantitative behaviour of ρu ′′ 1 ε Y , particularly towards the burned gas side of the flame-brush. Upon examination of Fig. 5, Eq. (7) has been found to demonstrate improved performance in capturing the general behaviour of ρu ′′ 1 ε Y for certain parts of the flame brush (e.g. Fig. 5a and b towards the burned gas side) but general overall performance of Eq. (7) can be considered to be comparable to − μ t =σ ð Þ∂ e ε Y =∂x 1 .

Modelling of T 2
It has been observed that the density variation term T 2 remains positive throughout the flamebrush in all cases. Furthermore, the term T 2 has been found to act as a leading order contributor in all cases considered. A previous study by Chakraborty et al. [45] proposed the following model for T 2 for perfectly premixed flames, which is consistent with the scaling arguments of Swaminathan and Bray [40], in the following manner: Þand A is a model parameter which is of the order of unity, but is dependent upon thermochemistry. In Eq. (8), ρ b is the burned gas density and S b is the characteristic burning velocity, which are evaluated as [46]: where ϕ min and ϕ max are minimum and maximum values of ϕ, p(ϕ) is the pdf of ϕ and ρ b(ϕ) is the burned gas density for an unstrained planar perfectly premixed flame at equivalence ratio ϕ. The function f(Ka L ) is given by f(Ka L ) = (1 + Ka L ) −1/2 according to Chakraborty et al. [45] where Ka L ≈ S b À Á −3=2 ε eδ b ð Þ 1=2 is the local Karlovitz number and δ b ¼ 2D 0 =S b is considered to be the characteristic flame thickness. The Karlovitz number dependence of Eq. (8) ensures that the contribution of heat release weakens with increasing Ka L as the broken reaction zones regime is approached. Detailed explanations for the formation of the model given by Eq. (8) have been presented elsewhere [6,[43][44][45], which interested readers are directed to, and will not be repeated here for the sake of brevity. However, it has been found that the model parameter A in Eq. (8), must account for the increase in the magnitude of T 2 with increasing droplet diameter. Accordingly, in the current study, a reasonable agreement has been found the normalised diameter of monodisperse droplets) was used as shown in Figs. 6a-f where the performance of Eq. (9) is shown in comparison to T 2 obtained from the DNS data.

Modelling of T 31
To close the behaviour of the scalar-turbulence interaction term T 3 , the components T 31 , T 32 and T 33 need to be modelled. Although for all cases considered here the components T 31 and T 33 are of negligible magnitude in comparison to the magnitude of component T 32 , this does not necessarily hold true for other choices of parameter values. Furthermore, since there are already existing models for T 31 and T 33 , it would be remiss to omit an assessment of these models. Several models have previously been proposed for T 31 . For turbulent premixed flames, the following model was proposed [47]: where C P1 is a model parameter [47]. Additionally, Mura et al. [48] proposed the following models for turbulent premixed high Damköhler number flames: (e) (f) where C PM is a model constant and n F ¼ −∇Y F = ∇Y F j j. Alternatively, Chakraborty and Swaminathan [49] proposed the following model to account for low Damköhler number premixed combustion: where is a local density-weighted Damköhler number and C 1 = 0.5 and Þ −2 are model parameters [49]. Malkeson and Chakraborty [6] extended this model for stratified combustion as: i s t h e B r a y n u m b e r , C * ! and ψ * = 2.0 are the model parameters [6]. It is worth noting that the model expressions given by Eqs. (10) (11) and (12) are in accordance with the non-reacting turbulent flow scaling according to Tennekes and Lumley [50], as used by Mantel and Borghi [47] and reacting flow scaling by Swaminathan and Bray [40], respectively. By contrast, the model expressions given by Eqs. (13) and (14) satisfy both reacting turbulent flow scaling by Mantel and Borghi [47] and reacting flow scaling by Swaminathan and Bray [39]. Detailed explanations of the formulation of the models given in Eqs. (10)- (14) have been presented elsewhere [6,[43][44][45][47][48][49], which interested readers are directed to, and are not repeated here for the sake of brevity. The variation of T 31 with c e across the flame-brush is shown for all cases in Figs. 7a-f along with the model predictions for Eqs. (10)- (14). It is evident that models given by Eqs. (10)- (12) do not capture the behaviour of T 31 but, given that they were not proposed in the context of droplet-laden mixture combustion, this is not unexpected. However, Figs. 7a-f demonstrate that Eqs. (13) and (14) can capture the behaviour of T 31 for the cases considered here with the best performance being shown by Eq. (13).

Modelling of T 32
A number of models have also been proposed for the contribution T 32 . For turbulent premixed flames, the following model was proposed by Mura and Borghi [51] which modified the model proposed by Mantel and Borghi [47]: where A e is a model parameter in the order of unity [51]. Mura et al. [52] subsequently proposed the following models for turbulent premixed high Damköhler number flames: where A M1 = 1, A M2 = 1, C MA = 0.6 and C MB = 1.6 are model parameters [52] and is a local Damköhler number.
(e) (f) Furthermore, Chakraborty and Swaminathan [49] proposed the following model to account for low Damköhler number premixed combustion: where C 3 = 1.5 and C 4 = 1.1(1 + Ka L ) −0.4 are the model parameters. It is worth noting that the model given by Eq. (15) and the first terms on the right hand side of Eqs. (16)-(18) account for the generation of dissipation rate due to the preferential collinear alignment of ∇Y F ′′ with the most compressive principal strain rate, which is qualitatively similar to the passive scalar mixing [9,49,53]. However, ∇Y F ′′ may preferentially align with the most extensive principal strain rate leading to destruction of e ε Y when the flame normal acceleration dominates over turbulent straining, and this strengthens with increasing Damköhler number [44,45,49,53]. This is accounted for by the negative contribution on the right hand side of the model expressions given by Eqs. (16)- (18). Moreover, Eq. (15) and the first term on the right hand side of Eqs. (16)- (18) are in accordance with the non-reacting turbulent flow scaling [50] used by Mantel and Borghi [47], whereas the second term on the right hand side of Eqs. (16)- (18) follows the reacting flow scaling by Swaminathan and Bray [39]. Detailed explanations of the formulation of the models given in Eqs. (15)- (18) have been presented elsewhere [6, 43-45, 47, 49, 52], and thus are not repeated here.
The variations of T 32 with c e across the flame-brush is shown for all cases in Figs. 8a-f along with the model predictions for Eqs. (15)- (18). It should be noted that Eqs. (15)- (17) for the considered cases can exhibit very similar values and as such their lines frequently coincide. It is evident that the models given by Eqs. (15)- (17) do not capture the behaviour of T 32 but, given that they were not proposed in the context of droplet-laden mixture combustion, this is not unexpected. However, Figs. 8a-f demonstrate that Eq. (18) can capture the general qualitative and quantitative behaviour of T 32 for the considered cases.

Modelling of T 33
Several models have also been proposed for the contribution T 33 . For turbulent premixed flames, the following model was proposed [47]: where C P2 is a model parameter in the order of unity [47]. Mura et al. [48] subsequently proposed the following models for turbulent premixed high Damköhler number flames: where A M1 = 1, A M2 = 1, C MA = 0.6 and C MB = 1.6 are model parameters [48] and Da L ¼ S b k e = δ b ε e ð Þ is a local Damköhler number. Furthermore, Chakraborty and Swaminathan [49] proposed the following model to account for low Damköhler number premixed combustion: where C T3 ¼ 1 þ a Ka −0:23 L À Á is a model parameter and . Equation (19) was originally proposed for passive scalar mixing and it follows the non-reacting turbulent flow scaling arguments [50] used by Mantel and Borghi [47], whereas the model expressions given by Eqs. (20) and (23) are in accordance with both the non-reacting turbulent flow scaling [47] and reacting flow scaling by Swaminathan and Bray [40]. Detailed explanations of the formulation of the models given in Eqs. (19)- (23) have been presented elsewhere [6,[43][44][45][47][48][49], and are not repeated here for the sake of brevity. The variations of T 33 with c e across the flame-brush is shown for all cases in Figs. 9a-f along with the model predictions for Eqs. (19)- (23). It is evident that Eqs. (21) and (22) do not capture the qualitative behaviour of T 33 for all cases considered but, given that these models were not proposed in the context of droplet-laden mixture combustion, this is not unexpected. However, Figs. 9a-f demonstrate that Eqs. (19), (20) and (23) can capture the general qualitative behaviour of T 33 for all considered cases with the best quantitative agreement being found for Eqs. (20) and (23).

Modelling of T
The combined contribution of the terms D 1 , T 4 , (−D 2 ) and f(D) can be written as [44,45]: where m Þis the local displacement speed. Equation (24) indicates that the net contribution of [D 1 + T 4 − D 2 + f(D)] signifies the effect due to flame normal propagation and flame curvature. As the molecular diffusion term D 1 is a closed term and it is often negligible in the context of RANS, it is often convenient to model the net contribution of [T 4 − D 2 + f(D)] rather than its individual components [6,[43][44][45]47]. It has been shown elsewhere [44,45] the net contribution of [T 4 − D 2 + f(D)] remains negative because of dominant contribution of the on the right hand of Eq.
(24). It has been observed that the net contribution of T 4 − D 2 + f(D) is negative across the flame-brush for all cases considered in the current study, which can be substantiated from Figs. 10a-f. This behaviour of the contribution of (T 4 − D 2 + f(D)) is dominated by the contribution of (−D 2 ), which is determined by the gradient of the fuel mass fraction (see Eq. (2viii)). Here, T 4 − D 2 + f(D) is modelled in the following manner according to previous studies [6,[43][44][45]47]: where β 2Y is a model parameter. The magnitude of (T 4 − D 2 + f(D)) increases with increasing droplet diameter a d reflecting the increases in the gradients of fuel mass fraction, which is accounted for by the model parameter β 2Y . Accordingly, satisfactory performance for Eq. (25) has been found when shown in Figs. 10a-f for all cases considered in the current study. 3.10 Modelling of (T 5 + T 6 ) To the best of the authors knowledge, there are currently no models in the existing literature for the terms arising due to droplet evaporation (i.e. T 5 and T 6 ). The contribution of (T 5 + T 6 ) has been found to be positive across the flame brush in all cases considered and, from a physical perspective, accounts for the contributions due to droplet evaporation. Therefore, from a modelling perspective, this term can be related to the mixture fraction dissipation rate ε e ξ , as the micro-mixing of mixture fraction variation induced by droplet evaporation is expected to affect the statistical behaviour of the evaporation contribution (T 5 + T 6 ) in the gaseous phase. It is worth considering that the magnitude of (T 5 + T 6 ) has qualitative and quantitative similarities to the magnitude of (T 4 − D 2 + f(D)), a similar form of model expression has been considered. It has been observed that the behaviour of the combined contribution of (T 5 + T 6 ) can be captured by the following expression: where β Y is a model parameter which provides reasonable agreement when The prediction of the model given by Eq. (26) along with the corresponding quantity obtained from DNS is shown in Figs. 11a-f for all cases considered in the current study, which show reasonable performance by the model.

Future Considerations
The current study considers a-priori DNS modelling of the fuel mass fraction dissipation rate f ε Y and the unclosed terms of its transport equation following a similar approach to several previous analyses [5, 6, 33, 39-45, 48-50, 54-58]. It should be noted that the strength of the a-priori modelling employed is to identify forms of models that can capture the behaviour of the quantities of interest, having access to all quantities from the DNS data that could be used for modelling. However, it is acknowledged that aposteriori analysis of the currently proposed models for the unclosed terms of the fuel mass fraction dissipation rate transport equation in the context of turbulent droplet-laden mixtures is required. Moreover, any such a-posteriori analyses of the models proposed in the current study must be carefully considered. For example, any a-posteriori analysis of these algebraic models to determine suitability faces additional issues as other quantities that these models are dependent upon the quantities, which need to modelled themselves (e.g. k e and ε e ). Therefore, determining the suitability, or not, of these models based upon such a-posteriori tests alone could be influenced by errors in modelling elsewhere which might lead to unsound conclusions. Furthermore, in actual RANS calculations the modelling and numerical errors interact in a complex non-intuitive manner. Thus, the results of a-posteriori assessment are likely to be problem dependent and code-specific so there is no reason to consider these findings as a definitive proof of the model performances. Therefore, both a-priori and a-posteriori analyses are necessary. This is beyond the scope of the current study and will form the basis of future investigations. Moreover, it is acknowledged that the effects of a detailed chemical mechanism have not yet been considered, as the current study employed a single step chemical reaction. The DNS simulations that have been considered in the current study are taken from a much larger database of simulations as part of a parametric investigation of the turbulent combustion of droplet-laden mixtures [9][10][11][12]. As such, the use of a detailed chemical mechanism for such parametric investigations where a range of simulation parameters have been considered is not feasible from the perspective of computational economy. However, it was previously demonstrated that the models proposed for the scalar dissipation rate and also for the unclosed terms of its transport equation based on apriori analysis of simple chemistry DNS data also remains valid in the presence of detailed chemistry and transport for turbulent premixed flames [57,58]. A similar outcome is expected also for turbulent spray combustion because the underlying physical mechanisms governing scalar dissipation rate and its transport are the same for both = .
= . (e) (f) simple and detailed chemical mechanisms. Nonetheless, future research in these directions will be necessary for a comprehensive assessment of the model performances.

Conclusions
The statistical behaviours of f ε Y and the unclosed terms of its transport equation have been analysed using three-dimensional DNS of statistically planar turbulent flames for which the fuel is supplied in the form of mono-disperse droplets for different a d and ϕ d . An algebraic closure based on presumed distribution of P e Y F ; ξ ð Þwhich was originally intended for high Damköhler number gaseous phase combustion does not adequately predict f ε Y obtained from DNS data. The behaviours of the unclosed terms of f ε Y transport equation have been analysed in the context of RANS simulations. It has been found that the density variation, evaporation and molecular dissipation contributions (i.e. T 2 , T 5 and −D 2 ) play significant roles in f ε Y transport. The suitability of the models previously proposed in the context of turbulent gaseous stratified flames have been assessed for the modelling of f ε Y transport in turbulent spray flames. Based on a-priori DNS analysis suitable model expressions have been identified for T 1 , T 2 , T 31 , T 32 , T 33 , [T 4 − D 2 + f(D)] and [T 5 + T 6 ], which have been shown to perform generally satisfactorily for all cases considered here. Further consideration of the modelling of f ε Y transport equation for spray flames is necessary as the current analysis deals with monodisperse droplets. It remains to be seen if the models for the scalar dissipation rate transport can be applicable to polydisperse systems if the Sauter mean diameter in the unburned gas is considered as being representative of such systems. Moreover, detailed chemical mechanism have not yet been considered, as the current study employed a modified single step chemical reaction. Thus, future research in these directions will be necessary for a comprehensive assessment of the model performances. Furthermore, the implementation of the proposed models in actual RANS simulations will be necessary for the purpose of a-posteriori assessment.