Modelling of Spray Flames with Doubly Conditional Moment Closure

Simulations of a pilot-stabilised flame in a uniformly dispersed ethanol spray are performed using a Doubly Conditional Moment Closure (DCMC) model. The DCMC equation for spray combustion is derived, using the mixture fraction and the reaction progress variable as conditioning variables, including droplet evaporation and differential diffusion terms. A set of closure sub-models is suggested to allow for a first, preliminary application of the DCMC model to the test case presented here. In particular, the DCMC model is used to provide complete closure for the Favre-averaged spray terms in the mean and variance equations of the conditioning variables and the present test case is used to assess the importance of each term. Comparison with experimental data shows a promising overall agreement, whilst differences are related to modelling choices.


Introduction
In a broad range of combustions devices, including most mobile applications such as aeroengines and IC-engines, fuel is supplied in liquid form. Modelling of spray flames is challenging, even when ignoring the difficulties related to dense sprays or the modelling of atomisation. In particular, the complex interplay of droplet evaporation, turbulent mixing and chemical reaction in the presence of large mixture inhomogeneities leads to a wide variety of spray combustion regimes and phenomena [1].
Many experimental and numerical studies have explored the propagation of flames through disperse sprays. Whilst flame propagation in a mist of very small droplets was similar to the case of a homogeneous mixture, larger droplets were found to have a positive effect on the burning velocity [2,3]. Yet, an inverse correlation of burning velocity with Sauter mean diameter was found above a certain droplet size [4]. In a numerical study of flame propagation in quiescent sprays, Neophytou and Mastorakos [5] showed that the effective equivalence ratio, as compared to the overall equivalence ratio, was an important parameter with respect to the burning velocity. Direct numerical simulations (DNS) showed that the flame propagation consisted of the successive ignition of flames engulfing individual droplets [6] and that premixed and non-premixed combustion modes co-exist in spray flames [7][8][9].
In comparison to laminar flame propagating in a droplet mist, turbulent spray flames have been studied less. In turbulent combustion, there is an additional need of modelling to take turbulent-flame interaction into account appropriately. A short overview of the modelling approaches for turbulent spray flames can be found in the review by Jenny et al. [1]. In order to capture finite kinetic effects, pyrolysis and pollutant formation, an advanced combustion model is required.
Conditional Moment Closure (CMC) is a statistical model for turbulent combustion, making use of a strong correlation between the reacting scalars and the conditioning variables. It then provides a very simple closure for the highly non-linear reaction source term, whilst its derivation only requires very light assumptions on the physics involved, such that its application is a priori not limited to a specific combustion mode or regime [10,11]. CMC was originally developed for non-premixed flames [10] but it has since been extended to spray combustion [12,13] and its applicability to premixed flames has also been demonstrated [14,15]. In all these models a single conditioning variable, either mixture fraction of reaction progress variable, is used. However, this was found insufficient in certain flames, e.g. in the presence of significant extinction, and Kronenburg [16] demonstrated the effectiveness of introducing a second conditioning variable to provide more accurate CMC closure of the reaction source term.
This wide range of successful applications makes CMC an attractive candidate for the modelling of flames that have both premixed and non-premixed features, in particular, by pursuing the strategy of double conditioning. The present work focuses on the development of a Doubly Conditional Moment Closure (DCMC) model for spray combustion. The objectives of this paper are to (i) present the DCMC equation for spray flames, (ii) provide closure by suggesting a set of sub-models, in order to (iii) apply the model in a first, preliminary test to a pilot stabilised flame in a uniformly dispersed ethanol spray and (iv) compare the simulation results with experimental data of this burner recently studied at the University of Cambridge [17]. This flame behaves as premixed due to pre-vaporisation but with the flame still generating its own vapour and hence provides a relevant test case.

