#### Abstract

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

#### 1. Introduction

Flame surface density (FSD) based reaction rate closure is one of the most well-established methods of turbulent premixed combustion modelling in the context of Reynolds averaged Navier Stokes (RANS) simulations [1–12]. In the FSD based formulation, the closure of reaction rate translates to the modelling of flame surface area to volume ratio [13]. The FSD, in RANS simulations, is either evaluated using an algebraic expression [1] or a modelled transport equation for FSD is solved alongside other modelled Favre averaged conservation equations [3–12]. However, algebraic closures may not be suitable when flame surface area generation and destruction are not in equilibrium and recently Richard et al. [14] demonstrated that algebraic closures may not be adequate for simulating in-cylinder processes in spark ignition piston engines. Therefore, it may be necessary to consider high-fidelity transport equation based closures of FSD for turbulent premixed flames. The flame surface area generation/destruction due to fluid straining plays a key role in the FSD transport and thus the statistical analysis and detailed assessments of the modelling of the tangential strain rate term in the FSD transport equation in the context of RANS simulations have been considered in the current work because the strain rate term remains a leading order contributor to the FSD transport [3–12].

The generalised FSD is defined as [15, 16]
where is the reaction progress variable and the quantities and denote the Reynolds averaged and Favre averaged values of a general quantity . The generalised FSD is often used in order to close the mean reaction rate of reaction progress variable in the following manner [15]:
where is the gas density, is the reaction progress variable diffusivity, is the displacement speed, and indicates a surface averaging operation of a general quantity [5, 15]. In the context of RANS simulations and therefore it is evident from (2) that an accurate prediction of is required for the closure of the mean reaction rate provided the statistical behaviour of the surface averaged density-weighted displacement speed (i.e., ) is satisfactorily accounted for. The exact transport equation for the generalised FSD (i.e., ) takes the following form [9, 11, 12, 17, 18]:
where and are the Favre mean and fluctuating velocity components in the th direction and is the local flame normal vector (which points towards the reactants). The reaction progress variable can be defined in terms of a suitable reactant mass fraction in such a manner that it increases monotonically from zero in unburned gases to unity in fully burned products. In (3), the first term in the left hand side is the transient term, whereas the second term represents the effects of mean advection. All the terms on the right hand side of (3) are unclosed and, hence, require modelling. The first term on the right hand side of (3) represents the effects of turbulent transport, the second term accounts for the flame area generation due to tangential strain rate (i.e., ), and the third term originates due to flame normal propagation and thus is commonly referred to as the propagation term. The final term on the right hand side of (3) (i.e., ) arises due to curvature (according to this convention a flame element, which is convex towards the reactants, has a positive curvature and *vice versa*) and thus is commonly referred to as the FSD curvature term. In the current study, only the *a priori* DNS assessment of the models of the tangential strain rate term of the generalised FSD transport equation will be considered.

The tangential strain rate term of the generalised FSD transport equation is often split in the following manner: Cant et al. [2] and Candel et al. [3] proposed models for both the resolved and the unresolved parts of the tangential strain rate contribution (i.e., and , resp.), which were subsequently assessed by Duclos et al. [4] and Prasad and Gore [8] based on RANS simulations. Veynante et al. [6] and Veynante et al. [7] assessed the performance of the models of and based on experimental data. The aforementioned studies [2–4, 6–8] were carried out in the corrugated flamelets regime, where the flame thickness remains smaller than the Kolmogorov length scale. Recently, Katragadda et al. [10] assessed the performances of the existing and models in both the corrugated flamelets and thin reaction zones regimes of premixed combustion based on direct numerical simulation (DNS) data, while accounting for the effects of Lewis number and heat release parameter . It was shown by Katragadda et al. [10] that both the dilatation rate and the relative alignment of with the fluid-dynamic strain rate significantly affect the behaviour of the strain rate term . Global Lewis number is known to have significant influences on both the dilatation rate and scalar gradient alignment statistics in turbulent premixed flames [19]. Additionally, it was shown in recent studies [19–21] that the Damköhler number significantly affects the relative strength of the dilatation rate in comparison to the turbulent straining and alignment characteristics of with the local principal strain rates. As the turbulent Reynolds number scales as [22] it is expected that the modelling of the strain rate term will also show some dependence.

