Effects of Body Forces on the Statistics of Flame Surface Density and Its Evolution in Statistically Planar Turbulent Premixed Flames

Body forces such as buoyancy and externally imposed pressure gradients are expected to have a strong influence on turbulent premixed combustion due to the considerable changes in density between the unburned and fully burned gases. The present work utilises Direct Numerical Simulation data of three-dimensional statistically planar turbulent premixed flames to study the influence of body forces on the statistical behaviour of the flame surface density (FSD) and its evolution within the flame brush. The analysis has been carried out for different turbulence intensities and normalised body force values (i.e., Froude numbers). A positive value of the body force signifies an unstable density stratification (i.e., body force is directed from the heavier unburned gas to the lighter burned gas) and vice versa. It is found that for a given set of turbulence parameters, flame wrinkling increases with an increase in body force magnitude in the unstable configuration. Furthermore, higher magnitudes of body force in the unstable density stratification configuration promote a gradient type transport of turbulent scalar and FSD fluxes, and this tendency weakens in the stable density stratification configuration where a counter-gradient type transport is promoted. The statistical behaviours of the different terms in the FSD transport equation and their closures in the context of Reynolds Averaged Navier–Stokes simulations have been analysed in detail. It has been demonstrated that the effects of body force on the FSD and the terms of its transport equation weakens with increasing turbulence intensity as a result of the diminishing relative strength of body force in comparison to the inertial force. The predictions of the existing models have been assessed with respect to the corresponding terms extracted from the explicitly averaged DNS data, and based on this evaluation, suitable modifications have been made to the existing models to incorporate the effects of body force (or Froude number).


Introduction
In many practical applications, turbulent premixed flames propagate within ducts, where the flame is subjected to external pressure gradients and body forces (e.g., gravity). As a considerable density change takes place between unburned and burned gases, buoyancy is expected to play a significant role in premixed turbulent flames. Therefore, external pressure gradients and buoyancy can have a significant impact on the various characteristics such as flame wrinkling and flame normal acceleration in turbulent premixed flames. Chomiak and Nisbet (1995) analysed the pressure gradient-density interactions in turbulent premixed flames and have proposed closures which incorporate these effects into the framework of k − model. It was indeed demonstrated by Veynante and Poinsot (1997) using two-dimensional direct numerical simulation (DNS) data that external pressure gradients and buoyancy can considerably affect the statistical behaviours of turbulent scalar flux and flame wrinkling in turbulent premixed combustion. Thus, it is expected that buoyancy and externally imposed pressure gradients will affect the statistical behaviour of the flame surface density (FSD), which provides the measure of the flame surface area to volume ratio Pope 1988), and is often used for the closure of the mean chemical reaction rate in turbulent premixed combustion (Poinsot and Veynante 2001). However, the effects of body force on the statistical behaviour of the FSD and its evolution have not been analysed in the existing literature, and this paper addresses this problem by using a three-dimensional direct numerical simulation (DNS) database of statistically planar turbulent premixed flames subjected to different strengths of body forces. The fine-grained FSD is defined as Σ � = |∇c|δ(c − c * ) (Pope 1988) where c * is the value of the reaction progress variable c that is used to denote the location of the flame. In order to eliminate the dependence of the FSD statistics on the choice of c * , a generalised FSD is defined as Σ gen = |∇c| (Boger et al. 1998), where the overbar refers to a Reynolds averaging or LES filtering operation, whichever is appropriate. The FSD can either be modelled using an algebraic expression (Boger et al. 1998;Cant and Bray 1988;Charlette et al. 2002;Knikker et al. 2002;Fureby 2005;Chakraborty and Klein 2008a;Ma et al. 2013;Keppeler et al. 2014;Klein et al. 2016) or modelled by solving a transport equation (Cant et al. 1990;Duclos et al. 1993;Veynante et al. 1996;Chakraborty and Cant, 2009aSellmann et al. 2017;Hawkes and Cant 2000;Hernandez-Perez et al. 2011;Reddy and Abraham 2012;Ma et al. 2014) alongside the other conservation equations. The latter approach is preferred in unsteady reacting flows where the production and destruction rates of Σ gen are not necessarily in equilibrium. Although this lack of equilibrium can be addressed by dynamic FSD models in the context of LES (Knikker et al. 2002), such an option does not exist for RANS and thus the transport equation based FSD closure plays an important role in RANS. Modelling of the various terms in the transport equation for FSD have been addressed in detail in several previous studies for RANS (Cant et al. 1990;Duclos et al. 1993;Veynante et al. 1996;Cant 2011, 2013;Sellmann et al. 2017) and LES 1 3 (Hawkes and Cant 2000;Cant 2007, 2009a;Hernandez-Perez et al. 2011;Reddy and Abraham 2012;Ma et al. 2014;Katragadda et al. 2011Katragadda et al. , 2014) but all of these analyses have been conducted in the absence of external pressure gradients and body forces. In this study a three-dimensional DNS database involving a number of statistically planar turbulent premixed flames subjected to different values of unburned gas turbulence intensity u � ∕S L (where u ′ is root-mean-square (rms) value of turbulent velocity fluctuation) and Froude number (i.e., a measure of inertial force to body force ratio) have been conducted. In this respect, the main objectives of the present study are as follows: 1. To investigate the effects of body force on flame-turbulence interaction and the statistical behaviour of the FSD transport for turbulent premixed flames. 2. To assess the performances of the existing closure models for the different terms of the FSD transport equation and propose suitable modifications, wherever necessary, to capture the effects of the body force.
The remainder of this paper is organised as follows. The next section describes the essential mathematical background which is followed by the discussion of the numerical implementation. The results are subsequently presented and discussed. The final section summarises the main findings and the conclusions that are drawn from the present study.

