Modelling of the Tangential Strain Rate Term in the Flame Surface Density Transport Equation in the Context of Reynolds Averaged Navier Stokes Simulations : A Direct Numerical Simulation Analysis

A direct numerical simulation (DNS) database of freely propagating statistically planar turbulent premixed flames with a range of different values of Karlovitz number Ka, turbulent Reynolds number Re t , heat release parameter τ, and global Lewis number Le has been used to assess the models of the tangential strain rate term in the generalised flame surface density (FSD) transport equation in the context of Reynolds averaged Navier Stokes (RANS) simulations.The tangential strain rate term has been split into contributions arising due to dilatation rate T D and flame normal strain rate (−T N ). Subsequently, T D and (−T N ) were split into their resolved (i.e., T D1 and (−T N1 )) and unresolved (T D2 and (−T N2 )) components. Detailed physical explanations have been provided for the observed behaviours of the components of the tangential strain rate term. This analysis gave way to the modelling of the unresolved dilatation rate and flame normal strain rate contributions. Models have been identified for T D2 and (−T N2 ) for RANS simulations, which are shown to perform satisfactorily in all cases considered, accounting for the variations in Ka, Re t , τ and Le. The performance of the newly proposed models for the FSD strain rate term have been found to be either comparable to or better than the existing models.


Introduction
Flame surface density (FSD) based reaction rate closure is one of the most well-established methods of turbulent premixed combustion modelling in the context of Reynolds averaged Navier Stokes (RANS) simulations [1][2][3][4][5][6][7][8][9][10][11][12].In the FSD based formulation, the closure of reaction rate translates to the modelling of flame surface area to volume ratio [13].The FSD, in RANS simulations, is either evaluated using an algebraic expression [1] or a modelled transport equation for FSD is solved alongside other modelled Favre averaged conservation equations [3][4][5][6][7][8][9][10][11][12].However, algebraic closures may not be suitable when flame surface area generation and destruction are not in equilibrium and recently Richard et al. [14] demonstrated that algebraic closures may not be adequate for simulating in-cylinder processes in spark ignition piston engines.Therefore, it may be necessary to consider high-fidelity transport equation based closures of FSD for turbulent premixed flames.The flame surface area generation/destruction due to fluid straining plays a key role in the FSD transport and thus the statistical analysis and detailed assessments of the modelling of the tangential strain rate term in the FSD transport equation in the context of RANS simulations have been considered in the current work because the strain rate term remains a leading order contributor to the FSD transport [3][4][5][6][7][8][9][10][11][12].
The generalised FSD Σ gen is defined as [15,16] where  is the reaction progress variable and the quantities  and Q = / denote the Reynolds averaged and Favre averaged values of a general quantity .The generalised FSD Σ gen is often used in order to close the mean reaction rate ω of reaction progress variable in the following manner [15]: where  is the gas density,  is the reaction progress variable diffusivity,   = (/)/|∇| is the displacement speed, and ()  = |∇|/|∇| indicates a surface averaging operation of a general quantity  [5,15].In the context of RANS simulations ω ≫ ∇ ⋅ (∇) and therefore it is evident from (2) that an accurate prediction of Σ gen is required for the closure of the mean reaction rate ω provided the statistical behaviour of the surface averaged density-weighted displacement speed (i.e., (  )  ) is satisfactorily accounted for.The exact transport equation for the generalised FSD (i.e., Σ gen = |∇|) takes the following form [9,11,12,17,18]: where ũ =   / and    =   − ũ are the Favre mean and fluctuating velocity components in the th direction and ⃗  = −∇/|∇| is the local flame normal vector (which points towards the reactants).The reaction progress variable  can be defined in terms of a suitable reactant mass fraction   in such a manner that it increases monotonically from zero in unburned gases to unity in fully burned products.In (3), the first term in the left hand side is the transient term, whereas the second term represents the effects of mean advection.All the terms on the right hand side of (3) are unclosed and, hence, require modelling.The first term on the right hand side of (3) represents the effects of turbulent transport, the second term accounts for the flame area generation due to tangential strain rate (i.e.,   = (  −     )  /  ), and the third term originates due to flame normal propagation and thus is commonly referred to as the propagation term.The final term on the right hand side of (3) (i.e., (  ∇ ⋅ ⃗ )  Σ gen ) arises due to curvature   = ∇ ⋅ ⃗ /2 (according to this convention a flame element, which is convex towards the reactants, has a positive curvature and vice versa) and thus is commonly referred to as the FSD curvature term.In the current study, only the a priori DNS assessment of the models of the tangential strain rate term of the generalised FSD transport equation will be considered.
The tangential strain rate term of the generalised FSD transport equation is often split in the following manner: Cant et al. [2] and Candel et al. [3] proposed models for both the resolved and the unresolved parts of the tangential strain rate contribution (i.e.,   and  UR , resp.), which were subsequently assessed by Duclos et al. [4] and Prasad and Gore [8] based on RANS simulations.Veynante et al. [6] and Veynante et al. [7] assessed the performance of the models of   and  UR based on experimental data.The aforementioned studies [2][3][4][6][7][8] were carried out in the corrugated flamelets regime, where the flame thickness remains smaller than the Kolmogorov length scale.Recently, Katragadda et al. [10] assessed the performances of the existing   and  UR models in both the corrugated flamelets and thin reaction zones regimes of premixed combustion based on direct numerical simulation (DNS) data, while accounting for the effects of Lewis number and heat release parameter .It was shown by Katragadda et al. [10] that both the dilatation rate   /  and the relative alignment of ∇ with the fluid-dynamic strain rate   = 0.5(  /  +   /  ) significantly affect the behaviour of the strain rate term (  )  Σ gen .Global Lewis number Le is known to have significant influences on both the dilatation rate and scalar gradient alignment statistics in turbulent premixed flames [19].Additionally, it was shown in recent studies [19][20][21] that the Damköhler number Da significantly affects the relative strength of the dilatation rate in comparison to the turbulent straining and alignment characteristics of ∇ with the local principal strain rates.As the turbulent Reynolds number scales as Re  ∼ Ka 2 Da 2 [22] it is expected that the modelling of the strain rate term will also show some Re  dependence.
The effects of Lewis number Le, turbulent Reynolds number Re  , Damköhler number Da, Karlovitz number Ka, and heat release parameter  on the statistical behaviour and RANS modelling of the tangential strain rate term in the FSD transport equation are yet to be addressed in detail and the present paper aims to address this gap in the existing literature.In order to investigate the effects of Le, Re  , Da, Ka, and  on the statistical behaviour and RANS modelling of the FSD tangential strain term, a database of three-dimensional compressible direct numerical simulations (DNS) of turbulent premixed flames has been considered.In this respect, the main objectives of this study are as follows: (1) to demonstrate and explain the statistical behaviours of the tangential strain rate term of the FSD transport equation and its components for a range of different values of Lewis numbers, turbulent Reynolds numbers, Karlovitz and Damköhler numbers, and heat release parameters spanning both the corrugated flamelets and thin reaction zones regimes of turbulent premixed combustion; (2) to identify models which satisfactorily capture the statistical behaviour of the tangential strain rate term in the FSD transport equation and its components for different values of Le, , and Re  and different combustion regimes in the context of RANS simulations.
The remainder of the paper will take the following form.
The next section will discuss the relevant mathematical background of the current study, which will then be followed by a brief discussion on the numerical implementation.Following this, the results will be presented and subsequently discussed.Finally, the main findings will be summarised and conclusions will be drawn.