The effects of Lewis number , turbulent Reynolds number , Damköhler number , Karlovitz number , and heat release parameter on the statistical behaviour and RANS modelling of the tangential strain rate term in the FSD transport equation are yet to be addressed in detail and the present paper aims to address this gap in the existing literature. In order to investigate the effects of , , , , and on the statistical behaviour and RANS modelling of the FSD tangential strain term, a database of three-dimensional compressible direct numerical simulations (DNS) of turbulent premixed flames has been considered. In this respect, the main objectives of this study are as follows:(1)to demonstrate and explain the statistical behaviours of the tangential strain rate term of the FSD transport equation and its components for a range of different values of Lewis numbers, turbulent Reynolds numbers, Karlovitz and Damköhler numbers, and heat release parameters spanning both the corrugated flamelets and thin reaction zones regimes of turbulent premixed combustion;(2)to identify models which satisfactorily capture the statistical behaviour of the tangential strain rate term in the FSD transport equation and its components for different values of , , and and different combustion regimes in the context of RANS simulations.

The remainder of the paper will take the following form. The next section will discuss the relevant mathematical background of the current study, which will then be followed by a brief discussion on the numerical implementation. Following this, the results will be presented and subsequently discussed. Finally, the main findings will be summarised and conclusions will be drawn.

#### 2. Mathematical Background

Direct numerical simulations (DNS) of turbulent combustion should, ideally, account for both the inherent three-dimensionality of turbulence and the detailed chemical structure of the flame. However, until recently combustion DNS studies were limited to be carried out either in two dimensions with detailed chemistry or in three dimensions with simple chemistry. Three-dimensional DNS of turbulent combustion is now possible but such simulations remain extremely computationally expensive [23] and thus are often not suitable for a detailed parametric analysis as carried out in this work. As such, a three-dimensional DNS database with simple chemistry has been considered here for the sake of computational economy. In the context of simple chemistry, the species field can be represented, uniquely, by a reaction progress variable which can be defined in terms of a suitable mass fraction (i.e., ) in the following manner: where subscripts 0 and were used to denote the values in the unburned and burned gases, respectively. In globally adiabatic, low Mach number unity Lewis number flames the reaction progress variable is the same as the nondimensional temperature: where is the instantaneous dimensional temperature, is the unburned gas temperature, and is the adiabatic flame temperature. However, for nonunity Lewis number flames the nondimensional temperature may locally assume superadiabatic values (i.e., ) [24, 25] even under globally adiabatic conditions, whereas is always bound within 0 and unity (i.e., ).

As discussed in Section 1, the strain rate term is often modelled by splitting it into the resolved and the unresolved components. For the purpose of evaluating , the quantity requires modelling. Cant et al. [2] modelled as where is the model expression for , which takes the following form according to Cant et al. [2]: Veynante et al. [7] modelled as where is the turbulent kinetic energy. Cant et al. [2] modelled for the unresolved part of the strain rate term in the following manner: where and are the unburned gas density and viscosity, respectively, and is the dissipation rate of . In the context of coherent-flamelet modelling (CFM) [4] is modelled as where is a model constant of the order of unity (i.e., ) and is the efficiency function proposed by Meneveau and Poinsot [26] which is a function of and with , , and being the thermal diffusivity in the unburned gas, unstrained laminar burning velocity, and local integral length scale, respectively.

The strain rate term can alternatively be decomposed in the following manner [10]: The terms and () represent contributions of dilatation rate and flame normal strain rate on the transport. The dilatation rate term can be split into the resolved and unresolved components in the following manner: From the expression above, it is evident that the resolved dilatation rate term can be closed if a suitable relationship between and can be found. The contribution of the normal strain rate term () can be split as The quantity () can be closed using the models for . This can be carried out as shown in (8) and (9).