Derivation of the DCMC equation
In the Conditional Moment Closure combustion model, transport equations are solved for the conditional moments Q α (Z c ; x, t) ≡ Y α (x, t)|Y c (x, t) = Z c of the reactive scalars Y α , which are conditionally averaged on a subset of the space of all reactive scalars Y c ⊂ Y [10]. The probability density function (pdf) of the conditioning variables is presumed and mean quantities of the reactive scalars are calculated by integrating the conditional moments with the pdf. In contrast to conventional CMC, in the present DCMC approach the reactive scalars are conditioned on two conditioning variables.
The DCMC model equation for spray flames is derived following the approach by Mortensen and Bilger [12] who used a separated flow model to incorporated the effects of spray evaporation. With a separated flow model the local instantaneous balance equations for a multi-phase flow are written by introducing a phase indicator function θ k (x, t), which is unity in the region occupied by phase k and zero everywhere else [18,19]; its governing equation is as follows.
Here Π k is the volumetric rate of phase change per unite volume. Then the continuity equation and the species transport for the phase k take the following form [12], In this work the subscript k may only refer either to the gas or to the liquid phase and is, therefore, omitted in further derivations. The present DCMC approach uses the mixture fraction ξ and the reaction progress variable c as conditioning variables. Whilst the mixture fraction is a passive scalar with respect to chemical reaction, it is produced through droplet evaporation. It is defined to be 0 in air and 1 in undiluted fuel vapour. For a pure liquid fuel the jump condition at the phase interface (5) for the mixture fraction gives ξV ξ = (1 − ξ)Π [12], which leads to the following transport equation of the mixture fraction, The progress variable is defined as a linear expression of one reactive scalar Y ψ , normalised by its un-reacted and equilibrium values, Y 0 ψ (ξ ) and Y Eq ψ (ξ ) respectively [20,21].
If Y ψ is governed by an equation of the same form as Eq. 3, then the instantaneous transport equation of the reaction progress variable c = c ψ is, where N ξ = D ξ ∇ξ · ∇ξ , N c = D c ∇c · ∇c and N ξc = D∇ξ · ∇c are the scalar dissipation rates. The first and second line of Eq. 8 are identical with the original result by Bray et al. [21], which was extended in the present work to include the spray source terms. Note that in the derivation of Eq. 8 it is assumed that D c = D ψ = D ξ . For the sake of brevity, the entire term in the second line of Eq. 8 is written as θρω * c . The presence of an evaporation source in the c-equation was discussed by Domingo et al. [7]. Using the jump conditions for Y ψ and ξ , the first term in the third line can be expressed as the product of the evaporation mass source ρΠ and a function, which can be evaluated for any ξ and c, since Y ψ (ξ, c) is known from Eq. 7. In this expression δ ψ is the mass fraction of species ψ in the liquid phase, which is equals 1 if c ψ was based on the fuel and 0 otherwise. By grouping the terms in the Eq. 8, the c-equation finally appears in a a very similar form compared to Eq. 3 with one term representing the apparent reaction rate and two separate source terms due to evaporation. The DCMC equation is the transport equation of the doubly conditional moments Q α ≡ Y α (x, t)|ξ(x, t) = η, c(x, t) = ζ . It is derived using the joint-pdf method assuming moderately high Reynolds number and invoking the primary closure hypothesis for the diffusive fluxes in conditional space [10]. This leads to, where ·|η, ζ denotes density-weighted conditional averaging and p is the density-weighted pdf, defined asρ p(η) = ρ|η, ζ p. Moreover,θ is the volume fraction of the gaseous phase, which is ≈ 1 in a dilute spray, and δ α signifies the species mass fraction in the liquid phase, which is 1 for the liquid fuel species and 0 otherwise. A similar transport equation can be derived for the conditional enthalpy Q h . This equation is presented in the Appendix and for completeness, the terms of molecular transport in the Q h -equation are also shown. The analogue term is neglected in Eq. 10 due to the high Reynolds number assumption.
In Eq. 10,ω * c represents the apparent reaction progress variable source term in a flow with non-uniform mixture fraction. Since Y ψ = Y ψ (ξ, c), according to the definition in Eq. 7, Y ψ = 0. Thus, the conditional apparent reaction rate is given as, Lewis number effects have been argued to be important for CMC of premixed flames [22]. Replacing the usual unity Lewis number assumption with the assumption of constant Le α for each species, adds a factor to the scalar dissipation rate terms in Eq. 10 and also adds a differential diffusion term, whose influence on major species might be small however [22]. We recall that in the derivation of the c-equation (8), it was assumed that Le c = Le ξ = Le ψ . This constraint can only be approximately valid in the general case of non-unity Lewis number and limits the choice of progress variable to major species with a Lewis number similar to Le ξ . If this limitation is respected and Le ψ is sufficiently close to Le ξ , the error due this assumption will be small compared to the effect of non-unity Le α which for minor species can be substantially different from Le ξ . Considering the effects of liquid fuel evaporation adds several terms to the DCMC equation (10), including doubly conditional terms equivalent to the terms first derived by Mortensen and Bilger [12]. Additional new terms appear to account for the effects of evaporation on the reaction progress variable. These terms contain the function C (η, ζ ), which is given by Eq. 9.
In the present derivation, the primary closure hypothesis is specifically applied to fluxes of the reactive scalar in conditional space only. This is in line with the original derivation [23], where it is clearly specified that closure is achieved by assuming a "first-order diffusion relation" for these fluxes. Source terms from evaporation or the reaction progress variable are not treated through this approximation, which leads to the appearance of terms containing the conditional correlations Y α Π |η, ζ and Y αω * c |η, ζ . In this way, Mortensen and Bilger [12] derived the conditional correlation term of mass fraction and evaporation rate using the joint-pdf method, whilst it seems that this term does not appear when the decomposition method is used. Indeed, the conditional correlation term of mass fraction and reaction progress variable source does not occur in the CMC equation for premixed flames derived by Mantel and Bilger [24] using the decomposition method and it has not been addressed in more recent work [22,25] either. The effect of the conditional correlation term of the evaporation rate has not been studied yet and it has been neglected so far [13,26]. Equation 10 represents the unclosed DCMC equation in its general form, whose derivation only requires very light modelling assumptions, viz. the primary closure hypothesis, moderate to high Reynolds number and constant Lewis numbers. Whilst the accuracy of first order closure for the doubly conditional reaction source term has been demonstrated [16], there is, in general, very little experience with the sub-models for DCMC. In some cases, however, it might be possible to generalise models used in conventional CMC or to adapt them from other combustion models that use a similar parametrisation, such as for instance mixture fraction-progress variable flamelet models.