Mathematical Background
Direct numerical simulations (DNS) of turbulent combustion should, ideally, account for both the inherent threedimensionality of turbulence and the detailed chemical structure of the flame.However, until recently combustion DNS studies were limited to be carried out either in two dimensions with detailed chemistry or in three dimensions with simple chemistry.Three-dimensional DNS of turbulent combustion is now possible but such simulations remain extremely computationally expensive [23] and thus are often not suitable for a detailed parametric analysis as carried out in this work.As such, a three-dimensional DNS database with simple chemistry has been considered here for the sake of computational economy.In the context of simple chemistry, the species field can be represented, uniquely, by a reaction progress variable  which can be defined in terms of a suitable mass fraction (i.e.,   ) in the following manner: where subscripts 0 and ∝ were used to denote the values in the unburned and burned gases, respectively.In globally adiabatic, low Mach number unity Lewis number flames the reaction progress variable  is the same as the nondimensional temperature: where T is the instantaneous dimensional temperature,  0 is the unburned gas temperature, and  ad is the adiabatic flame temperature.However, for nonunity Lewis number flames the nondimensional temperature  may locally assume superadiabatic values (i.e.,  > 1) [24,25] even under globally adiabatic conditions, whereas  is always bound within 0 and unity (i.e., 0 ≤  ≤ 1).
As discussed in Section 1, the strain rate term (  )  Σ gen is often modelled by splitting it into the resolved   and the unresolved  UR components.For the purpose of evaluating   , the quantity (    )  requires modelling.Cant et al. [2] modelled   as where   is the model expression for (    )  , which takes the following form according to Cant et al. [2]: Veynante et al. [7] modelled (    )  as where k = ũ     /2 is the turbulent kinetic energy.Cant et al. [2] modelled for the unresolved part of the strain rate term  UR in the following manner: where  0 and  0 are the unburned gas density and viscosity, respectively, and ε is the dissipation rate of k.In the context of coherent-flamelet modelling (CFM) [4]  UR is modelled as where  0 is a model constant of the order of unity (i.e.,  0 = 2.0) and Γ  is the efficiency function proposed by Meneveau and Poinsot [26] which is a function of   ⋅   / 0 and √ (2 k/3)/  with  0 ,   , and   being the thermal diffusivity in the unburned gas, unstrained laminar burning velocity, and local integral length scale, respectively.The strain rate term can alternatively be decomposed in the following manner [10]: The terms   and (−  ) represent contributions of dilatation rate and flame normal strain rate on the Σ gen transport.The dilatation rate term   can be split into the resolved and unresolved components in the following manner: From the expression above, it is evident that the resolved dilatation rate term  1 can be closed if a suitable relationship Mathematical Problems in Engineering between  and c can be found.The contribution of the normal strain rate term (−  ) can be split as The quantity (− 1 ) can be closed using the models for (    )  .This can be carried out as shown in (8) and (9).It is worth noting that the surface averaged quantities involving velocity components in (4) are modelled in terms of the gradients of unconditional Favre-averaged velocity components ũ  /  (see (7), ( 9), (10), and ( 11)).In actual RANS simulations the unconditional Favre-averaged velocities are readily obtained so it makes sense to model the surfaceaveraged terms involving velocity components using the Favre-averaged velocity components [2][3][4][5][6][7][8].The modelling of (− 2 ) and  2 (see ( 12) and ( 14)) will be discussed based on a priori DNS analysis using the current database in Section 4, where the surface averaged quantities involving velocity components are also modelled in terms of the gradients of unconditional Favre-averaged velocity components ũ  /  (see (15) and ( 17) later in this paper).The present formulation is based on the generalised FSD Σ gen , which is defined as [15,16]: Σ gen = |∇|.This indicates that the generalised FSD Σ gen is not dependent upon choice of  isosurface.Thus the difference between conditional and unconditional velocities does not play any role in the current analysis similar to several previous studies [2][3][4][5][6][7][8][9][10][11][12].