It is worth noting that the surface averaged quantities involving velocity components in (4) are modelled in terms of the gradients of unconditional Favre-averaged velocity components (see (7), (9), (10), and (11)). In actual RANS simulations the unconditional Favre-averaged velocities are readily obtained so it makes sense to model the surface-averaged terms involving velocity components using the Favre-averaged velocity components [2–8]. The modelling of () and (see (12) and (14)) will be discussed based on *a priori* DNS analysis using the current database in Section 4, where the surface averaged quantities involving velocity components are also modelled in terms of the gradients of unconditional Favre-averaged velocity components (see (15) and (17) later in this paper). The present formulation is based on the generalised FSD , which is defined as [15, 16]: . This indicates that the generalised FSD is not dependent upon choice of isosurface. Thus the difference between conditional and unconditional velocities does not play any role in the current analysis similar to several previous studies [2–12].

#### 3. Numerical Implementation

For the present study, the simulation parameters for the DNS database considered are shown in Table 1, namely, the initial values of root-mean-square turbulent velocity fluctuation normalised by the unstrained laminar burning velocity , integral length scale normalised by the thermal flame thickness , turbulent Reynolds number , heat release parameter , Lewis number , Damköhler number , and Karlovitz number . Standard values are chosen for Zel’dovich number , Prandtl number , and the ratio of specific heats (i.e., , , and ), where is the activation temperature. It should be noted that the simulation parameters have been chosen in such a manner that combustion takes place within the corrugated flamelets regime (thin reaction zones regime) [22] for case A (cases B–L). The turbulent flow conditions and heat release parameter are kept unaltered but the global Lewis number was modified from 0.34 to 1.2 in cases C–G. The cases B and F are identical in terms of turbulent flow conditions but only values are different between these cases. As the turbulent Reynolds number scales as [22], the variation of in cases H–L is brought about by modifying and independently of each other. In cases H, J, and L (I, J, and K) () is kept unaltered and () is altered to bring about the variation of . The range of considered here is comparable to that of previous laboratory-scale experiments (e.g., Kobayashi et al. [27]).

A domain of is taken for case A, which is discretised by a Cartesian grid of with uniform grid spacing in each direction [28]. In case A, inlet and outlet boundaries are specified in the mean direction of flame propagation (which is aligned with the negative -direction), whereas transverse boundaries are taken to be periodic. In cases B–G, a domain of size is discretised using a uniform grid of . In cases H–L, a domain of size is discretised using a uniform grid of . The domain boundaries in the direction of mean flame propagation in cases B–L are taken to be partially nonreflecting and the transverse boundaries are assumed to be periodic. In case A, a 6th order central-difference scheme has been used for spatial discretisation for the internal grid points in the direction of mean flame propagation, which gradually reduces to a one-sided 4th order scheme near nonperiodic boundaries, whereas a spectral method is used for spatial discretisation in the directions normal to the mean direction of flame propagation [28]. In cases B–L, a 10th central difference scheme is used for the internal grid points and the order of differentiation gradually reduces to a 2nd order one-sided scheme near nonperiodic boundaries. The time advancement for all viscous and diffusive terms in case A is carried out using an implicit solver, whereas the convection terms in case A and all the terms in cases B–L are time advanced with the help of a third order low storage Runge-Kutta method [29]. For all cases, the flame is initialised by a steady unstrained planar laminar flame solution and the turbulent fluctuating velocity field is initialised based on an incompressible homogeneous isotropic velocity distribution, which is generated using a standard pseudospectral method [30]. The grid spacing is determined by the resolution of the flame structure, and about 10 grid points are kept within for all cases considered here. In all cases flame-turbulence interaction takes place under decaying turbulence and simulations need to be carried out for , where is the initial eddy turn-over time and is the chemical time scale. The simulation in case A was run for about , whereas in cases B–G it was run for a time equivalent to . The simulation time remains either greater than (case A) or equal to (cases B–L) one chemical time scale. The simulation time remains either comparable to or greater than several previous studies [9, 31–35]. The turbulent kinetic energy and its dissipation rate in the unburned gas ahead of the flame were not varying significantly with time when statistics were extracted for all cases. Interested readers are referred to Chakraborty and Cant [11, 12] and Chakraborty et al. [36, 37] for further information on this database and the conditions under which statistics were extracted.