Closure for the DCMC equation
In order to apply the DCMC model to a preliminary test case, in the present work the modelling choices detailed below are made. The selection of sub-models was partly driven by their extremely limited availability. In this sense, the set of models proposed is a first suggestion and more work will be necessary in the future to improve and validate submodels for DCMC.
Unity Lewis number is assumed for simplicity in this first application of the model. For the transport in physical space the well-tested sub-models for conventional CMC can be easily adopted. The diffusion approximation is used for u Y α |η, ζ [10] and the conditional velocity is modelled as u|η, ζ = u, since gradients of ξ are small. For the Favre pdfs of the conditioning variables, β-pdfs are presumed and, assuming statistical independence, the joint-pdf is computed as p(η, ζ ) = p β (η; ξ, ξ 2 )p β (ζ ; c, c 2 ). This assumption is made for simplicity since accurate modelling of the joint-pdf in the present three-stream mixing problem plus droplet evaporation would be very complex and is not attempted here. The choice of sub-models for the doubly conditional scalar dissipation rates is very limited. In the present work, we follow the suggestions by Nguyen et al. [27]. For the scalar dissipation rate of the mixture fraction they assumed that N ξ is primarily imposed by the mixing with little dependence on chemistry and, thus, N ξ |η, ζ ≈ N ξ |η , which is modelled as a bell curve given by the Amplitude Mapping Closure (AMC) model [28]. Again following Nguyen et al. [27], N c |η, ζ is modelled as the product of two bell curves, b(η) centred on the stoichiometric mixture fraction η = ξ st and G(ζ ), centred on ζ = 0.5; for η ≥ 2ξ st it is zero. This model is a simple approximation of N c -values from premixed flamelet tabulation. The conditional cross-scalar dissipation rate is not considered, in line with the assumption of statistically independent conditioning variables in the modelling of the pdf. A dilute spray is assumed, i.e.θ ≈ 1 but the doubly conditional evaporation source terms are not included in the DCMC equation. Finally radiation and wall heat losses are neglect. Together with the unity Lewis number assumption and in the absence spray source terms, the transport equation of the conditionally averaged enthalpy Q h becomes trivial and does not need to be calculated. The conditional temperature Q T is then calculated from Q h and the conditional species mass fractions Q α .