Numerical Implementation
For the present study, the simulation parameters for the DNS database considered are shown in Table 1, namely, the initial values of root-mean-square turbulent velocity fluctuation normalised by the unstrained laminar burning velocity   /  , integral length scale normalised by the thermal flame thickness / th , turbulent Reynolds number Re  =  0   / 0 , heat release parameter  = ( ad −  0 )/ 0 , Lewis number Le, Damköhler number Da =   /   th , and Karlovitz number Ka = (  /  ) 3/2 (/ th ) −1/2 .Standard values are chosen for Zel' dovich number  =  ac ( ad −  0 )/ 2 ad , Prandtl number Pr, and the ratio of specific heats (i.e.,  = 6.0,Pr = 0.7, and  =   /  = 1.4),where  ac is the activation temperature.It should be noted that the simulation parameters have been chosen in such a manner that combustion takes place within the corrugated flamelets regime (thin reaction zones regime) [22] for case A (cases B-L).The turbulent flow conditions and heat release parameter are kept unaltered but the global Lewis number Le was modified from 0.34 to 1.2 in cases C-G.The cases B and F are identical in terms of turbulent flow conditions but only  values are different between these cases.As the turbulent Reynolds number Re  scales as Re  ∼ Da 2 Ka 2 [22], the variation of Re  in cases H-L is brought about by modifying Da and Ka independently of each other.In cases H, J, and L (I, J, and K) Da (Ka) is kept unaltered and Ka (Da) is altered to bring about the variation of Re  .The range of Re  considered here is comparable to that of  [28].In case A, inlet and outlet boundaries are specified in the mean direction of flame propagation (which is aligned with the negative  1 -direction), whereas transverse boundaries are taken to be periodic.In cases B-G, a domain of size 42.18 0 /  × 42.18 0 /  × 42.18 0 /  is discretised using a uniform grid of 230 × 230 × 230.In cases H-L, a domain of size 63.37 0 /  × 42.18 0 /  × 42.18 0 /  is discretised using a uniform grid of 345 × 230 × 230.The domain boundaries in the direction of mean flame propagation in cases B-L are taken to be partially nonreflecting and the transverse boundaries are assumed to be periodic.In case A, a 6th order central-difference scheme has been used for spatial discretisation for the internal grid points in the direction of mean flame propagation, which gradually reduces to a one-sided 4th order scheme near nonperiodic boundaries, whereas a spectral method is used for spatial discretisation in the directions normal to the mean direction of flame propagation [28].In cases B-L, a 10th central difference scheme is used for the internal grid points and the order of differentiation gradually reduces to a 2nd order one-sided scheme near nonperiodic boundaries.The time advancement for all viscous and diffusive terms in case A is carried out using an implicit solver, whereas the convection terms in case A and all the terms in cases B-L are time advanced with the help of a third order low storage Runge-Kutta method [29].For all cases, the flame is initialised by a steady unstrained planar laminar flame solution and the turbulent fluctuating velocity field is initialised based on an incompressible homogeneous isotropic velocity distribution, which is generated using a standard pseudospectral method [30].The grid spacing is determined by the resolution of the flame structure, and about 10 grid points are kept within  th for all cases considered here.In all cases flameturbulence interaction takes place under decaying turbulence and simulations need to be carried out for  sim ≥ Max(  ,   ), where   = /  is the initial eddy turn-over time and   =  th /  is the chemical time scale.The simulation in case A was run for about 4  , whereas in cases B-G it was run for a time equivalent to 3.34  .The simulation time remains either greater than (case A) or equal to (cases B-L) one chemical time scale.The simulation time remains either comparable to or greater than several previous studies [9,[31][32][33][34][35].The turbulent kinetic energy and its dissipation rate in the unburned gas ahead of the flame were not varying significantly with time when statistics were extracted for all cases.Interested readers are referred to Chakraborty and Cant [11,12] and Chakraborty et al. [36,37] for further information on this database and the conditions under which statistics were extracted.
The relevant quantities in all cases considered here have been Reynolds/Favre-averaged by ensemble averaging the relevant quantities in transverse directions ( 2 - 3 planes).The statistical convergence of the Reynolds/Favre-averaged quantities is examined by comparing the corresponding values obtained using half of the sample size from the distinct half of the domain in the transverse direction with those obtained based on full sample size.Both the qualitative and quantitative agreements between these sets of values are found to be satisfactory and more information can be obtained from [11].In the next section, only the results obtained based on full sample size will be presented.For statistically planar flames, c remains a unique function of the spatial coordinate in the direction of mean flame propagation (i.e.,  1 -direction); thus all the statistics in Section 4 will be presented as a function of c.