The relevant quantities in all cases considered here have been Reynolds/Favre-averaged by ensemble averaging the relevant quantities in transverse directions ( planes). The statistical convergence of the Reynolds/Favre-averaged quantities is examined by comparing the corresponding values obtained using half of the sample size from the distinct half of the domain in the transverse direction with those obtained based on full sample size. Both the qualitative and quantitative agreements between these sets of values are found to be satisfactory and more information can be obtained from [11]. In the next section, only the results obtained based on full sample size will be presented. For statistically planar flames, remains a unique function of the spatial coordinate in the direction of mean flame propagation (i.e., -direction); thus all the statistics in Section 4 will be presented as a function of .

#### 4. Results and Discussion

The contours of the reaction progress variable in the central plane when the statistics are extracted are shown in Figure 1 for cases A–L. In case A the contours of remain parallel to each other as a consequence of the corrugated flamelets regime combustion. In the corrugated flamelets regime the energetic turbulent eddies are unable to penetrate into the flame and wrinkles develop due to large-scale turbulent motion. However, in cases B–L the contours of are found to be distorted towards the unburned gas side (i.e., ) representing the preheat zone. This distortion of the preheat zone in cases B–L is characteristic of the thin reaction zones regime combustion. Additionally the contours of representing the reaction zone (i.e., ) remain parallel to each other for all cases. In cases B–L, where combustion takes place in the thin reaction zone regime, the energetic turbulent eddies penetrate into the flame, as the flame thickness remains greater than the Kolmogorov length scale though the reaction zone retains its quasi-laminar structure due to where and are the Kolmogorov length scale and the reaction zone thickness, respectively.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

**(h)**

**(i)**

**(j)**

**(k)**

**(l)**

The effects of Lewis number on the flame structure can be seen in Figures 1(c)–1(g), where the extent of the flame wrinkling increases as decreases. The reactants diffuse more rapidly into the reaction zone than the rate at which heat is conducted out for the flames with , which leads to simultaneous presence of high concentration of reactants and high temperature. As a result, the cases with exhibit greater burning rates than the unity Lewis number flames with the same turbulent flow conditions in the unburned gas. Just the opposite mechanism in terms of reactant and thermal diffusion takes place for flames, which reduces the flame area generation and burning rate in these flames in comparison to the corresponding unity Lewis number flames with the same unburned gas turbulence. This gives rise to increased burning rates and flame area generation with decreasing , which can be substantiated from Table 2 where the normalised turbulent flame speed and normalised flame surface area at the time when the statistics were extracted are presented. The values of have been evaluated by volume integrating the reaction rate using the expression where is the projected area of the flame in the direction of mean flame propagation, while the values of have been evaluated by volume integrating (i.e., ) under both turbulent and laminar conditions. The preheat zone () becomes increasingly distorted due the increased scale separation between and as a result of increasing Karlovitz number when Re_{t} increases for a given value of . Thus, the preheat zone is more susceptible to be penetrated by energetic turbulent eddies with increasing scale separation between and . It is clear from Figures 1(h)–1(l) that the flame wrinkling increases with increasing , which gives rise to high values of normalised turbulent flame speed and normalised flame surface area for high values of .

The variations of , , , , and () with for all cases are shown in Figures 2(a)–2(l), which show that , , , and assume positive values throughout the flame brush and the contribution of supersedes that of in all cases, in accordance with earlier findings [2–4, 6–8, 18, 38, 39]. As remains principally positive within the flame brush and is a positive semidefinite quantity (i.e., ), the dilatation rate term remains positive throughout the flame brush. It can be seen from Table 2 that the effects of heat release strengthen with decreasing so the magnitude of increases significantly with decreasing , which along with the increasing trend of with decreasing [11] leads to high magnitudes of for flames with small values of (see Figures 2(c)–2(g)).

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

**(h)**

**(i)**

**(j)**

**(k)**

**(l)**

