Modeling of Lewis number dependence of scalar dissipation rate transport for Large Eddy Simulations of turbulent premixed combustion

ABSTRACT The influences of differential diffusion of heat and mass on the Favre-filtered scalar dissipation rate (SDR) transport have been analyzed and modeled using a priori analysis of Direct Numerical Simulations (DNS) data of freely propagating statistically planar turbulent premixed flames with different values of global Lewis number, Le. The DNS data has been explicitly filtered using a Gaussian filter to obtain the unclosed terms of the Favre-filtered SDR transport equation, arising from turbulent transport (T1), density variation due to heat release (T2), strain rate contribution due to the alignment of scalar and velocity gradients (T3), correlation between the gradients of reaction rate and reaction progress variable (T4), molecular dissipation of SDR (−D2), and diffusivity gradients f(D). The statistical behaviors of these terms and their scaling estimates reported in a recent analysis have been utilized here to propose models for these unclosed terms in the context of Large Eddy Simulations (LES) and the performances of these models have been assessed using the values obtained from explicitly filtered DNS data. These newly proposed models are found to satisfactorily predict both the qualitative and quantitative behaviors of these unclosed terms for a range of filter widths Δ for all Le cases considered here.


Introduction
Lean premixed combustion has been identified as one of the possible ways to reduce pollutant emission from gasoline engines and industrial gas turbines [1]. Lean hydrogen and hydrogen-blended hydrocarbon combustion has the potential to attenuate pollutants and greenhouse gas emissions [2,3]. However, the flames with abundance of fast diffusing species such as hydrogen either in molecular or in atomic form give rise to a significant level of differential diffusion of heat and mass. The differential diffusion of heat and mass can be characterized by a nondimensional number known as the Lewis number Le, which is defined as the ratio of thermal diffusivity α T to mass diffusivity D (i.e., Le ¼ α T /D). In actual premixed combustion it is often not straightforward to assign a single global Lewis number in the presence of several species with different Lewis numbers. Often the Lewis number of the deficient species is considered to be the global Lewis number [4], whereas Law and Kwon [5] proposed a methodology of evaluating the effective Lewis number based on heat release measurements. More recently Dinkelacker et al. [6] proposed an algebraic expression for the effective Lewis number based on mole fractions of major species. A number of previous analyses concentrated on the effects of global Lewis number on different aspects of premixed combustion in isolation  and the same approach has been adopted here.
Modeling of the differential diffusion arising from non-unity global Lewis number remains pivotal to high-fidelity engineering simulations, which are likely to play important roles in the development of new-generation combustors using either hydrogen or hydrogen-blended fuels. Prediction of the micro-mixing rate of hot products and cold unburned gas plays a key role in the modeling of turbulent reacting flows and a quantity known as the scalar dissipation rate (SDR) characterizes this micro-mixing rate [29,[30][31][32]. Furthermore, the Favre-mean value of SDR of reaction progress variable c in premixed turbulent flames can be related to the mean reaction rate in the context of Reynolds Averaged Navier Stokes (RANS) simulations [23,[33][34][35]. The instantaneous SDR of reaction progress variable is defined as [23,[33][34][35][36][37][38] where D is the diffusivity of reaction progress variable c. Recent analyses have demonstrated [36][37][38] that the SDR-based reaction rate closure for RANS can also be used for the modeling of the filtered reaction rate _ w based on the Favre-filtered SDR of a reaction progress variable (i.e., Ñ c ¼ qDrc:rc=q) in the Nomenclature c reaction progress variable c m thermo-chemical parameter C P specific heat at constant pressure C V specific heat at constant volume C F model parameter C 3  reactant mass fraction in unburned gas Y R∝ reactant mass fraction in burned gas α T thermal diffusivity α T0 thermal diffusivity of the unburned gas β Zel'dovich number β 3 , β′ 3 model parameters γ ratio of specific heats (¼C P /C V ) γ 1  context of Large Eddy Simulation (LES) in the following manner when the filter size Δ remains greater than the thermal flame thickness δ th ¼ (T ad À T 0 )/Max| ∇ T| L (where T ad , T 0 , and T are the adiabatic flame, unburned gas, and instantaneous temperatures, respectively): where ρ is the density and q ¼ qQ=q is the Favre-filtered value of a quantity Q with the over-bar indicating an LES filtering operation. In Eq.
(2) f b (c) is the reacting mode probability density function (pdf) of c and the subscript "L" refers to the planar laminar flame conditions. By assuming f b (c) as a smooth function, regardless of the exact form, the numerical value of c m remains within a range of 0.7-0.9 for typical hydrocarbon-air mixtures [33]. The modeling of SDR not only is useful for the closure of filtered reaction rate but also plays a pivotal role in the closure of micro-mixing rate in the context of pdf methodology [30,[39][40][41]. For turbulent premixed flames, the Favre-filtered SDR Ñ c can be modeled either by using an algebraic expression in terms of the resolved quantities or by solving a modeled transport equation. A few recent analyses [36][37][38][39][40][41] have concentrated on the algebraic closure of SDR for turbulent premixed flames in the context of LES. Algebraic closures are suitable when an equilibrium is maintained between the generation and destruction rates of Ñ c , but this assumption may be rendered invalid under some conditions (e.g., low Damköhler number lean premixed combustion). A number of previous analyses [34,[42][43][44][45][46][47][48][49][50][51] concentrated on the modeling of SDR transport in turbulent premixed combustion in the context of RANS simulations. Interested readers are referred to Ref [34]. for a detailed review of the existing modeling methodologies for SDR transport in the context of RANS simulations. Recent advancements in highperformance computing have made LES of industrial flows more affordable than in the past, and LES is more successful in capturing unsteady flow features than RANS. However, relatively limited attention has been given to the investigation of SDR transport in the context of LES [52,53]. Recently, models for the unclosed terms of the SDR Ñ c transport equation for unity Lewis number flames in the context of LES have been proposed [53], but the differential diffusion effects due to non-unity Le were not addressed. A recent analysis [28] concentrated on the influences of global Le on the statistical behaviors of the unclosed terms of the Ñ c transport equation based on an order-of-magnitude approach, which successfully explained the effects of global Le and the filter width Δ dependences of the Favre-filtered SDR Ñ c and its transport. It has been found that Le has significant influences on both the qualitative and quantitative behaviors of the unclosed terms of the SDR Ñ c transport equation [22,28]; however, the modeling of Le effects on these unclosed terms is yet to be addressed, and the present analysis aims to address this gap in the existing literature. In this respect the main objectives of this paper are: i. To propose models for the unclosed terms of the SDR transport equation in such a manner that the performances of these models remain satisfactory for a range of Δ and Le. ii. To assess the performances of the newly proposed models with respect to explicitly filtered Direct Numerical Simulation (DNS) data. These objectives are addressed here by conducting a priori analysis using a DNS database of statistically planar turbulent premixed flames with a range of different values of Le (i.e., Le ¼ 0.34-1.2). The details related to mathematical background and numerical implementation are provided in the next section. This is followed by the presentation of the results and subsequent discussion. The main findings are summarized and conclusions are drawn in the final section of this paper.

Mathematical background & numerical implementation
Three-dimensional DNS simulations with detailed chemistry are now possible, but they remain extremely expensive and need several millions of CPU hours [54] for conducting extensive parametric variations and carrying out explicit filtering of DNS data using a range of filter widths Δ, as has been carried out in the current study. Thus, the chemical mechanism has been simplified here as a singlestep Arrhenius-type irreversible chemical reaction. Under this condition the species field is uniquely represented by a reaction progress variable c, which can be defined by using the mass fraction of a suitable reactant Y R as c ¼ (Y R0 À Y R )/(Y R0 À Y R∞ ), where subscripts 0 and ∞ denote the values in the unburned and burned gases, respectively. The transport equation of c can be used to derive a transport equation of Ñ c ¼ qDrc:rc=q, which takes the following form [28,53]: |ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl {zffl ffl ffl ffl ffl ffl ffl ffl ffl ffl } where u i is the ith component of the velocity vector. On the left-hand side of Eq. (3) the terms denote the transient effects and resolved advection of Ñ c , respectively. The term D 1 depicts the molecular diffusion of Ñ c , and the terms T 1 , T 2 , T 3 , T 4 , (À D 2 ), and f(D) are all unclosed and expressed as follows: where _ w is the reaction rate of c. The term T 1 represents the effects of sub-grid convection, whereas T 2 denotes the effects of density-variation due to heat release. The term T 3 is determined by the alignment of ∇c with local strain rate e ij ¼ 0.5(∂u i /∂x j þ ∂u i /∂x i ), and this term is commonly referred to as the scalar-turbulence interaction term. The term T 4 arises due to the correlation between r _ w and ∇c, whereas (À D 2 ) denotes the molecular dissipation of SDR; these terms will henceforth be referred to as the reaction rate term and the dissipation term, respectively. The term f(D) denotes the effects of D variation. A priori DNS modeling of the above-mentioned unclosed terms will be discussed in Section 3 of this paper.
For the present analysis, a DNS database of freely propagating turbulent premixed flames has been considered. The simulation domain is taken to be a cube of 24.1δ th � 24.1δ th � 24.1δ th , which is discretized using a uniform Cartesian grid of 230 � 230 � 230 points, ensuring about 10 grid points are kept within Min(δ L , δ th ), where δ L ¼ 1/(Max |∇c| L ) is an alternative flame thickness based on |∇c| and the values of δ L /δ th for cases A-E (with Le ¼ 0.34, 0.6, 0.8, 1.0, and 1.2) are provided in Table 1. The initial values of the normalized root-mean-square (rms) value of turbulent velocity u′/S L , integral length scale to thermal flame thickness ratio l/δ th , Damköhler number Da ¼ lS L /u′δ th , Karlovitz number Ka ¼ (u′/S L ) 3/2 (l/δ th ) À 1/2 , turbulent Reynolds number Re t ¼ ρ 0 u′l/μ 0 , and τ ¼ (T ad À T 0 )/T 0 are presented in Table 1 along with domain and grid sizes, where ρ 0 and μ 0 are the unburned gas density and viscosity, respectively, and S L is the unstrained laminar burning velocity. The flamelet assumption is likely to be valid for the values of u′/S L and l/δ th considered here, and all cases considered here represent the thin reaction zones regime combustion according to the regime diagram by Peters [55].
The simulations have been carried out using a well-known DNS code called SENGA [56]. For all cases the boundary conditions in the mean flame propagation direction are considered to be partially nonreflecting, whereas boundaries in the transverse directions are considered to be periodic. A 10th order central difference scheme is used for spatial differentiation for the internal grid points and the order of differentiation gradually drops to a one-sided second-order scheme at the non-periodic boundaries. A low-storage third-order Runge-Kutta method is used for explicit time advancement for all the governing equations. In all cases flame-turbulence interaction takes place under decaying turbulence, which necessitates the simulation time t sim ≥ max(t f , t c ), where t f ¼ l/u′ is the initial eddy turnover time and t c ¼ a T0 =S 2 L is the chemical time scale, with α T0 being the unburned gas thermal diffusivity. The simulations have been carried out for about 3.34t f ¼ 3.34l/u′, which amounts to approximately 1:75a T0 =S 2 L for all cases considered here. Several studies [12-15, 19, 57-61] with either similar or smaller simulation time have contributed significantly to the fundamental understanding and modeling of turbulent premixed combustion in the past. By the time the statistics were extracted, the value of u′/S L in the unburned reactants ahead of the flame had decayed by about 50%, whereas the value of l/δ th had increased by about 1.7 times, relative to their initial values. This database has been used in several previous analyses [20][21][22][23][24][25][26][27][28] and it was shown in Ref [23]. that the volume-integrated burning rate for the Le ¼ 1.0 and 1.2 flames reached quasi-steady state by the time statistics were extracted. However, the Le < 1.0 flames are thermo-diffusively unstable and thus the volumeintegrated burning rate increases with time for these cases [23]. The qualitative nature of the statistics was found to remain unchanged and the scaling estimates presented in the next section remain valid since t ¼ 2.0 l/u′ for all cases considered here.
The unclosed terms of the transport equation of Ñ c have been evaluated by explicitly filtering DNS data using a standard three-dimensional Gaussian filter [28,53,57,58,60]: GðrÞ ¼ ð6=pD 2 Þ 3=2 expðÀ 6r �r=D 2 Þ and the filtered values of a general quantity Q are given by the following integral: QðxÞ ¼ R Qðx ÀrÞGðrÞdr. In the next section, the results will be presented for Δ ranging from Δ ≈ 0.4δ th to Δ ≈ 2.8δ th . This range of filter widths is comparable to the range of Δ used in several previous a priori DNS analyses [57,58,60], and addresses a range of different length scales from Δ comparable to δ th ≈ 1.75δ z (δ z ¼ α T0 /S L is the Zel'dovich flame thickness) up to 2.8δ th ≈ 5.0δ z , where Δ is comparable to the integral length scale.

Results and discussion
The distributions of c on x 1 À x 2 mid-plane for Δ ¼ 0.8δ th , 1.6δ th , and 2.8δ th at t ¼1.75 t c for cases A-E are shown in Figure 1, which shows an increase in the extent of flame wrinkling with decreasing Le. The extent of flame wrinkling can be quantified in terms of the normalized turbulent flame surface area A T /A L , where the flame surface area is evaluated using the volume integration of the form: A ¼ ∫ V |∇c|dV with subscripts "T" and "L" denoting the turbulent and laminar flame values, respectively [28]. The values of A T /A L and the normalized turbulent burning velocity S T /S L (where Table 2, which demonstrates that both A T /A L and S T /S L increase significantly with decreasing Lewis number. The burning rate per unit area in turbulent flames increases (decreases) compared with the corresponding laminar value as a result of negative (positive) Markstein length [7][8][9][10] for the Le < 1 (Le > 1) flames. This, in turn, leads to Table 2). It can be further seen from Figure 1 that the flame brush thickens (i.e., the magnitude of rc decreases) and the extent of flame wrinkling decreases with increasing Δ as a result of the smearing of local information due to the convolution operation associated with LES filtering. As the SDR is related to the reaction rate, and the gradient of the reaction progress variable, the effects of Le on burning rate and Δ dependence of rc are expected to influence the statistical behavior of SDR Ñ c and its transport. The effects of Le and Δ on the statistical behavior of SDR Ñ c and its transport have been analyzed elsewhere [28] and the current analysis will only concentrate on the influences of global Lewis number on the modeling of SDR transport.
The normalized mean values of T 1 , T 2 , T 3 , T 4 , (À D 2 ), and f(D) conditional on c for cases A-E are shown in Figure 2 for Δ ≈ 0.4δ th and Δ ≈ 2.8δ th . Figure 2 shows that T 2 and (À D 2 ) act as source and sink, respectively, in all cases, which is consistent with previous findings [21,28]. The contribution of T 4 is positive for a major portion of the flame brush before becoming negative toward the burned gas side for Δ � δ th (e.g., Δ ≈ 0.4δ th ); however, for Δ > δ th (e.g., Δ ≈ 2.8δ th ) the contribution of T 4 remains a leading-order source throughout the flame brush. The term T 3 assumes positive values toward the unburned gas side of the flame brush before assuming mostly negative values for the major part of the flame brush in cases D and E, whereas T 3 is negative throughout the flame brush in cases A-C for all filter widths. The contribution of f(D) is negative (positive) toward the unburned (burned) gas side of the flame brush for all cases and for all filter widths. The magnitude of T 1 is negligible compared with T 2 , T 3 , T 4 , (À D 2 ), and f(D) for all Δ in all cases. It can be seen from Figure 2 that the magnitude of all the terms decrease with increasing Le and Δ, which is consistent with previous findings based on DNS data [22,28]. The observed behaviors of T 1 , T 2 , T 3 , T 4 , (À D 2 ), and f(D) in response to Le and Δ have recently been explained by Gao et al. [28] using a detailed scaling analysis, and the scaling estimates of the filtered SDR and the unclosed terms of the SDR transport equation are provided in Table 3. It is worth noting that m and n in Table 3 are positive numbers with magnitudes greater than unity, and the functions g(Le), φ(Le), φ 1 (Le), and Ψ (Le) increase with decreasing global Lewis number Le. It can be seen from the scaling estimates in Table 3 that the magnitudes of the terms T 1 , T 2 , T 3 , T 4 , (À D 2 ), and f(D) are expected to increase with decreasing filter width and global Lewis number. Interested readers are referred to [28] for further discussion on the derivation of the scaling estimates of T 1 , T 2 , T 3 , T 4 , (À D 2 ), and f(D), and only the modeling of these terms will be discussed in this paper in the following subsections.

Modeling of the turbulent transport term T 1
Equation (4i) indicates that the turbulent transport term T 1 could be satisfactorily closed if the subgrid flux of SDR (i.e., qu i N c À qũ iÑc ) is properly modeled. The sub-grid flux of SDR ðqu i N c À qũ iÑc Þ is often modeled using a gradient hypothesis as follows [34]: where D t is the sub-grid scale eddy diffusivity. It has been demonstrated earlier that the turbulent scalar flux of scalar gradients (e.g., flame surface density and SDR) may exhibit counter-gradient (gradient) transport for the flames when counter-gradient (gradient) transport is observed for ðqu i c À qũ ic Þ [20,22,23,62]. Thus, the modeling of T 1 needs to include both gradient and counter-gradient transports of ðqu i N c À qũ iÑc Þ.
The above expressions can be combined as ðqui cÀ qũicÞÑc Gao et al. [28] demonstrated that the unclosed term T 1 can be scaled in the following manner: where g(Le) is a function increasing with decreasing Le, which accounts for flame normal acceleration, S L is used to scale the sub-grid velocity fluctuations associated with sub-grid scalar gradients, and the sub-grid fluctuations of SDR are taken to scale with S L /δ L [28]. In Eq. (5ii), Da Δ ¼ ΔS L /u′ Δ δ th and Re Δ ¼ ρ 0 u′ Δ Δ/μ 0 are the sub-grid Damköhler number and sub-grid turbulent Reynolds number, respectively, with u 0 D ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi 2k sgs =3 p and k sgs ¼ ðqu i u i À qũ iũi Þ=2q being the subgrid turbulent velocity fluctuation and sub-grid kinetic energy, respectively. One obtains Da Δ Re Δ ∼ (Δ/δ th ) 2 using μ 0 ∼ ρ 0 S L δ th , which indicates that Da Δ Re Δ increases with increasing Δ. Alternatively, one obtains the following expression when the sub-grid velocity fluctuations are taken to scale with u′ Δ [28]: It is worth noting that the scaling estimate given by Eq. (5ii) (Eq. (5iii)) is more appropriate for counter-gradient (gradient) transport. Equations (5ii) and (5iii) can be combined to obtain the following scaling estimate, which is valid for both gradient and counter-gradient transport [28]: Gao et al. [53] have recently extended a RANS model proposed by Chakraborty and Swaminathan [51] for the purpose of modeling ðqu i N c À qũ iÑc Þ for the unity Lewis number flames in the context of LES in the following manner: where M i ¼ À qc=qx i Þ ð =jrcj is the ith component of the resolved flame normal vector for LES, Φ′ ¼ 0.7 is a model parameter, and the following values have been suggested for γ 1 , γ 2 , and C F [53]: The parameterization given by Eq. (6ii) ensures that γ 2 assumes an asymptotic value for large values of Re Δ (i.e., Re Δ → ∞ ). In Eq. (6i), the first term on the right-hand side principally accounts for the effects of flame normal acceleration due to heat release, whereas the last term on the righthand side of Eq. (6i) represents turbulent transport according to conventional gradient hypothesis. Moreover, the first and second terms on the right-hand side of Eq. (6i) are consistent with the scaling estimates given by Eqs. (5iv) and (5iii), respectively.
The predictions of J þ sg ¼ ðqu i N c À qũ iÑc ÞM i � d th =q 0 S 2 L according to Eq. (6i) with Φ′ ¼ 0.7 are compared to the corresponding quantity extracted from DNS data for Δ ≈ 0.4δ th , 1.6δ th , and 2.8δ th in Figure 3 for cases A-E. Figure 3 shows that even though Eq. (6i) predicts J þ sg in a reasonable manner in the cases with Le ≈ 1.0 (e.g., cases C-E), this model does not adequately capture the correct qualitative and quantitative behaviors of J þ sg for the flames with Le << 1.0 (i.e., cases A and B). The model given by Eqs. (6i) and (6ii) does not explicitly account for non-unity Lewis number effects, so it is not surprising that this model does not adequately capture the behavior of ðqu i N c À qũ iÑc Þ for Le << 1.0 flames where the nondimensional temperature T þ ¼ (T À T 0 )/(T ad À T 0 ) is significantly different from the reaction progress variable c, which alters the distribution of heat release and thermal expansion within the flame brush compared with the Le ≈ 1.0 flames. This behavior is mimicked here by introducing Le dependence of the model parameter Φ′ in the following manner: The predictions of the model given by Eq. (6i) with Φ′ according to Eq. (6iii) are also shown in Figure 3, which demonstrates that the model with new parameterization Φ′ ¼ 0.3(1 À Le) þ 0.7 predicts J þ sg satisfactorily for all filter widths in all cases considered here and the agreement between the predictions of Eq. (6i) and DNS data improves with increasing Δ (see Figure 3). The predictions of Eq. (6i) with Φ′ according to Eq. (6iii) become equal to the corresponding values obtained for Φ′ ¼ 0.7 for the Le ¼ 1.0 case, and these two predictions cannot be distinguished from each other for case D in Figure 3. It worth noting that the sub-grid flux of the reaction progress variable (i.e., qu i c À qũ ic ) requires modeling in LES, and the performances of the models for ðqu i N c À qũ iÑc Þ and the turbulent transport term T 1 depend on the modeling of ðqu i c À qũ ic Þ. The modeling of ðqu i c À qũ ic Þ is beyond the scope of the current analysis and interested readers are referred to recent investigations by Chakraborty and Cant [63] and Gao et al. [64] for further discussion on the modeling of turbulent scalar fluxes in premixed turbulent flames.

Modeling of the density variation term T 2
For unity Lewis number flames the gas density ρ can be expressed as ρ 0 /(1 þ τc) [33], which leads to an alternative expression for the density variation term T 2 as [22,47,48,51,53]: in the non-unity Lewis number flames because the equality between T þ and c no longer holds. Although T 2 ¼ 2ðqr �ũN c Þ does not strictly hold in non-unity Lewis number flames, the gas density can still be scaled as ρ ∼ ρ 0 /(1 þ τc); thus, the density variation term T 2 can be scaled for adiabatic flames with low Mach number as follows [28]: where m is a positive number greater than unity (i.e., m > 1). The resolved part of T 2 can be taken to scale as [28] where U ref is a velocity scale representing the Favre-filtered velocity components ũ i . It is worth noting that q ¼ q 0 =ð1 þ scÞ for unity Lewis number flames yields ðT 2 Þ res ¼ 2qDrc : rcðqũ i =qx i Þ; however, the expression q ¼ q 0 =ð1 þ scÞ does not strictly hold for non-unity Lewis number flames, but q and (T 2 ) res can still be scaled using ρ 0 =(1 þ sc) and 2qDrc : rcðqũ i =qx i Þ, respectively. The scaling estimates given by Eqs. (7i) and (7ii) demonstrate that T 2 remains of the order of q 0 sS 2 L =d 2 th irrespective of Δ. By contrast, the magnitude of (T 2 ) res remains comparable to q 0 S 2 L =d 2 th for U ref ∼ S L and Δ ≈ δ th ; however, the magnitude of (T 2 ) res is expected to decrease with increasing Δ. This suggests that the sub-grid component (T 2 ) sg ¼ T 2 À (T 2 ) res plays an increasingly important role with increasing Δ, which can be substantiated from Figure 4, where the variations of the mean values of T 2 and (T 2 ) sg ¼ T 2 À (T 2 ) res conditional on c are shown for cases A-E for Δ ≈ 0.4δ th , 1.6δ th , and 2.8δ th .
Gao et al. [53] recently proposed the following model T 2 for unity Le flames in the following manner: qq qx i |ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl {zffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl } ðT 2 Þ res þb T 2 sS L ½qÑ c À qDrc:rc� where Ka Δ ¼ (u′ Δ /S L ) 3/2 (Δ/δ th ) À 1/2 is local sub-grid Karlovitz number and b T 2 ¼ 2:7 is a model parameter. The first term on the right-hand side of Eq. (8) [22,28,47,48,51,52] as the combustion process is likely to show the attributes of the broken reaction zones regime [55] (where the effects of heat release are weak) for high values of Karlovitz number. The prediction of ) and (9) ( ) for Δ ≈ 0.4δ th (1st column), 1.6δ th (2nd column), and 2.8δ th (3rd column) in cases A-E (1st-5th row). All the terms are normalized with respect to q 0 S 2 L =d 2 th .
Eq. (8) is also shown in Figure 4 for cases A-E for Δ ≈ 0.4δ th , 1.6δ th , and 2.8δ th . A comparison between the predictions of Eq. (8) and the normalized T 2 extracted from explicitly filtered DNS data reveals that Eq. (8) satisfactorily predicts T 2 for a range of different filter widths for flames with Le ≈ 1.0 (e.g., cases C-E); however, this model significantly underpredicts the magnitude of T 2 for the Le << 1.0 cases (e.g., cases A and B). It can be seen from Eq. (7i) that the magnitude of T 2 is expected to increase with decreasing Le due to the strengthening of heat release effects as a result of the enhanced burning rate for the small values of Lewis number (see Table 2). As this effect is missing in Eq. (8), this model underpredicts the magnitude of T 2 for the Le << 1.0 cases (e.g., cases A and B), where the effects of enhanced heat release due to the differential diffusion of heat and mass are particularly strong.
Here the model given by Eq. (8) has been extended in order to account for the effects of Le in the following manner: qq qx i |ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl {zffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl } In Eq. (9ii) f T 2 ðLeÞ accounts for the strengthening of heat release effects with decreasing Le as suggested by the scaling estimate given by Eq. (7i). The parameter K � c is a thermo-chemical parameter, which provides information regarding the SDR-weighted dilatation rate r:ũ [34,36,65,66]. The thermo-chemical parameter K � c accounts for the correlation between r:ũ and ρN c within the flame front. It is possible to approximate f b (c) as f b (c) ¼ 1/|∇c| L [65,66], which enables one to evaluate K � c from laminar flame data. The thermo-chemical parameter K � c =s is also affected by Le and it is equal to 0.52, 0.67, 0.71, 0.78, and 0.79, respectively, for the Le ¼ 0.34, 0.6, 0.8, 1.0, and 1.2 flames considered here. The predictions of Eq. (9) are compared to the predictions of Eq. (8) and T 2 extracted from DNS data in Figure 4, which shows that Eq. (9) satisfactorily predicts the quantitative behavior of T 2 for a range of different values of Δ for flames with Le ranging from 0.34 to 1.2. Eq. (9) becomes exactly equal to Eq. (8) for the Le ¼ 1.0 case and thus the predictions of Eqs. (8) and (9) cannot be distinguished from each other for case D in Figure 4.

Modeling of the scalar turbulence interaction term T 3
The variations of the mean values of T 3 conditional on c are shown in Figure 5 for cases A-E at Δ ≈ 0.4δ th , 1.6δ th , and 2.8δ th . Figure 5 shows that T 3 assumes predominantly negative values throughout the flame brush for cases A-C, but this term exhibits weak positive values toward both the unburned gas sides of the flame brush before assuming mostly negative values for the major portion of the flame brush in cases D and E. The term T 3 can be expressed as follows [21,28,34,[45][46][47][48]: where e a , e β , and e γ are the most extensive, intermediate, and most compressive principal strain rates and their angles with∇c, respectively. Equation (10) suggests that a predominant collinear alignment of ∇c with e a (e γ ) leads to a negative (positive) value of T 3 . It was discussed elsewhere [21,28,34,[45][46][47][48] that ∇c predominantly aligns with e a , when the strain rate induced by the flame normal acceleration overcomes turbulent straining, whereas one obtains preferential alignment of ∇c with e γ when turbulent straining dominates over the strain rate due to flame normal acceleration. The flame normal acceleration strengthens with decreasing Le, and thus ∇c predominantly aligns with e a , for the Le << 1 flames (e.g., cases A and B), leading to negative values of T 3 [21,22,28]. By contrast, turbulent straining overcomes the flame normal acceleration on both ends of the flame brush for the Le ≈ 1.0 cases considered here (e.g., cases C-E), which leads to positive values of T 3 both on unburned and on burned ) and Eqs. (13i) and (13ii) ( ) for Δ ≈ 0.4δ th (1st column), Δ ≈ 1.6δ th (2nd column), and Δ ≈ 2.8δ th (3rd column) in cases A-E (1st-5th row). All the terms are normalized with respect to q 0 S 2 L =d 2 th .
gas sides of the flame brush for Δ ≈ 0.4δ th . However, the flame normal acceleration dominates over turbulent straining in the middle of the flame brush where the effects of heat release are strong even in the Le ≈ 1.0 cases considered here (e.g., cases C-E), which leads to negative values of T 3 for the major portion of the flame brush in these cases. The effects of ∇c alignment with e a , on T 3 can be scaled in the following manner [28]: The contribution of ∇c alignment with e γ on T 3 can be scaled as follows [28]: The Lewis number Le dependence in Eq. (11i) (with n > 1) accounts for the greater extent of ∇c alignment with e a for the flames with Le << 1.0. Gao et al. [28] proposed the following scaling estimate of the resolved part of T 3 : A comparison of Eqs. (11i)-(11iii) reveals that the contribution of (T 3 ) res to T 3 is expected to weaken with increasing Δ, and this behavior can indeed be seen from Figure 5, which shows that the magnitude of (T 3 ) res decreases with increasing Δ.
Gao et al. [53] utilized the scaling estimates given by Eqs. (11i) and (11ii) to propose a model for T 3 for Le ¼ 1.0 flames: where C 3 and C 4 are the model parameters and Da � L ¼ S L q 0 D=u 0 D qd th is the density-weighted local sub-grid Damköhler number. The symbol f T 3 is a bridging function in terms of ΔS L /α T0 , which ensures that (T 3 ) sg ≈ T 3 for Δ >> δ th and T 3 approaches (T 3 ) res when the flow is fully resolved: Gao et al. [53] proposed the following expressions for the model parameter C 3 , C 4 , and f T 3 : It is worth noting that the terms C 3 qðu 0 D =DÞÑ c and À C 4 q 0 sðS L =d th ÞÑ c are consistent with scaling estimates given by Eqs. (11i) and (11ii), respectively. However, a comparison between Eq. (11i) and À C 4 q 0 sðS L =d th ÞÑ c reveals that the increased alignment of ∇c with e α for small values of Le as a result of the strengthening of flame normal acceleration is not accounted for by the model given by Eq. (12ii). The effects of flame normal acceleration are expected to weaken with increasing Karlovitz number as the reacting flow field exhibits some attributes of passive scalar mixing for large values of Karlovitz number in the broken reaction zones regime [55]. This behavior is mimicked here by Ka Δ dependence of C 4 in Eq. (12iii).
The predictions of Eq. (12i) with the model parameters given by Eq. (12ii) are compared to T 3 extracted from DNS data in Figure 5, which shows that Eq. (12i) adequately captures the qualitative and quantitative behaviors of T 3 for the Le ≈ 1.0 cases considered here (e.g., cases C-E); however, this model has been found to underpredict the magnitude of the negative contribution of T 3 in the Le << 1.0 cases (e.g., cases A and B) for Δ > δ th . It has already been noted that the increased extent of scalar gradient destruction in the Le << 1.0 flames, due to the preferential alignment of ∇c with e α under strong actions of flame normal acceleration, is not addressed in the model given by Eq. (12i). Thus, this model underpredicts the negative contribution of T 3 for the flames with Le << 1.0. Here Eq. (12i) has been modified in the following manner to account for non-unity Lewis number effects: where CðLeÞ

Modeling of the combined reaction, dissipation, and diffusivity gradient contribution [T 4 À D 2 þ f(D)]
The variations of the mean values of [T 4 À D 2 þ f(D)] conditional on c are shown in Figure 6 for A-E for Δ ≈ 0.4δ th , 1.6δ th , and 2.8δ th . It can be seen from Figure 6 that [T 4 À D 2 þ f(D)] acts as a sink (source) term toward the burned (unburned) gas side of the flame brush for Δ ≈ 0.4δ th and Δ ≈ 1.6δ th ; however, the mean value of [T 4 À D 2 þ f(D)] conditional on c assumes predominantly negative values for Δ ≈ 2.8δ th . Table 3 shows that the order of magnitudes of T 4 , (À D 2 ), and f(D) remain comparable according to the scaling analysis by Gao et al. [28] and their magnitudes are expected to increase with decreasing Le. Furthermore, the scaling estimates of (T 4 ) res , (À D 2 ) res , and {f(D)} res in Table 3 suggest that their contributions are expected to weaken with increasing Δ, where (T 4 ) res , (À D 2 ) res , and {f(D)} res are the resolved components of T 4 , (À D 2 ), and f(D), which are given by Thus, the sub-grid components (T 4 ) sg ¼ T 4 À (T 4 ) res , (À D 2 ) sg ¼ À D 2 þ (D 2 ) res , and {f(D)} sg ¼ f(D) À {f(D)} res are expected to play major roles for Δ >> δ th . The aforementioned behaviors of the resolved and sub-grid components of T 4 , (À D 2 ), and f(D) can be confirmed from Figure 6. It can be seen from Table 3 that the magnitudes of (T 4 ) sg , (À D 2 ) sg , and {f(D)} sg remain of the order of q 0 S 2 L =d 2 th � qÑ 2 c for Δ >> δ th ; however, their magnitudes are expected to increase with decreasing Le, which can indeed be substantiated from Figure 6.
T 4 À D 2 þ f ðDÞ ¼ ðT 4 Þ res À ðD 2 Þ res þ ff ðDÞg res À ð1 À f TD Þb 3 ðc À c � Þq ½Ñ c ÀDrc:rc� 2 cð1 ÀcÞ ð15iÞ where b 3 ¼ 5:7; c � ¼ 1:0 À 0:83erf 0:5 DS L a T0 À 2:3 The involvement of ðc À c � Þ=½cð1 ÀcÞ� in Eq. (15i) is required for capturing the qualitative behavior of [T 4 À D 2 þ f(D)] across the flame brush, whereas f TD approaches unity for small values of Δ as the terms get fully resolved (i.e., lim D!0 ½T 4 À D 2 þ f ðDÞ� ¼ lim D!0 ½ðT 4 Þ res À ðD 2 Þ res þ ff ðDÞg res �). The transition from positive to negative contribution of [T 4 þ f(D) À D 2 ] with increasing Δ has been accounted for by c � . The predictions of Eq. (15i) are shown in Figure 6, which shows that this model captures both the qualitative and quantitative behaviors of [T 4 þ f(D) À D 2 ] for the Le ≈ 1.0 cases considered here (e.g., cases C-E); however, this model underpredicts the magnitude of [T 4 þ f(D) À D 2 ] significantly for the Le << 1.0 cases (e.g., cases A and B). It is worth noting that the model given by Eq. (15i) does not account for the increased magnitude of {T 4 À D 2 þ f(D)} sg for small values of Le (see Table 3); hence, perhaps it is not surprising that this model under-predicts the magnitude of [T 4 þ f(D) À D 2 ] for the flames with Le << 1.0 (e.g., cases A and B). The increased magnitude of [T 4 þ f(D) À D 2 ] for the small values of Le is accounted for by modifying Eq. (15i) in the following manner: where c � and f TD are considered according to Eq. (15ii). The predictions of Eq. (16) are shown in Figure 6, which demonstrates that Eq. (16) captures both the qualitative and quantitative behaviors of [T 4 þ f(D) À D 2 ] for a range of filter widths for different Le cases considered here. Equations (15i) and (16) become equal to each other for Le ¼ 1.0 and thus their predictions cannot be separated from each other in case D in Figure 6. It is worth noting that the combined contribution of the terms D 1 , T 4 , f(D), and (À D 2 ) can be expressed in the following manner if the SDR transport equation is derived based on the kinematic form of the progress variable transport equation (i.e., Dc/Dt ¼ S d |∇c|) [34,47]: where S d ¼ ½ _ w þ r � ðqDrcÞ�=ðq rc j jÞ and ñ ¼ À rc= rc j j are the flame displacement speed and local flame normal vector, respectively. Thus, Eq. (17) suggests that the net contribution of [T 4 À D 2 þ f(D)] originates due to flame normal propagation and flame curvature. This justifies the modeling of these terms together [21,34,44,47,48,53] because the molecular diffusion term D 1 is a closed term. Although Eq. (16) reasonably captures the qualitative behavior and the magnitude of [T 4 þ f(D) À D 2 ] for all cases considered here, the collective modeling of the terms T 4 , f(D), and (À D 2 ) may give rise to the loss of their individual physical significances. However, this is one of the first attempts to model the Lewis number effects on the SDR transport equation terms in the context of LES of premixed combustion and thus there is a scope for further improvement of this modeling in the future.

Implications of model implementation
The newly proposed models for the unclosed terms of the SDR Ñ c transport equation are summarized in Table 4 for the future potential users of these models. It is worth noting that the flamelet assumption is invoked while deriving these models, so they are expected to remain valid in the corrugated flamelets and thin reaction zones regimes of turbulent premixed combustion [55]. The scaling estimates in Table 3 indicate that the terms T 2 , T 3 , T 4 , (À D 2 ), and f(D) remain leading-order contributors to the SDR Ñ c transport and the magnitude of T 1 remains negligible compared with the terms T 2 , T 3 , T 4 , (À D 2 ), and f(D) irrespective of Damköhler and turbulent Reynolds numbers. This is consistent with the observations made from Figure 2. However, the turbulent transport term T 1 still needs to be modeled and included in the model implementation for LES for numerical stability.

Conclusions
The effects of global Lewis number Le on the modeling of the unclosed terms of the transport equation of Favre-filtered SDR Ñ c have been analyzed based on a priori analysis of a DNS database of freely propagating statistically planar turbulent premixed flames with Le ranging from 0.34 to 1.2. It has been found that Le has profound influence on the statistical behavior of the unclosed terms of Ñ c transport arising from turbulent transport T 1 , density variation due to heat release T 2 , alignment of scalar and velocity gradients T 3 , correlation between the gradients of reaction rate and reaction progress variable T 4 , molecular dissipation (À D 2 ), and diffusivity gradients f(D), and detailed physical explanations have been provided for the observed non-unity Lewis number effects. Recently proposed models for T 1 , T 2 , T 3 , T 4 , (À D 2 ), and f(D) for unity Lewis number flames have been extended here to account for the effects of Le based on the scaling estimates of these unclosed terms [28]. The newly proposed models have been found to satisfactorily predict the unclosed terms obtained from explicitly filtered DNS data for a range of Δ for different values of Le. However, it is still essential to implement these models into actual LES simulations for the purpose of a posteriori assessment. Moreover, these models need to be further validated based on detailed chemistry-based DNS simulations. Further validation of these models will form the basis of future investigations. Table 4. Summary of the proposed models for the unclosed terms of the SDR Ñ c transport equation (Eq. (3)) in this analysis.

Funding
The financial assistance of the Engineering and Physical Science (EPSRC) research council of the UK (EP/I028013/1) and computational support of N8 are gratefully acknowledged.