Flow field solver
The DCMC combustion model is coupled to the flow field solver. An Euler-Lagrangian approach is followed, where the gaseous phase is computed through an unsteady Reynoldsaveraged Navier-Stokes (RANS) simulation with the standard k-ε turbulence model and the parcels of liquid droplets are tracked as Lagrangian particles. The flow field solver integrates the momentum equation as well as the Favre-mean and variance equations of the conditioning variables, ξ and c. Note that (·) denotes the fluctuation around the Favre mean, i.e. Y = Y + Y , to distinguish it from the conditional fluctuation in the DCMC equation.
Here, the reaction progress variable is based on the mass fraction of carbon dioxide, The eddy diffusivity is calculated as D T = μ T /(ρSc T ), with a constant turbulent Schmidt number Sc T = 0.7; for the molecular diffusivity D = μ/(ρSc), with Sc = 0.7. In the case of the mixture fraction, the scalar dissipation rate is modelled as for a passive scalar, ( ε/ k) ξ 2 [29]. The scalar dissipation rate of c is closed using the model by Kolla et al. [30], where the model coefficients are set as described by Kolla and Swaminathan [31], notably β = 6.7 and K * c = 0.85τ . To account for non-uniform mixture τ , S 0 L and δ 0 L are evaluated at the local mean mixture fraction ξ [32]. For this purpose, the laminar flame speed and the thermal laminar flame thickness are pre-computed using the commercial software Cosilab and tabulated as functions of mixture fraction, S 0 L (ξ ) and δ 0 L (ξ ) respectively; τ ( ξ) is calculated based on the conditional temperature Q T (η, ζ ) with η = ξ and the Karlovitz number In this work an emphasis is put on the modelling of the spray combustion and we aim to provide closure for the complete set of evaporation terms that appear in the Favre-averaged transport equations of the conditioning variables. The mean evaporation source is computed by summing over the evaporated mass of all droplets in a CFD cell, where V represents the cell volume. The Favre-averaged transport equations contain several other evaporation source terms, which require modelling. In particular, for the mixture fraction variance equation, Giusti and Mastorakos [33] pointed out that both spray source terms could have a significant effect in the inner flame region. In single-conditional CMC, the terms ξΠ and ξ 2 Π are modelled either by summing ξ k s,iṁ d,i (k = 1 or 2) for all droplets in one cell [34] or, alternatively, by assuming that Π|η has the shape of a δ-function at the average surface mixture fraction ξ s [35]. In both cases, it is assumed that ξ s ≈ Y F,sat , calculated for the droplet temperature T d . This assumption is, however, only correct if the droplets evaporate upstream of the flame [13]. In general, F,Eq is the inverse function of the equilibrium fuel mass fraction Y F,Eq (ξ ). This reflects the case of a droplet evaporating in burned gases where some reaction products are present besides the fuel vapour at the droplet surface. In principle, this effect can be accounted for by integrating the spray source term in doubly conditional space. For a generic spray source term of the type FΠ, Fig. 1 shows a contour plot of the spray source term of the c-equation in conditional space, (C (η, ζ ) + ζ ), which can be directly evaluated using the conditional moments from CMC. Considering the conditional evaporation term shows that two necessary conditions for a source term of the progress variable are fulfilled: first, droplet evaporation in unburned mixture does not impact the progress variable since (C (η, ζ ) + ζ ) = 0 for ζ = 0 and, second, (C (η, ζ ) + ζ ) ≤ 1 for ζ = 1, which signifies c is automatically bounded at 1 with respect to this term.
In order to devise a model for Π|η, ζ it is recognised that the droplet temperature computed by the evaporation model fixes the fuel mass fraction -not the mixture fractionon the droplet surface. This means that the droplet surface corresponds to the isoline of Q F (η, ζ ) = Y F s in conditional space (represented by black lines in Fig. 1). Since evaporation only happens on the droplet surface, Π|η, ζ must be zero everywhere except on this isoline and FΠ can be calculated as a line integral. In Eq. 22, it was further assumed that Π(η, ζ ) is constant along this isoline. This is consistent with the evaporation model, wherė m d is computed based on the droplet temperature and the temperature far from the droplet surface, thus ignoring the exact temperature distribution in the near field.
Note that the model of Π|η, ζ presented above is only used to close the Favre-averaged evaporation source terms, but it is not used as a sub-model in the DCMC equation (see Section 2.2). The evaporation of the Lagrangian droplets was computed according to the model by Abramzon and Sirignano [36] with Stefan flow correction and non-unity Lewis number in the film; for the liquid droplets, the approximation of infinite conductivity is made.
where B M and B T are the Spalding mass and heat transfer numbers respectively. Sh and Nu are calculated using the Frössling correlation. The quantities ρ G , C pG , C pV , μ G , λ G , D G are evaluated at reference conditions for the gaseous boundary layer, computed according to the 1/3 rule for temperature and species; ρ L , C pL , latent heat L V and the fuel saturation vapour pressure p sat F are computed at the droplet temperature T d . The droplet experiences sphere drag and the impact of turbulence on the droplet motion is mimicked through stochastic dispersion for isotropic turbulence. No model for secondary break-up of droplets was used.

Chemistry
A detailed chemical mechanism for ethanol combustion [37] with 57 species and 383 reversible reactions is used. The mechanism performs well in predicting ignition delays and laminar flame speeds at ambient pressure when compared to experimental data. It has been successfully used in a CMC simulation by Giusti and Mastorakos [33].

