Turbulence Effects on the Statistical Behaviour and Modelling of Flame Surface Density and the Terms of Its Transport Equation in Turbulent Premixed Flames

The influence of the ratio of integral length scale to flame thickness on the statistical behaviours of flame surface density (FSD) and its transport has been analysed using a Direct Numerical Simulation database of three-dimensional statistically planar turbulent premixed flames for different turbulence intensities. It has been found that turbulent burning velocity based on volume-integration of reaction rate and flame surface area increase but the peak magnitudes of the FSD and the terms of the FSD transport term decrease with an increase in length scale ratio for a given turbulence intensity. The flame brush thickness and flame wrinkling increase with an increase in length scale ratio for all turbulence intensities. However, the qualitative behaviours of the unclosed terms in the FSD transport equation remain unaltered by the length scale ratio and in all cases the tangential strain rate term and the curvature term act as leading order source and sink, respectively. A decrease in length scale ratio for a given turbulence intensity leads to a decrease in Damköhler number and an increase in Karlovitz number. This has an implication on the alignment of reactive scalar gradient with local strain rate eigenvectors, which in turn increases positive contribution of the tangential strain rate term with a decrease in length scale ratio. Moreover, an increase in Karlovitz number increases the likelihood of negative contribution of the curvature term. Thus, the magnitude of the negative contribution of the FSD curvature term increases with a decrease in length scale ratio for a given turbulence intensity. The model for the tangential strain rate term, which explicitly considers the scalar gradient alignment with local principal strain rate eigenvectors, has been shown to be more successful than the models that do not account for the scalar gradient alignment characteristics. Moreover, the existing model for the curvature and propagation term needed modification to account for greater likelihood of negative values for higher Karlovitz number. However, the models for the unclosed flux of FSD and the mean reaction rate closure are not significantly affected by the length scale ratio.