Mathematical Background
In the present analysis, DNS of turbulent premixed combustion has been carried out in three dimensions with simple chemistry in the interest of a detailed parametric analysis because detailed chemistry simulations are unsuitable due to its expensive computational requirements. Hence, three-dimensional DNS with a single step Arrhenius type irreversible chemical reaction has been used in the present study. The reaction progress variable c is used to represent the species field and it is defined in terms of a suitable product mass fraction Y P according to c = Y P − Y P0 ∕ Y P∞ − Y P0 , where 0 and ∞ refer to the unburned and fully burned gases, respectively. Under adiabatic conditions, in the case of flames with unity Lewis number, the reaction progress variable c is identical to the non-dimensional temperature T = T − T 0 ∕ T ad − T 0 (Chakraborty and Cant 2011), where T is the instantaneous dimensional temperature, T 0 is the temperature of the unburned gas and T ad is the adiabatic flame temperature. Under the action of the body force, the momentum conservation equation in the ith direction takes the following form: where u i is the ith component of fluid velocity, p is the pressure, τ ij is the component of shear stress tensor and S i = ρΓ i is the source/sink term in the ith direction due to body force. The source term is assumed to take non-zero values only in the x 1 -direction, which is aligned with the mean direction of flame propagation. Thus, S 1 is expressed as S 1 = ρΓ = g * S 2 L ∕ Z where δ Z = α T0 ∕S L is the Zel'dovich flame thickness, S L is the unstrained laminar burning velocity, α T0 is the thermal diffusivity in the unburned gas, and g * stands for the inverse of Froude number-squared (i.e., g * = Γ Z ∕S 2 L = Fr −2 ). A positive (negative) value of g * indicates an externally imposed flow acceleration (deceleration). (1) Moreover, a positive (negative) value of g * is representative of an externally imposed adverse (favourable) pressure gradient. A positive value of g * also represents an unstable configuration where the heavier cold unburned reactants sit on top of the lighter hot burned products whereas a negative g * value is indicative of a stable configuration where the heavier cold unburned reactants sit underneath the lighter hot burned products. Considering an atmospheric stoichiometric CH 4 /air flame, a value of |g * | = 3.12 corresponds to a pressure gradient magnitude of 2000 Pa/m , which is consistent with the experimental parameters used by Shepherd et al. (1982). It is important to note that the buoyancy effects are usually negligible in turbulent premixed flames, except for large fires/explosions, but the methodology adopted in this paper and by Veynante and Poinsot (1997) provides a convenient numerical tool to analyse the effects of imposed pressure gradients. The g * values investigated in this study are chosen such that they induce pressure gradients similar to the experimental analysis by Shepherd et al. (1982), and these values are consistent with the values previously considered by Veynante and Poinsot (1997). The transport equation of the reaction progress variable c is given by: where ẇ and D are the reaction rate and mass diffusivity of c , respectively. Here q = q∕ and q �� = q −q are the Favre-average and Favre fluctuations of a general variable q , respectively. The quantity u ′′ j c ′′ is the turbulent scalar flux and needs closure, and interested readers are referred to Veynante and Poinsot (1997) for further discussion regarding the effects of pressure gradient on the modelling of turbulent scalar flux u ′′ j c ′′ . The first and second terms on the right-hand side can be modelled in the following manner (Boger et al. 1998;Cant 2007, 2011;Hawkes and Cant 2000): where S d = (1∕|∇c|)(Dc∕Dt) is the displacement speed, which signifies the speed with which the flame moves normal to itself with respect to the initially coincident material surface. In the context of RANS, the transport equation for the generalised FSD is given as (Cant et al. 1990;Duclos et al. 1993;Veynante et al. 1996;Cant 2011, 2013;Sellmann et al. 2017;Katragadda et al. 2011): where ⃗ N = −∇c∕|∇c| is the local flame normal vector, (Q) s is the surface average of a general quantity Q given by (Q) s = Q|∇c|∕|∇c| . The two terms on the left-hand side of Eq. 4 are the transient and mean advection terms, respectively, whereas the four terms on the right-hand side are referred to as the turbulent transport term T 1 , tangential strain rate term T 2 , propagation term T 3 and curvature term T 4 , respectively. The terms T 1 , T 2 , T 3 and T 4 are unclosed and need modelling. The performances of the existing closures of T 1 , T 2 , T 3 and T 4 will be assessed with respect to the corresponding terms extracted from DNS data in Sect. 4.

Numerical Implementation
A well-known compressible DNS code SENGA+ (Jenkins and Cant 1999) is used to generate the database in the present study. In SENGA+, all the spatial derivates for the internal grid points are evaluated using a 10th order central difference scheme and the order of accuracy gradually drops to a one-sided 2nd order scheme at the non-periodic boundaries.
The time advancement has been carried out using a low-storage 3rd order Runge-Kutta scheme (Wray 1990). The Lewis number of all species are considered to be unity and standard values are considered for the Prandtl number Pr ( Pr = 0.7) and Zel'dovich number β = T ac T ad − T 0 ∕T 2 ad (β = 6.0) , where T ac is the activation temperature. The heat release parameter τ = T ad − T 0 ∕T 0 is considered to be 4.5 . These values of and are representative of the stoichiometric methane-air flames preheated to T 0 = 415 K. The simulations have been carried out on a rectangular domain of size 70.2δ Z × 35.1δ Z × 35.1δ Z , which is discretised by a grid size of 400 × 200 × 200, ensuring that the grid spacing remains smaller than the Kolmogorov length scale and 10 grid points are accommodated within the thermal flame thickness The boundaries in the direction of mean flame propagation (i.e., x 1 -direction) are taken to be partially non-reflecting and the boundaries in the transverse directions are taken to be periodic. The boundary conditions are specified using the Navier-Stokes Characteristic Boundary Conditions (NSCBC) technique (Poinsot and Lele 1992). A divergence free, homogeneous, isotropic turbulence field generated using the pseudo-spectral method of Rogallo (1981) according to a Batchelor-Townsend kinetic energy spectrum (Batchelor and Townsend 1948) is used to initialise the turbulent velocity field. The flame and reacting scalar fields have been initialised by steady-state one-dimensional unstretched laminar premixed flame solution. The initial values of normalised root-mean-square turbulent velocity fluctuation u � ∕S L , integral length scale to the flame thickness ratio Table 1. For the values of u � ∕S L and l T ∕ th ( l T ∕ Z ) listed in Table 1, all flames considered in this analysis belong to the thin reaction zones regime of combustion (Peters 2000). The simulations have been conducted for g * = −3.12, −1.56, 0.0, 1.56, 3.12 for each set of turbulence parameters summarised in Table 1. Simulations have been carried out for 3.0−8.0 initial eddy turnover times for Set-A to Set-D, respectively, which is greater than 2δ th ∕S L for all the cases considered. The turbulent kinetic energy and dissipation rate evaluated over the unburned gas volume did not change significantly with time when the statistics are extracted. During the overall run time of the simulations, u ′ decayed by about 55% and l T increased by about 22% in comparison to their initial values. This simulation time remains comparable to several studies (Boger et al. 1998;Charlette et al. 2002;Hernandez-Perez et al. 2011;Hun and Huh 2008) which contributed to the FSD based modelling in the past. Moreover, the qualitative nature of the FSD transport statistics presented in this paper did not change since halfway through the simulation. The Reynolds/Favre-averaged values of a general quantity q are evaluated by ensemble averaging the relevant quantities in x 2 − x 3 planes at a given x 1 location and the Favre fluctuations are directly evaluated from the DNS data by subtracting the Favre averaged quantities from the instantaneous values following several previous analyses Cant 2011, 2013;Sellmann et al. 2017;Katragadda et al. 2011;Hun and Huh 2008;).