Numerical set-up
The computational fluid dynamics toolbox OpenFOAM-2.3.0 is used to solve for the Favreaveraged flow field variables, u, ξ , ξ 2 , c and c 2 and to track the Lagrangian droplets. The coupling of the flow field solver and the DCMC model is effectuated as shown in a schematic in Ref. [38]. In addition toρ and T , the DCMC model also returns Y F , ω * c , c ω * c and the Favre-averaged evaporation terms (22) to the flow field solver. The motion and evaporation of the Lagrangian droplets are solved at the at beginning of each time step.
Applying the closure models and simplifying assumptions presented in Section 2.2, the DCMC equation (10) reduces to, Eq. 25 is integrated with a finite volume method whilst employing an operator splitting procedure. The fractional step approach is well established in CMC modelling and the significance of operator splitting errors in spray flames was investigated by Wright et al. [39]. The operator splitting approach allows for higher computational efficiency, since the stiff chemistry term is decoupled from the non-stiff convective terms. First, the transport in physical space by advection and turbulent diffusion, as well as the dilatation term (terms 2, 3 and 4 in the first line of Eq. 25 respectively) are computed. An upwind scheme is used for the advection term. Second, transport in conditional space (lines 2 and 3 in Eq. 25) is integrated for every DCMC cell. Finally, the chemical reaction source is solved independently for every (η, ζ )-node in every DCMC cell. For the second and the third sub-step the solver VODPK is used. This operator splitting procedure automatically assures that Q ψ conserves its linear dependence on ζ according to the definition of the progress variable. A 1D DCMC grid was used to discretise the physical domain with 35 DCMC cells along the burner axis for 0 < z/d < 5. The transfer of data from the fine CFD mesh to the coarse DCMC grid is achieved by integrating the conditional scalar dissipation rates over the volume of the DCMC cell [40].
The conditional space, D = {(η, ζ ) ∈ R 2 : 0 ≤ η ≤ 1, 0 ≤ ζ ≤ 1}, is discretised with 51 η-nodes, clustered around the stoichiometric mixture fraction η = ξ st (for the reaction of ethanol with air ξ st ≈ 0.1006) and 41 ζ -nodes, which are more closely spaced at ζ = 1. The derivatives in conditional space are discretised with second-order finite differences apart from ∂/∂ζ , which is computed with an upwind scheme. Dirichlet boundary conditions are set on all four sides of the conditional domain; mixing line and equilibrium condition for ζ = 0 and 1 respectively, air at η = 0 and fuel vapour at the boiling point at η = 1. The equilibrium condition was approximated by solving single-conditioned, nonpremixed "0D-CMC" equation, similar to the non-premixed flamelet equation, with a very low scalar dissipation rate N ξ,max = 1 1/s, compared to the critical scalar dissipation rate of approximately 367 1/s [41]. The initial condition for the doubly conditional moments in the DCMC cells is computed by solving the 0D-DCMC equation, that is to say the DCMC equation (25) with prescribed, fixed scalar dissipation rates and for a spatially homogeneous case, i.e. without the terms that represent transport in physical space, until steady state is reached. The solutions of this equation are similar to the results presented by Nguyen et al. [27]. As DCMC initial condition a very weakly strained solution of the 0D-DCMC equation was used; notably, N ξ,max = 1 1/s and N c,max = 200 1/s. The set of conditional moments used as inlet boundary condition is identical with the initial condition. Hence the conditional moments in the DCMC cell located at z/d = 0, shown later in Figs. 4 and 5, are representative of the DCMC initial condition.