It can be seen from Figure 2 that the normal strain rate term () exhibits marked differences from one case to another. The term () assumes negative value throughout the flame brush for cases A, C, D, E, H, I, J, and K, whereas this term assumes a positive (negative) value towards the unburned (burned) gas side of the flame brush in cases B, F, G, and L. The behavior of is determined by the alignment of with the local principal strain rates where , , and are the most extensive, intermediate, and compressive principal strain rates, respectively, and , , and are the angles between and , , and , respectively. For passive scalars, the scalar gradient is known to align with the most compressive principal strain rate [28, 40, 41] but recent studies [19–21, 42] have shown that may align with the most extensive principal strain rate , where the strain rate induced by flame normal acceleration overcomes the effects of turbulent straining.

The strain rate induced by flame normal acceleration can be scaled as for flames, whereas the turbulent straining can be scaled as , which gives rise to [20] with being a function which accounts for the weakening of the effects of flame normal acceleration with increasing as the broken reaction zones regime combustion is approached. Alternatively, the turbulent straining can be scaled as following Tennekes and Lumley [43] where is the Taylor microscale, which leads to . The scaling () indicates that may align with to yield a positive (negative) contribution of () for large flames. By contrast, in cases B, F, G, and L, remains small (i.e., ), so that the effects of overcome the effects of , giving rise to the predominant alignment of with for the major portion of the flame brush, except in the heat releasing zone where the effects of overcome the effects of . This is reflected in the predominantly negative (positive) value of () towards the unburned gas side of the flame brush, and positive (negative) values of () are obtained towards the burned gas side of the flame brush in cases B, F, G, and L. Chakraborty et al. [19] demonstrated that strengthens with decreasing and this effect is particularly strong for the flames due to high values of chemical heat release. As a result, can be scaled as () for nonunity Lewis number flames, where the function increases with decreasing [19]. This suggests that may become sufficiently strong to overcome the effects of for small values of , even for flames. Strong in the , , and flames overcomes the effects of and induces a preferential alignment of with . This gives rise to positive (negative) values of () for the major portion of the flame brush. As the effects of weaken with an increase in , the extent of the negative contribution of () progressively decreases with increasing , which can be substantiated from Figures 2(c)–2(g). Figures 2(a)–2(l) show that the magnitude of () remains smaller than the magnitude of for all cases, which yields a positive contribution of throughout the flame brush even when the contribution of () remains negative. In all cases, remains of the same order of magnitude as that of (i.e., or ) and thus the effects of are partially nullified by the effects induced by , which give rise to a smaller magnitude of the normal strain rate contribution () than the magnitude of the dilatation rate term .

It is worth noting that the variations of some of the strain rate terms with in Figure 2 and subsequent figures show wavy behaviour although the variation of with for all cases has been found to be reasonably smooth. However, the variation of with exhibits weak undulations for some cases (e.g., cases A, C, K, and L), which originate due to large extent of wrinkling in these flames. These undulations in variations with give rise to wavy behaviour of (= for statistically planar flames) when a 10th order central difference scheme is used for spatial differentiation. This also contributes to the wavy behaviour of the variations of the strain rate terms with . Similar wavy behaviour has been seen in several previous DNS studies [9, 10, 20, 21, 28, 42]. This is not due to a lack of statistical convergence but is a consequence of the high levels of wrinkling in these flames and high order of finite difference scheme used for postprocessing purpose. Moreover, this waviness is prevalent only for case A, which was not simulated by the present authors but taken from a widely used well-respected DNS database (see [28]). This waviness is also observed in several previous analyses [20, 21, 28, 42, 44], where case A was used.