Flame-Turbulence Interaction and Statistical Behaviour of FSD
The instantaneous views of the iso-surfaces of reaction progress variable c for g * = −3.12, 0.0 and 3.12 for all sets of turbulence parameters considered are shown in Fig. 1 when the statistics were extracted. The corresponding contour plots of the reaction progress variable in the x 1 − x 2 mid plane of the simulation domain are shown in Fig. 2. It can be seen from Fig. 1 that the value of g * influences the flame morphology and the largescale flame wrinkling in the direction of mean flame propagation (i.e., x 1 -direction) is promoted with increasing g * , possibly due to Rayleigh-Taylor type instability. Figures 1 and 2 indicate that the flame wrinkling increases with increasing g * and u � ∕S L , which is also reflected in the values of normalised flame surface area A T ∕A L and normalised turbulent flame speed S T ∕S L shown in Fig. 3. Here, the turbulent flame speed and flame surface area have been evaluated using the volume integrals S T = ∫ Vẇ dV∕ρ 0 A P and A T = ∫ V |∇c|dV , respectively, where A P is the projected flame surface area in the direction of mean flame propagation. The subscripts T and L refer to the turbulent and laminar flame conditions, respectively. The increased flame wrinkling with increasing values of g * is consistent with the previous two-dimensional findings by Veynante and Poinsot (1997) and the theory proposed by Libby (1989).
The variation of the generalized FSD Σ gen = |∇c| , the resolved FSD | | ∇c | | and the wrinkling factor Ξ = Σ gen ∕ | | ∇c | | with the Favre averaged reaction progress variable c are shown in Fig. 4. It can be seen from Fig. 4 that smaller values of Σ gen × δ Z are obtained for higher g * and this behaviour is most prominent for small initial values of u � ∕S L (e.g., for Set-A and Set-B). However, the magnitudes of Σ gen × δ Z do not change with the variation of g * for large values of u � ∕S L (e.g., for Set-C and Set-D). The corresponding variation of the resolved part of the generalised FSD (i.e., | | ∇c | | ) shows that | | ∇c | | × δ Z decreases considerably with increasing g * . Since the inverse of the maximum value of | | ∇c | | is a measure of the turbulent flame brush thickness (i.e., δ T ∼ 1∕max | | ∇c | | ), a decrease in | | ∇c | | indicates a broadening of the flame brush with increasing g * , which is consistent with the contours of c shown in Fig. 2 and previous results from two-dimensional DNS data (Veynante and Poinsot 1997). It can also be observed from Fig. 4 that the wrinkling factor Ξ increases with increasing g * but this increase is insufficient to overcome the effects of increased flame brush thickening, thereby resulting in a reduction in Σ gen magnitude with increasing g * . It is also evident from Fig. 2 that the iso-surfaces of c are more wrinkled on the unburned gas side representing the preheat zone of the flame, which is valid for the thin reaction zone regime of combustion (Peters 2000).

Statistical Behaviour of FSD Transport
The variations of normalised unclosed terms of the FSD transport equation { T 1 , T 2 , T 3 and T 4 } × 2 Z ∕S L for g * = −3.12 , 0.0 and 3.12 are shown in Fig. 5 for Set-A to Set-D. The variations of these terms for the remaining g * values (i.e., g * = −1.56 and 1.56 ) follow the monotonic trend that can be derived from g * = −3.12 , 0.0 and 3.12 and thus are not shown here for the sake of conciseness. The following observations can be made from • The magnitude of the turbulent transport term T 1 decreases as the value of g * increases for a given value of u � ∕S L . The magnitude of T 1 also decreases with increasing initial u � ∕S L . • The FSD tangential strain rate term T 2 acts as a source term in all cases and its magnitude increases with increasing g * for small values of u � ∕S L (e.g., Set-A). For high values of turbulence intensity u � ∕S L , the magnitude of T 2 on the reactant side of the flame brush increases with increasing g * but the magnitude of T 2 on the product side remains mostly insensitive to the variation of g * . • The FSD propagation term generally assumes positive values towards the unburned gas side and negative values towards the burned gas side of the flame brush. The magnitude of the propagation term T 3 decreases with increasing g * for a given value of u � ∕S L . The magnitude of T 3 also decreases with increasing initial turbulence intensity u � ∕S L . The case with g * = −3.12 for Set-A remains weakly turbulent with a small extent of flame wrinkling (see Figs. 1, 2). For a steady laminar 1D planar premixed flame, the contributions of T 1 , T 2 and T 4 remain identically zero, whereas T 3 assumes positive value towards the unburned gas side of the flame brush before assuming negative values towards the burned gas side of the flame brush. The contribution of −C 1 = − ũ k Σ gen ∕ x k is balanced by T 3 in a steady laminar 1D planar premixed flame. A similar situation holds in the case g * = −3.12 for Set-A (not shown here). The departure from the expected behaviour for a steady 1D laminar premixed flame increases with increasing initial turbulence intensity u � ∕S L due to the strengthening of turbulence.