Experimental test case
The DCMC model is applied to a piloted ethanol spray flame recently studied experimentally by Kariuki and Mastorakos [17]. The presence of droplet evaporation and, thus, the existence of mixture inhomogeneities, combined with substantial pre-vaporisation and premixing makes this flame a suitable candidate to test the present DCMC model.
In the experimental set-up, the ethanol spray is injected into the main flow of air ≈ 38 cm upstream of the nozzle with a diameter d = 42 mm. The flame is stabilised at the end of the nozzle through an annular pilot of diameter ≈ 51 mm and width ≈ 6 mm. The pilot is a open premixed methane-air flame with an approximately stoichiometric fuel-air ratio and a cold volume flow rate of 33.5 L/min. In the present work, we consider two reacting cases with overall equivalence ratios (liquid fuel to main air flow) φ = 0.62 and 0.82 (corresponds to ξ ≈ 0.065 and 0.084 respectively) and a non-reacting case. In all three cases the air-flow rate is 235 L/min and the liquid fuel mass injected in the cold case is the same as in the richer flame. Hence, the flames have the character of a pilot-stabilised turbulent premixed jet flames with a bulk velocity U b ≈ 3.04 m/s, calculated using the area of the nozzle.
The three-dimensional computational domain stretches −1 < z/d < 19 along the burner axis, where z = 0 is the position of the nozzle outlet and the pilot and the main flow inlet is retracted by one nozzle diameter; the diameter of the entire domain is 33d. Turbulence levels are set according to Laser Doppler Anemometry (LDA) measurements to u /U b ≈ 0.18 and estimating L T ≈ d/3. Since the degree of pre-vaporisation of the ethanol spray was not measured in the experiment, the mixture fraction at the inlet had to be estimated. According to a separate RANS simulation of the upstream part of the burner corresponding to the richer flame and the cold case, the inlet boundary condition for the mixture fraction was set to ξ = 0.04 and ξ 2 = 0. Even though the injected amount of liquid fuel is smaller in the leaner case, the same boundary condition is also used for the simulation of the flame with φ = 0.62. This will allow us to explore the sensitivity of the simulation results to the level of pre-vaporisation and premixing in the upstream region of a spray flame; in the rich flame approximately half the fuel is pre-vaporised compared to two thirds in the leaner flame. For the main flow, the inlet boundary condition for the reaction progress variable upstream of the nozzle is set to c = 0. For the purpose of the present paper whose aim is to discuss the model and its application, the uncertainty related to the inlet boundary condition is considered satisfactory.
The pilot is an open flame, which is not retracted relative to the main flow. Hence, dilatation in the pilot flame does not lead to a significant increase of the mean axial velocity but to an increase in width of the hot pilot stream as it is expected for a usual triangular flame. Following a single mixture fraction based approach in the simulations, the pilot is modelled as a premixed stoichiometric ethanol-air flame. The annular pilot flame itself is not resolved but instead modelled as a uniform inlet boundary condition with a laminar inflow of burnt gases, i.e. ξ = ξ st , ξ 2 = 0 and c = 1, with a mass flow rate corresponding to the experimental configuration. In the experimental rig the annulus that stabilises the pilot is very narrow and the pilot flow is broadened by the dilatation across the flame. If the width of the original annulus was used for the uniform inlet of hot gas, whilst the pilot mass flow rate was kept equal to its value in the experiments, the axial velocity of the pilot flow would be greatly over-estimated. Instead, the surface area of the pilot inlet is increased to assure a realistic axial velocity for the pilot flow. For this purpose, the surface of the pilot inlet is three times the surface area of the annulus in the rig.
A small, laminar co-flow of 0.1 m/s is set around the burner nozzle, where ξ = 0 and c = 1. Walls are considered adiabatic.  Figure 2 shows the radial profiles of the mean axial velocity of the gas phase for cold flow and the two flames considered here. The results are compared to Phase Doppler Anemometry (PDA) measurements of the mean axial velocity of the droplets in the range 0 ≤ r/d ≤ 0.5, not including the pilot stream, by Kariuki and Mastorakos [17]. In the cold case, the effect of the pilot stream is small and the flow spreads like a turbulent jet with a radial profile similar to a Gaussian bell curve. In contrast, for both flames studied, the axial mean velocity is almost constant for r/d < 0.5. This feature is well reproduced by the simulation. Next the droplet size information from the PDA measurements [17] is compared to simulation results. simulations, the droplet volume distributions at the inlet were set equal to the experimental ones from the axial position z = 6 mm. In general, a reasonable agreement between the experimental volume distributions and the simulation results is achieved. In both flames, the droplet volume distribution flattens in the range of small droplets (d < 30 μm) whilst it increases for larger droplets (40 μm < d < 70 μm). This can be explained by (i) a shorter heating-up period and thus quicker evaporation of smaller droplets and (ii) a faster decrease in diameter for smaller droplets when their temperature is nearly constant and d 2 is known to decrease approximately linearly in time. The phenomenon is slightly over-predicted in the simulations but can also be observed in the PDA data, in particular, for the richer flame. In Fig. 4 the conditional moments of CO 2 , OH, CH 2 O and temperature are shown for two different DCMC cells of the richer flame. Due to the definition of the progress variable based on the mass fraction of CO 2 , the conditional moment Q CO 2 is fixed for all DCMC cells equal to its shape defined by Eq. 7. The DCMC cell at z/d = 0 is located at the exit of the nozzle very close to the inlet and the scalar dissipation rates upstream of the flame are low. Hence, the structure of a weakly strained flame is mostly influenced by the advective term of the DCMC equation. The conditional moments shown here are almost identical with the DCMC inlet boundary condition and the DCMC initial condition. For η = ξ st the temperature rises quickly for 0 ≤ ζ < 0.5 and then flattens out for higher ζ , slowly approaching the temperature at equilibrium condition. In contrast, the DCMC cell at z/d = 1 contains the flame and thus experiences increased scalar dissipation rates. In particular, a high N c diffuses reactants in conditional space. Consequently, Q T rises almost linearly from ζ = 0 to 1 at η = ξ st . In the region further downstream of the flame N c = 0 and N ξ decays. As a result, the conditional moments experience less straining due to the scalar dissipation rates and slowly relax back to the structure of the weakly strained flame.

