A Priori Direct Numerical Simulation Modelling of the Curvature Term of the Flame Surface Density Transport Equation for Nonunity Lewis Number Flames in the Context of Large Eddy Simulations

A Direct Numerical Simulation (DNS) database of freely propagating statistically planar turbulent premixed flames with Lewis numbers Le ranging from 0.34 to 1.2 has been used to analyse the statistical behaviours of the curvature term of the generalised Flame surface Density (FSD) transport equation, in the context of the Large Eddy Simulation (LES). Lewis number is shown to have significant influences on the statistical behaviours of the resolved and sub-grid parts of the FSD curvature term. It has been found that the existing models for the sub-grid curvature term Csg do not capture the qualitative behaviour of this term extracted from the DNS database for flames with Le << 1. The existing models of Csg only predict negative values, whereas the sub-grid curvature term is shown to assume positive values within the flame brush for the Le = 0.34 and 0.6 flames. Here the sub-grid curvature terms arising from combined reaction and normal diffusion and tangential diffusion components of displacement speed are individually modelled, and the new model of the sub-grid curvature term has been found to capture Csg extracted from DNS data satisfactorily for all the different Lewis number flames considered here for a wide range of filter widths.


Introduction
Flame Surface Density (FSD) based reaction rate closure is well established in the context of Reynolds Averaged Navier-Stokes (RANS) simulations of turbulent premixed flames [1,2].The increased affordability of high performance computing has made Large Eddy Simulation (LES) an alternative simulation tool, where the large-scale physical processes are resolved, but modelling is still required for the subgrid quantities.The FSD-based reaction rate closure has recently been successfully extended for the purpose of LES [3][4][5][6][7][8][9][10][11][12][13][14].In LES simulation of premixed combustion, a Favre-filtered reaction progress variable transport equation is solved alongside other filtered conservation equations.The reaction progress variable is defined as c = (Y R0 − Y R )/(Y R0 − Y R∞ ), where Y R is the mass fraction of a suitable reactant and the subscripts 0 and ∞ denote the values in the fully unburned and burned gases, respectively.The generalised FSD is defined as Σ gen = |∇c| [3][4][5][6][7][8][9][10][11][12][13][14], where the overbar indicates the LES filtering process.The Favre-filtered reaction progress variable transport equation takes the following form: where Q = ρQ/ρ indicates the Favre filtered value of a general variable Q, u j is the velocity component in the jth direction, ρ is the density, D is the molecular diffusivity, and ẇ is the filtered reaction rate.The first two terms on right hand side of (1) denote the filtered molecular diffusion and reaction rates, respectively, and their combined contribution can be modelled using Σ gen in the following manner: where (Q) s = Q|∇c|/Σ gen indicates the surface-weighted filtered value of a general quantity Q and S d = Dc/Dt/|∇c| is the displacement speed, which denotes the speed at which a given c isosurface moves normal to itself with respect to an initially coincident material surface.The generalised FSD Σ gen is an unclosed quantity and is closed either by using an algebraic expression or by solving a modelled transport equation alongside other conservation equations.The algebraic closure is valid when the generation rate of flame surface area remains in equilibrium with its destruction rate, but this assumption is rendered invalid under unsteady conditions (e.g., combustion instabilities).
Under unsteady conditions, it is often advantageous to solve a modelled transport equation of Σ gen .The exact transport equation for the generalised FSD Σ gen is given as [1, 4-7, 9, 10, 12]: where N i = −(∂c/∂x i )/|∇c| is the ith component of flame normal vector.The terms on the left hand side of (3) denote transient and mean advection effects, respectively.The first three terms on the right hand side of (3) denote the effects of subgrid convection, flame surface area generation due to fluid-dynamic straining, and flame normal propagation, respectively.The last term of (3) describes the production/destruction of Σ gen due to flame curvature κ m = (∂N i /∂x i )/2 and thus referred to as the FSD curvature term [4-7, 9, 10, 12].It has been found in several previous studies [5-7, 9, 14] that the FSD curvature term remains a leading order contributor to the FSD transport for both unity and nonunity Lewis number turbulent premixed combustion.As the curvature term remains a leading order contributor to the FSD transport, the modelling of (S d ∇ • N) s Σ gen is crucial for the transport equation-based FSD closure.The statistical behaviour of (S d ∇ • N) s Σ gen is significantly affected by curvature dependence on S d [9,10,12].Earlier a priori Direct Numerical Simulation (DNS) analyses [9,10,12] showed that existing models for the subgrid curvature term C sg do not adequately capture the qualitative behaviour of this term obtained from DNS data.Moreover, the model parameters for the existing subgrid curvature term C sg models are found to be strong functions of the LES filter width Δ [9,10,12].To date, most existing FSD-based models have been proposed for unity Lewis number flames where the differential diffusion of heat and mass has been ignored.The Lewis number is defined as the ratio of thermal diffusivity α T to mass diffusivity D (i.e., Le = α T /D).
The effects of Le on the statistical behaviour of the FSD curvature term (S d ∇ • N) s Σ gen are yet to be analysed in detail, and this paper aims to bridge this gap in the existing literature.It is worth noting that, in a premixed flame, different species have different values of Lewis number.Thus, specifying a global Lewis number Le characterising the whole combustion process is not straightforward.The Lewis number of the deficient reactant is often considered to be the characteristic Le of the combustion process in question [15,16].Moreover, several previous studies [16][17][18][19][20][21][22][23][24][25][26][27][28][29] analysed the effects of differential diffusion of heat and mass by modifying the characteristic Lewis number in isolation, and the same procedure has been adopted here.In the present study, a simplified chemistry-based DNS database of statistically planar turbulent premixed flames with global Lewis numbers ranging from 0.34 to 1.2 has been considered to analyse the statistical behaviour of the FSD curvature term (S d ∇ • N) s Σ gen in the context of LES.In this context, the main objectives of this study are as follows: (1) to analyse the statistical behaviours of the subgrid FSD curvature term in the context of LES, for flames with different values of Lewis number; (2) to propose models for different components of the subgrid FSD curvature terms and assess their performances in comparison to the corresponding quantities extracted from DNS data.
The rest of the paper will be organised as follows.The necessary mathematical background will be provided in the next section.This will be followed by a brief description of the numerical implementation related to the DNS database.Following this, results will be presented and subsequently discussed.The main findings will be summarised, and conclusions will be drawn in the final section.