Results and Discussion
The contours of the reaction progress variable  in the central  1 - 2 plane when the statistics are extracted are shown in Figure 1 for cases A-L.In case A the contours of  remain parallel to each other as a consequence of the corrugated flamelets regime combustion.In the corrugated flamelets regime the energetic turbulent eddies are unable to penetrate into the flame and wrinkles develop due to largescale turbulent motion.However, in cases B-L the contours of c are found to be distorted towards the unburned gas side (i.e.,  < 0.5) representing the preheat zone.This distortion of the preheat zone in cases B-L is characteristic of the thin reaction zones regime combustion.Additionally the contours of  representing the reaction zone (i.e., 0.7 ≤  ≤ 0.9) remain parallel to each other for all cases.In cases B-L, where combustion takes place in the thin reaction zone regime, the energetic turbulent eddies penetrate into the flame, as the flame thickness remains greater than the Kolmogorov length scale though the reaction zone retains its quasi-laminar structure due to  >   where  and   are the Kolmogorov length scale and the reaction zone thickness, respectively.
The effects of Lewis number on the flame structure can be seen in Figures 1(c)-1(g), where the extent of the flame wrinkling increases as Le decreases.The reactants diffuse more rapidly into the reaction zone than the rate at which heat  throughout the flame brush.It can be seen from Table 2 that the effects of heat release strengthen with decreasing Le so the magnitude of ∇ ⋅ ⃗  increases significantly with decreasing Le, which along with the increasing trend of |∇| with decreasing Le [11] leads to high magnitudes of   for flames with small values of Le (see Figures 2(c)-2(g)).
It can be seen from Figure 2 that the normal strain rate term (−  ) exhibits marked differences from one case to another.The term (−  ) assumes negative value throughout the flame brush for cases A, C, D, E, H, I, J, and K, whereas this term assumes a positive (negative) value towards the unburned (burned) gas side of the flame brush in cases B, F, G, and L. The behavior of   = (      /  |∇|) = (  cos 2  +   cos 2  +   cos 2 )|∇| is determined by the alignment of ∇ with the local principal strain rates where   ,   , and   are the most extensive, intermediate, and compressive principal strain rates, respectively, and , , and  are the angles between |∇| and   ,   , and   , respectively.For passive scalars, the scalar gradient is known to align with the most compressive principal strain rate   [28,40,41] but recent studies [19][20][21]42]   where the strain rate induced by flame normal acceleration overcomes the effects of turbulent straining.The strain rate induced by flame normal acceleration can be scaled as  chem ∼ (Ka)  / th for Le = 1.0 flames, whereas the turbulent straining can be scaled as  turb ∼   /, which gives rise to  chem / turb ∼ (Ka)  /   th ∼  ⋅ (Ka)Da [20] with (Ka) being a function which accounts for the weakening of the effects of flame normal acceleration with increasing Ka as the broken reaction zones regime combustion is approached.Alternatively, the turbulent straining  turb can be scaled as  turb ∼   / following Tennekes and Lumley [43] where  is the Taylor microscale, which leads to  chem / turb ∼  ⋅ (Ka)Da/Re 1/2  .The scaling  chem /  turb ∼ (Ka)  /   th ∼ ⋅(Ka)Da ( chem / turb ∼ ⋅(Ka) Da/Re 1/2  ) indicates that ∇ may align with   to yield a positive (negative) contribution of   (−  ) for large Da flames.By contrast, in cases B, F, G, and L, Da remains small (i.e., Da < 1), so that the effects of  turb overcome the effects of  chem , giving rise to the predominant alignment of ∇ with   for the major portion of the flame brush, except in the heat releasing zone where the effects of  chem overcome the effects of  turb .This is reflected in the predominantly negative (positive) value of   (−  ) towards the unburned gas side of the flame brush, and positive (negative) values of   (−  ) are obtained towards the burned gas side of the flame brush in cases B, F, G, and L. Chakraborty et al. [19] demonstrated that  chem strengthens with decreasing Le and this effect is particularly strong for the Le ≪ 1 flames due to high values of chemical heat release.As a result,  chem / turb can be scaled as  chem / turb ∼  ⋅ (Le, Ka)⋅Da ( chem / turb ∼ (Le, Ka)Da/Re 1/2  ) for nonunity Lewis number flames, where the function  increases with decreasing Le [19].This suggests that  chem may become sufficiently strong to overcome the effects of  turb for small values of Le, even for Da < 1 flames.Strong  chem in the Le = 0.34, 0.6, and 0.8 flames overcomes the effects of  turb and induces a preferential alignment of ∇ with   .This gives rise to positive (negative) values of   (−  ) for the major portion of the flame brush.As the effects of  chem weaken with an increase in Le, the extent of the negative contribution of (−  ) progressively decreases with increasing Le, which can be substantiated from Figures 2(c)-2(g).Figures 2(a)-2(l) show that the magnitude of (−  ) remains smaller than the magnitude of   for all cases, which yields a positive contribution of (  )  Σ gen throughout the flame brush even when the contribution of (−  ) remains negative.In all cases,  turb remains of the same order of magnitude as that of  chem (i.e.,  chem / turb ∼ (Le, Ka) ⋅  ⋅ Da ∼  (1) or  chem / turb ∼ (Le, Ka) ⋅ Da/Re 1/2  ∼ (1)) and thus the effects of  turb are partially nullified by the effects induced by  chem , which give rise to a smaller magnitude of the normal strain rate contribution (−  ) than the magnitude of the dilatation rate term   .
It is worth noting that the variations of some of the strain rate terms with c in Figure 2 and subsequent figures show wavy behaviour although the variation of  with c for all cases has been found to be reasonably smooth.However, the variation of  with c exhibits weak undulations for some cases (e.g., cases A, C, K, and L), which originate due to large extent of wrinkling in these flames.These undulations in  variations with c give rise to wavy behaviour of |∇| (=|/ 1 | for statistically planar flames) when a 10th order central difference scheme is used for spatial differentiation.This also contributes to the wavy behaviour of the variations of the strain rate terms with c.Similar wavy behaviour has been seen in several previous DNS studies [9,10,20,21,28,42].This is not due to a lack of statistical convergence but is a consequence of the high levels of wrinkling in these flames and high order of finite difference scheme used for postprocessing purpose.Moreover, this waviness is prevalent only for case A, which was not simulated by the present authors but taken from a widely used well-respected DNS database (see [28]).This waviness is also observed in several previous analyses [20,21,28,42,44], where case A was used.
Katragadda et al. [10] modelled the tangential strain rate term based on the separate modelling of the strain rate components   and (−  ).The term  1 = ũ  /  |∇| can be closed if a suitable relationship can be established between c and .According to Bray-Moss-Libby (BML) analysis [45], this can be achieved by  = (1 + )c/(1 + c) + () for unity Lewis number globally adiabatic flames, where () accounts for the contribution of the reacting mixture.The contribution of () remains negligible for high Da flames but it might be nonnegligible for low Da combustion.Katragadda et al. [10] proposed an alternative empirical expression  = (1 +  1.5 Le −0.26 )c/(1 +  1.5 Le −0.26 c) where  = c 2 /c(1 − c) is the segregation factor and it was shown in [10][11][12] that (1 +  1.5 Le −0.26 )c/(1 +  1.5 Le −0.26 c) adequately captures  for low Da flames considered in this analysis and thus predictions of  = (1+ 1.5 Le −0.26 )c/(1+ 1.5 Le −0.26 c) are not shown here for the sake of conciseness.Moreover, (1 +  1.5 Le −0.26 )c/(1 +  1.5 Le −0.26 c) approaches the BML expression  = (1+)c/(1+ c) for unity Lewis number flames for Da ≫ 1 combustion under which  = c 2 /c(1 − c) approaches unity.Thus  = (1 +  1.5 Le −0.26 )c/(1 +  1.5 Le −0.26 c) enables one to close  1 = (ũ  /  )|∇|.The unresolved part of the dilatation rate term  2 is taken to model as [10] where the model parameters  2 = 1.8,  = 0.35,  1 = 1.8Le − 1.5, and  2 = 1.845.In (15) dilatation rate is scaled as ∇ ⋅ ⃗  ∼  ⋅   / th Le  2 as the strength of the dilatation rate increases with decreasing Le, and (1 − c)  1 is used to adequately capture the variation of  2 with c obtained from DNS where the exponent  1 is likely to be a function of Le, as the variation of  2 with c is skewed with a peak towards the unburned gas side for Le = 0.34 and 0.6 flames, whereas the peak value of  2 is attained close to the middle of the flame brush (i.e., c = 0.4) for the flames with Le close to unity.According to (15), the contribution of  2 vanishes when the flame is fully resolved (i.e., Σ gen = |∇|) and the net contribution of ( 1 +  2 ) becomes identical to (ũ  /  )|∇|.The contribution of  2 is similar to the dilatation rate contribution   = 2(∇ ⋅ ⃗ )∇  ⋅ ∇  in the scalar dissipation rate ε = ∇  ⋅ ∇  / transport equation [42].The term   = 2(∇ ⋅ ⃗ )∇  ⋅ ∇  is modeled as [42] In ( 16),   is taken to be   =   /(1 + Ka  ) 1/2 with   being a constant of the order of unity and Ka  = (ε th ) 1/2 / 3/2  is the local Karlovitz number.The local Karlovitz number Ka  dependence of   ensures that the strength of the dilatation rate contribution to the scalar dissipation rate ε transport diminishes with increasing Ka  as combustion starts to exhibit the attributes of the broken reaction zones regime [26].Following the same procedure, the local Karlovitz number Ka  dependence is used in (15) in the form of 1/(1 + Ka  )  .However, this local Karlovitz number dependence is one of the several possibilities and any other function of Ka  , which accurately predicts the quantitative behaviour of  2 and its diminishing strength with increasing Karlovitz number as the broken reaction zones regime is approached, can in principle be used for parameterising the local Ka  dependence of  2 .The numerical values of  2 , ,  1 , and  2 are optimised based on a least-squares method.
The predictions of ( 15) are shown in Figures 3(a)-3(l) along with  2 obtained from the DNS data for cases A-L, respectively.The model performs satisfactorily for a range of different values of Le, Re  , and  spanning both the corrugated flamelets and thin reaction zones regimes of premixed turbulent combustion.
The variations of (− 1 ) and (− 2 ) are shown in Figure 4 with c for cases A-L.The predictions for (− 1 ) according to the (    )  models proposed by Cant et al. [2] and Veynante et al. [7] are also shown in Figure 4 for cases A-L.The predictions of (− 1 ) according to Cant et al. [2] and Veynante et al. [7] models (denoted as TN1CPB and TN1V models, resp.) are shown in Figure 4, which shows that the performances of both of the models are comparable, and both satisfactorily predict (− 1 ) for all cases considered in the current analysis.However, the prediction of the TN1V model is closer to the DNS data than the TN1CPB model.The model given by ( 8) overpredicts the unresolved part of (    )  (i.e., (    )  − (  )  (  )  =   (1 − (  )  (  )  )) due to the assumption of isotropy of the unresolved flame normal fluctuations.However, experimental data [6,7] suggested that the assumption of isotropy of the unresolved flame normal fluctuations does not hold, and this result has been verified here by the overprediction of (− 1 ) by the TN1CPB model.It should be noted that (    )  reverts to     when the flame is fully resolved according to the TN1CPB model, but this condition is not satisfied by the model given by (9).
It is evident from (−  ) = −(  cos 2  +   cos 2  +   cos 2 )|∇| that a preferential alignment between ∇ and   (  ) leads to a negative (positive) contribution of (−  ).As the magnitude of (− 2 ) remains either greater than or comparable to (− 1 ) for all cases, it can be expected that ∇ alignment with local principal strain rates affects not only (−  ) but also (− 2 ).
The alignment of ∇ with local principal strain rates is determined by the competition between the strain rate due to turbulent motion  turb and the strain rate  chem induced by flame normal acceleration.The scalar gradient ∇ aligns with   (  ) when  chem ( turb ) dominates over  turb ( chem ).According to Tennekes and Lumley [43]  turb can be scaled as   / ∼ (  /)Re 1/2  ∼ (ε/ k)Re 1/2  , whereas  chem can be scaled as (Le, Ka)  / th .Thus, the effects of  turb on the term (− 2 ) through the preferential alignment of ∇ with   can be scaled as  1 (ε/ k)Σ gen , where the model parameter  1 is expected to be a function of local turbulent Reynolds number Re  =  0 k2 / 0 ε.By contrast, the effects of ∇ alignment with   on the term (− 2 ) due to the action of  chem can be scaled as − 2 ⋅ (Le)  Σ gen / th , where (Le) is a function of Le which increases with decreasing Le due to strengthening of flame normal acceleration as a result of augmented burning rate, and the model parameter  2 is expected to decrease with increasing Ka  because the effects of  chem are likely to weaken progressively with increasing Ka  .Combining the above scaling estimations, the following model for (− 2 ) can be proposed: where Da  = k  /ε th is the local Damköhler number and  1 is given as while the following two parameterisations are proposed for the model parameter  2 : where  in (19i) is given as and (Le) from (19i) and (19ii) is defined as follows: The parameterisations of  1 given by (18) Figure 4: Variations of (− 1 ) (solid line) and (− 2 ) (dash line) with c across the flame brush for cases ((a)-(l)) A-L along with the predictions of TN1CPB given by (8) (star line), TN1V given by (9) (box line), and the newly proposed models for (− 2 ) given by ( 17) using  2 given by (19i) (triangle line) and  2 given by (19ii) (cross line).All the terms in Figure 4 are normalised with respect to   / Figure 5: Variations of (  )  Σ gen (solid line) with c across the flame brush for cases ((a)-(l)) A-L along with the predictions of CPB (dash line), CFM (star line), and the newly proposed model given by ( 22) using  2 given by (19i) (square line) and  2 given by (19ii) (triangle line).All the terms in Figure 5 are normalised with respect to   / 2 th .
of  chem on the term (− 2 ) weaken with increasing Ka  , whereas (Le) according to (21) It should be noted that the involvement of (Σ gen − |∇|) in the first term on the right hand side of ( 22) and [1 − (  )  (  )  ] in  2 ensures that (  )  Σ gen becomes equal to (  −     )  /  |∇| when the flow becomes fully resolved (i.e., (  )  (  )  = 1.0).The predictions of (  )  Σ gen according to (22) based on the parameterisations given by (19i) and (19ii) are shown in Figures 5(a)-5(l) for cases A-L, respectively, along with the predictions of the models proposed by Cant et al. [2] (i.e., CPB model) and CFM methodology [4,6,7] where   and  UR are modeled using ( 8) and (10) in the context of the CPB model and by ( 9) and (11) in the context of the CFM model, respectively.The CFM model is found to show a peak at c ≈ 0.3 which is quantitatively similar to cases E-G and J-L, while in cases H and I the peak value of (  )  Σ gen occurs towards the products (i.e., c ≈ 0.6) and in cases C and D the peak value of (  )  Σ gen occurs towards the reactants (i.e., c ≈ 0.2).Although the CFM model performs adequately for cases F, G, and J, it fails to capture the behavior of (  )  Σ gen in cases H, I, K, and L. From Figure 5, it is evident that the model CPB can be seen to predict the correct qualitative behavior for most of the cases but in cases H and I it fails to capture the correct magnitude of (  )  Σ gen .
The new model performs well for all cases considered here, and the difference between the two parameterisations for  2 can be seen to make very little difference.The newly proposed model given by ( 22) with  2 given by (19i) and (19ii) also captures the behaviours of (  )  Σ gen satisfactorily for a range of different values of Le, Re  , and  spanning both the corrugated flamelets and thin reaction zones regimes of premixed turbulent combustion.

