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

Abstract

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 $\overset{\sim }{\upvarepsilon_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 $\overset{\sim }{\upvarepsilon_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 $\overset{\sim }{\upvarepsilon_Y}$ transport equation (i.e. the turbulent transport term T₁, the density variation term T₂, the scalar turbulence interaction term T₃, the reaction rate term T₄, the evaporation contribution terms T₅ and T₆, and the dissipation rate term −D₂) 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₁, T₂, T₃₁, T₃₂, T₃₃, [T₄ − D₂ + f(D)] and [T₅ + T₆], which have been shown to perform generally satisfactorily for all cases considered here.

1 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. $\overset{\sim }{\varepsilon_Y}=\overline{\uprho D\nabla {Y}_F^{\prime \prime}\cdotp \nabla {Y}_F^{\prime \prime }}/\overline{\uprho}$ [1,2,3], where $\overline{\mathrm{q}}$, $\overset{\sim }{\mathrm{q}}=\overline{\uprho \mathrm{q}}/\overline{\uprho}$ and ${\mathrm{q}}^{\prime \prime }=\mathrm{q}-\overset{\sim }{\mathrm{q}}$ 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 $\overset{\sim }{\varepsilon_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 $\overset{\sim }{\varepsilon_Y}$ is yet to be addressed in detail. Furthermore, the validity of existing closures of $\overset{\sim }{\varepsilon_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 $\overset{\sim }{\varepsilon_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 $\overset{\sim }{\varepsilon_Y}$ and its transport can be analysed in detail. The main objectives of this study are:

• (a) To analyse the statistical behaviours of $\overset{\sim }{\varepsilon_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 $\overset{\sim }{\varepsilon_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.

2 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\left({\phi}_g\right)}$ 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\left({\phi}_g\right)}/{\left\{{S}_{b\left({\phi}_g\right)}\right\}}_{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\left({\phi}_g\right)}/{\left\{{S}_{b\left({\phi}_g\right)}\right\}}_{max}$ and non-dimensional adiabatic flame temperature ${\theta}_{\left({\phi}_g\right)}=\left({T}_{ad\left({\phi}_g\right)}-{T}_0\right)/\left({T}_{ad}\left({\phi}_{g=1}\right)-{T}_0\right)$ (with ${T}_{ad\left({\phi}_g\right)}$ and T₀ are the adiabatic flame temperature at gaseous equivalence ratio of ${\phi}_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\left({\phi}_g\right)}/{\left\{{S}_{b\left({\phi}_g\right)}\right\}}_{max}$ and ${\theta}_{\left({\phi}_g\right)}$ 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, ${\overrightarrow{x}}_d$, velocity, ${\overrightarrow{u}}_d$, diameter, a_d and temperature, T_d. The transport equations of ${\overrightarrow{x}}_d$,${\overrightarrow{u}}_d$, a_d and T_d are given as [9,10,11,12, 17,18,19,20,21,22,23]:

$$ \frac{d{\overrightarrow{x}}_d}{dt}={\overrightarrow{u}}_d;\frac{d{\overrightarrow{u}}_d}{dt}=\frac{\overrightarrow{u}\left({\overrightarrow{x}}_d,t\right)-{\overrightarrow{u}}_d}{\tau_d^p};\frac{d{a}_d^2}{dt}=-\frac{a_d^2}{\tau_d^u}\ \mathrm{and}\ \frac{d{T}_d}{dt}=-\frac{\hat{T}\left({\overrightarrow{x}}_d,t\right)-{T}_d-{B}_d{L}_v/{C}_P^g}{\tau_d^T} $$

(1i)

where $\hat{T}$ is the instantaneous dimensional temperature, L_v is the latent heat of vaporization, and ${\tau}_d^p$,${\tau}_d^u$ and ${\tau}_d^T$ are relaxation/decay timescales for droplet velocity, diameter and temperature respectively, which are defined as [17,18,19,20,21,22]:

$$ {\tau}_d^p=\frac{\rho_d{a}_d^2}{18{C}_u\mu };{\tau}_d^u=\frac{\rho_d{a}_d^2}{4\mu}\frac{Sc}{S{h}_c}\frac{1}{\mathit{\ln}\left(1+{B}_d\right)}\ \mathrm{and}\ {\tau}_d^T=\frac{\rho_d{a}_d^2}{6\mu}\frac{\mathit{\Pr}}{N{u}_c}\frac{B_d}{\mathit{\ln}\left(1+{B}_d\right)}\frac{C_p^L}{C_p^g} $$

(1ii)

where ρ_d is the droplet density, ${C}_p^L$ is the specific heat for the liquid phase, ${C}_p^g$ 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]:

$$ {C}_u=1+\frac{1}{6}{\mathit{\operatorname{Re}}}_d^{2/3} $$

(1iii)

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]:

$$ R{e}_d=\frac{\rho \left|\overrightarrow{u}\left({\overrightarrow{x}}_d,t\right)-{\overrightarrow{u}}_d\right|{a}_d}{\mu };{B}_d=\frac{Y_F^S-{Y}_F\left({\overrightarrow{x}}_d,t\right)}{1-{Y}_F^s}\ \mathrm{and}\ S{h}_c=N{u}_c=2+\frac{0.555R{e}_d Sc}{{\left(1.232+R{e}_d{Sc}^{4/3}\right)}^{1/2}} $$

(1iv)

where ${Y}_F^s$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}_F^s$, is used to evaluate the Spalding number B_d, which leads to:

$$ {p}_F^s={p}_{ref}\mathit{\exp}\left(\frac{L_v}{R_0}\left[\frac{1}{T_{ref}^s}-\frac{1}{T_d^s}\right]\right);{Y}_F^s={\left(1+\frac{W_{air}}{W_F}\left[\frac{p\left({\overrightarrow{x}}_d,t\right)}{p_F^s}-1\right]\right)}^{-1} $$

(1v)

where ${T}_{ref}^s$ is the boiling point of the fuel at pressure p_ref, R₀ is the gas constant, ${T}_d^s$ is assumed to be T_d, and W_air and W_F are the molecular weights of air and fuel, respectively.

The droplet and gaseous phases are coupled in the gaseous transport equations [9,10,11,12, 17,18,19,20,21,22,23]:

$$ \frac{\partial \rho \psi}{\partial t}+\frac{\partial \rho {u}_j\psi }{\partial {x}_j}=\frac{\partial }{\partial {x}_j}\left({\varGamma}_{\psi}\frac{\partial {\psi}_1}{\partial {x}_j}\right)+{\overset{.}{\omega}}_{\psi }+{\overset{.}{S}}_g+{\overset{.}{S}}_{\psi } $$

(1vi)

where ψ = {1, u_j, e, Y_F, Y_O} for the conservation of mass, momentum, energy and mass fractions respectively, ${\psi}_1=\left\{1,{u}_j,\hat{T},{Y}_F,{Y}_O\right\}$ for ψ = {1, u_j, e, Y_F, Y_O}, and Γ_ψ = μ/σ_ψ and λ for ψ = {u_j, Y_F, Y_O} and ψ = e respectively, with u_j, μ, λ and σ_ψ being the velocity component in the j^th direction, dynamic viscosity, thermal conductivity and an appropriate Schmidt number for ψ, respectively. The term ${\overset{.}{\omega}}_{\psi }$ arises due to chemical reaction rate and ${\overset{.}{S}}_g$ is an appropriate gaseous phase source term. The droplet source term arising from evaporation, ${\overset{.}{S}}_{\psi }=-1/V{\sum}_dd\left({m}_d{\psi}_d\right)/ dt$, is interpolated from the droplet’s sub-grid position to the 8 surrounding nodes, where V is the cell volume, ${m}_d={\rho}_d\left(1/6\right)\pi {a}_d^3$ is the droplet mass and the summation is carried out over all droplets in the vicinity of each node [9,10,11,12, 17,18,19,20,21,22,23].

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₇H₁₆, 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]: c = [(1 − ξ)Y_O∞ − Y_O]/[(1 − ξ)Y_O∞ −  max (0, [ξ_st − ξ]/ξ_st)Y_O∞] 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 Y_Ou = (1 − ξ)Y_O∞ and Y_Ob =  max (0, [ξ_st − ξ]/ξ_st)Y_O∞ (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 $\overset{\sim }{\varepsilon_Y}$ [6]:

(2i)

(2ii)

$$ {T}_2=-2D\overline{\frac{\left[\ {\overset{.}{\omega}}_F+\nabla \cdotp \left(\rho D\nabla {Y}_F\right)\right]}{\rho}\frac{\partial {Y}_F}{\partial {x}_k}\frac{\partial \rho }{\partial {x}_k}}+2\frac{\overset{\sim }{D}}{\overline{\rho}}\frac{\partial \overset{\sim }{Y_F}}{\partial {x}_k}\frac{\partial \overline{\rho}}{\partial {x}_k}\left[\overline{{\overset{.}{\omega}}_F+\nabla .\left(\rho D\nabla {Y}_F\right)}-\frac{\partial \left(\overline{\rho {u}_l^{\prime \prime }{Y}_F^{\prime \prime }}\right)}{\partial {x}_l}\right] $$

(2iii)

(2iv)

$$ {T}_4=2\overline{D\frac{\partial {\overset{.}{\omega}}_F}{\partial {x}_k}\frac{\partial {Y}_F}{\partial {x}_k}}-2\overset{\sim }{D}\frac{\partial \overline{{\overset{.}{\omega}}_F}}{\partial {x}_k}\frac{\partial \overset{\sim }{Y_F}}{\partial {x}_k} $$

(2v)

$$ {T}_5=2\overline{D\frac{\partial {\overset{.}{S}}_E}{\partial {x}_k}\frac{\partial {Y}_F}{\partial {x}_k}}-2\overset{\sim }{D}\frac{\partial {\overset{.}{S}}_{ER}}{\partial {x}_k}\frac{\partial \overset{\sim }{Y_F}}{\partial {x}_k} $$

(2vi)

$$ {T}_6=\overline{D\nabla {\mathrm{Y}}_F\cdotp \nabla {Y}_F\Gamma}-\overset{\sim }{D}\nabla \overset{\sim }{{\mathrm{Y}}_F}\cdotp \nabla \overset{\sim }{Y_F}\overset{\sim }{\Gamma} $$

(2vii)

$$ {D}_2=2\overline{\rho {D}^2\frac{\partial^2{Y}_F^{\prime \prime }}{\partial {x}_k\partial {x}_i}\frac{\partial^2{Y}_F^{\prime \prime }}{\partial {x}_k\partial {x}_i}} $$

(2viii)

$$ f(D)=\overline{2D\frac{\partial {Y}_F}{\partial {x}_k}\frac{\partial \left(\rho D\right)}{\partial {x}_k}\frac{\partial^2{Y}_F}{\partial {x}_j\partial {x}_j}}+\overline{2D\frac{\partial {Y}_F}{\partial {x}_k}\frac{\partial^2\left(\rho D\right)}{\partial {x}_j\partial {x}_k}\frac{\partial {Y}_F}{\partial {x}_j}}-\overline{\frac{\partial }{\partial {x}_j}\left(\rho {N}_Y\frac{\partial D}{\partial {x}_j}\right)}-\overline{2\rho D\frac{\partial D}{\partial {x}_j}\frac{\partial }{\partial {x}_j}\left(\frac{\partial {Y}_F}{\partial {x}_k}\frac{\partial {Y}_F}{\partial {x}_k}\right)}+\overline{\rho \left(\frac{\partial {Y}_F}{\partial {x}_k}\frac{\partial {Y}_F}{\partial {x}_k}\right)\left[\frac{\partial D}{\partial t}+{u}_j\frac{\partial D}{\partial {x}_j}\right]} $$

(2ix)

where Γ is the source term in the mass conservation equation due to evaporation, ${\overset{.}{S}}_{ER}=\overline{\Gamma}\left(1-\overset{\sim }{Y_F}\right)$ and ${\overset{.}{S}}_E=\Gamma \left(1-{Y}_F\right)$. The term T₁ is the turbulent transport term, T₂ arises due to density variation, T₃ originates due to the alignment of the scalar gradient with the fluid-dynamic strain rates and essentially signifies generation/destruction of the scalar gradient by the velocity gradients, T₄ arises due to the chemical reaction rate, T₅ and T₆ arise due to droplet evaporation, whereas the term D₂ arises due to molecular dissipation. In Eq. (2), the terms of T₁, T₂, T₃₁, T₃₂, T₃₃, T₄, T₅, T₆ and (−D₂) 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.35{D}_0/{S}_{b\left({\phi}_g=1\right)}\times 42.17{D}_0/{S}_{b\left({\phi}_g=1\right)}\times 42.17{D}_0/{S}_{b\left({\phi}_g=1\right)}$ has been considered, where D₀ and ${S}_{b\left({\phi}_g=1\right)}$ are the unburned gas diffusivity and the unstrained laminar burning velocity of the stoichiometric mixture, respectively. For the present thermo-chemistry ${D}_0/{S}_{b\left({\phi}_g=1\right)}\approx 0.625{\updelta}_{\mathrm{th}}$ where ${\delta}_{th}=\left({T}_{ad\left({\phi}_g=1\right)}-{T}_0\right)/\mathit{\max}{\left(\left|\nabla \hat{T}\right|\right)}_L$ 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\le x{S}_{b\left({\phi}_g=1\right)}/{D}_0\le 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₀ = 300K, which leads to $\tau =\left({T}_{ad\left({\phi}_g=1\right)}-{T}_0\right)/{T}_0=6.54$ where ${T}_{ad\left({\phi}_g=1\right)}$ 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.53{D}_0/{S}_{b\left({\phi}_g=1\right)}$ from the point in the laminar flame at which $\hat{\mathrm{T}}=400\mathrm{K}$. It should be noted that the droplet diameter a_d has been non-dimensionalised using the well-defined 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\le x{S}_{b\left({\phi}_g=1\right)}/{D}_0\le 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}^{\prime }/{S}_{b\left({\phi}_g=1\right)}=7.5$ and non-dimensional longitudinal integral length-scale L₁₁/δ_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}^{\prime }/{S}_{b\left({\phi}_g=1\right)}=7.5$. All simulations have been carried out until t_final = max(3t_turb, 4t_chem), where t_turb = L₁₁/u^′ and ${\mathrm{t}}_{\mathrm{chem}}={D}_0/{S}_{b\left({\upphi}_g=1\right)}^2$ 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 volume-integrated 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. $\overline{\mathrm{q}}$ and $\overset{\sim }{\mathrm{q}}$) 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 $\overset{\sim }{c}$ and the scalar dissipation rate $\overset{\sim }{\varepsilon_Y}$ and the terms of its transport equation are unique functions of the x₁-direction, 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 $\overset{\sim }{c}=f\left({x}_1\right)$.