Mathematical Background
The curvature term of the FSD (S d ∇ • N) s Σ gen is often decomposed in the following manner [4-7, 9, 10, 12]: where C mean and C sg are the resolved and subgrid components of the FSD curvature term, respectively.The resolved curvature term C mean can be expressed in three different manners [5,9,10,12]: where M i = −(∂c/∂x i )/|∇c| is the ith component of the resolved flame normal vector.It was demonstrated by Chakraborty and Cant [10,12] that (5a) provides the best option for the resolved curvature term C mean , as it gives rise to the smallest magnitude of C sg among all the possibilities shown in (5a)-(5c).Equation (5a) was found to perform the best among the three possibilities shown in (5a), (5b), (5c) for this database.This is advantageous from the perspective of efficient modelling of the FSD curvature term (S d ∇ • N) s Σ gen as most of the modelling uncertainty is associated with C sg .Moreover, (5a) has also been used for the modelling of (S d ∇ • N) s Σ gen in previous LES simulations [5][6][7]13].For the present analysis (5a), (i.e., C mean = (S d ) s [∂(N i ) s /∂x i ]Σ gen ) will be considered for the resolved curvature term C mean .
It is often useful to decompose the flame displacement speed S d = (Dc/Dt)/|∇c| = [ ẇ + ∇ • (ρD∇c)]/ρ|∇c| in the following manner for the purpose of modelling the FSD curvature term [9-12, 30, 31]: where S r and S n are the reaction and normal diffusion components of displacement speed and S t is the tangential diffusion component of displacement speed.The following expression for C sg can be obtained using (6a)-(6d) and (5a) (i.e., C mean = (S d ) s [∂(N i ) s /∂x i ]Σ gen ): where dependences of (S r + S n ) and |∇c| significantly influence the statistical behaviour of C sg1 .Equation (8b) suggests that C sg2 is expected to assume negative values throughout the flame brush.
Hawkes and Cant [6,7] modified a version of the Coherent Flamelet Model (CFM) by Candel et al. [2] for the purpose of LES as: where α N = 1 − (N k ) s (N k ) s is a resolution parameter which vanishes when the flow is fully resolved and β 1 is a model parameter.Hawkes [5] discussed a possibility of modifying a RANS model proposed by Cant et al. [1] for the purpose of LES as: where A = 10.0, u Δ = 2 k/3 is the subgrid turbulent velocity fluctuation, k = (ρu i u i − ρ u i u i )/2ρ is the subgrid kinetic energy, and β 2 is a model parameter.Another model of C sg was proposed by Charlette et al. [4]: where β 3 is a model parameter.The models given by ( 9)-( 11) (henceforth will be referred to as CSGCFM, CSGCPB, and CSGCHAR, resp.)ensure that C sg vanishes when the flow is fully resolved (i.e., (N k ) s (N k ) s = 1.0 and Σ gen = |∇c|).
A priori DNS assessment of the CSGCFM, CSGCPB, and CSGCHAR models and the modelling of C sg1 and C sg2 will be addressed in Section 4 of this paper.