Katragadda et al. [10] modelled the tangential strain rate term based on the separate modelling of the strain rate components and (). The term can be closed if a suitable relationship can be established between and . According to Bray-Moss-Libby (BML) analysis [45], this can be achieved by for unity Lewis number globally adiabatic flames, where accounts for the contribution of the reacting mixture. The contribution of remains negligible for high flames but it might be nonnegligible for low combustion. Katragadda et al. [10] proposed an alternative empirical expression where is the segregation factor and it was shown in [10–12] that adequately captures for low flames considered in this analysis and thus predictions of are not shown here for the sake of conciseness. Moreover, approaches the BML expression for unity Lewis number flames for combustion under which approaches unity. Thus enables one to close . The unresolved part of the dilatation rate term is taken to model as [10] where the model parameters , , , and . In (15) dilatation rate is scaled as as the strength of the dilatation rate increases with decreasing , and is used to adequately capture the variation of with obtained from DNS where the exponent is likely to be a function of , as the variation of with is skewed with a peak towards the unburned gas side for and flames, whereas the peak value of is attained close to the middle of the flame brush (i.e., ) for the flames with close to unity. According to (15), the contribution of vanishes when the flame is fully resolved (i.e., ) and the net contribution of () becomes identical to . The contribution of is similar to the dilatation rate contribution in the scalar dissipation rate transport equation [42]. The term is modeled as [42] In (16), is taken to be with being a constant of the order of unity and is the local Karlovitz number. The local Karlovitz number dependence of ensures that the strength of the dilatation rate contribution to the scalar dissipation rate transport diminishes with increasing as combustion starts to exhibit the attributes of the broken reaction zones regime [26]. Following the same procedure, the local Karlovitz number dependence is used in (15) in the form of . However, this local Karlovitz number dependence is one of the several possibilities and any other function of , which accurately predicts the quantitative behaviour of and its diminishing strength with increasing Karlovitz number as the broken reaction zones regime is approached, can in principle be used for parameterising the local dependence of . The numerical values of , , , and are optimised based on a least-squares method.

The predictions of (15) are shown in Figures 3(a)–3(l) along with obtained from the DNS data for cases A–L, respectively. The model performs satisfactorily for a range of different values of , , and spanning both the corrugated flamelets and thin reaction zones regimes of premixed turbulent combustion.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

**(h)**

**(i)**

**(j)**

**(k)**

**(l)**

The variations of () and () are shown in Figure 4 with for cases A–L. The predictions for () according to the models proposed by Cant et al. [2] and Veynante et al. [7] are also shown in Figure 4 for cases A–L. The predictions of () according to Cant et al. [2] and Veynante et al. [7] models (denoted as TN1CPB and TN1V models, resp.) are shown in Figure 4, which shows that the performances of both of the models are comparable, and both satisfactorily predict () for all cases considered in the current analysis. However, the prediction of the TN1V model is closer to the DNS data than the TN1CPB model. The model given by (8) overpredicts the unresolved part of (i.e., ) due to the assumption of isotropy of the unresolved flame normal fluctuations. However, experimental data [6, 7] suggested that the assumption of isotropy of the unresolved flame normal fluctuations does not hold, and this result has been verified here by the overprediction of () by the TN1CPB model. It should be noted that reverts to when the flame is fully resolved according to the TN1CPB model, but this condition is not satisfied by the model given by (9).

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

**(h)**

**(i)**

**(j)**

**(k)**

**(l)**

It is evident from that a preferential alignment between and () leads to a negative (positive) contribution of . As the magnitude of remains either greater than or comparable to for all cases, it can be expected that alignment with local principal strain rates affects not only but also . The alignment of with local principal strain rates is determined by the competition between the strain rate due to turbulent motion and the strain rate induced by flame normal acceleration. The scalar gradient aligns with () when () dominates over (). According to Tennekes and Lumley [43] can be scaled as , whereas can be scaled as . Thus, the effects of on the term () through the preferential alignment of with can be scaled as , where the model parameter is expected to be a function of local turbulent Reynolds number . By contrast, the effects of alignment with on the term () due to the action of can be scaled as , where is a function of which increases with decreasing due to strengthening of flame normal acceleration as a result of augmented burning rate, and the model parameter is expected to decrease with increasing because the effects of are likely to weaken progressively with increasing . Combining the above scaling estimations, the following model for () can be proposed: where is the local Damköhler number and is given as while the following two parameterisations are proposed for the model parameter :where in is given as and from and is defined as follows: The parameterisations of given by (18) include turbulent Reynolds number dependence but this model parameter needs to be turbulent Reynolds number independent for high values of (i.e., ). Equation (18) ensures that assumes an asymptotic value for . Moreover, according to and ensures that the effects of on the term () weaken with increasing , whereas according to (21) accounts for strengthening of flame normal acceleration with decreasing . However, these local Karlovitz number dependences according to and are only two possibilities and any other function of , which accurately predicts the quantitative behaviour of and the diminishing strength of on the term with increasing Karlovitz number , can in principle be used for . Similarly, any other function which satisfies the expected behaviour of () in response to () can in principle be used to parameterise () instead of (18) (see (21)). The presence of in and ensures that the second term on the right hand side of (17) vanishes when the flow becomes fully resolved (i.e., ). The numerical values of , , , and are optimised based on a least-squares method.