3 Results & Discussion

3.1 Flame Behaviours

Figures 1a, b, and c present the instantaneous distributions of normalised fuel mass fraction Y_F/Y_Fst, mixture fraction ξ and non-dimensional temperature $T=\left(\hat{T}-{T}_0\right)/\left({T}_{\left( ad\left({\phi}_g=1\right)\right)}-{T}_0\right)$ 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 Da (i.e. $Da={L}_{11}{S}_{b\left({\phi}_g=1\right)}^2/\left({u}^{\prime }{D}_0\right)<1$ [9, 10].

Fig. 1

Instantaneous fields of a normalized fuel mass fraction, Y_F/Y_Fst, b normalised temperature T with c-isolines (left to right c = 0.1,0.3,0.5,0.7,0.9) in white and c mixture fraction ξ fields at the central x − z plane at t = 4.0t_chem for case a_d/δ_th = 0.08, ϕ_d = 1.0. Droplets are shown by black dots (not to scale). The stoichiometric mixture fraction is ξ_st = 0.0621

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 $D{a}_{\phi_g=1}={L}_{11}{S}_{b\left({\phi}_g=1\right)}^2/{u}^{\prime }{D}_0$ is 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, $\mathrm{G}=3\left(1+0.276{\mathit{\operatorname{Re}}}_d^{1/2}{Sc}^{1/3}\right) Le{N}^{2/3}\left({a}_d/{s}_d\right)$ (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.

3.2 Algebraic Modelling of $\overset{\sim }{\varepsilon_Y}$

The variations of the fuel mass fraction dissipation rate $\overset{\sim }{\varepsilon_Y}$ with $\overset{\sim }{c}$ across the flame-brush are shown in Figs. 2a-f. The scalar dissipation rate of fuel mass fraction, $\overset{\sim }{\varepsilon_Y}$, assumes non-zero values close to $\overset{\sim }{c}=0$, which decrease before rising and attaining a maximum roughly at the middle of the flame-brush for all cases. The magnitude of $\overset{\sim }{\varepsilon_Y}$ then decreases as $\overset{\sim }{c}$ approaches unity. A discussion of the variation of $\overset{\sim }{\varepsilon_Y}$ across the flame brush has been provided in the following sections.

Fig. 2

Variation of ${\overset{\sim }{\varepsilon}}_Y\times {D}_0/{Y}_{Fst}^2{S}_{b\left({\phi}_g=1\right)}^2$ and predictions of ${C}_{YLR}\left(\overset{\sim }{\varepsilon }/\overset{\sim }{k}\right)\overset{\sim }{Y_F^{\prime \prime 2}}$, Eqs. (4) and (5) [ , , , ] with $\overset{\sim }{c}$ for (a) a_d/δ_th = 0.06, ϕ_d = 1.0; (b) a_d/δ_th = 0.06, ϕ_d = 1.7; for (c) a_d/δ_th = 0.08, ϕ_d = 1.0; (d) a_d/δ_th = 0.08, ϕ_d = 1.7; (e) a_d/δ_th = 0.10, ϕ_d = 1.0; (f) a_d/δ_th = 0.10, ϕ_d = 1.7

The fuel mass fraction dissipation rate $\overset{\sim }{\varepsilon_Y}$ can be modelled using an algebraic expression following the linear relaxation approach as $\overset{\sim }{\varepsilon_Y}={C}_{YLR}\left(\overset{\sim }{\varepsilon }/\overset{\sim }{k}\right)\overset{\sim }{Y_F^{\prime \prime 2}}$ where C_YLR is a model parameter, $\overset{\sim }{k}$ is the turbulent kinetic energy and $\overset{\sim }{\varepsilon }$ is the dissipation rate of turbulent kinetic energy. Alternatively, Mura et al. [3] proposed algebraic expressions for $\overset{\sim }{\varepsilon_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. $\overset{\sim }{P}\left({Y}_F,\xi \right)=\rho P\left({Y}_F,\xi \right)/\overline{\rho}$) for turbulent stratified gaseous mixture combustion as:

$$ \overset{\sim }{P}\left({Y}_F,\xi \right)={\lambda}_w\overset{\sim }{P}\left(\xi |{Y}_{max}\right)\delta \left({Y}_F-{Y}_{max}\left(\xi \right)\right)+\left(1-{\lambda}_w\right)\overset{\sim }{P}\left(\xi |{Y}_{min}\right)\delta \left({Y}_F-{Y}_{min}\left(\xi \right)\right)+O\left(1/ Da\right) $$

(3)

where $\overset{\sim }{P}\left(\xi |{Y}_F\right)$ 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 $\overset{\sim }{P}\left(\xi |{Y}_{max}\right)=\overset{\sim }{P}\left(\xi |{Y}_{min}\right)$ [3] such that ${\uplambda}_{\mathrm{w}}=\left(\overset{\sim }{Y_F}-{\overset{\sim }{Y}}_{min}\right)/\left({\overset{\sim }{Y}}_{max}-{\overset{\sim }{Y}}_{min}\right)$. Based upon Eq. (3), Mura et al. [3] proposed a model for $\overset{\sim }{\varepsilon_Y}$ in the following manner:

$$ {\displaystyle \begin{array}{l}\overset{\sim }{\varepsilon_Y}=\left(-\overline{\rho D}\frac{\partial {\overset{\sim }{Y}}_F}{\partial {x}_k}\frac{\partial {\overset{\sim }{Y}}_F}{\partial {x}_k}-\frac{{\overline{\overset{.}{\omega}}}_F}{2}\left({\overset{\sim }{Y}}_{max}-{\overset{\sim }{Y}}_F+{\overset{\sim }{Y}}_{min}\right)+\left[\overline{{\overset{.}{\omega}}_F{Y}_F}-{\overline{\overset{.}{\omega}}}_F{\overset{\sim }{Y}}_F\right]\right)\frac{1}{\overline{\rho}}\\ {}+\left(\frac{{\overset{\sim }{Y}}_F-{\overset{\sim }{Y}}_{min}}{{\overset{\sim }{Y}}_{max}-{\overset{\sim }{Y}}_{min}}+\frac{{\overset{\sim }{Y}}_{max}-{\overset{\sim }{Y}}_F}{{\overset{\sim }{Y}}_{max}-{\overset{\sim }{Y}}_{min}}\overset{\sim }{A^2}\right)\times {\overset{\sim }{\epsilon}}_{\xi}\end{array}} $$

(4)

Mura et al. [3] further proposed the following model for $\overset{\sim }{\varepsilon_Y}$:

$$ {\displaystyle \begin{array}{l}\overset{\sim }{\varepsilon_Y}=S\left(-\overline{\rho D}\frac{\partial {\overset{\sim }{Y}}_F}{\partial {x}_k}\frac{\partial {\overset{\sim }{Y}}_F}{\partial {x}_k}-\frac{{\overline{\overset{.}{\omega}}}_F}{2}\left({\overset{\sim }{Y}}_{max}-{\overset{\sim }{Y}}_F+{\overset{\sim }{Y}}_{min}\right)+\left[\overline{{\overset{.}{\omega}}_F{Y}_F}-{\overline{\overset{.}{\omega}}}_F{\overset{\sim }{Y}}_F\right]\right)\frac{1}{\overline{\rho}}\\ {}+S\left(\frac{{\overset{\sim }{Y}}_F-{\overset{\sim }{Y}}_{min}}{{\overset{\sim }{Y}}_{max}-{\overset{\sim }{Y}}_{min}}+\frac{{\overset{\sim }{Y}}_{max}-{\overset{\sim }{Y}}_F}{{\overset{\sim }{Y}}_{max}-{\overset{\sim }{Y}}_{min}}\overset{\sim }{A^2}\right)\times \overset{\sim }{\varepsilon_{\xi }}+\left(1-{S}_{Mod}\right){C}_Y\frac{\overset{\sim }{\varepsilon }}{\overset{\sim }{k}}\overset{\sim }{{Y_F^{\prime \prime}}^2}\end{array}} $$

(5)

where C_Y is a model parameter and S is a segregation factor that is defined as [3]:

$$ S=\left\{\overset{\sim }{{Y_F^{\prime \prime}}^2}-\left[\left({\overset{\sim }{Y}}_F-{\overset{\sim }{Y}}_{min}\right)/\left({\overset{\sim }{Y}}_{max}-{\overset{\sim }{Y}}_{min}\right)+\left\{\left({\overset{\sim }{Y}}_{max}-{\overset{\sim }{Y}}_F\right)/\left({\overset{\sim }{Y}}_{max}-{\overset{\sim }{Y}}_{min}\right)\right\}{\overset{\sim }{A}}^2\right]\overset{\sim }{\xi^{\prime \prime 2}}\right\}/\left[\left({\overset{\sim }{Y}}_{max}-{\overset{\sim }{Y}}_F\right)\left({\overset{\sim }{Y}}_F-{\overset{\sim }{Y}}_{min}\right)\right] $$

(6)

The performances of the models given by $\overset{\sim }{\varepsilon_Y}={C}_{YLR}\left(\overset{\sim }{\varepsilon }/\overset{\sim }{k}\right)\overset{\sim }{Y_F^{\prime \prime 2}}$ (where C_YLR = 1), Eq. (4) and Eq. (5) (where C_Y = 1 and S is replaced by a modified segregation factor such that S_Mod =  max (0, S) because S according to Eq. (6) may locally assume negative values in droplet-flame 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 $\overset{\sim }{\varepsilon_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 $\overset{\sim }{\varepsilon_Y}$.

3.3 Statistical Behaviours of the Unclosed Terms of the $\overset{\sim }{\varepsilon_Y}$Transport Equation

The variations of the unclosed terms T₁, T₂, T₃, T₄, T₅ (here multiplied by 0.05 for visualisation purposes), T₆, −D₂ (here multiplied by 0.05 for visualisation purposes) and f(D) with $\overset{\sim }{c}$ across the flame-brush for all cases considered are shown in Figs. 3a-f. The turbulent transport term T₁ is shown to exhibit large negative values towards the unburned gas side of the flame-brush (i.e. $\overset{\sim }{c}\approx 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₂ 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₂ increases with increasing ϕ_d and increasing a_d/δ_th. The scalar turbulence interaction term T₃ is shown to exhibit large positive values towards the unburned gas side of the flame-brush (i.e. $\overset{\sim }{c}\approx 0$) before reducing in magnitude, but remaining mainly positive across the flame-brush. The relative importance of T₃ 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₃ 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₃₂ is likely to dominate other components of T₃ for large and moderate values of Damköhler number which can be observed from Figs. 4a-f where the contributions to the T₃ term (i.e. T₃₁, T₃₂ and T₃₃) across the flame-brush for all cases considered in the current study are shown. Figures 4a-f show that T₃₁ exhibits large positive values and T₃₃ large negative values towards the unburned gas side of the flame-brush (i.e. $\overset{\sim }{c}\approx 0$). However, in most cases considered here both T₃₁ and T₃₃ assume negligible values across the entire flame brush (i.e. $\overset{\sim }{c}>0$), with the exception of ϕ_d = 1.7 for which T₃₃ 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₃ and its components, particularly T₃₁ and T₃₂, are influenced directly by the inner products between the velocity and the scalar gradients. Therefore, the behaviour of T₃ 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₃₁, T₃₂ and T₃₃ in turbulent stratified flames under globally fuel-lean conditions [6].

Fig. 3

Variation of (T₁, T₂, T₃, T₄, T₅ × 0.05, T₆, −D₂ × 0.05 and f(D) [, , , , , , , ] with $\overset{\sim }{\mathrm{c}}$ for (a) a_d/δ_th = 0.06, ϕ_d = 1.0; (b) a_d/δ_th = 0.06, ϕ_d = 1.7; (c) a_d/δ_th = 0.08, ϕ_d = 1.0; (d) a_d/δ_th = 0.08, ϕ_d = 1.7; (e) a_d/δ_th = 0.10, ϕ_d = 1.0; (f) a_d/δ_th = 0.10, ϕ_d = 1.7. All terms are normalised by ${D}_0^2/{\rho}_0{Y}_{Fst}^2{S}_{b\left({\phi}_g=1\right)}^4$

Fig. 4

Variation of T₃, T₃₁, T₃₂ and T₃₃ [, , ,] with $\overset{\sim }{c}$ for (a) a_d/δ_th = 0.06, ϕ_d = 1.0; (b) a_d/δ_th = 0.06, ϕ_d = 1.7; (c) a_d/δ_th = 0.08, ϕ_d = 1.0; (d) a_d/δ_th = 0.08, ϕ_d = 1.7; (e) a_d/δ_th = 0.10, ϕ_d = 1.0; (f) a_d/δ_th = 0.10, ϕ_d = 1.7. All terms are normalised by ${D}_0^2/{\rho}_0{Y}_{Fst}^2{S}_{b\left({\phi}_g=1\right)}^4$

The reaction rate contribution term T₄ is shown to assume positive values across the flame-brush 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₅ 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. $\overset{\sim }{c}\approx 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₆ has been found to be positive across the flame-brush, but small in magnitude in comparison to T₅ with a peak value in the flame-brush shown to shift towards the burned gas side of the flame-brush increasing a_d/δ_th. Furthermore, the magnitudes of T₅ and T₆ 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₅ and T₆, and more noticeably of T₅, close to $\overset{\sim }{c}=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 ($\overset{\sim }{c}=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₅ and T₆ drop) as the droplets traverse the wide region of unburnt gas in which the value of $\overset{\sim }{c}$ 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}_F^s$ (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}_F^s$ more quickly for smaller droplets. Consequently, the gradient of the decrease in T₅ and T₆ for small $\overset{\sim }{c}$ is noticeably greater for smaller droplets than for larger ones. The molecular dissipation term (−D₂) 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. $\overset{\sim }{c}\approx 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₂) 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₅, T₆ and (−D₂) with increasing a_d/δ_th (ϕ_d) for a given value of ϕ_d (a_d/δ_th) is consistent with increasing values of $\overset{\sim }{\varepsilon_Y}$ under these conditions. The physical reasons discussed earlier for explaining the increasing trend of $\overset{\sim }{\varepsilon_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₅, T₆ and (−D₂).

It should be noted that the scalar dissipation rate transport of fuel mass fraction $\overset{\sim }{\varepsilon_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 $\overset{\sim }{\varepsilon_Y}$ transport in the purely gaseous and stratified flames.

3.4 Modelling of T₁

According to the scaling arguments of Swaminathan and Bray [40], the modelling of T₁ translates to the modelling of T₁₁ for large turbulent Reynolds numbers Re_t. Therefore, in the context of statistically planar flames, the turbulent transport term becomes equal to ${T}_1=-\partial \left(\overline{\rho {u}_1^{\prime \prime }{\varepsilon}_Y}\right)/\partial {x}_j$. The quantity $\left(\overline{\rho {u}_1^{\prime \prime }{\varepsilon}_Y}\right)$ is often modelled using the gradient hypothesis as $\overline{\rho {u}_1^{\prime \prime }{\varepsilon}_Y}=-\left({\mu}_t/\sigma \right)\partial \overset{\sim }{\varepsilon_Y}/\partial {x}_1$ where ${\mu}_t=0.09\overline{\rho}{\overset{\sim }{k}}^2/\overset{\sim }{\varepsilon }$ is the eddy viscosity and σ is an appropriate turbulent Schmidt number. An alternative model for $\overline{\rho {u}_1^{\prime \prime }{\varepsilon}_Y}$, capable of addressing both gradient and counter-gradient behaviour was proposed for flames with varying model parameter ψ_m in the following manner [6]:

$$ \overline{\rho {u}_1^{\prime \prime }{\varepsilon}_Y}=\left({\psi}_m-{\overset{\sim }{Y}}_F\right)\overline{\rho {u}_1^{\prime \prime }{Y}_F^{\prime \prime }}\times \overline{\rho}\overset{\sim }{\varepsilon_Y}/\overline{\rho {Y}_F^{\prime \prime 2}} $$

(7)

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 $-\left({\mu}_t/\sigma \right)\partial \overset{\sim }{\varepsilon_Y}/\partial {x}_1$ (where σ = 1) and Eq. (7) (where ψ_m = 0.25 × Y_Fst) compared to $\overline{\rho {u}_1^{\prime \prime }{\varepsilon}_Y}$ obtained from DNS is shown in Figs. 5a-f. It is evident from Fig. 5a-f that $-\left({\mu}_t/\sigma \right)\partial \overset{\sim }{\varepsilon_Y}/\partial {x}_1$ captures neither the qualitative nor the quantitative behaviour of $\overline{\rho {u}_1^{\prime \prime }{\varepsilon}_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 $\overline{\rho {u}_1^{\prime \prime }{\varepsilon}_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 $-\left({\mu}_t/\sigma \right)\partial \overset{\sim }{\varepsilon_Y}/\partial {x}_1$.

Fig. 5

Variation of $\overline{\rho {u}_1^{\hbox{'}\hbox{'}}{\varepsilon}_Y}$ and predictions of $-\left({\mu}_t/\sigma \right)\partial \overset{\sim }{\varepsilon_Y}/\partial {x}_1$ and Eq. (7) [ , , ] with $\overset{\sim }{c}$ for (a) a_d/δ_th = 0.06, ϕ_d = 1.0; (b) a_d/δ_th = 0.06, ϕ_d = 1.7; (c) a_d/δ_th = 0.08, ϕ_d = 1.0; (d) a_d/δ_th = 0.08, ϕ_d = 1.7; (e) a_d/δ_th = 0.10, ϕ_d = 1.0; (f) a_d/δ_th = 0.10, ϕ_d = 1.7. All terms are normalised by ${D}_0/{\rho}_0{Y}_{Fst}^2{S_b^3}_{\left({\phi}_g=1\right)}$

3.5 Modelling of T₂

It has been observed that the density variation term T₂ remains positive throughout the flame-brush in all cases. Furthermore, the term T₂ 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₂ for perfectly premixed flames, which is consistent with the scaling arguments of Swaminathan and Bray [40], in the following manner:

$$ {T}_2=2{\alpha}_{T_2}f\left({Ka}_L\right){\overline{S}}_b\overline{\rho\ }\overset{\sim }{\varepsilon_Y}/{\delta}_b $$

(8)

where ${\alpha}_{T_2}\approx A\left({\rho}_0/{\overline{\rho}}_b-1\right)$ and A is a model parameter which is of the order of unity, but is dependent upon thermochemistry. In Eq. (8), ${\overline{\rho}}_b$ is the burned gas density and ${\overline{S}}_b$ is the characteristic burning velocity, which are evaluated as [46]:

$$ {\overline{\rho}}_b={\int}_{\phi_{min}}^{\phi_{max}}{\rho}_{b\left(\phi \right)}\mathrm{p}\left(\phi \right)\mathrm{d}\phi\ \mathrm{and}\ {\overline{S}}_b={\int}_{\phi_{min}}^{\phi_{max}}{S}_{b\left(\phi \right)}\mathrm{p}\left(\phi \right)\mathrm{d}\phi $$

(9)

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\approx {\left({\overline{S}}_b\right)}^{-3/2}{\left(\overset{\sim }{\varepsilon }{\delta}_b\right)}^{1/2}$ is the local Karlovitz number and ${\delta}_b=2{D}_0/{\overline{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₂ with increasing droplet diameter. Accordingly, in the current study, a reasonable agreement has been found when $A=15{\left[0.01+{\left({a}_d/{\delta}_{th}\right)}_0\right]}^{0.25}/{\left[{\overline{S}}_b/{S}_{b\left({\phi}_g=1\right)}\right]}^{1.5}$ (where (a_d/δ_th)₀ is the initial value of 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₂ obtained from the DNS data.

Fig. 6

Variation of T₂ and predictions of Eq. (8) [ , ] with $\overset{\sim }{c}$ for (a) a_d/δ_th = 0.06, ϕ_d = 1.0; (b) a_d/δ_th = 0.06, ϕ_d = 1.7; (c) a_d/δ_th = 0.08, ϕ_d = 1.0; (d) a_d/δ_th = 0.08, ϕ_d = 1.7; (e) a_d/δ_th = 0.10, ϕ_d = 1.0; (f) a_d/δ_th = 0.10, ϕ_d = 1.7. All terms are normalised by ${D}_0^2/{\rho}_0{Y}_{Fst}^2{S}_{b\left({\phi}_g=1\right)}^4$

3.6 Modelling of T₃₁

To close the behaviour of the scalar-turbulence interaction term T₃, the components T₃₁, T₃₂ and T₃₃ need to be modelled. Although for all cases considered here the components T₃₁ and T₃₃ are of negligible magnitude in comparison to the magnitude of component T₃₂, this does not necessarily hold true for other choices of parameter values. Furthermore, since there are already existing models for T₃₁ and T₃₃, it would be remiss to omit an assessment of these models.

Several models have previously been proposed for T₃₁. For turbulent premixed flames, the following model was proposed [47]:

$$ {T}_{31}=-{C}_{P1}\overline{\rho}\left(\overset{\sim }{\varepsilon }/\overset{\sim }{k}\right)\overset{\sim }{u_j^{\prime \prime }{Y}_F^{\prime \prime }}\left(\partial {\overset{\sim }{Y}}_F/\partial {x}_j\right) $$

(10)

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:

$$ {T}_{31}=-{C}_{PM}\overline{\rho}\left({\overset{\sim }{\varepsilon}}_Y/\overset{\sim }{Y_F^{\prime \prime 2}}\right)\overset{\sim }{u_j^{\prime \prime }{Y}_F^{\prime \prime }}\left(\partial {\overset{\sim }{Y}}_F/\partial {x}_j\right) $$

(11)

$$ {T}_{31}=\left({\rho}_0/{\overline{\rho}}_b-1\right){\overline{S}}_b\overline{\rho}{\overset{\sim }{\varepsilon}}_Y\left\langle {\overrightarrow{n}}_F\cdotp {\overrightarrow{x}}_j\right\rangle \left(\partial {\overset{\sim }{Y}}_F/\partial {x}_j\right) $$

(12)

where C_PM is a model constant and ${\overrightarrow{n}}_F=-\nabla {\mathrm{Y}}_F/\left|\nabla {Y}_F\right|$. Alternatively, Chakraborty and Swaminathan [49] proposed the following model to account for low Damköhler number premixed combustion:

$$ {T}_{31}=-\left[{C}_1+{C}_2{Da}_L^{\ast}\right]\overline{\rho}\left(\overset{\sim }{\varepsilon }/\overset{\sim }{k}\right)\overset{\sim }{u_j^{\prime \prime }{Y}_F^{\prime \prime }}\left(\partial {\overset{\sim }{Y}}_F/\partial {x}_j\right) $$

(13)

where ${Da}_L^{\ast }={\overline{S}}_b{\rho}_0\overset{\sim }{k}/\left({\delta}_b\overline{\rho}\overset{\sim }{\varepsilon}\right)$ is a local density-weighted Damköhler number and C₁ = 0.5 and ${C}_2=1.3{Ka}_L^2/{\left(1+{Ka}_L\right)}^{-2}$ are model parameters [49]. Malkeson and Chakraborty [6] extended this model for stratified combustion as:

$$ {T}_{31}=-\left[{C}_1^{\ast }+{C}_2^{\ast }{Da}_L^{\ast}\right]\overline{\rho}\left(\overset{\sim }{\varepsilon }/\overset{\sim }{k}\right)\overset{\sim }{u_j^{\prime \prime }{Y}_F^{\prime \prime }}\left(\partial {\overset{\sim }{Y}}_F/\partial {x}_j\right)\left({\psi}^{\ast }-{N}_B\right) $$

(14)

where ${N}_B=\left({\rho}_0/{\overline{\rho}}_b-1\right){\overline{S}}_b/\sqrt{2\overset{\sim }{k}/3}$ is the Bray number, ${C}_1^{\ast }=0.7$, ${C}_2^{\ast }=4.5\left[1-\mathit{\operatorname{erf}}\left({Ka}_L{\left({\overline{S}}_b/{S}_{b\left({\phi}_g=1\right)}\right)}^2\right)/2.5\right]$ 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₃₁ with $\overset{\sim }{c}$ 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₃₁ 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₃₁ for the cases considered here with the best performance being shown by Eq. (13).

Fig. 7

Variation of T₃₁ and predictions of Eqs. (10)–(14) [, , , , , ] with $\overset{\sim }{c}$ for (a) a_d/δ_th = 0.06, ϕ_d = 1.0; (b) a_d/δ_th = 0.06, ϕ_d = 1.7; (c) a_d/δ_th = 0.08, ϕ_d = 1.0; (d) a_d/δ_th = 0.08, ϕ_d = 1.7; (e) a_d/δ_th = 0.10, ϕ_d = 1.0; (f) a_d/δ_th = 0.10, ϕ_d = 1.7. All terms are normalised by ${D}_0^2/{\rho}_0{Y}_{Fst}^2{S}_{b\left({\phi}_g=1\right)}^4$

3.7 Modelling of T₃₂

A number of models have also been proposed for the contribution T₃₂. For turbulent premixed flames, the following model was proposed by Mura and Borghi [51] which modified the model proposed by Mantel and Borghi [47]:

$$ {T}_{32}={A}_e\overline{\rho}\left(\overset{\sim }{\varepsilon }/\overset{\sim }{k}\right)\overset{\sim }{\varepsilon_Y} $$

(15)

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:

$$ {T}_{32}=\overline{\rho}\left[{A}_{M1}\left(\overset{\sim }{\varepsilon }/\overset{\sim }{k}\right)-2{C}_{MA}\left({\rho}_0/{\overline{\rho}}_b-1\right){Da}_L\overset{\sim }{\varepsilon_Y}\right]\overset{\sim }{\varepsilon_Y} $$

(16)

$$ {T}_{32}=\overline{\rho}\left[{A}_{M2}\left(\overset{\sim }{\varepsilon }/\overset{\sim }{k}\right)-2{C}_{MB}\ln \left({\rho}_0/{\overline{\rho}}_b\right){Da}_L\overset{\sim }{\varepsilon_Y}\right]\overset{\sim }{\varepsilon_Y} $$

(17)

where A_M1 = 1, A_M2 = 1, C_MA = 0.6 and C_MB = 1.6 are model parameters [52] and ${Da}_L={\overline{S}}_b\overset{\sim }{k}/\left({\delta}_b\overset{\sim }{\varepsilon}\right)$ is a local Damköhler number.

Furthermore, Chakraborty and Swaminathan [49] proposed the following model to account for low Damköhler number premixed combustion:

$$ {T}_{32}=\overline{\rho}\left(\overset{\sim }{\varepsilon }/\overset{\sim }{k}\right)\left[{C}_3-{C}_4\left({\rho}_0/{\overline{\rho}}_b-1\right){Da}_L^{\ast}\right]\overset{\sim }{\varepsilon_Y} $$

(18)

where C₃ = 1.5 and C₄ = 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 $\overset{\sim }{\varepsilon_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,44,45, 47, 49, 52], and thus are not repeated here.

The variations of T₃₂ with $\overset{\sim }{c}$ 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₃₂ 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₃₂ for the considered cases.

Fig. 8

Variation of T₃₂ and predictions of Eqs. (15)–(18) [, , , , ] with $\overset{\sim }{c}$ for (a) a_d/δ_th = 0.06, ϕ_d = 1.0; (b) a_d/δ_th = 0.06, ϕ_d = 1.7; (c) a_d/δ_th = 0.08, ϕ_d = 1.0; (d) a_d/δ_th = 0.08, ϕ_d = 1.7; (e) a_d/δ_th = 0.10, ϕ_d = 1.0; (f) a_d/δ_th = 0.10, ϕ_d = 1.7. All terms are normalised by ${D}_0^2/{\rho}_0{Y}_{Fst}^2{S}_{b\left({\phi}_g=1\right)}^4$

3.8 Modelling of T₃₃

Several models have also been proposed for the contribution T₃₃. For turbulent premixed flames, the following model was proposed [47]:

$$ {T}_{33}=-{C}_{P2}\overline{\rho}\tilde{\varepsilon_Y}\left(\tilde{u_{\mathrm{J}}^{{\prime\prime}}\kern0.33em {u}_k^{{\prime\prime} }}/\tilde{k}\right)\partial {\tilde{u}}_j/\partial {x}_k $$

(19)

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:

$$ {T}_{33}=-\left(2/3\right)\overline{\rho}\tilde{\varepsilon_Y}\partial {\tilde{u}}_j/\partial {x}_j $$

(20)

$$ {T}_{33}=-2\overline{\rho}{\tilde{\varepsilon_Y}}^2\left[\left(\overline{\rho {u}_{\mathrm{J}}^{{\prime\prime} }{Y}_F^{{\prime\prime} }}\cdot \overline{\rho {u}_k^{{\prime\prime} }{Y}_F^{{\prime\prime} }}\right)/\left(\tilde{\varepsilon}\overline{\rho {Y_F^{{\prime\prime}}}^2}\cdot \overline{\rho {Y}_F^{{\prime\prime} 2}}\right)\right]\partial {\tilde{u}}_j/\partial {x}_k $$

(21)

$$ {T}_{33}=-2\overline{\rho}{\tilde{\varepsilon_Y}}^2{\left({\rho}_0/{\overline{\rho}}_b-1\right)}^2{{\overline{S}}_b}^2\left(\left\langle {\overrightarrow{n}}_f\cdot {\overrightarrow{x}}_j\right\rangle \left\langle {\overrightarrow{n}}_f\cdot {\overrightarrow{x}}_k\right\rangle /\tilde{\varepsilon}\right)\partial {\tilde{u}}_j/\partial {x}_k $$

(22)

where A_M1 = 1, A_M2 = 1, C_MA = 0.6 and C_MB = 1.6 are model parameters [48] and ${Da}_L={\overline{S}}_b\overset{\sim }{k}/\left({\delta}_b\overset{\sim }{\varepsilon}\right)$ is a local Damköhler number. Furthermore, Chakraborty and Swaminathan [49] proposed the following model to account for low Damköhler number premixed combustion:

$$ {T}_{33}=-{C}_{T3}\overline{\rho}\tilde{\varepsilon_Y}\left[{n}_i{n}_j+{\delta}_{ij}\left(1-{n}_k{n}_k\right)/3\right]\partial {\tilde{u}}_i/\partial {x}_j $$

(23)

where ${C}_{T3}=\left(1+a\ {Ka}_L^{-0.23}\right)$ is a model parameter and ${n}_i=\left(-\partial \sqrt{\overset{\sim }{Y_F^{\prime \prime 2}}}/\partial {x}_i\right).\sqrt{\rho_0{D}_0/\left(\overline{\rho}\overset{\sim }{\varepsilon_Y}\right)}$. 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₃₃ with $\overset{\sim }{c}$ 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₃₃ 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₃₃ for all considered cases with the best quantitative agreement being found for Eqs. (20) and (23).

Fig. 9

Variation of T₃₃ and predictions of Eqs. (19)–(23) [, , , , ,] with $\overset{\sim }{c}$ for (a) a_d/δ_th = 0.06, ϕ_d = 1.0; (b) a_d/δ_th = 0.06, ϕ_d = 1.7; (c) a_d/δ_th = 0.08, ϕ_d = 1.0; (d) a_d/δ_th = 0.08, ϕ_d = 1.7; (e) a_d/δ_th = 0.10, ϕ_d = 1.0; (f) a_d/δ_th = 0.10, ϕ_d = 1.7. All terms are normalised by ${D}_0^2/{\rho}_0{Y}_{Fst}^2{S}_{b\left({\phi}_g=1\right)}^4$

3.9 Modelling of T₄ − D₂ + f(D)

The combined contribution of the terms D₁, T₄, (−D₂) and f(D) can be written as [44, 45]:

$$ {\displaystyle \begin{array}{c}{D}_1+{T}_4-{D}_2+f(D)\approx \overline{-}-\overline{2D\nabla \cdot \left(\rho {S}_d\overrightarrow{n}|\nabla {Y}_F|\right)\mid \nabla {Y}_F\mid }+2\tilde{D}\nabla \cdot \left(\overline{\rho {S}_d\mid \nabla {Y}_F\mid}\overrightarrow{m}\right)\mid \overline{\nabla {Y}_F}\mid \\ {}+2\overline{\rho D{S}_{LS}\nabla \cdot \overrightarrow{n}{\left|\nabla {Y}_F\right|}^2}-2\tilde{D}\overline{\rho {S}_{LS}}\nabla \cdot \overrightarrow{m}\mid \overline{\nabla {Y}_F}\mid \\ {}-2\overline{\rho {D}^2{\left(\nabla \cdot \overrightarrow{n}\right)}^2{\left|\nabla {Y}_F\right|}^2}+2\tilde{D}\;\overline{\rho D\left(\nabla \cdot \overrightarrow{n}\right)\mid \nabla {Y}_F\mid}\nabla \cdot \overrightarrow{m}\mid \overline{\nabla {Y}_F}\mid \end{array}} $$

(24)

where $\overrightarrow{m}=\nabla \overset{\sim }{Y_F}/\left|\overline{\nabla {\mathrm{Y}}_{\mathrm{F}}}\right|$ is the resolved flame normal, $\overrightarrow{n}=\nabla {Y}_F/\mid \nabla {Y}_F\mid$ is the local flamelet normal, ${S}_d={S}_{LS}-D\nabla .\overrightarrow{n}=-{\left|\nabla {Y}_F\right|}^{-1}\left(D{Y}_F/ Dt\right)$ is the local displacement speed. Equation (24) indicates that the net contribution of [D₁ + T₄ − D₂ + f(D)] signifies the effect due to flame normal propagation and flame curvature. As the molecular diffusion term D₁ 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₄ − D₂ + f(D)] rather than its individual components [6, 43,44,45, 47]. It has been shown elsewhere [44, 45] the net contribution of [T₄ − D₂ + f(D)] remains negative because of dominant contribution of the term $-\left\{2\overline{\rho {D}^2{\left(\nabla \cdotp \overrightarrow{n}\right)}^2{\left|\nabla {Y}_F\right|}^2}-2\overset{\sim }{D}\overline{\rho D\left(\nabla \cdotp \overrightarrow{n}\right)\left|\nabla {Y}_F\right|}\nabla \cdotp \overrightarrow{m}\left|\overline{\nabla {Y}_F}\right|\right\}$ on the right hand of Eq. (24). It has been observed that the net contribution of T₄ − D₂ + 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₄ − D₂ + f(D)) is dominated by the contribution of (−D₂), which is determined by the gradient of the fuel mass fraction (see Eq. (2viii)). Here, T₄ − D₂ + f(D) is modelled in the following manner according to previous studies [6, 43,44,45, 47]:

$$ \left({T}_4-{D}_2+f(D)\right)=-{\beta}_{2Y}\overline{\rho}{\tilde{\varepsilon}}_Y^2\tilde{{Y_F^{{\prime\prime}}}^2} $$

(25)

where β_2Y is a model parameter. The magnitude of (T₄ − D₂ + 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 ${\beta}_{2Y}=5\left[1+2{\left[{\left({a}_d/{\delta}_{th}\right)}_0\right]}^{0.1}{\left[{\overline{S}}_b/{S}_{b\left({\phi}_g=1\right)}\right]}^{-2}\right]$as shown in Figs. 10a-f for all cases considered in the current study.

Fig. 10

Variation of [T₄ − D₂ + f(D)] and predictions of Eq. (25) [, ] with $\overset{\sim }{c}$ for (a) a_d/δ_th = 0.06, ϕ_d = 1.0; (b) a_d/δ_th = 0.06, ϕ_d = 1.7; (c) a_d/δ_th = 0.08, ϕ_d = 1.0; (d) a_d/δ_th = 0.08, ϕ_d = 1.7; (e) a_d/δ_th = 0.10, ϕ_d = 1.0; (f) a_d/δ_th = 0.10, ϕ_d = 1.7. All terms are normalised by ${D}_0^2/{\rho}_0{Y}_{Fst}^2{S}_{b\left({\phi}_g=1\right)}^4$

3.10 Modelling of (T₅ + T₆)

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₅ and T₆). The contribution of (T₅ + T₆) 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 ${\overset{\sim }{\varepsilon}}_{\xi }$, as the micro-mixing of mixture fraction variation induced by droplet evaporation is expected to affect the statistical behaviour of the evaporation contribution (T₅ + T₆) in the gaseous phase. It is worth considering that the magnitude of (T₅ + T₆) has qualitative and quantitative similarities to the magnitude of (T₄ − D₂ + f(D)), a similar form of model expression has been considered. It has been observed that the behaviour of the combined contribution of (T₅ + T₆) can be captured by the following expression:

$$ \left({T}_5+{T}_6\right)={\beta}_Y\overline{\rho}\;{\tilde{\varepsilon}}_{\xi}^2/\overset{\sim }{\xi^{\prime \prime 2}} $$

(26)

where β_Y is a model parameter which provides reasonable agreement when ${\beta}_Y=7.5\left[1+5{\left[{\left({a}_d/{\delta}_{th}\right)}_0\right]}^{0.1}{\left[{\overline{S}}_b/{S}_{b\left({\phi}_g=1\right)}\right]}^{-1.5}\right]$. 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.

Fig. 11

Variation of [T₅ + T₆] and predictions of Eq. (26) [ , ] with $\overset{\sim }{\mathrm{c}}$ for (a) a_d/δ_th = 0.06, ϕ_d = 1.0; (b) a_d/δ_th = 0.06, ϕ_d = 1.7; (c) a_d/δ_th = 0.08, ϕ_d = 1.0; (d) a_d/δ_th = 0.08, ϕ_d = 1.7; (e) a_d/δ_th = 0.10, ϕ_d = 1.0; (f) a_d/δ_th = 0.10, ϕ_d = 1.7. All terms are normalised by ${D}_0^2/{\rho}_0{Y}_{Fst}^2{S}_{b\left({\phi}_g=1\right)}^4$

3.11 Future Considerations

The current study considers a-priori DNS modelling of the fuel mass fraction dissipation rate $\overset{\sim }{\upvarepsilon_Y}$ and the unclosed terms of its transport equation following a similar approach to several previous analyses [5, 6, 33, 39,40,41,42,43,44,45, 48,49,50, 54,55,56,57,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 a-posteriori 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. $\overset{\sim }{k}$ and $\overset{\sim }{\varepsilon }$). 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 a-priori 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 simple and detailed chemical mechanisms. Nonetheless, future research in these directions will be necessary for a comprehensive assessment of the model performances.

4 Conclusions

The statistical behaviours of $\overset{\sim }{\upvarepsilon_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 $\overset{\sim }{P}\left({Y}_F,\xi \right)$ which was originally intended for high Damköhler number gaseous phase combustion does not adequately predict $\overset{\sim }{\upvarepsilon_Y}$ obtained from DNS data. The behaviours of the unclosed terms of $\overset{\sim }{\upvarepsilon_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₂, T₅ and −D₂) play significant roles in $\overset{\sim }{\upvarepsilon_Y}$ transport. The suitability of the models previously proposed in the context of turbulent gaseous stratified flames have been assessed for the modelling of $\overset{\sim }{\upvarepsilon_Y}$ transport in turbulent spray flames. Based on a-priori DNS analysis suitable model expressions have been identified for T₁, T₂, T₃₁, T₃₂, T₃₃, [T₄ − D₂ + f(D)] and [T₅ + T₆], which have been shown to perform generally satisfactorily for all cases considered here. Further consideration of the modelling of $\overset{\sim }{\upvarepsilon_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.