Results and Discussion
The conditional apparent reaction rates computed in the two DCMC cells previously discussed are shown in Fig. 5. They are compared toω c (η, ζ ) from unstrained premixed flamelets, which were computed using the commercial software Cosilab. In the weakly strained case at z/d = 0 the conditional reaction rate reaches significantly higher values than for the DCMC cell located at z/d = 1. This shows how the conditional reaction rate adjusts as the scalar dissipation rates increase. Moreover, in the strained flame ω * c |η, ζ also takes negative values. This is due to the contribution of the non-premixed term, which is also shown in Fig. 5.
Next the shape of the flames is assessed. Figure 6 shows sample averaged OH planar laser-induced fluorescence (OH-PLIF) images by Kariuki and Mastorakos [17], in comparison to Y OH from the present simulations. Comparisons are qualitative and no scale is shown for the mean OH PLIF intensity. On top of the simulation data c-isolines are plotted to emphasis the position and width of the flame brush. They can be directly compared to isolines of ensemble-averaged progress variable based on experimental data, even at a quantitative level. For this purpose, the instantaneous field of progress variable was determined experimentally by tracking the flame front through a thresholding procedure applied to the OH PLIF images [17]. For the richer flame with an overall equivalence ratio of φ = 0.82 (Fig. 6, bottom), the mean flame brush is relatively broad, occupying the range 1 < z/d < 2.5 along the burner axis, where the mean progress variable increases from 0.1 to 0.9. The simulation predicts a slightly longer flame but a thinner flame brush, stretching over the range 2.1 < z/d < 2.8. Moreover, Y OH in the burnt gases, downstream of the flame is lower than in the pilot stream, whilst this trend is not observed in the OH PLIF signal intensity. The simulation results also show that a significant proportion of fuel is not burned and the presence of this fuel in the hot combustion products leads to the production of CH 2 O downstream of the flame. The leaner flame (Fig. 6, top) with an overall equivalence ratio of φ = 0.62 is longer and experiments showed that the mean flame brush is present in the range 2.5 < z/d < 3.5 on the burner axis. This flame length is well predicted by the simulation. Moreover, the simulation also predicts a lower Y OH in the main flow compared to the pilot stream. Indeed, this behaviour is also found in the OH PLIF signal of the leaner flame. However, this feature is less pronounced in this case and the simulations, unexpectedly, predict that the combustion products of the leaner flame contain more Y OH than in the case of the richer flame. Since the pre-vaporised fraction of fuel is different in both cases, this feature can be directly related to the droplet terms in the transport equations of the conditioning variables, which will be discussed next.
In the following, some information on the relative magnitude of the various terms in the transport equations of ξ , ξ 2 , c, c 2 and the DCMC equation is given. Figure 7 shows the field of the mean mixture fraction and its variance for the richer flame with φ = 0.82. The mean mixture fractions set at the inlet, one nozzle diameter upstream of the burner exit, is below the lean flammability limit of an ethanol air mixture (ξ lean ≈ 0.05) to ξ . Upstream of the flame the mean temperature is low and the evaporation rate Π is low. Nevertheless ξ increases along the burner axis, first slowly to ξ ≈ 0.048 at z/d = 1 and then faster, reaching 0.59 at z/d = 2 in the preheat zone of the flame. Even though the evaporation rate is highest in the range of 0.5 < c < 0.9, some droplets still exist downstream of the flame brush and Π takes significant values until z/d = 3.5. Droplet evaporation is also the dominant source in the mixture fraction variance equation (14) and the production term due to mean mixture fraction gradients is negligible in the flame investigated in this work.
In contrast to Π , the mixture fraction variance source term is very small upstream of the c = 0.1 isoline because the droplets evaporating in this region have a low T d and thus ξ s is not much larger than ξ . Consequently, ξ 2 rises significantly in the region of the mean flame brush. In particular, this increase of ξ 2 is not directly counterbalanced by the scalar dissipation rate term, since a model for passive scalar mixing [29], was used due to the lack of alternatives. Even though the scalar dissipation rate is globally of the same order of magnitude as the droplet source of ξ 2 , we note that ξ 2 is locally over-predicted. In particular, this is the case in the core of the main flow for 2.5 < z/d < 3.5, where the production of variance outweighs its destruction (Fig. 7). This high level of mixture fraction variance in the post flame region of the richer flame explains low levels of Y OH in the main flow of the richer flame, compared to the pilot stream (Fig. 6). For high ξ 2 a broad β-pdf is presumed and the probability of finding flammable mixture that can react to form OH is low. This connection is also demonstrated by the unexpected fact that the simulations show higher Y OH for the leaner flame. Since the mass fraction of prevaporised fuel is ξ = 0.04 for both flames, in the lean case two thirds of fuel mass are fully premixed. Hence, evaporation produces lower levels of ξ 2 in the lean case leading to higher Y OH . For the same reasonρ is over-predicted leading to an underprediction of u z at downstream locations.
Fields of c and c 2 are shown in Fig. 8. In addition to the mean reaction source term, the cequation (8) also contains a droplet source term, which is computed as detailed in Eq. 22. As previously discussed, this source term is zero as long as the droplets evaporate in unburned mixture. On a global scaleṠ c is also two orders of magnitude lower than the apparent reaction rate ω * c and, thus, it only has a small effect on the shape of the flame brush. The same applies to the total evaporation variance source termṠ c 2 compared to c ω * c , such that the effect of evaporation on the c 2 -equation is negligible in the present case. In contrast, thė S c plays an important role in the region downstream of the flame, where it counter-balances a large fraction of the decrease in c, otherwise caused by the evaporative mass source acting on the mean density. For this purpose, a simplified model for this term can be proposed by comparingṠ c to c Π . Figure 9 shows that instead of using Eq. 22, a simplified model, can be used for this purpose.  This leaves the reaction source c ω * c and the production due to mean progress variable gradients as the main source terms in the c 2 -equation. Note that c ω * c locally also takes the role of an important variance sink term as the reaction reaches completion. On the side of the sink terms Fig. 8 shows a comparison of the mean scalar dissipation rate ε c , as computed in the present work using the model by Kolla et al. [30], and the model for a passive scalar [29] applied to the reaction progress variable. The latter takes the largest values in the thinnest part of the reacting shear layer that forms between the main flow and the pilot stream, but takes much smaller values than ε c in the region where the flame closes. Indeed, in the present simulations it was not possible to stabilise a flame using 2( ε/ k) c 2 to model the scalar dissipation rate, such that the model by Kolla et al. [30] was employed.