Conclusions
A DNS database of freely propagating statistically planar turbulent premixed flames with a range of different values of turbulent Reynolds number Re  , heat release parameter, and global Lewis number Le spanning both the corrugated flamelets and thin reaction zones regimes of premixed combustion has been used to assess the modelling of the tangential strain rate term (  )  Σ gen of the generalised FSD transport equation.In order to propose a model for (  )  Σ gen , which will be valid for both the corrugated flamelets and thin reaction zones regimes, the statistical behaviours of the strain rate contributions to the FSD transport due to the dilatation rate and flame normal strain rate have been analysed separately.
It has been shown that the dilatation rate contribution   to the generalised FSD transport equation strengthens with decreasing Lewis number.The contribution of normal strain rate term (−  ) assumes predominantly negative values for cases with Le ≪ 1.0 and Da ≫ 1 where the effects of straining due to flame normal acceleration  chem dominate over turbulent straining  turb , while in the Le ≈ 1.0 cases with Da < 1 the normal strain rate term (−  ) assumes positive (negative) values towards the reactants (products) side of the flame brush.This analysis gave rise to the explicit modelling of the unresolved dilatation rate and flame normal strain rate contributions.The predictions of the new model have been compared with the predictions of the existing models (i.e., CPB and CFM models) and the performance of the newly developed model has been found to be either comparable to or better than the existing models for the FSD strain rate term.
The present analysis has been carried out using a simple chemistry based DNS database with moderate range of values of turbulent Reynolds number Re  .Although the model parameters have been found to attain asymptotic values for Re  ≥ 50, further validations will be required based on experimental and three-dimensional detailed chemistry DNS data at higher values of Re  .Finally a posteriori assessment of the new models needs to be carried based on actual RANS simulations in a configuration for which well-documented experimental data is available, which will form the basis of future communications.