Introduction
The non-linear temperature and species mass fraction dependences of the chemical reaction rate make the closure of mean reaction rate a challenging task in the modelling of turbulent reacting flows. One of the ways that this challenge can be bypassed in turbulent premixed combustion modelling is the closure of flame surface area to volume ratio (Pope 1988; according to the Flame Surface Density (FSD) based methodology (Pope 1988;Poinsot and Veynante 2001) under flamelet assumption, which assumes that the consumption rate of reactants per unit flame surface area remains identical to that for the laminar flame. However, the FSD is an unclosed quantity and several attempts have been made to model this quantity based on algebraic (Bray 1980;Bray 1988, 1989;Abu-Orf and Cant 2000;Boger et al. 1998;Charlette et al. 2002;Knikker et al. 2002;Fureby 2005;Chakraborty and Klein 2008;Ma et al. 2013;Keppeler et al. 2014;Klein et al. 2016;Ahmed et al. 2021;Rasool et al. 2022) and transport equation (Cant et al. 1990;Duclos et al. 1993;Trouvé and Poinsot 1994;Veynante et al. 1996;Cant 2000, 2001;Chakraborty and Cant 2007, 2009aHun and Huh 2008;Hernandez-Perez et al. 2011;Katragadda et al. 2011Katragadda et al. , 2014aReddy and Abraham 2012;Ma et al. 2014;Sellmann et al. 2017;Papapostolou et al. 2019;Luca et al. 2019;Keil et al. 2020;Varma et al. 2022a;Berger et al. 2022) closures in the context of Reynolds Averaged Navier-Stokes (RANS) (Bray 1980;Bray 1988, 1989;Cant et al. 1990;Duclos et al. 1993;Trouvé and Poinsot 1994;Veynante et al. 1996;Abu-Orf and Cant 2000;Hun and Huh 2008;Cant 2011, 2013;Ahmed et al. 2021;Katragadda et al. 2011Katragadda et al. , 2014bSellmann et al. 2017;Papapostolou et al. 2019;Luca et al. 2019;Rasool et al. 2022;Varma et al. 2022a;Berger et al. 2022) and Large Eddy Simulations (LES) (Boger et al. 1998;Charlette et al. 2002;Knikker et al. 2002;Fureby 2005;Chakraborty and Klein 2008;Cant 2007, 2009a; Hun and Huh 2008;Hernandez-Perez et al. 2011;Reddy and Abraham 2012;Ma et al. 2013Ma et al. , 2014Katragadda et al. 2014a;Keppeler et al. 2014;Klein et al. 2016;Keil et al. 2020). The algebraic closures are designed for the situations where an equilibrium is maintained between the generation and destruction rates of flame surface area (Bray 1980;Bray 1988, 1989;Abu-Orf and Cant 2000;Boger et al. 1998;Charlette et al. 2002;Knikker et al. 2002;Fureby 2005;Chakraborty and Klein 2008;Ma et al. 2013;Keppeler et al. 2014;Klein et al. 2016;Ahmed et al. 2021;Rasool et al. 2022). However, this equilibrium may not be maintained in flame instabilities and under this condition a modelled transport equation for the FSD might be necessary (Cant et al. 1990;Duclos et al. 1993;Trouvé and Poinsot 1994;Veynante et al. 1996;Cant 2000, 2001;Chakraborty and Cant 2007, 2009aHun and Huh 2008;Hernandez-Perez et al. 2011;Katragadda et al. 2011Katragadda et al. , 2014aReddy and Abraham 2012;Ma et al. 2014;Sellmann et al. 2017;Papapostolou et al. 2019;Luca et al. 2019;Keil et al. 2020;Varma et al. 2022a;Berger et al. 2022). The advancements of computational power have enabled significant improvements in the FSD modelling using Direct Numerical Simulations (DNS) data (Trouvé and Poinsot 1994;Boger et al. 1998;Charlette et al. 2002;Chakraborty and Cant 2007, 2009aHun and Huh 2008;Chakraborty and Klein 2008;Katragadda et al. 2011Katragadda et al. , 2014aReddy and Abraham 2012;Klein et al. 2016;Ahmed et al. 2021;Rasool et al. 2022;Papapostolou et al. 2019;Luca et al. 2019;Keil et al. 2020;Varma et al. 2022a;Berger et al. 2022). Among these studies, there have been analyses where the effects of Lewis number (Trouvé and Poinsot 1994;Hun and Huh 2008;Chakraborty and Cant 2011;Katragadda et al. 2011Katragadda et al. , 2014aBerger et al. 2022), turbulent Reynolds number (Chakraborty and Cant 2013;Luca et al. 2019), buoyancy (or Froude number) (Varma 1 3 et al. 2022a) and thermodynamic pressure (Keil et al. 2020) on the statistical behaviour of the FSD and its transport were addressed. A number of recent experimental (Kim et al., 2022) and computational (Yu et al. 2017;Varma et al. 2022a;Song et al. 2021Song et al. , 2022Trivedi et al. 2022) analyses revealed that the flame surface area and turbulent burning velocity for a given turbulence intensity increase with an increase in integral length scale to flame thickness ratio. Moreover, the analysis by Varma et al. (2022a) further demonstrated that the volume integrated values of the strain rate and curvature contributions to the FSD transport increase with increasing integral length scale to flame thickness ratio. Thus, it is expected that the ratio of integral length scale to flame thickness is likely to have an influence on the statistical behaviours of the FSD and the terms of its transport equation. However, to date, the effects of integral length scale to flame thickness ratio on the statistical behaviour of the FSD and the terms of its transport equation for different turbulence intensities are yet to be analysed in detail. This gap in the existing literature has been addressed in this paper and in this respect, the main objectives of this paper are: (a) to demonstrate the effects of the length scale separation on the magnitudes and statistical behaviours of FSD, wrinkling factor and the unclosed terms of the FSD transport equation, and (b) to indicate modelling implications of the different values of turbulence intensity and length scale ratio on the closure of different terms of the FSD transport equation and identify the appropriate model expressions.
The rest of the paper will be organised in the following manner. The mathematical background and numerical implementation relevant to this analysis are presented in Sects. 2 and 3, respectively. This will be followed by the presentation of results and their discussion. The main findings of the present analysis are summarised, and conclusions are drawn in Sect. 5.

Mathematical Background
The reactive scalar field in premixed combustion is often represented by reaction progress variable c , which can be defined based on a suitable species mass fraction Y in the following manner: where subscripts R and P represent values in the fully unburned reactants and fully burned products, respectively. According to Eq. (1) c increases from 0 from the unburned reactants to 1.0 in the fully burned products. The transport equation for c takes the following form: where u j is the jth component of velocity, is density, D is the reaction progress variable diffusivity and ̇=̇∕(Y ,P − Y ,R ) is the reaction rate of the progress variable with ̇ being the reaction rate of species . On Reynolds averaging Eq. (2) one obtains: where q,q = q∕ and q �� = q −q are the Reynolds averaged, Favre averaged and Favre fluctuations of a general quantity q , respectively. All the terms on the right-hand side of Eq. (3) are unclosed and need modelling. The closure of the last term needs modelling of x j turbulent scalar flux ( u �� j c �� ) and the first two terms are modelled in the following manner in the context of FSD based methodology (Trouvé and Poinsot 1994;Boger et al. 1998): where Σ gen = |∇c| is the generalised FSD (Boger et al. 1998;Poinsot and Veynante 2001), S d = (Dc∕Dt)∕|∇c| is the displacement speed and (q) s = q|∇c|∕|∇c| is the surface averaged value of a general variable q (Trouvé and Poinsot 1994;Boger et al. 1998). The transport equation of Σ gen is given by (Pope 1988; where �� ⃗ N = −∇c∕|∇c| is the local flame normal vector. In Eq. (5), T 1 is the turbulent transport term. The term T 2 is commonly referred to as the strain rate term because it originates due to the tangential strain rate a T = ij − N i N j u i ∕ x j , whereas T 3 is referred to as the FSD propagation term and T 4 is known as the FSD curvature term because it arises as a result of flame curvature m = 0.5∇ ⋅ �� ⃗ N = 0.5( N i ∕ x i ) . All the terms on the right-hand side of Eq. (5) (i.e., T 1 , T 2 , T 3 and T 4 ) are unclosed and therefore they need to be modelled. The modelling of T 1 , T 2 , T 3 and T 4 will be discussed in detail in Sect. 4 of this paper.

Numerical Implementation
DNS simulations of statistically planar premixed flames for different turbulence intensities for two different integral length scale to flame thickness ratios have been considered for this analysis. These simulations have been conducted using a well-known compressible DNS code SENGA+ (Jenkins and Cant 1999). In SENGA+ , the conservation equations of mass, momentum, energy, and species are solved in non-dimensional form. In SENGA+ , all first and second-order derivatives for the internal grid points are evaluated using a 10th order central difference scheme but the order of accuracy gradually drops to a one-sided 2nd order scheme at the non-periodic boundaries (Jenkins and Cant 1999). An explicit 3rd order Runge-Kutta scheme (Wray 1990) is used for the timeadvancement of all governing equations. All the simulations consider inlet and outlet boundaries in the direction of mean flame propagation, whereas transverse boundaries are considered to be periodic. The outflow boundary is considered to be partially non-reflecting and specified according to the well-known Navier-Stokes Characteristic Boundary Conditions (NSCBC) technique (Poinsot and Lele 1992). Here the mean inlet velocity U mean was gradually altered to match the turbulent burning velocity in order to ensure that the flame remains stationary in the statistical sense within the computational domain. A well-known pseudo-spectral method (Rogallo 1981) is used to initialise the fluctuating velocity field by a homogeneous isotropic incompressible distribution following the Batchelor-Townsend spectrum (Batchelor and Townsend 1948) for prescribed values of root-mean-square turbulent velocity u ′ and the integral length scale of turbulence l . The reacting scalar field is initialised by a unstretched laminar flame solution.
A modified bandwidth filtered forcing method (Klein et al. 2017) in physical space is used for the current analysis, which yields a forcing term proportional to (1 − c) in the momentum equation such that the forcing is employed for the unburned gas ahead of the flame. This forcing scheme maintains both the required values of root-mean-square velocity u ′ and the integral length scale of turbulence l d in the unburned gas (i.e., c < 0.001 ). Table 1 lists the simulation domain size, the uniform Cartesian grid for discretisation along with the inlet values of root-mean-square turbulent velocity fluctuation normalised by the unstrained laminar burning velocity u � ∕S L , integral length scale to the Zel'dovich flame thickness ratio l∕ z , Damköhler number Da = lS L ∕u � z , Karlovitz number Ka = u � ∕S L 3∕2 l∕ z −1∕2 and heat release parameter = (T ad − T 0 )∕T 0 where the Zel'dovich flame thickness z is defined as: z = T0 ∕S L with T0 and S L being the thermal diffusivity and unstrained laminar burning velocity, respectively. The locations of these cases on Borghi-Peters diagram (Peters 2000) are shown in Fig. 1, which shows that the cases considered here span from the wrinkled flamelet regime to the high Karlovitz number thin reaction zones regime. The grid spacing used here ensures at least 10 grid points within th = (T ad − T 0 )∕max|∇T| L with T, T 0 and T ad being the dimensional temperature, unburned gas temperature and the adiabatic flame temperature, respectively. This implies that the steepest gradients of temperature and reaction progress variable have been resolved using at least 10 grid points. Moreover, the reaction zone for the present thermochemistry (i.e., 0.4 < c < 0.95 ) has at least 10 grid points for the current database. Moreover, at least 1.5 grid points are kept within the Kolmogorov length scale for the simulation parameters considered here. The values of l∕ z = 2.625 and 5.25 correspond to l∕ th = 1.5 and 3.0 for the present thermochemistry. The chemical mechanism is simplified by a single-step Arrhenius mechanism representing stoichiometric methaneair flame in this analysis in the interests of the computational economy for a detailed parametric analysis. The Lewis number of all the species is taken to be unity and the gaseous mixture is taken to obey the ideal gas law. Standard values are considered for Prandtl number (i.e., Pr = 0.7 ), Zel'dovich number (i.e., = T ac T ad − T 0 ∕T 2 ad = 6.0 with T ac being the activation temperature) and the ratio of specific heats (i.e., = 1.4 ). It was demonstrated that the flame propagation and reactive scalar statistics obtained from simple chemistry DNS of stoichiometric methane-air premixed flames remain in good qualitative agreement with the corresponding detailed chemistry simulations (Keil et al. 2021a, b). The quantitative differences remain comparable to the uncertainties associated with the different definitions of reaction progress variable in detailed chemistry simulations (Keil et al. 2021b). All the simulations listed in Table 1 have been conducted until the desired values of turbulent kinetic energy and integral length scale have been obtained and also the turbulent burning velocity S T and flame surface area A T reach statistically stationary states. The readers are directed to Varma et al. (2021) to find out statistically steady nature of volume integrated FSD transport for the database considered in this work. This simulation time remains greater than the through-pass time (i.e., t sim > L x ∕U mean ) and at least 10 eddy turn over times (i.e., t sim > 10l∕u � ) for all cases.
In order to obtain Reynolds/Favre averaged quantities, the quantities of interest are averaged in the transverse directions normal to the mean flame propagation (i.e., homogeneous directions) and also in time over 2 chemical times (i.e., 2t c = 2 z ∕S L ) following previous analyses (Trouvé and Poinsot 1994;Veynante et al. 1997;Katragadda et al. 2011Katragadda et al. , 2014aCant 2011, 2013;Sellmann et al. 2017;Varma et al. 2022a). For statistically planar flames, the Favre-averaged reaction progress variable c is a function of the coordinate in the direction of mean flame propagation (i.e., x-direction). As all the terms in the FSD transport equation and c are functions of x , the results presented in this paper are not dependent on the averaging over bins of c . Thus, all the quantities in Sect. 4 will be presented as a function of c for the sake of generality. These results do not change appreciably (i.e., the maximum variation is of the order of 2%) if the distinct half of the domain and half of the time duration of averaging are considered instead of the full sample size and time series, as used in the results section of this paper.

Flame-Turbulence Interaction and Statistical Behaviour of the FSD
The reaction progress variable c isosurfaces when the flames have reached statistically stationary state are shown in Fig. 2. It can be seen from Fig. 2 that the extent of flame wrinkling increases with increasing u � ∕S L for a given value of l∕ z , as expected. Figure 2 further shows that the extent of flame wrinkling also increases with increasing l∕ z for a given value of u � ∕S L . The observations from Fig. 2 can be substantiated from the variations of the values of normalised flame surface area A T ∕A L and normalised turbulent burning velocity S T ∕S L , which are shown Fig. 3. Here, A T and S T are defined as A T = ∫ V |∇c|dV  Figure 3 shows that A T ∕A L increases with increasing l∕ z for a given value of u � ∕S L , which is also reflected with an increase in S T ∕S L with an increase in  Fig. 4 (first column) that the peak value of Σ gen mostly increases (mildly) with increasing u � ∕S L , whereas the peak value of Σ gen for smaller values of l∕ z is significantly greater than the higher value of l∕ z . In order to understand this behaviour of Σ gen = Ξ × |∇c| , it is worthwhile to consider the variations of |∇c| and Ξ with c , which are also shown in Fig. 4 (second and third columns, respectively). Here, wrinkling factor Ξ is defined as: Ξ = Σ gen ∕|∇c| (Weller et al. 1998;Charlette et al. 2002). It can be seen from Fig. 4 (second column) that the peak value of |∇c| is smaller for higher values of l∕ z . The turbulent flame brush thickness can be quantified as b ∼ 1∕max|∇c| , and thus a smaller peak value of |∇c| for higher values of l∕ z is representative of a thicker flame brush because of greater extent of flame wrinkling (see Fig. 2). The greater extent of flame wrinkling for higher values of l∕ z is reflected in the higher values of wrinkling factor Ξ , as shown in Fig. 4 (third column). However, the decrease in |∇c| with an increase in l∕ z

Variations of the Unclosed Terms in the FSD Transport Equation
The challenges of algebraic closures of Σ gen in the context of RANS were discussed elsewhere Rasool et al. 2022) and qualitatively similar results have been obtained for the present database so will not be discussed further. Instead, the present analysis will focus on the effects of u � ∕S L and l∕ z on the unclosed terms of the FSD Σ gen transport equation. The variations of the normalised unclosed terms of the FSD transport equation (i.e., T 1 , T 2 , T 3 , T 4 × 2 z ∕S L ) with c are shown in Fig. 5 for the cases considered here. It can be seen from Fig. 5 that the qualitative behaviours of the terms of the FSD transport equation are not influenced by l∕ z but the magnitudes of the leading order source and sink contributions by the strain rate term T 2 and the curvature term T 4 , respectively are smaller for higher values of l∕ z . Moreover, the magnitudes of the leading order source and sink contributions of T 2 and T 4 , respectively increase with increasing u � ∕S L for a given value of l∕ z . It can be seen from Fig. 5 that the strain rate term T 2 acts as a leading order source for all cases, where T 4 assumes positive values towards the unburned gas side of the flame brush and becomes negative towards the burned gas side for small turbulence intensities (e.g., u � ∕S L = 1.0 and 3.0) but for high turbulence intensities (e.g., u � ∕S L = 7.5 and 10.0) the curvature term T 4 predominantly assumes negative values throughout the flame brush. The propagation term T 3 and the turbulent transport term T 1 assume locally positive and negative values within the flame brush but the magnitude of their relative contributions in comparison to the magnitudes of T 2 and T 4 diminish with increasing u � ∕S L . These behaviours are qualitatively similar to previous studies for both simple (Trouvé and Poinsot 1994;Cant 2011, 2013;Varma et al. 2022a) and detailed (Luca et al. 2019;Papaposolou et al. 2019;Berger et al. 2022) chemistry and also with the scaling arguments by Chakraborty and Cant (2013). The modelling of these unclosed terms will be addressed next in this paper.

Modelling of the Turbulent Transport Term T 1
Equation 5 indicates that the closure of the turbulent transport term T 1 depends on the modelling of the turbulent flux of the generalised FSD u i s −ũ i Σ gen . The conventional way to close u i s −ũ i Σ gen is to employ a classical gradient hypothesis (Poinsot and Veynante 2001;Cant 2000, 2001): ∕ being the turbulent kinetic energy and its dissipation rate, respectively (Poinsot and Veynante 2001;Cant 2000, 2001). The only nonzero component of the FSD flux is u 1 s −ũ 1 Σ gen for statistically planar flames propagating in the negative x 1 − direction. The variations of u 1 s −ũ 1 Σ gen × z ∕S L and ( Σ gen ∕ x 1 ) × 2 z with c for all cases considered here are shown in Figs. 6 and 7, respectively. A comparison between Figs. 6 and 7 reveals that u 1 s −ũ 1 Σ gen and ( Σ gen ∕ x 1 ) show different signs for a major part of the flame brush and this tendency is particularly prevalent for small values of turbulence intensity. This suggests u 1 s −ũ 1 Σ gen exhib- its counter-gradient transport for small values of u � ∕S L and therefore the gradient hypothesis based model given by Eq. (6) is not suitable for the purpose of modelling u 1 s −ũ 1 Σ gen and therefore its prediction is not shown in Fig. 6.
Previously a model for the turbulent flux of generalised FSD was proposed by Cant (2009b, 2011) in the following manner, which can predict both gradient and counter-gradient behaviour: This model accounts for the fact the FSD flux u 1 s −ũ i Σ gen exhibits a counter-gradient (gradient) behaviour when u ′ ′ i c ′ ′ shows a counter-gradient (gradient) transport. It can be seen from Fig. 6 that the model given by Eq. (7) provides satisfactory prediction of u 1 s −ũ i Σ gen for all cases considered here. The competition between turbulent velocity fluctuation √ 2k∕3 and the velocity jump due to heat release S L determines the nature of turbulent transport (Veynante et al. 1997;Chakraborty and Cant 2009b-d) and thus l∕ z is unlikely have any influence on the closure of u 1 s −ũ i Σ gen . Although Eq. (7) captures the qualitative behaviour of u 1 s −ũ i Σ gen , there is a scope for improvement for quantitative agreement. As the magnitude of T 1 remains small in comparison to the leading order contributions of T 2 and T 4 and thus the modelling inaccuracies of T 1 using Eq. (7) is not expected to affect the closure of the Σ gen transport equation. Moreover, it is worth noting that the successful performance of Eq. (7) depends on the successful closures of u ′′ 1 c ′′ and c ′′ 2 , which have been addressed elsewhere (Veynante et al. 1997;Chakraborty and Cant 2009b-d;Varma et al. 2022b) and will not be addressed further in this paper.

Modelling of the Turbulent Strain Rate Term T 2
The tangential strain rate term T 2 is traditionally modelled by splitting it in the following manner (Cant et al. 1990;Poinsot and Veynante 2001): where S R and S UR are the resolved and unresolved components of the strain rate term, respectively. It can be appreciated from Eq. (8) that the closure of S R depends on the closures of the orientation factors N i N j s , which is modelled by Cant et al (1990) (Bray et al. 1985), which is expected to be negligible for large values of Damköhler number (i.e., Da ≫ 1 ). However, PDFs of c cannot be approximated by a bi-modal distribution in the thin reaction zones regime flames with Da < 1 . Therefore, Katragadda et al. (2011) proposed a modification to the BML expression for c for Da < 1 in the following manner: c = (1 + .g a )c∕ 1 + .g ac where g =c ��2 ∕[c 1 −c ] is the segregation factor, and a = 1.5 is a model parameter. An alternative expression for c was proposed by Bray et al. (2011) in the following manner: c =c + c �� 2 ∕ 0 . It was demonstrated elsewhere Rasool et al. 2022) that c models proposed by Katragadda et al. (2011) and Bray et al. (2011) exhibit comparable performance and satisfactorily predicts c extracted from DNS data. The same conclusion holds also for this study and thus the predictions of c = (1 + .g a )c∕ 1 + .g ac and c =c + c �� 2 ∕ 0 are not explicitly shown here for the sake of conciseness. Another alternative model for N i N j s was proposed by Veynante et al. (1996), which is given as: The models by Veynante et al. (1996) and Cant et al. (1990) will henceforth be referred to as the VPDM and MCPB models, respectively. The predictions of these models are compared to the variation of S R with c across the flame brush, as obtained from DNS data, in Fig. 8. Figure 8 shows that the predictions of the MCPB model mostly captures S R behaviour obtained from DNS data, although it slightly underpredicts the magnitude of S R especially for small values of u � ∕S L (e.g., u � ∕S L = 1.0 and 3.0). The VPDM model, on the other hand, overpredicts S R values for u � ∕S L = 3.0 and 5.0 towards the burned gas side for l∕ z = 2.625 but a satisfactory agreement between DNS data and VPDM model is obtained for l∕ z = 5.25.
The discrepancy between model predictions and DNS data in the modelling of S R originates due to the uncertainty of closure of the orientation factors N i N j s . Cant et al. (1990) assume isotropic distribution of the fluctuating part of the flame normal components while deriving the MCPB model. However, experimental analysis by Veynante et al. (1996) and DNS data by Chakraborty and Cant (2006) revealed that the assumption of the isotropic fluctuation of the flame normal vector does not remain valid especially for small and moderate turbulence intensities. This is responsible for the underprediction of S R by the MCPB model. Although the VPDM model does not assume isotropic nature of fluctuating part of the flame normal components, but the modelling inaccuracies involved in estimating N i N j s leads to the disagreement between the VPDM model and DNS data in Fig. 8. The level of isotropy of flame normal component fluctuations is expected to hold better for smaller values of l∕ z (e.g., l∕ z = 2.625 ) and thus the VPDM model performs relatively better for higher values of l∕ z (e.g., l∕ z = 5.25 ) where the aforementioned isotropy effects are expected to be weak. The unresolved component S UR is modelled by the models proposed by Cant et al. (1990) and , which will henceforth be referred to as the SCPB and SCFM models, respectively. According to Cant et al. (1990) S UR is modelled as S UR = 0.28 √ε ∕ν 0 Σ gen , where ν 0 is the kinematic viscosity of the unburned gas. By contrast, the SCFM model considers S UR = a 0 Γ k ε∕k Σ gen , where a 0 = 2.0 is a model parameter and Γ k is the efficiency function (as proposed by Meneveau and Poinsot (1991)) which is a function of l t ∕δ z and √ 2k∕3∕S L and l t = C kk 3∕2 ∕ε is the local integral length scale where C k = 1.5 (Sreenivasan 1984). The predictions of the SCPB and SCFM models are shown along with the S UR extracted from the DNS data in Fig. 9. It can be seen from Fig. 9 that both SCPB and SCFM models significantly overpredict S UR for all cases. The extent of this overprediction is greater for the SCFM model than the SCPB model for small values of turbulence intensity (e.g., u � ∕S L = 1.0 and 3.0) but the predictions of these models become comparable for large values of turbulence intensity (e.g., u � ∕S L = 10.0 ). Moreover, S UR assumes locally small negative values for u � ∕S L = 1.0 for both l∕ z = 2.625 and 5.25 whereas both SCPB and SCFM models can only predict positive values.
It is worth noting that both SCPB and SCFM models consider only the turbulent time scale to model S UR . The SCPB model considers the Kolmogorov timescale (i.e., √ ν 0 ∕ε ), whereas the SCFM model accounts for large scale turbulent timescale (i.e., k ∕ε ) but the unsatisfactory performance of both the SCPB and SCFM models indicate that turbulent timescale alone is not sufficient to model S UR and T 2 . The term S UR can be expressed as: S UR = (e sin 2 + e sin 2 + e sin 2 )( (where e , e and e are the most extensive, intermediate and the most compressive principal fluctuating strain rate tensor (i.e., e ij = 0.5 ( u �� i ∕ x j + u �� j ∕ x i ) ) and , and are angles between ∇c and eigenvectors associated with e , e and e , respectively. The SCPB and SCFM models implicitly assume the collinear alignment of ∇c with the eigenvector corresponding to e , which is valid only for passive scalar mixing. For high values of Damköhler number premixed flames, ∇c may locally align with the eigenvector corresponding to e (Chakraborty and Swaminathan 2007), which is ignored in these models. Therefore, both SCPB and SCFM models cannot predict the negative values of S UR for small values of u � ∕S L where Damköhler number remains high. In order to avoid the aforementioned limitations of the SCPB and SCFM models, the tangential strain rate term T 2 can alternatively be decomposed as (Katragadda et al. 2011;Sellmann et al. 2017;Varma et al. 2022a): Here, the terms D FSD and N FSD denote the effects of dilatation rate and normal strain rate, respectively. The dilatation rate term can be split as D FSD = D1 FSD + D2 FSD where D1 FSD and D2 FSD are the mean and fluctuating components, respectively, which are defined as (Katragadda et al. 2011;Sellmann et al. 2017;Varma et al. 2022a): The term D2 FSD is modelled as (Katragadda et al. 2011;Sellmann et al. 2017;Varma et al. 2022a): , B 1 and ζ are the model parameters. The model expressions proposed by Varma et al. (2022a) for the components of D FSD and N FSD involving fluctuating quantities previously are summarised in Table 2. In comparison to the model proposed by Varma et al. (2022a), the model parameter B 1 was modified as a function of turbulent Reynolds number Re L (i.e., Re L = 0k 2 ∕ ̃ 0 ) (Katragadda et al. 2011;Sellmann et al. (9) Varma et al. 2022a). The following expressions of B 1 and ζ are proposed here so that the model given by Eq. (11) can be applied for a range of different values of u � ∕S L and l∕ z : It was demonstrated by Chakraborty and Cant (2013) that most model parameters related to the FSD transport reach asymptotic values for Re L ≥ 100 . However, it is important to recognise that Eq. (12) is one of the several possibilities of parameterising B 1 such that it assumes an asymptotic value (= 2.55) for Re L → ∞ and provides the desired values for moderate values of Re L corresponding to the DNS data. It is admitted that the parameterisation of B 1 in Eq. (12) (and some of the other model parameters in the following discussion) involves a degree of empiricism but it is not unusual for turbulence and turbulent combustion modelling. Even if these parameterisations are considered to be purely empirical fits to the DNS data, such fits have a value similar to any new experimental or DNS result that shows an effect, which is yet to be analysed. Moreover, these parameterisations are proposed in such a manner that they assume asymptotic values for Re L → ∞.
The predictions of D FSD according to Eqs. (10) and (11) (with the model parameter given by Eq. (12)) are compared to the corresponding term extracted from the DNS data in Fig. 10. The predictions using the model proposed by Varma et al. (2022a) are also shown in Fig. 10 for the sake of completeness. Figure 10 shows the model given by Eq. (11) (12) B 1 = 1.8 + 0.75erf [ Re L ∕60 − 1] and ζ = −0.3 Table 2 Models of the fluctuating contributions of D FSD and N FSD (i.e., D2 FSD and N2 FSD ) proposed by Varma et al. (2022a) where g * = Γ Z ∕S 2 L with Γ being the body force in the direction of mean flame propagation
The normal strain rate contribution N FSD can also be split into its respective mean and fluctuating components, i.e., N FSD = N1 FSD + N2 FSD , where N1 FSD and N2 FSD are defined as (Katragadda et al. 2011;Sellmann et al. 2017;Varma et al. 2022a): The predictions of N1 FSD with N i N j s evaluated according to the MCPB and VPDM models are compared to the corresponding term extracted from the DNS data in Fig. 11. It can be seen from Fig. 11 that both MCPB and VPDM models reasonably capture the behaviour of N1 FSD extracted from DNS data but their predictions become comparable for high turbulence intensities (e.g., u � ∕S L = 7.5 and 10.0 cases) for both values of l∕ z considered in this analysis. However, the VPDM model shows better agreement with DNS data than the MCPB model for u � ∕S L = 1.0 case for l∕ z = 2.625 , whereas just the opposite behaviour is observed for u � ∕S L = 3.0 case for l∕ z = 2.625 . For the of l∕ z values considered in this analysis, the MCPB model tends to overpredict the magnitude of the negative values of N1 FSD , which is consistent with the underprediction of the magnitude of S R by the MCPB model shown earlier (see Fig. 8). The observations made from Figs. 8 and 11 reveal that the differences in model predictions for S R and N1 FSD become less sensitive for higher values of l∕ z .
In order to consider the modelling of the term N2 FSD it is worthwhile to express N FSD and T 2 as (Katragadda et al. 2011;Sellmann et al. 2017;Varma et al. 2022a): Here e , e and e are the most extensive, intermediate and the most compressive principal fluctuating strain rate tensor (i.e., 0.5( u �� i ∕ x j + u �� j ∕ x i ) ) and , and are the angles between ∇c with the eigenvectors associated with e , e and e , respectively. It has been discussed and demonstrated elsewhere (Chakraborty and Swaminathan 2007;) that ∇c preferentially collinearly aligns with the eigenvectors associated with e (i.e., high probability of |cos | ≈ 1.0 ) under the dominance of the strain rate induced by flame normal acceleration arising from heat release over turbulent straining. By contrast, ∇c demonstrates preferential collinear alignment with the eigenvectors associated with e (i.e., high probability of |cos | ≈ 1.0 ) as a result of the dominance of turbulent straining over the strain rate induced by flame normal acceleration arising from heat release (Chakraborty and Swaminathan 2007;. The high probability of obtaining |cos | ≈ 1.0 is obtained for Da ≫ 1 (Chakraborty and Swaminathan 2007;, which leads to negative values of N2 FSD (and in turn promotes negative values of N FSD and T 2 ). By contrast, high probability of obtaining |cos | ≈ 1.0 is obtained for Da < 1 , which gives rise to positive values of N2 FSD (and in turn promotes positive values of N FSD and T 2 ). For a given value of l∕ z , an increase in u � ∕S L leads to a reduction in Da and therefore, the likelihood of obtaining positive values of N2 FSD increases with increasing turbulence intensity, which can be substantiated from Fig. 12 where the variations of N2 FSD × 2 z ∕S L with c for all cases considered here are shown. It can indeed be seen from Fig. 12 that N2 FSD assumes negative values throughout N2 FSD = −(e cos 2 + e cos 2 + e cos 2 )|∇c| the flame brush for small values of turbulence intensity (e.g., u � ∕S L = 1.0 and 3.0), whereas for large turbulence intensities (e.g., u � ∕S L = 7.5 and 10.0) N2 FSD assumes positive values towards the unburned gas side of the flame brush where the effects of heat release are weak but this term becomes negative from the middle to the burned gas side of the flame brush where the effects of strain rate arising from flame normal acceleration remain strong. It is also worth noting that for a given value of u � ∕S L , an increase in l∕ z leads to an increase in Damköhler number Da and thus the extent of collinear alignment of ∇c with the eigenvector associated with e ( e ) is greater (smaller) in the l∕ z = 5.25 cases than in the l∕ z = 2.625 cases. This is reflected in the higher (smaller) of the propensity negative (positive) contribution of N2 FSD for higher l∕ z case for a given value of u � ∕S L , which can also be confirmed from Fig. 12.
The foregoing discussion suggests that the competition between the strain rate induced by flame normal acceleration and turbulent straining, and the alignment characteristics of ∇c need to be addressed to model N FSD . This is addressed by the following model expression in previous studies (Katragadda et al. 2011;Sellmann et al. 2017;Varma et al. 2022a): where Da L =kS L ∕εδ th is the local Damköhler number. For small values of Damköhler number ∇c preferentially aligns with the eigenvector associated with e γ due to the dominance of turbulent straining over the effects of heat release, which is modelled as ̃ ∕k C 1 Σ gen in Eq. (15). For high Damköhler numbers, heat release effects dominate over the turbulent straining, which leads to the preferential alignment of ∇c with the eigenvector associated with e α . This aspect is accounted for by − C 2 Σ gen S L ∕δ th in Eq. (15). In the present analysis, C 1 and C 2 have been revised in comparison to the previous expressions proposed by Varma et al. (2022a) in the following manner to account for l∕ z dependences which have been embedded in Re L : It can be seen from Fig. 12 that the predictions of the model given by Eq. (15) with model parameters according to Eqs. (16) and (17) remain in satisfactory agreement with N2 FSD extracted from DNS data including its negative value for small values of turbulence intensity (e.g., u � ∕S L = 1.0 and 3.0) and the transition from positive to negative values of N2 FSD for large turbulence intensities (e.g., u � ∕S L = 7.5 and 10.0). It can also be observed from Fig. 12 that the modifications to the closures proposed by Varma et al. (2022a) have resulted in significantly improved predictions of the N2 FSD term.
A combination of the MCPB model for N1 FSD with Eqs. (10), (11) and (15-17) for D1 FSD , D2 FSD and N2 FSD respectively yield the following model for the strain rate term T 2 : . Also shown are the predictions using the closures proposed by Varma et al. (2022a) 1 3 The predictions of Eq. (18) are compared to the predictions of the CPB (MCPB + SCPB) and CFM (VPDM + SCFM) models in Fig. 13 along with the T 2 term extracted from DNS data and the predictions of T 2 according to the closures proposed by Varma et al. (2022a). Figures 13 shows that Eq. (18) agrees better with the DNS data than the CPB and CFM models for u � ∕S L ≥ 5.0 cases for both values of l∕ z but Eq. (18) overpredicts T 2 for l∕ z = 5.25 cases for u � ∕S L = 1.0 and 3.0, whereas an underprediction is obtained for the l∕ z = 2.625 case with u � ∕S L = 3.0 . Despite these discrepancies between the predictions of Eq. (18) and DNS data, the performance of Eq. (18) remains better than the CPB and CFM models. It is worth noting that Eq. (18) takes care about the competition between the strain rates due to turbulence and flame normal acceleration through the model for N2 FSD (see Eq. (15)) in response to the variations of Damköhler number, which is missing in the CPB and CFM models. Thus, Eq. (18) performs better than the SPB and SCFM models for a broad range of values of u � ∕S L and l∕ z . Figure 13 also shows that Eq. (18) gives better predictions of the T 2 term in comparison to the closures proposed by Varma et al. (2022a) especially for higher u � ∕S L values for both length scales. This can be attributed to the use of forced turbulence in the present study compared to decaying turbulence used by Varma et al. (2022a), which results in higher turbulence levels at the time of extraction of statistics in the present study for all initial u � ∕S L cases compared to the analysis by Varma et al. (2022a).

Modelling of the Combined Propagation and Curvature Terms ( T 3 + T 4 )
The curvature term T 4 can be expressed as: T 4 = 2 S r + S n m s Σ gen − 4 D 2 m s Σ gen (Chakraborty and Cant 2011;Katragadda et al. 2014b) where reaction and normal diffusion components of displacement speed are defined as: S r =̇∕ |∇c| and S n = �� ⃗ N • ( D �� ⃗ N • ∇c)∕ |∇c| (Peters et al. 1998;Echekki and Chen 1999). For high values of Karlovitz number (i.e., Ka ≫ 1 ) the negative contribution of {−4 D 2 m s Σ gen } becomes the dominant term, whereas the contribution of 2 S r + S n m s Σ gen remains dominant for small and moderate values of Ka (i.e., Ka < 1 and Ka ∼ 1.0 ). The qualitative behaviour of 2 S r + S n m s Σ gen is determined by the behaviour of ( m ) s because S r + S n remains predominantly positive and Σ gen is a positive semi-definite quantity (i.e., Σ gen ≥ 0 ). It was demonstrated earlier Cant 2011, 2013;Katragadda et al. 2014b) that ( m ) s assumes positive (negative) values towards the unburned (burned) gas side of the flame brush and the same qualitative behaviour has been observed here but not shown here for brevity. The aforementioned behaviour of ( m ) s leads to positive (negative) values of 2 S r + S n m s Σ gen towards the unburned (burned) gas side of the flame brush. A decrease in l∕ z for a given value of u � ∕S L gives rise to an increase in Ka . Thus, the negative contribution of {−4 D 2 m s Σ gen } is more likely to dominate over 2 S r + S n m s Σ gen for smaller values of l∕ z (due to an increase in Karlovitz number with a decrease in l∕ z ) according to the scaling argument by Peters (2000) and it can indeed be seen from Fig. 5 that the (18) Fig. 13 Variation of T 2 × 2 z ∕S L with c for l∕ z = 2.625 (left) and 5.25 (right) along with the predictions of the CPB, CFM models and the model given by Eq. (18) for u � ∕S L = 1.0 , 3.0, 5.0, 7.5 and 10.0 (1st-5th row). Also shown are the predictions using the closures given by Varma et al. (2022a) magnitude of the negative contribution of T 4 is higher for smaller values of l∕ z . Moreover, Karlovitz number increases with increasing u � ∕S L , which also acts to increase the range of curvature values which gives rise to the increase in the magnitude of the negative contribution of {−4 D 2 m s Σ gen } . Accordingly it can be seen from Fig. 14 that ( T 3 + T 4 ) assumes positive values towards the unburned gas and negative values towards the burned gas side of the flame brush for u � ∕S L ≤ 5.0 but the net contribution of these terms assumes predominantly negative values throughout the flame brush for u � ∕S L = 7.5 and 10.0. The curvature and propagation terms are usually modelled together due to their dependence on the displacement speed S d but the behaviour of ( T 3 + T 4 ) is principally determined by the curvature term as the magnitude of the propagation term T 3 remains much smaller than T 4 . The term ( T 3 + T 4 ) is usually modelled together as (Hun and Huh 2008;Cant 2011, 2013;Sellmann et al. 2017;Varma et al. 2022a): The last term on the right-hand side of Eq. (19) denotes the fluctuating component of the combined propagation and curvature term where 0 and c cp are model parameters, and α N represents an orientation factor (resolution factor) in terms of RANS (LES) which is expressed as α N = 1 − N k s N k s and thus becomes zero under fully resolved condition causing the last term to vanish (Fureby 2005;Hawkes and Cant 2000;Jenkins and Cant 1999;Katragadda et al. 2014b). Previous analyses (Hun and Huh 2008;Cant 2011, 2013;Sellmann et al. 2017;Varma et al. 2022a Fig. 14, which shows that this model captures the behaviour of ( T 3 + T 4 ) throughout the flame brush for all cases considered here.

Mean Reaction Rate Closure Using Generalised FSD
The variations of ̇ and ̇+ ∇.( D∇c) with c for all cases considered here are shown in Fig. 15. It can be seen from Fig. 15 that there is a good level of agreement between ̇ and ̇+ ∇.( D∇c) , which shows that the contribution of the mean molecular diffusion rate in statistically planar flames is much smaller than the mean reaction rate and hence 1 3 ≈ S d s Σ gen holds in the context of RANS simulations. For unity Lewis number flames, the surface-averaged density-weighted displacement speed is usually modelled as S d s ≈ 0 S L Cant 2000, 2001) although this approximation may not be valid, especially for the thin reaction zones regime flames (Sabelnikov et al. 2017). Thus, the model for mean reaction rate may be expressed as ̇= 0 S L Σ gen . It can further be seen from Fig. 15 that 0 S L Σ gen satisfactorily captures the variation of ̇ across the flame brush but 0 S L Σ gen underpredicts the magnitude of ̇ towards the product side of the flame brush for all cases.
It is worth noting that the quantity ̇+ ∇.( D∇c) , by definition, equals to S d s Σ gen and thus a comparison between S d s Σ gen and 0 S L Σ gen in Fig. 15 indicates that the magnitude of S d s remains close to 0 S L for the cases considered here. Thus, the discrepancy between ̇+ ∇.( D∇c) (i.e., ̇ due to ̇≫ ∇.( D∇c) ) and 0 S L Σ gen originates due to the approximation of S d s ≈ 0 S L and a better representation of S d s is likely to provide an improved prediction of the mean reaction rate. the discrepancy between 0 S L Σ gen and ̇ originates due to the approximation of S d s ≈ 0 S L and a better representation of S d s is likely to provide an improved prediction of the mean reaction rate.

Conclusions
The effects of both turbulence intensity u � ∕S L and integral length scale to flame thickness ratio l∕ z on the statistical behaviours of FSD and its transport have been evaluated using a DNS database of three-dimensional statistically planar turbulent premixed flames. In contrast to several previous analyses (Trouvé and Poinsot 1994;Cant 2000, 2001;Chakraborty and Cant 2007, 2009aHun and Huh 2008;Katragadda et al. 2011Katragadda et al. , 2014aPapapostolou et al. 2019;Varma et al. 2022a), the current study assesses the performances of the models for different unclosed terms in the FSD transport equation in response to the variations of both u � ∕S L and l∕ z . The main findings are summarised below: • It has been found that the peak magnitudes of the FSD and the terms of the FSD transport term decrease with an increase in length scale ratio. • The turbulent burning velocity and flame surface area assume greater values for the higher values of length scale ratio l∕ z for a given value of u � ∕S L . Both turbulent burning velocity and flame surface area increase with increasing u � ∕S L before the bending behaviour at high values of u � ∕S L . • The flame brush thickness and flame wrinkling increase with an increase in l∕ z for a given value of u � ∕S L . However, the qualitative behaviours of the unclosed terms in the FSD transport equation remain unchanged by the length scale ratio and in all cases the tangential strain rate term T 2 acts as a leading order source and the curvature term T 4 remains the major sink term towards the burned gas of the flame brush. • A decrease in l∕ z for a given value of u � ∕S L leads to a decrease in Damköhler number and an increase in Karlovitz number. This has an implication on the alignment of reactive scalar gradient with local strain rate eigenvectors. Therefore, the existing models for the tangential strain rate term, which do not explicitly account for the alignment of reactive scalar gradient with local strain rate eigenvectors in response to the changes in Damköhler number, have been found to be inadequate. An alternative model for the tangential strain rate term T 2 , which explicitly considers the alignment of ∇c with local strain rate eigendirections and the competition between turbulent straining and chemical timescales, has been shown to be more successful in capturing the behaviour obtained from DNS data. • It has been found that an existing model for the turbulent flux of the FSD, which accounts for both gradient and counter-gradient behaviour, does not need any major modifications due to l∕ z variation. The nature of turbulent transport is determined by relative strengths of turbulent velocity fluctuation and velocity jump across the flame, which is not directly affected by l∕ z . • According to scaling argument by Peters (2000), the curvature term is expected to show greater likelihood of negative values for higher Karlovitz number and accordingly larger magnitudes of the curvature term are obtained for smaller (higher) value of l∕ z ( u � ∕S L ) for a given value of u � ∕S L ( l∕ z ). As a result, the combined propagation and curvature term T 3 + T 4 showed predominantly negative contribution for high Karlovitz number values, which needed modification to an existing model. • The FSD based mean reaction rate closure mostly works satisfactorily irrespective the value of l∕ z but underpredicts the magnitude of ̇ for all cases towards the product side of the flame brush due to inaccuracies involved in modelling surface-averaged densityweighted displacement speed.
It is worth noting that the current analysis has been conducted with simple chemistry for the sake of a detailed parametric analysis. Although previous analyses (Papapostolou et al. 2019;Luca et al. 2019) revealed no major differences in the qualitative behaviours of the unclosed terms of the FSD transport equation between simple and detailed chemistry, future analyses need to include the effects of detailed chemistry for comprehensive validation of the models proposed in this analysis. The turbulent Reynolds number Re t = 0 k 2 ∕ 0 ranges from 12.5 (for case A2) to 250 (case E1) for the DNS database considered here. It is worth noting that most laboratory scale experiments have Re t values of the order of 100. The range of u � ∕S L and l∕ z considered here is comparable to previous studies by other authors (Trouvé and Poinsot 1994;Boger et al. 1998;Charlette et al. 2002;Hun and Huh 2008;Reddy and Abraham 2012) which contributed to the FSD transport modelling. The FSD/wrinkling models which have been either proposed or recommended based on a priori DNS analysis (Boger et al. 1998;Charlette et al. 2002;Fureby 2005;Chakraborty and Cant 2009a) for moderate turbulent Reynolds numbers have been found to work well based on a posteriori analyses (Ma et al. 2013(Ma et al. , 2014Klein et al. 2016). However, further analyses for higher values of Re t will be needed to confirm the current findings. Finally, these models need to be implemented in actual RANS simulations of an experimental configuration for the purpose of a posteriori assessment of the closures of the unclosed terms.