Conclusions
The CMC framework was used to build a new model for spray flames, which may feature characteristics of both premixed and non-premixed combustion. The DCMC equation for spray flames, conditioned on mixture fraction and reaction progress variable, was presented in this work. In the derivation, spray evaporation terms in the transport equations of the reactive scalar and both conditioning variables, including c, were explicitly considered, which added new terms to the DCMC equation.
In a preliminary application, the model was used to simulate the behaviour of a piloted ethanol spray flame with significant pre-vaporisation and strong premixing at two different, lean conditions. For this purpose, a set of sub-models was proposed to provide closure for the DCMC model and an operator splitting procedure was suggested to solve the model equation. The velocity field and the droplet distributions showed good agreement with experimental data and the flame shape prediction was promising in revealing the experimental trend due to overall equivalence ratio.
Using the conditional moments of the reactive scalars, available from DCMC, the unclosed spray terms in the Favre-averaged transport equations of the conditioning variables were modelled. Droplet evaporation was the dominant source term in the mixture fraction variance equation and had a large effect on the result, which had also been reported in previous studies. At the same time, closure of the mixture fraction scalar dissipation rate with a passive scalar mixing model seemed to be inadequate and main differences between simulation and experiment could be related to this imbalance of variance source and sink terms. The spray terms in the reaction progress variables equations were small compared to the reaction source terms and, thus, had little effect on the shape and width of the mean flame brush. This suggests that the effect of droplet evaporation can be neglected in the c 2equation. However, the evaporation source term should be included in the mean progress variable equation, but it can be approximated with a simple model.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.