Fig. 5
Variations of T 1 , T 2 , T 3 and T 4 with c for g * = −3.12 (first column), 0.0 (second column) and 3.12 (third column) for Set-A, Set-B, Set-C and Set-D (rows 1-4), respectively 1 3 • The FSD curvature term T 4 is weakly positive towards the reactant side but acts as a sink (negative) over the major part of the flame brush. The positive contribution of T 4 weakens with increasing initial turbulence intensity u � ∕S L . • The magnitudes of FSD transport and propagation terms (i.e., T 1 and T 3 ) relative to the FSD tangential strain rate and curvature terms (i.e., T 2 and T 4 ) decrease with increasing u � ∕S L . This behaviour is consistent with previous findings by Chakraborty and Cant (2013) and their scaling arguments.
It is worthwhile to consider the ratio of body forces to the inertial forces due to turbulent velocity fluctuations as: where Fr T = u � ∕ √ Γ 1 l T is the turbulent Froude number. This suggests that an increase in u � ∕S L for a given set of values of l T ∕ Z and g * leads to weakening of the body force effects and thus the FSD and the terms of its transport equation do not change appreciably with the variation of g * for large values of u � ∕S L . The observed behaviours of the terms T 1 , T 2 , T 3 and T 4 from Fig. 4 can be explained further by considering these terms individually and their closures.

Modelling of the Turbulent Transport Term T 1
It can be seen from Eq. 4 that the modelling of the turbulent transport term T 1 depends on the closure of the turbulent flux of the generalised FSD given by u i s −ũ i Σ gen (Cant et al. 1990;Duclos et al. 1993;Veynante et al. 1996;Chakraborty and Cant 2009aSellmann et al. 2017;Hawkes and Cant 2000;Hernandez-Perez et al. 2011;Reddy and Abraham 2012;Ma et al. 2014). This closure is often achieved using a classical gradient hypothesis (Cant et al. 1990;Hawkes and Cant 2000;Reddy and Abraham 2012;Ma et al. 2014) according to which one obtains: where ν t = C μk 2 ∕̃ is the kinematic eddy viscosity (with C μ = 0.09 ) and Sc Σ is a suitable turbulent Schmidt number (Cant et al. 1990;Duclos et al. 1993) ∕ being the turbulent kinetic energy and its dissipation rate, respectively and evaluated directly from the DNS data. However, previous studies Cant 2011, 2013;Sellmann et al. 2017;) demonstrated that u i s −ũ i Σ gen might exhibit counter-gradient (gradient) behaviour when the turbulent flux of reaction progress variable u ′′ 1 c ′′ is counter-gradient (gradient) in nature. The variations of the only non-zero normalised turbulent scalar flux component u �� 1 c �� ∕ 0 S L for statistically planar flames with the Favre averaged reaction progress variable c are shown in Fig. 6 for all cases considered here. As c∕ x 1 is positive in the present configuration, a positive (negative) value of u �� 1 c �� ∕ 0 S L indicates countergradient (gradient) type transport. It is evident from Fig. 6 that counter-gradient transport is observed over the majority of the flame brush for all cases. However, a gradient transport is promoted for positive g * values (i.e., g * = 1.56 and 3.12 ), towards the reactant side of the flame brush where the effects of heat release are weak and this tendency weakens as the value of g * decreases, which is consistent with the previous findings (Chomiak and Nisbet 1995;Veynante and Poinsot 1997). This behaviour can (5) u i s −ũ i Σ gen = − ν t ∕Sc Σ Σ gen ∕ x i 1 3 be explained using a hydrostatic approximation (Veynante and Poinsot 1997) under which the momentum equation reduces to p∕ x 1 = g * S 2 L ∕ Z . Thus, a positive (negative) value of g * denotes an adverse (favourable) pressure gradient across the simulation domain, leading to gradient (counter-gradient) transport (Veynante and Poinsot 1997).
The variation of the only non-zero component of the FSD flux u 1 s −ũ 1 Σ gen for statistically planar flames with c are shown in Fig. 7 for the cases considered in the present analysis. It can be seen from Fig. 4 that Σ gen attains peak value at the middle of the flame brush (i.e., at c = c Σ with 0 < c Σ < 1.0 ), and therefore Σ gen ∕ x 1 assumes positive values for 0 <c < c Σ and negative values for c Σ <c < 1.0 because c∕ x 1 assumes positive values within the flame brush, in statistically planar flames. This indicates that Σ gen ∕ x 1 is expected to be positive (negative) towards the reactants (products) side of the flame brush. Therefore, u 1 s −ũ 1 Σ gen is expected to assume negative (positive) values towards the reactants (products) side of the flame brush according to Eq. 5. However, the behaviour of u 1 s −ũ 1 Σ gen obtained from DNS in Fig. 7 reveals that the FSD flux assumes positive (negative) values towards the reactants (products) side of the flame brush for small values of u � ∕S L and also for negative values of g * . This suggests that counter-gradient behaviour is observed for u 1 s −ũ 1 Σ gen for small values of u � ∕S L and also for negative values of g * . Thus, the gradient hypothesis-based model given by Eq. 5 is not suitable for the purpose of modelling u 1 s −ũ 1 Σ gen and therefore its prediction is not shown in Fig. 7. It can indeed be dis-  Fig. 7 that the qualitative behaviour of u 1 s −ũ 1 Σ gen in Set-C and Set-D are different from those in Set-A and Set-B for a given value of g * . Moreover, the qualitative behaviours of u 1 s −ũ 1 Σ gen also change with the variation of g * for a given turbulence intensity u � ∕S L . This difference in behaviour originates from the fact that u 1 s −ũ 1 Σ gen shows counter-gradient behaviour for small values of u � ∕S L (e.g., Set-A and Set-B), whereas the gradient behaviour is prevalent towards the reactants side of the flame brush for large values of u � ∕S L (e.g., Set-C and Set-D) similar to the behaviour of turbulent scalar flux u ′′ 1 c ′′ (see Fig. 6). As shown in Fig. 6, the propensity of counter-gradient behaviour of turbulent scalar flux u ′′ 1 c ′′ increases with decreasing g * and a similar behaviour is observed for u 1 s −ũ 1 Σ gen with counter-gradient behaviour for small values of g * and a gradient type behaviour on the reactant side of the flame brush is observed for high positive values of g * . Fig. 7 Variation of u i s −ũ i Σ gen × Z ∕S L (solid line represents DNS data and dash-dotted line represents prediction of Eq. 6) with c for Set-A to Set-D (1st -4th column) for g * = −3.12, −1.56, 0.0, 1.56 and 3.12 (1st-5th row) A model for the turbulent flux of generalised FSD was proposed in previous analyses Cant 2011, 2013;Sellmann et al. 2017) in the following manner, which can predict both gradient and counter-gradient behaviour: It can be seen from Fig. 7 that Eq. 6 satisfactorily captures the qualitative and quantitative behaviour of u 1 s −ũ 1 Σ gen irrespective of the values of u � ∕S L and g * . However, it is worth noting that the successful performance of Eq. 6 depends on the satisfactory closure of u ′′ 1 c ′′ and c ′′2 . The closures of u ′′ 1 c ′′ and c ′′2 have been addressed elsewhere Chakraborty and Cant 2009b, c, d;Papapostolopu et al. 2019; Chakraborty and Swaminathan 2011) and will not be addressed further in this paper.