In Figure 4, the predictions based on (17) using according to and are compared with () obtained from DNS. Figure 4 shows that (17) performs satisfactorily for the range of and considered here. Using (see (8)) and combining (15) and (17) yield the following expression: It should be noted that the involvement of in the first term on the right hand side of (22) and in ensures that becomes equal to when the flow becomes fully resolved (i.e., ). The predictions of according to (22) based on the parameterisations given by and are shown in Figures 5(a)–5(l) for cases A–L, respectively, along with the predictions of the models proposed by Cant et al. [2] (i.e., CPB model) and CFM methodology [4, 6, 7] where and are modeled using (8) and (10) in the context of the CPB model and by (9) and (11) in the context of the CFM model, respectively. The CFM model is found to show a peak at which is quantitatively similar to cases E–G and J–L, while in cases H and I the peak value of occurs towards the products (i.e., ) and in cases C and D the peak value of occurs towards the reactants (i.e., ). Although the CFM model performs adequately for cases F, G, and J, it fails to capture the behavior of in cases H, I, K, and L. From Figure 5, it is evident that the model CPB can be seen to predict the correct qualitative behavior for most of the cases but in cases H and I it fails to capture the correct magnitude of . The new model performs well for all cases considered here, and the difference between the two parameterisations for can be seen to make very little difference. The newly proposed model given by (22) with given by and also captures the behaviours of satisfactorily for a range of different values of , , and spanning both the corrugated flamelets and thin reaction zones regimes of premixed turbulent combustion.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

**(h)**

**(i)**

**(j)**

**(k)**

**(l)**

#### 5. Conclusions

A DNS database of freely propagating statistically planar turbulent premixed flames with a range of different values of turbulent Reynolds number , heat release parameter, and global Lewis number spanning both the corrugated flamelets and thin reaction zones regimes of premixed combustion has been used to assess the modelling of the tangential strain rate term of the generalised FSD transport equation. In order to propose a model for , which will be valid for both the corrugated flamelets and thin reaction zones regimes, the statistical behaviours of the strain rate contributions to the FSD transport due to the dilatation rate and flame normal strain rate have been analysed separately.

It has been shown that the dilatation rate contribution to the generalised FSD transport equation strengthens with decreasing Lewis number. The contribution of normal strain rate term () assumes predominantly negative values for cases with and where the effects of straining due to flame normal acceleration dominate over turbulent straining , while in the cases with the normal strain rate term () assumes positive (negative) values towards the reactants (products) side of the flame brush. This analysis gave rise to the explicit modelling of the unresolved dilatation rate and flame normal strain rate contributions. The predictions of the new model have been compared with the predictions of the existing models (i.e., CPB and CFM models) and the performance of the newly developed model has been found to be either comparable to or better than the existing models for the FSD strain rate term.

The present analysis has been carried out using a simple chemistry based DNS database with moderate range of values of turbulent Reynolds number . Although the model parameters have been found to attain asymptotic values for , further validations will be required based on experimental and three-dimensional detailed chemistry DNS data at higher values of . Finally *a posteriori* assessment of the new models needs to be carried based on actual RANS simulations in a configuration for which well-documented experimental data is available, which will form the basis of future communications.

#### Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

#### Acknowledgment

The financial support of EPSRC is gratefully acknowledged.