Numerical Implementation
In principle combustion, DNS should account for both three dimensionality of turbulence and detailed chemical mechanism.However, until recently, most combustion DNS studies were carried out either in two dimensions with detailed chemistry or in three dimensions with simplified chemistry due to the limitation of computer storage capacity.Although it is now possible to carry out three-dimensional DNS with detailed chemistry, they remain extremely expensive (e.g., millions of CPU hours and thousands of processors [32]) and the cost of an extensive parametric analysis based on threedimensional detailed chemistry-based DNS often becomes prohibitive.As the present analysis concentrates on an extensive parametric variation in terms of Lewis number, the chemical mechanism is simplified here by an Arrhenius-type irreversible single-step chemical reaction (i.e., Reactants → Products) following several previous studies [1][2][3][4][5][6][7][8][9][10][11][12]14].It has been found that the strain rate and curvature dependences of S d and |∇c| obtained from three-dimensional simplified chemistry DNS [25-27, 33, 34] are found to be qualitatively similar to the corresponding behaviours obtained from detailed chemistry-based DNS simulations [16,30,31,35].As the statistical behaviours of the FSD curvature term (S d ∇ • N) s Σ gen are strongly dependent on the curvature dependences of S d and |∇c|, the results for this analysis are expected to be valid even for detailed chemistry based simulations at least in a qualitative sense without much loss of generality.Several studies [3][4][5][6][7][9][10][11][12] have concentrated on a priori DNS modelling of FSD based on simplified chemistry in the past and the same approach has been adopted here.