Figure 2 :
Figure 2: Variations of the values of (  )  Σ gen (solid line),   (dash line),  UR (star line),   (box line), and (−  ) (triangle line) with c across the flame brush for cases ((a)-(l)) A-L, where all the terms are normalised with respect to   / 2 th .

Figure 3 :
Figure 3: Variations of  2 (solid line) and (15) (dash line) with c across the flame brush for cases ((a)-(l)) A-L, where all the terms are normalised with respect to   / 2 th .

Table 1 :
Initial values of the simulation parameters corresponding to the DNS database considered here.0 /  × 131.95 0 /  × 131.95 0 /  is taken for case A, which is discretised by a Cartesian grid of 261 × 128 × 128 with uniform grid spacing in each direction [27]ious laboratory-scale experiments (e.g., Kobayashi et al.[27]).A domain of 105

Table 2 :
The values of normalised turbulent flame speed   /  and normalised flame surface area   /  when statistics were extracted.fortheflameswithLe<1,whichleadsto simultaneous presence of high concentration of reactants and high temperature.As a result, the cases with Le < 1 exhibit greater burning rates than the unity Lewis number flames with the same turbulent flow conditions in the unburned gas.Just the opposite mechanism in terms of reactant and thermal diffusion takes place for Le > 1 flames, which reduces the flame area generation and burning rate in these flames in comparison to the corresponding unity Lewis number flames with the same unburned gas turbulence.This gives rise to increased burning rates and flame area generation with decreasing Le, which can be substantiated from Table2where the normalised turbulent flame speed   /  and normalised flame surface area   /  at the time when the statistics were extracted are presented.The values of   /  have been evaluated by volume integrating the reaction rate ẇ using the expression   = (1/ 0   ) ∫  ẇ  where   is the projected area of the flame in the direction of mean flame propagation, while the values of   /  have been evaluated by volume integrating |∇| (i.e., ∫  |∇|) under both turbulent and laminar conditions.The preheat zone ( < 0.5) becomes increasingly distorted due the increased scale separation between  th and  as a result of increasing Karlovitz number Ka ∼ Re 1/2  /Da ∼  2 th / 2 when Re  increases for a given value of Da.Thus, the preheat zone is more susceptible to be penetrated by energetic turbulent eddies with increasing scale separation between  th and .It is clear from Figures1(h)-1(l) that the flame wrinkling increases with increasing   /  ∼ Re 1/4  Ka 1/2 ∼ Re 1/2  /Da 1/2 , which gives rise to high values of normalised turbulent flame speed   /  and normalised flame surface area   /  for high values of   /  .The variations of (  )  Σ gen ,   ,  UR ,   , and (−  ) with c for all cases are shown in Figures2(a)-2(l), which show that (  )  Σ gen ,   ,  UR , and   assume positive values throughout the flame brush and the contribution of  UR supersedes that of   in all cases, in accordance with earlier findings[2-4,  6-8, 18, 38, 39].As ∇ ⋅ ⃗  remains principally positive within the flame brush and Σ gen is a positive semidefinite quantity (i.e., Σ gen ≥ 0), the dilatation rate term   remains positive have shown that ∇ may align with the most extensive principal strain rate   , (18)ude turbulent Reynolds number dependence but this model parameter needs to be turbulent Reynolds number independent for high values of Re  (i.e., Re  → ∝).Equation(18)ensures that  1 assumes an asymptotic value for Re  → ∝.Moreover,  2 according to (19i) and (19ii) ensures that the effects (15)unts for strengthening of flame normal acceleration with decreasing Le.However, these local Karlovitz number dependences according to (19i) and (19ii) are only two possibilities and any other function of Ka  , which accurately predicts the quantitative behaviour of (− 2 ) and the diminishing strength of  chem on the term (− 2 ) with increasing Karlovitz number Ka  , can in principle be used for  2 .Similarly, any other function which satisfies the expected behaviour of  1 () in response to Re  (Le) can in principle be used to parameterise  1 () instead of(18)(see(21)).The presence of [1 − (  )  (  )  ] in (19i) and (19ii) ensures that the second term on the right hand side of (17) vanishes when the flow becomes fully resolved (i.e., (  )  (  )  = 1.0).The numerical values of  1 ,  2 , , and  are optimised based on a least-squares method.In Figure4, the predictions based on (17) using  2 according to (19i) and (19ii) are compared with (− 2 ) obtained from DNS. Figure4shows that (17) performs satisfactorily for the range of Re  and Le considered here.Using  1 = −  (ũ  /  )Σ gen (see(8)) and combining(15)and (17) yield the following expression:(  )  Σ gen = 2  2 (1 + Ka  )    (Σ gen − |∇|)  th (1 − c) 1 Le  2 + ε k ( 1 −  2 Da  ) Σ gen + ũ    (  |∇| −   Σ gen ) .