Modelling of the tangential strain rate term T 2
For the purpose of understanding the statistical behaviour of the tangential strain rate term T 2 , it has been decomposed in several previous studies (Cant et al. 1990;Duclos et al. 1993;Veynante et al. 1996;Cant 2011, 2013;Sellmann et al. 2017;Katragadda et al. 2011) in the following manner: Here S R and S UR are the resolved and unresolved flow contributions of the strain rate term, respectively. The modelling of S R depends on the appropriate modelling of the orientation factors N i N j s . Cant et al. (1990) modelled N i N j s as N i N j s = N i s N j s + δ ij ∕3 1 − N k s N k s where N i s = − c∕ x i ∕Σ gen is the surface averaged flame normal vector. The evaluation of N i s depends on the evaluation of c from c , which is usually achieved using the Bray-Moss-Libby (BML) analysis (Bray et al. 1985), according to which c = (1 + τ)c∕(1 + τc) + O(γ) , where O( ) represents the contribution of the burning mixture. The BML analysis implicitly assumes bi-modal probability density function (PDF) of c with impulses at c = 0 and c = 1.0 (Bray et al. 1985), which is likely to be obtained for large values of Damköhler number (i.e., Da ≫ 1 ) but the bi-modal PDF distribution is unlikely to be realised in the thin reaction zones regime flames with Da < 1 . Katragadda et al. (2011) proposed a modification to the BML expression for c in the case of flames with Da < 1 , where the bi-modal PDF of c is unlikely to be realised and this expression takes the following form for unity Lewis number flames: where g = � c ��2 ∕[c(1 −c)] is the segregation factor, and a = 1.5 is a model parameter (Katragadda et al. 2011). Equation 8 has been found to capture the interrelation between c and c for all cases considered in this work, which is not explicitly shown here for the sake of conciseness.
1 3 Veynante et al. (1996) proposed an alternative model for N i N j s , which is given as: The models by Cant et al. (1990) and Veynante et al. (1996) will henceforth be referred to as the MCPB and VPDM models, respectively. The variation of S R with c across the flame brush from DNS data is compared with predictions of the MCPB and VPDM models in Fig. 8. It can be seen from Fig. 8 that the predictions of the MCPB model agree better with the DNS data for most of the cases, although it slightly underpredicts the magnitude of S R for all cases considered. The VPDM model, on the other hand, highly overpredicts S R for cases with small initial u � ∕S L values (Set-A and Set-B) and negative g * values. However, the agreement of the prediction of the VPDM model with DNS data improves with increases in g * and the initial u � ∕S L (e.g., the VPDM model performs better in Set-C and Set-D than in Set-A and Set-B).
The MCPB model assumes isotropic fluctuations of surface averaged flame normal components (Cant et al. 1990), whereas the VPDM model considers anisotropic fluctuations of  1.56, 0.0, 1.56 and 3.12 (1st-5th row) surface averaged flame normal components (Veynante et al. 1996). Moreover, there is another major difference between these models. The quantity N i N j s approaches N i s N j s in the limit of a laminar premixed flame (where N k s N k s ≈ 1.0 ), which is satisfied by the MCPB model but not by the VPDM model. The case with g * = −3.12 for Set-A exhibits the least turbulent behaviour and shows weak flame wrinkling (see Figs. 1, 2) and thus in this case the MCPB model shows a good agreement with S R obtained from DNS data but the VPDM model does not adequately capture the statistical behaviour of the orientation factors N i N j s . The departure from the laminar premixed flame condition increases with increasing g * and/or u � ∕S L and thus the agreement between the VPDM model and DNS data for S R improves with increases in g * and/or u � ∕S L .
The two most widely used models for the unresolved component S UR are the ones proposed by Cant et al. (1990) and , which are referred to as the SCPB and SCFM models, respectively. Cant et al. (1990) modelled S UR as S UR = 0.28 √̃∕ ν 0 Σ gen , where ν 0 is the kinematic viscosity of the unburned gas. According to the SCFM model, S UR = a 0 Γ k ̃∕k Σ gen , where a 0 = 2.0 is a model constant 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 the constant C k is taken to be C k = 1.5 (Sreenivasan 1984). The predictions of the SCPB and SCFM models are compared with the S UR extracted from the DNS data in Fig. 9, which shows neither the SCPB nor SCFM model predicts the variation of S UR satisfactorily. Both models predict only positive values of S UR , although negative values of S UR are obtained from DNS data for small values of u � ∕S L and g * . Moreover, both SCPB and SCFM models tend to over-predict the magnitude of the unresolved strain rate term S UR by a large extent for most cases considered here. However, for the initial u � ∕S L = 3.0 case (Set-A), the SCFM model provides satisfactory prediction for positive g * values.
In order to explain the negative values of S UR , the tangential strain rate term T 2 can alternatively be decomposed as Cant 2011, 2013;Sellmann et al. 2017;Katragadda et al. 2011): Here the terms D FSD and N FSD denote the effects of dilatation rate and normal strain rate, respectively. Katragadda et al. (2011) andSellmann et al. (2017) further split the dilatation rate term as D FSD = D1 FSD + D2 FSD where D1 FSD and D2 FSD are the resolved and unresolved components, respectively, and these are defined as: Katragadda et al. (2011) andSellmann et al. (2017) proposed the following model expression for D2 FSD ∶ (9) , B 1 and ζ = −0.3 are the model parameters. The model parameter B 1 has been modified in this analysis in comparison to the previous studies (Sellmann et al. 2017;Katragadda et al. 2011) to account for the effects of the local turbulent Reynolds number Re L (i.e., Re L = 0 � k 2 ∕ ̃0 ) and the Froude number (i.e., Fr = 1∕ √ g * ) in the following manner: where b 1 and b 2 are functions of g * as shown below: The term D FSD extracted from the DNS data and the combined predictions of D1 FSD and D2 FSD according to the closures given by Eqs. 10 and 11-13 are compared in Fig. 10. It can be seen from Fig. 10 that the magnitude of D FSD × δ 2 Z ∕S L decreases with increasing g * and u � ∕S L , which necessitates the inclusion of g * and Re L dependences in Eq. 13. It can also be seen that the models given by Eqs. 10, 11-13 satisfactorily predict the behaviour of the D FSD term for all the different g * values and for all sets of initial u � ∕S L values considered in this analysis. However, some over-predictions are obtained towards the reactant side of the flame brush for positive g * values for Set-A.
Similar to D FSD , the normal strain rate contribution N FSD can also be split into its respective resolved and unresolved flow contributions (Sellmann et al. 2017;Katragadda et al. 2011), i.e., N FSD = N1 FSD + N2 FSD , where N1 FSD and N2 FSD are defined as: Fig. 10 Variation of D FSD × δ 2 Z ∕S L (solid line represents DNS data and dash-dotted line represents the prediction of the closure model) with c for Set-A to Set-D (1st-4th column) for g * = −3.12, −1.56, 0.0, 1.56 and 3.12 (1st-5th row)