International Journal of Chemical Engineering
A compressible three-dimensional DNS code SENGA [36] was used for the simulations where the conservation equations of mass, momentum, energy, and species are solved in nondimensional form.A cubic domain of each side equal to 24.1δ th is considered for the present DNS database where δ th is the thermal flame thickness, which is defined as δ th = (T ad − T 0 )/ Max |∇ T| L , and the subscript L refers to quantities in an unstrained planar laminar flame with T ad , T 0 , and T being the adiabatic flame, unburned gas, and instantaneous gas temperatures, respectively.The computational domain was discretised using a Cartesian grid of 230 × 230 × 230 with equal grid spacing in each direction.The grid spacing Δx is determined based on the flame resolution, and about 10 grid points are kept within the thermal flame thickness δ th for all the cases considered here.This grid spacing Δx corresponds to 0.73η, where η is the Kolmogorov length scale.The boundaries in the mean flame propagation were taken to be partially nonreflecting and were implemented using the Navier-Stokes Characteristic Boundary Conditions (NSCBC) technique [37].The boundary conditions in the transverse direction were taken to be periodic.The spatial derivatives for the internal grid points were evaluated using a tenth-order central differencing scheme, and the order of differentiation gradually decreases to a one-sided 2nd order scheme at the partially nonreflecting boundaries.The time advancement was carried out using an explicit low storage third-order Runge-Kutta scheme [38].
For the current DNS database, the turbulent velocity field was initialised using a pseudospectral method [39] following the Batchelor-Townsend turbulent kinetic energy spectrum [40].The flame is initialised using a steady planar unstrained laminar flame solution.The initial values of normalised root mean square (rms) turbulent velocity fluctuation u /S L , integral length scale to thermal flame thickness ratio l/δ th , heat release parameter τ = (T ad − T 0 )/T 0 , Damköhler number Da = lS L /u δ th , and Karlovitz number Ka = (u /S L ) 3/2 (l/δ th ) −1/2 are listed in Table 1.According to Peters [41], all the cases considered here can be taken to represent the thin reaction zone regime combustion, as Ka remains greater that unity.Standard values are considered for Prandtl number Pr and the Zel'dovich number β (i.e., Pr = 0.7 and β = T ac (T ad − T 0 )T 2 ad = 6.0,where T ac is the activation temperature).
Under decaying turbulence, DNS simulations should be carried out for a simulation time t sim ≥ max(t f , t c ) [42], where t f = l/u is the initial eddy turn over time and t c = δ th /S L is the chemical time scale.For this database, the statistics were extracted after about three eddy turn over times (i.e., 3t f = 3l/u ), which corresponded to one chemical time scale (i.e., t c = δ th /S L ).This simulation time remains small but comparable to several studies [3,24,28,[43][44][45][46][47] which contributed significantly to the FSD-based modelling in the past.The statistics presented in this paper did not change significantly since halfway through the simulation (i.e., 1.5t f = 1.5l/u ).The value of u /S L in the unburned gas ahead of the flame had decayed by 50% of its initial value when the statistics were extracted.By this time, the normalised integral value l/δ th had increased to around 1.7 times of its initial value.The values of u /S L and l/δ th at the time statistics were extracted are also representative of the thin reaction zones regime combustion [41].This DNS database was used extensively earlier for the purpose of RANS modelling [27,28,48,49], and the interested readers are referred to these papers for further details.
The DNS data was explicitly LES filtered using a Gaussian filter kernel in physical space for the purpose of a priori analysis.The filtered quantity Q( x, t), is given by where G( r) is the Gaussian filter kernel, which is defined in the following manner: The filtered quantities of interest were extracted for filter widths Δ ranging from 0.4δ th to 2.4δ th in steps of 0.4δ th .These filter sizes are comparable to the range of Δ used in several previous studies [3,4,[9][10][11][12]14] for a priori DNS analysis and span a useful range of length scales (i.e., from Δ ≈ 0.4δ th , where the flame is partially resolved, up to 2.4δ th , where the flame becomes fully unresolved and Δ is comparable to the integral length scale).

Results and Discussion
The instantaneous isosurfaces of c ranging from 0.01 to 0.99 at t c = δ th /S L are shown in Figure 1, which indicates that the flame wrinkling increases with decreasing Lewis number and this tendency is particularly prevalent for the Le 1 flames due to thermodiffusive instabilities [17][18][19][20][21][22][23][24][25][26][27][28][29].The unburned reactants diffuse into the reaction zone at a faster rate than the rate at which heat diffuses out in the Le < 1 flames.This gives rise to simultaneous presence of high temperature and reactant concentration in the reaction zone for the Le < 1 flames, which in turn leads to greater burning rate and flame surface area generation in comparison to the unity Lewis number flame.By contrast, heat diffuses faster than the diffusion rate of reactants into the reaction zone in the case of Le > 1, which reduces the burning rate and the rate of flame area generation in comparison to the unity Lewis number flame.The increase in burning rate and flame area generation with decreasing Lewis number can be substantiated by the values of normalised turbulent flame speed S T /S L and normalised flame surface area A T /A L which are presented in Table 2.The values of S T /S L have been evaluated by volume integrating the reaction rate ẇ using the expression S T = (1/ρ 0 A P ) V ẇdV , where A P is the projected area of the flame in the direction of mean flame propagation, while the values of A T /A L have been evaluated by volume integrating |∇c| (i.e., V |∇c|dV ) under both turbulent and laminar conditions.Table 2 shows that both S T /S L and A T /A L increase with decreasing Lewis number, and this effect is particularly prevalent in the flames with Le < 1 due to the presence of thermodiffusive instabilities [17][18][19][20][21][22][23][24][25][26][27][28][29].The increase in flame wrinkling with decreasing Lewis number is also visually evident from the c isosurfaces presented in Figure 1.
The variations of (S d ∇ • N) s Σ gen , C mean , and C sg conditionally averaged in bins of c isosurfaces for cases (a)-(e) are shown in Figure 2 for filter widths Δ = 8Δ m ≈ 0.8δ th and Δ = 24Δ m ≈ 2.4δ th , where Δ m is the DNS grid size.It is evident from Figure 2 that Le significantly affects the statistical behaviours of the curvature terms.The filter widths Δ = 8Δ m ≈ 0.8δ th , and Δ = 24Δ m ≈ 2.4δ th span a useful range of length scales (i.e., from Δ ≈ 0.8δ th , where the flame is partially resolved, up to 2.4δ th where the flame becomes fully unresolved and Δ is comparable to the integral length scale).In the Le 1 flames (e.g., cases (a) and (b)), the FSD curvature term (S d ∇ • N) s Σ gen behaves as a source term for the major part of the flame brush before assuming negative values towards the burned gas side for Δ = 8Δ m ≈ 0.8δ th .For Δ = 24Δ m ≈ 2.4δ th , the FSD curvature term (S d ∇ • N) s Σ gen acts as a source (sink) term towards the unburned (burned) gas side of the flame brush in the Le 1 flames.In the case of Le ≈ 1.0 flames (i.e., cases (c)-(e)) the curvature term (S d ∇ • N) s Σ gen behaves as a sink type term throughout the flame brush for all filter widths.It can be seen from Figure 2 9)-( 11)).This suggests that new models for C sg are warranted to account for the influences of nonunity Lewis number (i.e., Le / = 1.0).In order to be able to model the subgrid curvature term C sg , the decomposition prescribed in (8a)-(8b) has been used here.The variations of ( m ) s Σ gen , C sg1 and C sg2 conditionally averaged in bins of c isosurfaces for cases (a)-(e) are shown in Figure 3 for filter widths Δ = 8Δ m ≈ 0.8δ th and Δ = 24Δ m ≈ 2.4δ th .It is evident from Figure 3 that C sg2 remains negative throughout the flame brush for all cases and follows the qualitative behaviour of (−4(Dκ 2 m ) s Σ gen ).A comparison between ((S r + S n )∇ • N) s Σ gen and −4(Dκ 2 m ) s Σ gen reveals that −4(Dκ 2 m ) s Σ gen remains the major contributor to (S d ∇ • N) s Σ gen for all the flames at all values of Δ, which is consistent with the expected behaviour in the thin reaction zones regime [41].The contribution of ((S r + S n )∇ • N) s Σ gen remains significant for the Le < 1 cases (i. to the magnitude of −4(Dκ 2 m ) s Σ gen in the Le = 1.0 and 1.2 flames (i.e., cases (d) and (e)).Figure 3 demonstrates that C sg1 remains close to the magnitude of ((S r + S n )∇ • N) s Σ gen for all Δ for the Le = 1.0 flame (i.e., case (d)), indicating that (S r + S n ) s ∂(N i ) s /∂x i Σ gen does not play a major role in capturing the behaviour of ((S r + S n )∇ • N) s Σ gen .However, there is a significant difference in magnitudes of ((S r + S n )∇ • N) s Σ gen and C sg1 for small values of Δ (i.e., Δ < δ th ) in the nonunity Lewis number flames (i.e., cases (a)-(c) and (e)), which indicates that (S r + S n ) s ∂(N i ) s /∂x i Σ gen plays a key role for small values of filter width in these flames.For large values of filter width (i.e., Δ δ th ) C sg1 remains the major contributor to ((S r + S n )∇ • N) s Σ gen for all cases considered here, indicating that (S r + S n ) s ∂(N i ) s /∂x i Σ gen plays progressively less important role for increasing values of Δ. Figure 3 shows that there is a significant difference between −4(Dκ 2 m ) s Σ gen and C sg2 for all cases for small values of Δ, and the difference between these quantities decreases with increasing Δ.As most of the contribution of −4(Dκ 2 m ) s Σ gen remains unresolved for large values of Δ, the subgrid curvature term C sg2 remains the major contributor to −4(Dκ 2 m ) s Σ gen , indicating that (−(D∂N i /∂x i ) s ∂(N i ) s /∂x i Σ gen ) plays progressively less important role for increasing values of Δ where the flame is fully unresolved.However, the contribution of (−(D∂N i /∂x i ) s ∂(N i ) s /∂x i Σ gen ) remains significant for small values of Δ, where the flame is partially resolved.Figure 3 further shows that the order of magnitudes of both C sg1 and C sg2 remains comparable and thus accurate modelling of C sg1 and C sg2 is necessary for precise predictions of C sg .
As the range of κ m values obtained on a flame surface increases with increasing flame wrinkling, the magnitude of −4(Dκ 2 m ) s increases with decreasing Le, which in turn leads to increasing magnitude of −4(Dκ 2 m ) s Σ gen and C sg2 (see Figure 3).The positive contribution of C sg1 overcomes the negative contribution of C sg2 towards the unburned gas side of the flame brush for the Le = 0.34 and 0.6 flames (i.e., cases (a) and (b)) and yields a net positive contribution of C sg towards the reactant side of the flame brush (see Figure 2).
The statistical behaviours of ((S r + S n )∇ • N) s Σ gen and C sg1 depend on the nature of the curvature κ m = ∇ • N/2 dependences of (S r + S n ) and |∇c|, and the variation of (κ m ) s across the flame brush.The correlation coefficients for (S r + S n ) − κ m and |∇c| − κ m for five different c isosurfaces across the flame brush for all the cases are shown in Figures 4(a) and 4(b), respectively.For all cases, S t = −2Dκ m remains negatively correlated with κ m with a correlation coefficient close to (−1.0).However, Figures 4(a) and 4(b) demonstrate that Le significantly affects the curvature κ m dependences of (S r + S n ) and |∇c|.It can be seen from Figures 4(a) and 4(b) that (S r +S n ) and |∇c| remain positively (negatively) correlated with κ m for the Le < 1.0 (Le > 1.0) flames, whereas both (S r + S n ) and |∇c| show weak curvature dependences in the unity Lewis number flame.The positive correlations between (S r + S n ) and κ m , and between |∇c| and κ m strengthen with decreasing Le for the Le < 1 flames.The physical explanations for the observed influences of Lewis number on the curvature dependence of (S r + S n ) and |∇c| have been discussed elsewhere [25][26][27] and thus will not be discussed in this paper.The variations of (κ m ) s conditionally averaged in bins of c isosurfaces for cases A-E are shown in Figures 4(c) and 4(d) for filter widths Δ = 8Δ m ≈ 0.8δ th and Δ = 24Δ m ≈ 2.4δ th , respectively.It is evident from Figures 4(c) and 4(d) that (κ m ) s assumes positive (negative) values towards the unburned (burned) gas side of the flame brush.For small values of Δ, the surface-weighted filtered value of curvature (κ m ) s approaches to κ m (i.e., lim Δ → 0 (κ m ) s = κ m |∇c|/|∇c| = κ m ) and thus the ensemble averaged value of (κ m ) s remains small for small values of filter width as the ensemble averaged value of κ m remains negligible for statistically planar flames.The difference between the ensemble averaged values of (κ m ) s and κ m increases with increasing filter width Δ, as flame wrinkling increasingly takes place at the subgrid level.For the Le = 1.0 flame (i.e., case D), the combination of positive (negative) value of (κ m of ((S r + S n )∇ • N) s Σ gen and C sg1 towards the unburned (burned) gas side of the flame brush for all cases considered here, including the nonunity Lewis number flames where the curvature dependences of (S r + S n ) and |∇c| are particularly strong.
The dependences of (S r + S n ) s and Σ gen on 0.5 × ∂(N i ) s /∂x i Σ gen are likely to capture some of κ m dependences of (S r + S n ) and |∇c| at small values of filter widths Δ (i.e., Δ < δ th ), where the flame is partially resolved.This effect is particularly prevalent in the nonunity Lewis number flames where both (S r + S n ) and |∇c| are strongly correlated with curvature κ m even though the flames are statistically planar in nature.As a result of this, the contribution of (S r + S n ) s ∂(N i ) s /∂x i Σ gen remains close to that of ((S r + S n )∇ • N) s Σ gen for small filter widths (i.e., Δ < δ th ) for the non-unity Lewis number flames, which is reflected in the small contribution of C sg1 (see Δ = 0.8δ th variations in Figures 3(a)-3(c) and Figure 3(e)).The correlation between the resolved quantities (e.g., dependences of (S r + S n ) s and Σ gen on 0.5 × ∂(N i ) s /∂x i Σ gen ) weakens with increasing filter width Δ due to smearing of local information.Moreover, physical processes take place increasingly at the subgrid level for Δ δ th , and thus (S r + S n ) s ∂(N i ) s /∂x i Σ gen does not capture the behaviour of ((S r + S n )∇ • N) s Σ gen for large filter widths in all cases considered here, including the nonunity Lewis number flames where the curvature dependences of (S r +S n ) and |∇c| are particularly strong.This leads to C sg1 ≈ ((S r + S n )∇ • N) s Σ gen for Δ δ th in all cases considered here (see Δ = 2.4δ th variations in Figures 3(f  C sg1 overcomes the negative contribution of C sg2 towards the unburned gas side of the flame brush in the Le = 0.34 and 0.6 flames, which lead to positive value of C sg = C sg1 + C sg2 towards the unburned gas side for all values of Δ in these cases (see Figure 2).By contrast, negative values of C sg2 overcome the positive contributions of C sg1 towards the unburned gas side of the flame brush in the Le = 0.8, 1.0, and 1.2 flames, which lead to negative values of C sg = C sg1 + C sg2 throughout the flame brush in these cases (see Figure 2).The subgrid fluctuations of the surface-weighted contributions of (S r + S n ) and ∇ • N are scaled here using S L and (Σ gen − |∇c|), respectively, to propose the following model for C sg1 :

Correlation coefficients
where β 4 , c * , a Σ , and m are the model parameters.The function is used to capture the correct qualitative behaviour of C sg1 throughout the flame brush.In a compressible, LES simulation c is readily available and c needs to be extracted from c.The methodology of extracting c from c in the context of LES was discussed elsewhere [9,10,12] and will not be discussed in detail in this paper.The model parameter c * ensures that the transition from positive to negative value of C sg1 takes place at the correct location within the flame brush.The quantity (Σ gen − |∇c|) vanishes when the flow is fully resolved (i.e., lim , and thus C sg1 becomes exactly equal to zero when the flow is fully resolved (i.e., Δ → 0) according to (14).It has been found that m = 1.85 enables (14)  these local dependences also appear in the resolved scale but their strength diminishes with increasing Δ due to filtering operation.As the resolved and subgrid curvature terms are closely related [9,10,12], the qualitative behaviour of C sg1 is also affected by the curvature dependences of displacement speed components and scalar gradient at the resolved scale, which leads to the variation of the optimum values of a Σ , β 4 , and c * with Le and Δ.The model parameter β 4 needs to be deceased for decreasing values of Σ gen for satisfactory prediction of (14).The prediction of ( 14) ensemble averaged on c isosurfaces is compared with the ensemble averaged values of C sg1 in Figure 5 for all cases considered here for the optimum values of β 4 , c * , and a Σ for Δ = 0.8δ th and Δ = 2.4δ th when m is taken to be m = 1.85.The optimum values of β 4 , c * , and a Σ are estimated by calibrating the prediction of ( 14) with respect to the ensemble averaged values of C sg1 obtained from DNS data and the variation of the optimum values of β 4 /Σ gen , c * , and a Σ with Δ/δ th for all cases are shown in Figure 6.The optimum values of β 4 /Σ gen , c * , and a Σ are parameterised here as where ; where ; Figure 5 shows that ( 14) satisfactorily predicts C sg1 when m is taken to be m = 1.85, and the optimum values of β 4 , c * , and a Σ are used.According to the parameterisation given by (15a)-(15e), β 4 increases with decreasing Le, as the effects of chemical reaction strengthen with decreasing values of Lewis number (see Table 2).Moreover, β 4 /Σ gen , c * , and a Σ approach to asymptotic values for large values of Δ and turbulent Reynolds number based on LES filter width Re Δ = 4ρ 0 2 k/3Δ/μ 0 , where ρ 0 and μ 0 are the unburned gas density and viscosity, respectively.
Here, the contribution of (Dκ  sub-grid fluctuations of D are taken to scale with S L /Σ gen .
The above relations are utilised here to propose a model for C sg2 in the following manner: where Ξ Δ = Σ gen /|∇c| is the wrinkling factor [8,11,43,50,51], β 5 and n are the model parameters, and c(1 − c) is used to capture the correct qualitative behaviour of C sg2 .The subgrid curvature term C sg2 vanishes when the flow is fully resolved according to ( 16), (i.e., lim It has been found that (16) satisfactorily captures the behaviour of C sg2 throughout the flame brush for n = 1.0 in all cases when a suitable value of β 5 is used.The variation of the global mean optimum values of β 5 with Δ/δ th is shown in Figure 6 for all cases considered here.The optimum values of β 5 have been parameterised here in the following manner: where  (17d) The predictions of ( 16) ensemble averaged on c isosurfaces are compared with ensemble averaged values of C sg2 in Figure 5 for all cases at Δ = 0.8δ th and Δ = 2.4δ th , which show that (16) satisfactorily predicts the statistical behaviour of C sg2 when n is taken to be n = 1.0 and the optimum value of β 5 is used.According to (17a)-(17d), β 5 approaches to asymptotic values for large values of Δ and Re Δ = 4ρ 0 2 k/3Δ/μ 0 .Equations ( 13) and (15a)-(15e) can be combined to propose a model for C sg in the following manner: The above model will henceforth be referred to CSGNEW model in this paper.Equation ( 18) allows for a positive contribution of C sg through the contribution of which is absent in the CSGCAND, CSGCANT, and CSGCHAR models.The predictions of the CSGCAND, CSGCANT, CSGCHAR, and CSGNEW models for Δ = 0.8δ th and Δ = 2.4δ th are compared with C sg obtained from DNS in Figure 7 for the optimum values of β 1 , β 2 , β 3 , and β 5 .The optimum values of β 1 , β 2 , and β 3 are estimated by calibrating the models based on the ensemble averaged value of C sg obtained from DNS data.The variations of the optimum values of β 1 , β 2 , and β 3 with Δ for all cases are also shown in Figure 6.It is evident from Figure 6 that β 1 , β 2 , β 3 , and β 5 remain greater than unity for all cases.This is found to be consistent with the realisability analysis by Hawkes and Cant [52]. Figure 6 further demonstrates that the optimum values of β 1 , β 2 , and β 3 change appreciably with increasing Δ, which is consistent with earlier findings [9,10,12].Moreover, optimum values of β 1 , β 2 , and β 3 for a given Δ are affected by Le (see Figure 6).It is worth noting that parameterisation of the optimum values of β 1 , β 2 and β 3 also yields complex relations similar to (15a)-(15e) and (17a)-(17d).However, such parameterisation is not presented here because the CSGCAND, CSGCANT, and CSGCHAR models do not capture the qualitative behaviour of C sg for the Le = 0.34 and 0.6 flames.
It can further be seen from Figure 7 that the CSGCHAR model tends to overpredict the negative values of C sg towards the unburned gas side in cases C-E (Le = 0.8, 1.0 and 1.2 flames), and this behaviour becomes more prominent with increasing filter size.It is clear from Figure 7 that for Δ = 24Δ m = 2.4δ th , the CSGCHAR model predicts the maximum magnitude of C sg near the middle of the flame whereas the actual maximum magnitude of C sg is attained slightly towards the burned gas side.The CSGCAND and CSGCANT models give comparable performance for optimum values of β 1 and β 2 in cases C-E.However, the CSGCAND and CSGCANT models do not satisfactorily capture the qualitative behaviour of C sg and underpredict (overpredict) the magnitude of C sg towards the burned gas side (middle) of the flame brush in the Le = 0.8, 1.0 and 1.2 flames.Figure 7 demonstrates that the CSGNEW model captures the qualitative behaviour of C sg in a better manner than the CSGCAND and CSGCANT models and the quantitative agreement between C sg and the CSGNEW model remains better than the CSGCAND, CSGCANT, and CSGCHAR models in all cases for all values of Δ when optimum values of β 4 , β 5 , a Σ , and c * are used.

Conclusions
The LES modelling of the curvature term (S d ∇ • N) s Σ gen of the generalised FSD Σ gen transport equation has been addressed here using a simplified chemistry-based DNS database of freely propagating statistically planar turbulent International Journal of Chemical Engineering premixed flames with Lewis number Le ranging from 0.34 to 1.2.The statistical behaviours of the subgrid curvature term C sg for a range of different values of Δ have been analysed in terms of its contributions C sg1 and C sg2 arising from the combined reaction and normal diffusion component and tangential diffusion components of displacement speed (i.e., (S r + S n ) and S t = −2Dκ m ), respectively.The Lewis number is shown to have significant influences on the statistical behaviours of the resolved and subgrid components of the FSD curvature term.Detailed physical explanations have been provided for the observed filter size and Lewis number dependences of the different components of (S d ∇ • N) s Σ gen .Models have been identified for individual components of the subgrid curvature term (i.e., C sg1 and C sg2 ), and the performances of these models have been compared to the corresponding quantities extracted from DNS data.It has been found that the new models for C sg1 and C sg2 satisfactorily capture the statistical behaviours of the corresponding terms extracted from DNS data.It has been found that the existing models for the subgrid curvature term C sg do not satisfactorily capture the qualitative behaviour of the corresponding quantity extracted from DNS data for all the flames considered here.This problem is particularly prevalent for flames with small values of Lewis number (e.g., Le = 0.34 and 0.6) where C sg locally assumes positive values, whereas the existing models can only predict negative values of C sg .The performance of the newly proposed model for C sg has been found to be better than the existing models, and it has been shown to capture positive contributions of C sg for the Le 1 flames.The present analysis has been carried out using a DNS database with moderate value of Re t in the absence of the effects of detailed chemistry and transport.As simplified chemistry-based DNS qualitatively captures the curvature κ m = ∇ • N/2 and strain rate dependences of S d and |∇c| obtained from chemistry based simulations, it can be expected that the statistical behaviours of the curvature term (S d ∇ • N) s Σ gen presented in this paper will be valid at least in a qualitative sense in the context of detailed chemistry.However, the quantitative values of the model parameters (i.e., β 4 , β 5 , a Σ , and c * ) may need to be altered in the presence of detailed chemistry.Thus, three-dimensional DNS data with detailed chemistry and experimental data at higher values of Re t will be necessary for more comprehensive modelling of (S d ∇ • N) s Σ gen and C sg in the context of LES.Moreover, the newly proposed models need to be implemented in LES simulations for the purpose of a posteriori assessments.However, this is kept beyond the scope of this paper.Several previous studies [3-7, 9-12, 43-49] concentrated purely on the model development based on a priori analysis of DNS data and the same approach has been adopted here.Implementation of the newly developed models in LES simulations will form the basis of future investigations.
that C mean acts as a source (sink) term for cases (a)-(b) ((c)-(e)).The magnitude of C mean (C sg ) decreases

(
increases) with increasing Δ in all cases, and for large filter widths (S d ∇ • N) s Σ gen is principally made up of C sg .The LES filtering is a convolution process, and the weighted averaging involved in the filtering process leads to a decrease in the magnitude of C mean with increasing filter width Δ.The flow becomes increasingly unresolved with increasing filter width Δ, and this is reflected in the rise in C sg magnitude with increasing filter width Δ.The resolved curvature term C mean = (S d ) s ∂(N i ) s /∂x i Σ gen can be seen to capture the behaviour of the curvature term (S d ∇ • N) s Σ gen , well at small filter widths (i.e., Δ ≤ δ th ) for flames with Le ≈ 1.0 (i.e., cases (c)-(e)).However, the magnitude of C mean decreases with increasing Δ and it does not capture the behaviour of the FSD curvature term (S d ∇ • N) s Σ gen for the Le 1.0 flames (i.e., cases (a) and (b)).The subgrid curvature term, C sg follows the qualitative behaviour of the FSD curvature term (S d ∇ • N) s Σ gen for all filter widths.The subgrid curvature term C sg almost entirely makes up the FSD curvature term (S d ∇ • N) s Σ gen for Δ δ th , and this is especially true for the Le 1.0 cases (i.e., cases (a) and (b)).It can further be observed from Figure 2 that C sg assumes positive values towards the unburned gas side of the flame brush in the Le 1 flames (e.g., cases (a) and (b)), whereas the existing models for C sg allow for only negative values (see (
) s and weak (S r +S n )−κ m and |∇c|−κ m correlations gives rise to positive (negative) values of the ensemble averaged values of ((S r + S n )∇ • N) s Σ gen and C sg1 towards the unburned (burned) gas side of the flame brush for all values of Δ.The predominant positive (S r + S n ) − κ m and |∇c| − κ m correlations give rise to positive values of the ensemble averaged values of ((S r + S n )∇.N) s Σ gen and C sg1 throughout the flame brush for small values of Δ in the Le = 0.34, 0.6, and 0.8 flames.By contrast, negative (S r + S n ) − κ m and |∇c| − κ m correlations (see Figures 4(a) and 4(b)) give rise to negative values of the ensemble averaged values of ((S r + S n )∇ • N) s Σ gen and C sg1 throughout the flame brush for small values of Δ in the Le = 1.2 flame.These local dependences are progressively smeared with increasing Δ because of the convolution operation associated with LES filtering process, and this leads to positive (negative) values )-3(j)).It can be seen from Figure3that the positive contribution of International Journal of Chemical Engineering

Table 1 :
Initial values of the simulation parameters and nondimensional numbers relevant to the DNS database.Case Le u /S L l/δ L l/δ th δ th /η τ

Table 2 :
The effects of Lewis number on normalised turbulent flame speed S T /S L and normalised flame surface area A T /A L after 3.0 initial eddy turn over times.