3
The predictions of N1 FSD according to the CPB and VPDM models are compared to the corresponding term extracted from the DNS data in Fig. 11. The magnitude of N1 FSD × 2 Z ∕S L is found to decrease with increasing values of g * and initial u � ∕S L . It is evident from Fig. 11 that overall, for all sets of initial u � ∕S L and g * values, the CPB model captures the quantitative and qualitative behaviour of the N1 FSD component more accurately than the VPDM model. However, the CPB model tends to overpredict the magnitude of the N1 FSD component whereas the VPDM model underpredicts (overpredicts) it for negative (positive) g * values. It can also be seen that the agreement between the predictions of the VPDM model and DNS data improve significantly as the values of initial u � ∕S L and g * increases. This behaviour is consistent with the S R predictions by the MCPB and VPDM Fig. 11 Variation of N1 FSD × δ 2 Z ∕S L (solid line represents DNS data, dotted line represents the predictions of the CPB model and dashed line represents the predictions of the VPDM model) with c for Set-A to Set-D (1st-4th column) for g * = −3.12, −1.56, 0.0, 1.56 and 3.12 (1st-5th row) 1 3 models shown in Fig. 8 and the differences in these two model predictions have already been explained in the context of S R predictions and the same explanation is valid in the context of the N1 FSD .
The behaviour of the unresolved component of the normal strain rate term N2 FSD depends on the alignment of ∇c with the local principal fluctuating strain rates. The terms N FSD and S UR can be expressed as (Katragadda et al. 2011): where 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 and the eigenvectors associated with e , e and e , respectively. The reaction progress variable gradient ∇c preferentially collinearly aligns with the eigenvectors associated with e (i.e., high probability of | cos | ≈ 1.0 ) when the strain rate induced by flame normal acceleration arising from heat release overcomes turbulent straining (Chakraborty and Swaminathan 2007;. By contrast, ∇c preferentially collinearly aligns with the eigenvectors associated with e (i.e., high probability of when turbulent straining dominates 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 high values of Da (Chakraborty and Swaminathan 2007;, which leads to negative values of T 2 and N FSD (and its unresolved component N2 FSD ). This behaviour is promoted further for small values of g * because turbulence straining is expected to weaken with decreasing g * especially for negative values of g * . By the same token, the high probability of obtaining | cos | ≈ 1.0 is obtained for low values of Da (Chakraborty and Swaminathan 2007;, which leads to positive values of S UR and N FSD (and its unresolved component N2 FSD ). This behaviour strengthens further for high positive values of g * because turbulence straining strengthens with increasing g * especially for positive values of g * . The aforementioned 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 in the modelling of T 2 . The SCPB and SCFM models consider the Kolmogorov time scale and large-scale turbulent time scale, respectively, which are both turbulent flow time scales, and thus ignores the chemical time scale associated with the strain rate induced by flame normal acceleration. In order to address this deficiency, Katragadda et al. (2011) andSellmann et al. (2017) proposed the following model for N2 FSD (for unity Lewis number): where Da L =kS L ∕̃δ th is the local Damköhler number. For small Damköhler number flames ∇c tends to preferentially align with the eigenvector associated with e γ due to the dominance of turbulent straining over the effects of heat release and this is scaled as ̃∕k C 1 Σ gen in Eq. 16. Similarly, for flames with high Damköhler numbers, the effects of heat release tend to dominate over the turbulent straining, leading to the alignment of ∇c with the eigenvector associated with e α , which is scaled as − C 2 Σ gen S L ∕δ th in Eq. 16. In the (15i) N FSD = − e cos 2 + e cos 2 + e cos 2 |∇c| (15ii) S UR = e sin 2 + e sin 2 + e sin 2 |∇c| (16) N2 FSD = ̃∕k C 1 − C 2 Da L Σ gen present analysis, the model constants C 1 and C 2 have been modified in order to incorporate the effects of local turbulent Reynolds number Re L and Froude number Fr on the turbulent straining. The model parameter C 1 can be expressed as: Similarly, the effects of Re L and Fr = 1∕ √ g * are included into the strain rate due to heat release using the model parameter B 3 as given below: The predictions of Eqs. 16-20 are compared to N2 FSD extracted from DNS data in Fig. 12, which shows that the model satisfactorily captures the qualitative behaviour of the unresolved normal strain rate term N2 FSD , although there is an over-prediction (underprediction) of the magnitude of N2 FSD towards the reactant side of the flame brush for positive (negative) g * values. It is also evident that the predictions of the model improve with decreasing Damköhler number (Set-A to Set-D).
Combining the MCPB model for N1 FSD with Eqs. 10, 11 and 16-20 for D1 FSD , D2 FSD and N2 FSD respectively, the strain rate term T 2 may be written as: The predictions of Eq. 21, along with the predictions of the CPB (MCPB + SCPB) and CFM (VPDM + SCFM) models are shown in Fig. 13 along with the T 2 term extracted from DNS data. It can be clearly observed that the predictions of Eq. 21 agree better with the DNS data than the CPB and CFM models. However, Eq. 21 tends to underpredict the magnitude of T 2 towards the product side, mainly for the cases with positive g * values.

Modelling of the Combined Propagation and Curvature Terms ( T 3 + T 4 )
The propagation and curvature terms are often modelled together due to their dependence on the displacement speed S d . The combined ( T 3 + T 4 ) term is usually modelled together as B 3 = P∕erf Re L + 1.0 ∕a 3 P = 18.9 − 17.0erf (g * + 1.44) for g * ≥ 0.0 P = 1.67 + 0.98erf (g * + 1.36) for g * ≤ 0.0 a 3 = 16.16 − 12.56erf (g * + 1.79) Duclos et al. 1993;Cant 2011, 2013;Sellmann et al. 2017;Hun and Huh 2008): The third term on the right-hand side of Eq. 22 denotes the unresolved 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 evaluated as α N = 1 − N k s N k s and thus becomes zero under fully resolved condition causing the entire unresolved term to vanish Cant 2011, 2013;Sellmann et al. 2017). Cant (2011, 2013) found that 0 = 8.0 and c cp = 0.35 was successful in capturing the behaviour of ( T 3 + T 4 ) across the flame Variation of N2 FSD × δ 2 Z ∕S L (solid line represents DNS data, dash-dotted line represents the predictions of the model given by  with c for Set-A to Set-D (1st-4th column) for g * = −3.12, −1.56, 0.0, 1.56 and 3.12 (1st-5th row) 1 3 brush for unity Lewis number flames. The variation of the combined contribution of ( T 3 + T 4 ) across the flame brush extracted from DNS data is presented in Fig. 14 along with the predictions of Eq. 22. It can be observed from Fig. 14 that the ( T 3 + T 4 ) term is positive towards the unburned side of the flame brush but becomes negative towards the burned gas side, which is consistent with previous analyses Cant 2011, 2013;Sellmann et al. 2017;Hun and Huh 2008). The magnitude of the ( T 3 + T 4 ) term decreases with increases in g * and initial u � ∕S L values. The decrease in magnitude of ( T 3 + T 4 ) with increasing g * is more pronounced for cases with lower initial u � ∕S L (i.e., Set-A and Set-B). It is also evident from Fig. 14 that the model given by Eq. 22 accurately predicts the qualitative behaviour of the combined propagation and Fig. 13 Variation of T 2 × δ 2 Z ∕S L (solid line represents DNS data, dotted line represents the predictions of the CPB model, dashed line represents the predictions of the CFM model and line with circles represents the predictions given by Eq. 21) with c for Set-A to Set-D (1st-4th column) for g * = −3.12, −1.56, 0.0, 1.56 and 3.12 (1st-5th row) curvature term for all cases considered in the present study irrespective of the values of g * and initial u � ∕S L .

Mean Reaction Rate Closure Using Generalised FSD
According to Boger et al. (1998), one obtains: where the two terms on the left-hand side represent the contributions of the reaction rate and molecular diffusion rate, respectively. In the context of RANS simulations, the contribution of the mean molecular diffusion rate in statistically planar flames is much smaller than the mean reaction rate (i.e., ẇ ≫ ∇.(ρD∇c) ) and hence ẇ ≈ ρS d s Σ gen . For unity Lewis number flames, the surface-averaged density-weighted displacement speed is usually (23) w + ∇.(ρD∇c) = ρS d s Σ gen Fig. 14 Variation of the combined propagation and curvature term (T 3 + T 4 ) × 2 Z ∕S L (solid line represents the DNS data and dash-dotted line represents the predictions of Eq. 22) with c for Set-A to Set-D (1st-4th column) for g * = −3.12, −1.56, 0.0, 1.56 and 3.12 (1st-5th row) modelled as ρS d s ≈ 0 S L (Boger et al. 1998;Cant et al. 1990;Hawkes and Cant 2000;Hernandez-Perez et al. 2011) although this approximation may not be valid, especially for the thin reaction zones regime flames Cant 2007, 2009a;Sabelnikov et al. 2017). Thus, the model for mean reaction rate may be expressed as ẇ = ρ 0 S L Σ gen and this model will henceforth be referred to as the RRFSD model.
The predictions of the RRFSD model are shown in Fig. 15 alongside the mean reaction rate ẇ and the combined contribution of the mean values of reaction rate and molecular diffusion rate ẇ + ∇.( D∇c) extracted from the DNS data. It is evident from Fig. 15 that ∇.(ρD∇c) indeed remains negligible when compared to ẇ in most cases (thus ẇ ≈ ρS d s Σ gen remains valid) but the contribution of ∇.(ρD∇c) cannot be ignored for negative values of g * for small turbulence intensities (e.g., g * = −3.12 and − 1.56 for Set-A). The effects of turbulence weaken with increasingly negative values of g * (see Figs. 1,   Fig. 15 Variation of ẇ × δ Z ∕ρ 0 S L (solid line represents the DNS data, line with circles represents combined contribution of mean reaction rate and molecular diffusion rate, and dotted line represents the predictions of the RRFSD model) and ẇ + ∇.( D∇c) × δ Z ∕ρ 0 S L (line with circles) with c for Set-A to Set-D (1st-4th column) for g * = −3.12, −1.56, 0.0, 1.56 and 3.12 (1st-5th row) 2, 3), and this tendency is particularly strong for small values of u � ∕S L e.g., g * = −3.12 and − 1.56 for Set-A) and thus in these cases ẇ ≫ ∇.( D∇c) do not remain strictly valid. Figure 15 shows that the RRFSD model mostly satisfactorily predicts ẇ extracted from DNS data but under-predicts the magnitude of ẇ for all cases considered, especially towards the product side of the flame brush. This under-prediction is particularly prevalent for negative values of g * and the agreement of the prediction of the RRFSD model and DNS data improves with increasing g * and initial u � ∕S L values.
The variation of I 0 =ẇ∕ 0 S L Σ gen with c for the cases considered here are shown in Fig. 16, which shows that the stretch factor remains close to unity throughout the flame brush for most cases considered here except for small turbulence intensities with negative values of g * (e.g., Set-A and Set-B for g * = −3.12 and − 1.56). It is evident from Fig. 16 that the discrepancy between the RRFSD model and ẇ originates due to the approximation of ρS d s ≈ 0 S L for most cases and a better representation of ρS d s is likely to provide an improved prediction of the mean reaction rate. Moreover, the assumption ẇ ≫ ∇.( D∇c) Fig. 16 Variation of stretch factor I 0 = ẇ∕ 0 S L Σ gen with c (between 0.1 and 0.9) for Set-A to Set-D (1st-4th column) for g * = −3.12, −1.56, 0.0, 1.56 and 3.12 (1st-5th row) is rendered invalid for small values of u � ∕S L in the case of negative g * values (e.g., Set-A with g * = −3.12 and 1.56).This also contributes to the disagreement between the RRFSD model and ẇ for small values of u � ∕S L in the case of negative g * values.

Conclusions
The effects of body force on the statistical behaviours of FSD and its transport have been evaluated using a DNS database of three-dimensional statistically planar turbulent premixed flames with unity Lewis number for different sets of initial u � ∕S L values. A range of different values of normalised body forces (or Froude numbers) have been considered representing both unstable and stable configurations in terms of density stratification for each set of initial turbulence parameters. The body force is directed from the heavier unburned gas to the lighter burned gas side in the case of unstable density stratification, whereas the body force is directed from the lighter burned gas side to the heavier unburned gas side of the flame in the case of stable density stratification. The unstable (stable) density stratification configuration can alternatively be taken to represent the situations where the flame is subjected to an externally imposed adverse (favourable) pressure gradient. It has been found that an unstable density stratification promotes flame wrinkling and a gradient type transport, especially towards the reactant side of the flame brush where heat release effects are weak. By contrast, the flame wrinkling weakens, and a counter-gradient type transport is promoted for a stable density stratification. The effects of body forces on the behaviour of the various terms of the FSD transport equation have been analysed in detail and it is observed that the tangential strain rate term T 2 acts as a source and the curvature term T 4 remains the major sink term towards the burned gas of the flame brush for all cases considered here. The effects of body force on the FSD and the terms of its transport equation weakens with increasing turbulence intensity as the relative strength of body force in comparison to the inertial force weakens with an increase in u � ∕S L . The magnitudes of T 2 and T 4 generally increase with the strengthening of body force from the stable to the unstable density stratification. The performance of the existing models for the tangential strain rate term has been found to be inadequate especially for small values of turbulence intensity under stable density stratification. Based on a-priori analysis of DNS data, an alternative model for the tangential strain rate term T 2 has been proposed by considering the alignment of ∇c with local principal strain rates and the competition between turbulent straining and chemical timescales. In this new modelling methodology of T 2 , the effects of turbulent Reynolds number and Froude number have also been explicitly accounted for. It has also been shown that the existing models for the turbulent flux of the FSD and the combined propagation and curvature term T 3 + T 4 do not need any modifications due to the presence of the body force. Moreover, the FSD based mean reaction rate closure mostly works satisfactorily irrespective of the strength of the body force but under-predicts the magnitude of ẇ for all cases towards the product side of the flame brush. This under-prediction is particularly prevalent for stable density stratification configuration, but this situation improves for the unstable density stratification configuration and also with increasing turbulence intensity. This study has been conducted with simple chemistry in the interest of a detailed parametric analysis. Although previous analyses (Chakraborty and Klein 2008b;Chakraborty et al. 2008) revealed no major differences in the qualitative behaviours of the |∇c| transport equation between simple and detailed chemistry, future analyses need to include the effects 1 3 of detailed chemistry and transport in order to conduct further validation of the models proposed in this analysis.
Finally, it is important to note that a degree of empiricism is involved in Eqs. 12, 13, 17-20 but they satisfy all the asymptotic requirements in terms of Re L and g * . Thus, further validation of these models based on a-posteriori assessment will be necessary. It is important to note that the unclosed quantities, such as k and ̃ , also act as input parameters to the models for the unclosed terms of the FSD transport equation. Thus, turbulence modelling (e.g., the models for k and ̃ ) also interacts with the models of the FSD transport equation in actual RANS simulations. In addition to this, the numerical scheme used in the RANS code is expected to influence the simulation predictions. Thus, it is not straightforward to assess the performances of the terms of the FSD transport equation in isolation. Thus, a model, which fails based on a-priori analysis, is unlikely to perform well in actual RANS simulation but the models, which exhibit promising performance based on a-priori DNS analysis, need to be assessed whether they perform satisfactorily while interacting with turbulence models and numerical schemes in actual RANS simulations. This necessitates a detailed a-posteriori assessment of the newly developed models based on actual RANS simulations for a configuration with well-documented experimental data.