Sub-grid Reaction Progress Variable Variance Closure in Turbulent Premixed Flames

Abstract

The statistical behaviour and modelling of the sub-grid variance of reaction progress variable have been analysed based on a priori analysis of Direct Numerical Simulation (DNS) data of freely propagating statistically planar turbulent premixed flames with different turbulence intensities. It has been found that an algebraic expression, which can be derived based on a presumed bi-modal sub-grid distribution of reaction progress variable with impulses at unburned reactants and fully burned products, is inadequate for the purpose of prediction of sub-grid variance. An algebraic model, which is often used for modelling sub-grid variance of passive scalars, has been found to significantly overpredict the sub-grid variance of reaction progress variable for the cases considered here. The statistical behaviours of the terms of the sub-grid variance have been analysed in detail and explained based on scaling arguments. It has been found that reaction rate and molecular dissipation contributions play leading order roles in the transport of sub-grid reaction progress variable variance and they remain in approximate equilibrium for large filter widths. Suitable model expressions have been identified for the sub-grid flux of variance, reaction rate contribution and scalar dissipation rate based on a priori analysis of DNS data.

1 Introduction

The closure of sub-grid scalar variance plays a crucial role in the modelling of micro-mixing in the context of large eddy simulations (LES) (Pierce and Moin 1998; Jimenez et al. 2001). The knowledge of sub-grid scale (SGS) variance of reaction progress is often necessary to construct the sub-grid probability density function (PDF) of reaction progress variable $c$ in the context of flamelet (Langella and Swaminathan 2016; Langella et al. 2016) and Linear Eddy Modelling (Ranjan et al. 2016) based modelling methodologies. The SGS variance of reaction progress variable $c$ is defined as (Pierce and Moin 1998; Jimenez et al. 2001; Langella and Swaminathan 2016; Langella et al. 2016; Ranjan et al. 2016):

$$ \sigma_{v}^{2} = \widetilde{{c^{2} }} - \tilde{c}^{2} $$

(1)

where $\tilde{q} = \overline{\rho q} /\bar{\rho }$ and $\bar{q}$ are the Favre-filtered and LES-filtered values of a general variable $q$, respectively and $\rho$ is the gas density. Based on a presumed bi-modal PDF of $c$ with peaks at 0.0 and 1.0, the SGS variance of reaction progress variable $c$ can be expressed as (Bray et al. 1985): $\sigma_{v}^{2} = \tilde{c}\left( {1 - \tilde{c}} \right)$ but the sub-grid PDF of $c$ is unlikely to be bi-modal in practice (Chakraborty and Cant 2007, 2009; Dunstan et al. 2013; Gao et al. 2014; Ma et al. 2014). An alternative algebraic expression is often used for the evaluation of SGS variance of passive scalars in the following manner (Pierce and Moin 1998; Jimenez et al. 2001; Ranjan et al. 2016):

$$ \sigma_{v}^{2} = C\Delta^{2} \left| {\nabla \tilde{c}} \right|^{2} $$

(2)

where $C = 0.5$ is the model constant (Ranjan et al. 2016) and $\Delta$ is the LES filter width. It is worth noting that the evolution of SGS variance is expected to be influenced by both chemical and turbulent processes (Chakraborty and Swaminathan 2011; Lai and Chakraborty 2016) but this is not explicitly accounted for in Eq. 2. Thus, it is often necessary to solve a modelled transport equation of $\sigma_{v}^{2}$. The exact transport equation of $\sigma_{v}^{2}$ contains unclosed terms, which need closures. The LES modelling of Flame Surface Density (FSD) (Chakraborty and Cant 2007, 2009; Boger et al. 1998; Chakraborty and Klein 2008; Katragadda et al. 2012a, b; Reddy and Abraham 2012; Klein et al. 2016; Lai et al. 2017; Klein and Chakraborty 2018), SGS scalar flux (Tullis and Cant 2003; Lecocq et al. 2010; Gao et al. 2015a, b; Klein et al. 2016, 2018) and scalar dissipation rate (SDR) (Dunstan et al. 2013; Ma et al. 2014; Gao et al. 2014a, b, 2015a, b, 2016; Langella et al. 2015; Gao and Chakraborty 2016) have been addressed based on a priori analysis of Direct Numerical Simulation (DNS) data. However, relatively limited effort (Nilsson et al. 2019) has been directed to a priori analysis of the closure of SGS scalar variance in premixed turbulent flames. Recently, Nilsson et al. (2019) analysed the transport characteristics of sub-grid scale variance of reaction progress variable for high Karlovitz number flames, which mostly nominally belonged to the broken reaction zones regime (Peters 2000). The statistical behaviours of SGS variance and its transport characteristics for moderate values of Karlovitz numbers traditionally associated with the flamelet regime of combustion, and for a range of different Damköhler numbers are yet to be analysed in detail based on a priori DNS analysis in the existing literature. In this respect, it is worth noting that the Damköhler number range analysed by Nilsson et al. (2019) was much smaller than the values for which flamelet based SGS variance closure is often needed. Furthermore, some of the thermal expansion related effects (e.g. counter-gradient transport) are likely to be weak for high Karlovitz numbers, which might not be the case for moderate values of Karlovitz number. Thus, it is necessary to analyse the statistical behaviours of SGS variance and its transport characteristics for Damköhler and Karlovitz number values traditionally associated with the flamelet regime of combustion. This paper addresses this void in the existing literature by analyzing the statistical behaviours of the SGS variance transport using a three-dimensional DNS database of statistically planar turbulent premixed flames representing the flamelet regime of combustion for a range of different turbulence intensities. The DNS data has been explicitly filtered to extract the unclosed terms of the SGS variance transport equation so that the performances of the existing models can be assessed and based on this exercise, either the most appropriate model expression can be identified or new models can be proposed. In this respect, the main objectives of this analysis are:

• 1. To analyse the statistical behaviours of the different terms of the SGS variance transport equation for a range of different LES filter widths for Damköhler and Karlovitz number values representative of the flamelet regime of combustion.

• 2. To identify the most suitable model expressions for the unclosed terms in the SGS variance transport equation for the flamelet combustion regime based on a detailed a priori analysis.

The mathematical background and numerical implementation pertaining to this work are presented in the next section. This will be followed by the presentation of the results and their discussion. The main findings will be summarized and conclusions will be drawn in the final section of this paper.

2 Mathematical Background and Numerical Implementation

The statistics of the transport of SGS variance of reaction progress variable have been analysed here for a range of different turbulence intensities and filter widths for the purpose of the identification of the appropriate model expressions for the unclosed terms, which are expected to perform satisfactorily for a range of different conditions. This necessitates a computationally extensive parametric analysis and therefore a single-step Arrhenius type chemical mechanism has been considered for the purpose of computational economy. In this context, the reaction progress variable $c$ is defined based on the reactant mass fraction $Y_{R}$ as: $c = \left( {Y_{R0} - Y_{R} } \right)/\left( {Y_{R0} - Y_{R\infty } } \right)$ where the subscripts 0 and $\infty$ refer to the values in the unburned reactants and fully burned products, respectively. According to this definition, $c$ increases monotonically from 0.0 in the unburned gas to 1.0 in fully burned products. The transport equations of instantaneous reaction progress variable $c$ and Favre-filtered reaction progress variable $\tilde{c}$ are given by:

$$ \begin{aligned} & \rho \frac{\partial c}{\partial t} + \rho u_{j} \frac{\partial c}{{\partial x_{j} }} = \dot{w} + \frac{\partial }{{\partial x_{j} }}\left( {\rho D\frac{\partial c}{{\partial x_{j} }}} \right) \\ & {\text{and}} \\ & \bar{\rho }\frac{{\partial \tilde{c}}}{\partial t} + \bar{\rho }\tilde{u}_{j} \frac{{\partial \tilde{c}}}{{\partial x_{j} }} = \bar{\dot{w}} + \frac{\partial }{{\partial x_{j} }}\left( {\rho D\frac{{\partial \tilde{c}}}{{\partial x_{j} }}} \right) - \frac{{\partial \left[ {\overline{{\rho u_{j} c}} - \bar{\rho }\tilde{u}_{j} \tilde{c}} \right]}}{{\partial x_{j} }} \\ \end{aligned} $$

(3)

where $u_{j} , \dot{w}$ and $D$ are the jth component of velocity, reaction rate of progress variable and reaction progress variable diffusivity, respectively. It is worth noting that Eq. 3 is valid only for premixed combustion (i.e. for constant values of equivalence ratio) and for partially-premixed/stratified mixture combustion, there will be additional terms due to mixture inhomogeneity on the right-hand side of Eq. 3. As $c$ is defined based on a species mass fraction in this analysis, Eq. 3 remains valid for non-unity Lewis number and non-adiabatic conditions. Henceforth the present analysis will focus on the closure of SGS variance $\sigma_{v}^{2}$ and the unclosed terms of its transport equation for unity Lewis number globally adiabatic flames for purely premixed combustion.

Multiplying the transport equation of $c$ by $2c$ and subsequent filtering yields a transport equation for $\widetilde{{c^{2} }}$. Similarly, multiplying the transport equation of $\tilde{c}$ by $2\tilde{c}$ yields a transport equation for $\tilde{c} ^{2}$. Subtracting the transport equation of $\tilde{c} ^{2}$ from the transport equation of $\widetilde{{c^{2} }}$ yields the transport equation of SGS variance $\sigma_{v}^{2}$ as:

$$ \begin{aligned} & \frac{{\partial \left( {\bar{\rho }\sigma_{v}^{2} } \right)}}{\partial t} + \frac{{\partial \left( {\bar{\rho }\tilde{u}_{j} \sigma_{v}^{2} } \right)}}{{\partial x_{j} }} = \underbrace {{ - \frac{\partial }{{\partial x_{j} }}\left[ {\overline{{\rho u_{j} c^{2} }} - 2\left[ {\overline{{\rho u_{j} c}} - \bar{\rho }\tilde{u}_{j} \tilde{c}} \right]\tilde{c} - \bar{\rho }\tilde{u}_{j} \widetilde{{c^{2} }}} \right]}}_{{T_{1} }} \\ & \underbrace {{ - 2\left[ {\overline{{\rho u_{j} c}} - \bar{\rho }\tilde{u}_{j} \tilde{c}} \right]\frac{{\partial \tilde{c}}}{{\partial x_{j} }}}}_{{T_{2} }} + \underbrace {{2\left[ {\overline{{\dot{w}c}} - \overline{{\dot{w}}} \tilde{c}} \right]}}_{{T_{3} }} + \underbrace {{\frac{\partial }{{\partial x_{j} }}\left( {\rho D\frac{{\partial \sigma_{v}^{2} }}{{\partial x_{j} }}} \right)}}_{{T_{4} }}\underbrace {{ - 2\bar{\rho }\widetilde{{\epsilon_{c} }}}}_{{D_{v} }} \\ \end{aligned} $$

(4)

where $\widetilde{{\epsilon_{c} }} = \overline{\rho D\nabla c \cdot \nabla c} /\bar{\rho } - \widetilde{D}\nabla \tilde{c} \cdot \nabla \tilde{c}$ is the sub-grid scalar dissipation rate (SDR), whereas $\widetilde{{N_{c} }} = \overline{\rho D\nabla c \cdot \nabla c} /\bar{\rho }$ will henceforth be referred to as the Favre-filtered SDR in this paper. The term $T_{1}$ represents turbulent transport of SGS variance, whereas the term $T_{2}$ represents generation/destruction of SGS variance due to resolved-scale scalar gradients. The term $T_{3}$ represents generation/destruction of SGS variance due to chemical reaction rate, whereas $T_{4}$ depicts molecular diffusion of $\sigma_{v}^{2}$. The last term $D_{v} = - 2\bar{\rho }\widetilde{{\epsilon_{c} }}$ represents molecular dissipation of SGS variance due to SDR. In LES, the SGS flux of reaction progress variable $\left[ {\overline{{\rho u_{j} c}} - \bar{\rho }\tilde{u}_{j} \tilde{c}} \right]$ needs closure to solve the transport equation of $\tilde{c}$, and thus the term $T_{2}$ can be considered closed. Interested readers are referred to Refs. (Gao et al. 2015a, b; Klein et al. 2016, 2018) for the discussion on the closures of SGS scalar flux in LES of turbulent premixed flames. The molecular diffusion term $T_{4}$ is also a closed term and thus the closure of Eq. 4 depends on the modelling of $T_{1} ,T_{3}$ and $D_{v}$. The terms on the right hand side of Eq. 4 have been extracted to analyse their statistical behaviours by explicitly filtering the DNS data in this paper. Furthermore, the terms $T_{1} ,T_{3}$ and $D_{v}$ extracted from DNS data have been utilised for the purpose of assessment of their respective models.

A three-dimensional compressible DNS database of freely propagating statistically planar flames with different turbulence intensities has been considered for this analysis. A well-known DNS code SENGA (Jenkins and Cant 1999) has been used for generating the database. In SENGA, the conservation equations of mass, momentum, energy and species have been solved in non-dimensional form and spatial discretization and explicit time advancements have been achieved by using high-order finite-difference (10th order at the internal grid points but the order of accuracy drops gradually to 2nd order one-sided scheme at non-periodic boundaries) and Runge–Kutta (3rd order low-storage) schemes. The simulation domain is taken to be cube a of $\left( {45.7\delta_{th} } \right)^{3}$ where $\delta_{th} = \left( {T_{ad} - T_{0} } \right)/\hbox{max} \left| {\nabla T} \right|_{L}$ is the thermal flame thickness with $T_{0}$ and $T_{ad}$ being the unburned gas temperature and adiabatic flame temperature, respectively and the subscript ‘L’ refers to unstrained laminar premixed flame quantities. The computational domain is discretized by a Cartesisan grid of $512^{3}$ points with uniform grid spacing. This resolution ensures about 11 grid points within $\delta_{th}$ and the Kolmogorov length scale $\eta$ remains at least 1.5 times the grid spacing for the cases considered here. The domain boundaries in the direction of mean flame propagation (i.e. $x_{1}$-direction) are taken to be partially non-reflecting and the transverse directions are taken to be periodic. Standard Navier–Stokes Characteristic Boundary Conditions (NSCBC) have been utilised to specify the non-periodic boundary conditions. The turbulent velocity fluctuations have been specified by homogeneous isotropic divergence-free distributions created by a standard pseudo-spectral method (Rogallo 1981) following a model turbulent kinetic energy spectrum (Pope 2000). The reacting scalar field has been initialized by a steady-state planar unstrained laminar flame solution. All the species are considered to be perfect gases and standard values have been taken for Prandtl number (i.e. ${ \Pr } = 0.7$) and ratio of specific heats (i.e. $\gamma = 1.4$). The heat release parameter $\tau = \left( {T_{ad} - T_{0} } \right)/T_{0}$ and Zel’dovich number $\beta = T_{ac} \left( {T_{ad} - T_{0} } \right)/T_{ad}^{2}$ are taken to be 3.0 and 6.0, respectively where $T_{ac}$ is the activation temperature. These values are representative of atmospheric stoichiometric methane-air combustion when the mixture is preheated to 590 K. In all cases, the flame-turbulence interaction takes place under decaying turbulence. The initial values of normalised root-mean square velocity $u^{{\prime }} /S_{L}$, integral length scale of turbulence to flame thickness ratio $l/\delta_{th}$, Damköhler number ${\text{Da}} = lS_{L} /u^{{\prime }} \delta_{th}$ and Karlovitz number ${\text{Ka}} = \left( {u^{{\prime }} /S_{L} } \right)^{3/2} \left( {l/\delta_{th} } \right)^{ - 1/2}$ are listed in Table 1 where $S_{L}$ is the unstrained laminar burning velocity. All these cases nominally represent the thin reaction zones regime combustion according to the regime diagram by Peters (Peters 2000). All the simulations have been conducted for $t_{sim} = \hbox{max} (2.0l/u^{{\prime }} ,\delta_{th} /S_{L} )$. By this time turbulence intensity $u^{{\prime }} /S_{L}$ decayed by 36%, 38%, 41% and 53% of its initial value for cases A–D respectively. This simulation time remains comparable to several previous studies (Boger et al. 1998; Reddy and Abraham 2012; Tullis and Cant 2003; Lecocq et al. 2010; Han and Huh 2008; Veynante et al. 1997), which contributed significantly to the simulation and modelling of turbulent premixed combustion in the past. The turbulent kinetic energy and dissipation rate were not varying rapidly with time when the statistics were extracted.

Table 1 Initial values of turbulence parameters for the cases considered here

For the purpose of this analysis, the DNS data has been explicitly LES filtered using a Gaussian filter kernel $G\left( \varvec{r} \right)$ so that the LES filtered values of a general quantity $Q$ can be calculated as follows (Chakraborty and Cant 2007, 2009; Dunstan et al. 2013; Gao et al. 2014; Ma et al. 2014; Boger et al. 1998; Chakraborty and Klein 2008; Katragadda et al. 2012a, b; Reddy and Abraham 2012; Klein et al. 2016; Lai et al. 2017; Klein and Chakraborty 2018; Tullis and Cant 2003; Lecocq et al. 2010; Gao et al. 2014, 2015a, b, c, d, 2016; Klein et al. 2016, 2018; Langella et al. 2015; Gao and Chakraborty 2016):

$$ \overline{{Q\left( \varvec{x} \right)}} = \int Q\left( {\varvec{x} - \varvec{r}} \right)G\left( \varvec{r} \right)\varvec{dr},\varvec{ }\quad G\left( \varvec{r} \right) = \left( {6/\pi\Delta ^{2} } \right)^{3/2} \exp \left( { - 6 \varvec{r} \cdot \varvec{r} /\Delta ^{2} } \right) $$

(5)

In this paper, results will be presented from $\Delta \approx 0.8 \delta_{th}$ where the flame is partially resolved, up to $\Delta \approx 4.6\delta_{th}$ where the flame becomes fully unresolved and $\Delta$ becomes equal to the integral length scale $l$.

3 Results and Discussion

The instantaneous views of $c = 0.1,0.3, 0.5,0.7$ and 0.9 isosurfaces for cases A–D are shown in Fig. 1, which demonstrates that the extent of flame wrinkling increases with increasing $u^{\prime}/S_{L}$. It can be discerned from Fig. 1 that $c$-isosurfaces are mostly parallel to each other in case A, whereas the $c = 0.1$ isosurface is more distorted than the $c = 0.9$ isosurface in cases C and D. The value of ${\text{Ka}}\sim\delta_{th}^{2} /\eta^{2}$ increases from case A to D (see Table 1) and thus the length scale separation between $\delta_{th}$ and $\eta$ increases from case A to D. Thus, energetic turbulent eddies are more likely to perturb the preheat zone for higher values of Karlovitz number.

Fig. 1

Views of $c = 0.1,0.3, 0.5,0.7$ and 0.9 isosurfaces for cases A–D at $t_{sim}$

The variations of mean values of $\sigma_{v}^{2}$ conditional upon bins of $\tilde{c}$ for cases A–D for different normalised filter widths $\Delta /\delta_{th}$ are shown in Fig. 2. It can be seen from Fig. 2 that $\sigma_{v}^{2}$ increases with increasing $\Delta$ before approaching an asymptotic value which is considerably smaller than the maximum possible value $\tilde{c}\left( {1 - \tilde{c}} \right)$. For small filter widths both $\widetilde{{c^{2} }}$ and $\tilde{c}^{2}$ approach $c^{2}$ ($lim_{\Delta \to 0} \widetilde{{c^{2} }} = c^{2}$ and $lim_{\Delta \to 0} \tilde{c}^{2} = c^{2}$) and thus $\sigma_{v}^{2} = \widetilde{{c^{2} }} - \tilde{c}^{2}$ vanishes when the LES filter width approaches the DNS grid size. For a presumed bi-modal sub-grid PDF with impulses at 0.0 and 1.0 one obtains: $\sigma_{v}^{2} = \tilde{c}\left( {1 - \tilde{c}} \right) + O\left( {{\text{Da}}_{\Delta }^{ - 1} } \right)$ (Bray et al. 1985) where ${\text{Da}}_{\Delta } = \Delta S_{\text{L}} /u_{\Delta }^{'} \delta_{th}$ is the sub-grid Damköhler number with $u_{\Delta }^{'} = \sqrt {\left( {\widetilde{{u_{i} u_{i} }} - \tilde{u}_{i} \tilde{u}_{i} } \right)/3} = \sqrt {2k_{sgs} /3}$ and $k_{sgs} = \left( {\widetilde{{u_{i} u_{i} }} - \tilde{u}_{i} \tilde{u}_{i} } \right)/2$ being the SGS velocity fluctuation and the SGS turbulent kinetic energy, respectively. The term $O\left( {{\text{Da}}_{\Delta }^{ - 1} } \right)$ can be ignored only for ${\text{Da}}_{\Delta } \gg 1$ and under that condition the maximum possible value of the SGS variance (i.e. $\sigma_{v}^{2} = \tilde{c}\left( {1 - \tilde{c}} \right)$) is obtained. The variation of mean values of ${\text{Da}}_{\Delta }$ conditional upon bins of $\tilde{c}$ for cases A–D for different normalised filter widths $\Delta /\delta_{th}$ are shown in Fig. 3, which shows that ${\text{Da}}_{\Delta }$ increases with increasing $\Delta$ but ${\text{Da}}_{\Delta }$ does not assume large (i.e. ${\text{Da}}_{\Delta }$ ≫ 1) values even for the largest value of $\Delta$ considered here. Thus, the $O\left( {{\text{Da}}_{\Delta }^{ - 1} } \right)$ contribution is expected to play a significant role in the present flames for the range of $\Delta$ analysed here.

Fig. 2

Variations of mean values of $\sigma_{v}^{2}$ conditional upon bins of $\tilde{c}$ along with the predictions of Eq. 2 and $\tilde{c}\left( {1 - \tilde{c}} \right)$ for different $\Delta /\delta_{th}$ in cases A–D

Fig. 3

Variations of mean values of $Da_{\Delta }$ conditional upon bins of $\tilde{c}$ for different $\Delta /\delta_{th}$ in cases A–D

The extent of the deviation of $\sigma_{v}^{2}$ from $\tilde{c}\left( {1 - \tilde{c}} \right)$ provides a measure of the departure from the bi-modal distribution with impulses at $c = 0.0$ and 1.0, and the results in Fig. 2 suggest that the sub-grid PDF cannot be assumed to be bi-modal for all cases considered here irrespective of the filter width, which is consistent with previous findings (Chakraborty and Cant 2007, 2009; Dunstan et al. 2013; Gao et al. 2014; Ma et al. 2014).

The predictions of Eq. 2 are also shown in Fig. 2, which indicates that this model overpredicts $\sigma_{v}^{2}$ significantly and predicts unphysical $\sigma_{v}^{2} > \tilde{c}\left( {1 - \tilde{c}} \right)$ for $C = 0.5$ for large filter widths. The quantitative agreement between the algebraic model prediction (i.e. Eq. 2) and $\sigma_{v}^{2}$ extracted from DNS data improves only slightly from case A–D. This suggests that the algebraic closure given by Eq. 2 might not be sufficient for the prediction of $\sigma_{v}^{2}$ under all conditions. Thus, it might be worthwhile to solve a modelled transport equation of SGS variance $\sigma_{v}^{2}$.

The variations of the normalised mean values of the terms on the right hand side of the SGS variance transport equation (i.e. $\{ T_{1} ,T_{2} ,T_{3} ,T_{4}$ and $D_{v} \} \times \delta_{th} /\rho_{0} S_{L}$ with $\rho_{0}$ being the unburned gas density) conditional upon bins of $\tilde{c}$ for $\Delta /\delta_{th} = 0.8$ and 4.6 are shown in Fig. 4 for cases A–D. It can be seen from Fig. 4 that $D_{v}$ acts as a leading order sink term for all cases irrespective of the filter width. By contrast, the reaction rate contribution $T_{3}$ acts as a leading order source term apart from the negative contribution towards the burned gas side of the flame brush for all cases for all filter widths considered here. The mean contribution of the resolved scalar gradient term $T_{2}$ remains negative for cases A and B for all filter widths which is indicative of counter-gradient transport (i.e. $\left[ {\overline{{\rho u_{j} c}} - \bar{\rho }\tilde{u}_{j} \tilde{c}} \right]\partial \tilde{c}/\partial x_{j} > 0$), whereas this mean contribution is positive in cases C and D for all filter widths due to gradient transport (i.e. $\left[ {\overline{{\rho u_{j} c}} - \bar{\rho }\tilde{u}_{j} \tilde{c}} \right]\partial \tilde{c}/\partial x_{j} < 0$) behaviour in these cases. This is consistent with the fact that counter-gradient transport is likely when $u^{\prime} < \tau S_{L}$ (e.g. cases A and B when the statistics are extracted), whereas the gradient transport is expected to be predominant for $u^{\prime} > \tau S_{L}$ (e.g. in cases C and D when statistics are extracted) (Lecocq et al. 2010; Gao et al. 2015a, b; Klein et al. 2016, 2018; Veynante et al. 1997). The magnitude of $T_{2}$ remains small in comparison to the magnitudes of $T_{3}$ and $D_{v}$. The mean contribution of the molecular diffusion term $T_{4}$ remains positive on both unburned and burned gas sides of the flame brush with a negative dip in the middle for all cases irrespective of $\Delta$, whereas the mean behaviour of the turbulent transport term $T_{1}$ assumes both positive and negative values within the flame brush and its qualitative behaviour changes from one case to another. The magnitudes of mean contributions of $T_{1}$ $(T_{1}$ and $T_{4} )$ remain smaller than the magnitudes of $T_{3}$ and $D_{v}$ for all cases for $\Delta \le \delta_{th}$ ($\Delta \gg \delta_{th}$). It can further be seen from Fig. 4 that $T_{3}$ and $D_{v}$ remain in approximate equilibrium for $\Delta \gg \delta_{th}$ but this equilibrium is not maintained for small filter sizes. The aforementioned behaviour can be explained in the following manner using scaling relations:

Fig. 4

Variations of mean values of $\{ T_{1} ,T_{2} ,T_{3} ,T_{4}$ and $D_{v} \} \times \delta_{th} /\rho_{0} S_{L}$ conditional upon bins of $\tilde{c}$ for $\Delta /\delta_{th} = 0.8$ and 4.6 in cases A–D

$$ \begin{aligned} & T_{1} \sim\frac{{\rho_{0} u_{\Delta }^{'} }}{\Delta }\sim\left( {\frac{{\rho_{0} S_{L} }}{{\delta_{th} }}} \right){\text{Da}}_{\Delta }^{ - 1} \;\;{\text{for}}\;\;u^{\prime} \gg \tau S_{L} \;\;{\text{or}}\;\;T_{1} \sim\frac{{\rho_{0} S_{L} }}{\Delta }\sim\left( {\frac{{\rho_{0} S_{L} }}{{\delta_{th} }}} \right){\text{Re}}_{\Delta }^{ - 0.5} {\text{Da}}_{\Delta }^{ - 0.5} \;\;{\text{for}}\;\;u^{\prime} \ll \tau S_{L} \\ & T_{2} \sim\frac{{\rho_{0} u_{\Delta }^{'} }}{\Delta }\sim\left( {\frac{{\rho_{0} S_{L} }}{{\delta_{th} }}} \right){\text{Da}}_{\Delta }^{ - 1} \;\;{\text{for}}\;\;u^{\prime} \gg \tau S_{L} \;\;{\text{or}}\;\;T_{2} \sim\frac{{\rho_{0} S_{L} }}{\Delta }\sim\left( {\frac{{\rho_{0} S_{L} }}{{\delta_{th} }}} \right){\text{Re}}_{\Delta }^{ - 0.5} {\text{Da}}_{\Delta }^{ - 0.5} \;\;{\text{for}}\;\;u^{\prime} \ll \tau S_{L} \\ & T_{3} \sim\frac{{\rho_{0} S_{L} }}{{\delta_{th} }} \\ & T_{4} \sim\frac{{\rho_{0} S_{L} \delta_{th} }}{{\Delta ^{2} }}\sim\left( {\frac{{\rho_{0} S_{L} }}{{\delta_{th} }}} \right){\text{Re}}_{\Delta }^{ - 1} {\text{Da}}_{\Delta }^{ - 1} \\ & D_{v} = \left( {\frac{{\rho_{0} S_{L} }}{{\delta_{th} }}} \right) \\ \end{aligned} $$

(6)

Here, the Favre-filtered velocity $\tilde{u}_{j}$ is taken to scale with the mean velocity $U_{mean}$ ahead of the flame, whereas sub-grid fluctuations of $\dot{w}$ and $\nabla c$ are taken to scale with $\rho_{0} S_{L} /\delta_{th}$ and $1/\delta_{th}$ respectively (Gao et al. 2014; Swaminathan and Bray 2005). Moreover, the length scale associated with the gradients of resolved/filtered quantities are scaled with respect to $\Delta$. Furthermore, the SGS scalar flux $\left[ {\overline{{\rho u_{j} c}} - \bar{\rho }\tilde{u}_{j} \tilde{c}} \right]$ is scaled with respect to $\rho_{0} S_{L}$ for counter-gradient transport where $u^{\prime} \ll \tau S_{L}$, whereas the SGS scalar flux is taken to scale with $\rho_{0} u_{\Delta }^{ '}$ in the case of gradient transport where $u^{\prime} \gg \tau S_{L}$ (Gao et al. 2014; Swaminathan and Bray 2005). The SGS flux of variance $\left[ {\overline{{\rho u_{j} c^{2} }} - 2\left[ {\overline{{\rho u_{j} c}} - \bar{\rho }\tilde{u}_{j} \tilde{c}} \right]\tilde{c} - \bar{\rho }\tilde{u}_{j} \widetilde{{c^{2} }}} \right]$ is also scaled in the same manner as that of the SGS scalar flux $\left[ {\overline{{\rho u_{j} c}} - \bar{\rho }\tilde{u}_{j} \tilde{c}} \right]$. It is worth noting that in Eq. 6, $Re_{\Delta } = \rho_{0} u_{\Delta }^{'}\Delta /\mu_{0}$ is the sub-grid Reynolds number with $\mu_{0}$ being the unburned gas viscosity. The term ${\text{Re}}_{\Delta } {\text{Da}}_{\Delta }$ scales as ${\text{Re}}_{\Delta } {\text{Da}}_{\Delta } \sim\left( {\Delta /\delta_{th} } \right)^{2}$ and it can be seen from Fig. 3 that ${\text{Da}}_{\Delta }$ increases with increasing $\Delta$, which is also consistent with previous analyses (Dunstan et al. 2013; Gao et al. 2014; Ma et al. 2014). According to inertial-range scaling ${\text{Da}}_{\Delta }$ scales as: ${\text{Da}}_{\Delta } \sim\left( {u_{\Delta }^{'} /S_{L} } \right)^{2} (S_{L}^{3} /\delta_{th} \epsilon_{\Delta } )$ and as $(S_{L}^{3} /\delta_{th} \epsilon_{\Delta } )$ is a constant in the inertial range, ${\text{Da}}_{\Delta }$ is expected to increase with increasing $\Delta$. Although an inertial range is not obtained for this database, the findings from Fig. 3 are also qualitatively consistent with the inertial scaling discussed above. The scaling relations in Eq. 6 indicate that $T_{3}$ and $D_{v}$ are expected to play leading roles for all cases, whereas relative magnitudes of $T_{1} ,T_{2} , T_{4}$ and $C_{v}$ are expected to be increasingly small in comparison to those of $T_{3}$ and $D_{v}$ with increasing $\Delta$. This is found to be qualitatively consistent with the observations made from Fig. 4. As $T_{2}$ and $T_{4}$ are closed terms, the closure of $\sigma_{v}^{2}$ depends on the modelling of $T_{1} , T_{3}$ and $D_{v}$. These findings are also consistent with the findings in a recent analysis by Nilsson et al. (2019), where a scaling analysis with respect to the Karlovitz number came to the same conclusions. That analysis used a detailed chemistry database to extract the different terms of the transport equation for different filter widths and they were shown to scale in the same way as found in this analysis.

It can be seen from Eq. 4 that the closure of $T_{1}$ depends on the appropriate closure of the SGS flux of variance $F_{j} = \left[ {\overline{{\rho u_{j} c^{2} }} - 2\left[ {\overline{{\rho u_{j} c}} - \bar{\rho }\tilde{u}_{j} \tilde{c}} \right]\tilde{c} - \bar{\rho }\tilde{u}_{j} \widetilde{{c^{2} }}} \right]$. The SGS flux of variance is often modelled using the gradient hypothesis (Langella and Swaminathan 2016; Langella et al. 2016):

$$ F_{j} = \left[ {\overline{{\rho u_{j} c^{2} }} - 2\left[ {\overline{{\rho u_{j} c}} - \bar{\rho }\tilde{u}_{j} \tilde{c}} \right]\tilde{c} - \bar{\rho }\tilde{u}_{j} \widetilde{{c^{2} }}} \right] = - \frac{{\bar{\rho }}}{{{\text{Sc}}_{t} }}\left( {C_{s} \Delta } \right)^{2} \sqrt {2\tilde{S}_{ij} \tilde{S}_{ij} } \frac{{\partial \sigma_{v}^{2} }}{{\partial x_{j} }} $$

(7)

where $C_{s} = 0.18$ is the Smagorinsky constant, $\tilde{S}_{ij} = 0.5\left( {\partial \tilde{u}_{i} /\partial x_{j} + \partial \tilde{u}_{j} /\partial x_{i} } \right)$ is the resolved strain rate and ${\text{Sc}}_{t}$ is the turbulent Schmidt number. The mean values of the components $F_{1} /\rho_{0} S_{L}$ and $F_{2} /\rho_{0} S_{L}$ conditional upon bins of $\tilde{c}$ are shown along with the prediction of Eq. 7 (for ${\text{Sc}}_{t} = 1.0$) in Fig. 5 for cases A–D for $\Delta /\delta_{th} = 0.8$ and 4.6. The behaviour of $F_{3}$ is not explicitly shown because $F_{2}$ and $F_{3}$ are statistically similar in this configuration. It can be seen from Fig. 5 that Eq. 7 does not capture the qualitative behaviour of $F_{1}$ for both filter widths, and it predicts the opposite sign to the value obtained from DNS in cases A and B, which is indicative of the counter-gradient behaviour of the SGS variance flux. The predictions of $F_{2}$ by Eq. 7 in cases A and B are in better agreement with DNS data than for $F_{1}$ because the effects of heat release and flow normal acceleration are stronger in the direction of mean flame propagation than in the transverse direction. However, Eq. 7 predicts both qualitative and quantitative behaviours of $F_{1}$ and $F_{2}$ satisfactorily for both filter widths in cases C and D, which indicates the predominance of gradient transport in these cases. These findings are consistent with the expectation that the effects of flame normal acceleration may overcome the influences of turbulent fluctuation to yield counter-gradient transport for small values of $u^{\prime}/S_{L}$, whereas gradient transport is observed for high values of $u^{\prime}/S_{L}$ where turbulence effects overwhelm the influence of flame normal acceleration (Lecocq et al. 2010; Gao et al. 2015a, b; Klein et al. 2016; Veynante et al. 1997). Thus, it will be desirable to have a model, which is capable of predicting both gradient and counter-gradient behaviours of $F_{j}$. The satisfactory performance of Clark’s model (1979) for SGS scalar flux (Gao et al. 2015a, b; Klein et al. 2016) prompts an alternative model as:

Fig. 5

Variations of mean values of $F_{1} /\rho_{0} S_{L}$ and $F_{2} /\rho_{0} S_{L}$ conditional upon bins of $\tilde{c}$ for $\Delta /\delta_{th} = 0.8$ and 4.6 in cases A–D along with the predictions of Eq. 7 (for $Sc_{t} = 1.0$), 8, 9

$$ F_{j} = \left[ {\overline{{\rho u_{j} c^{2} }} - 2\left[ {\overline{{\rho u_{j} c}} - \bar{\rho }\tilde{u}_{j} \tilde{c}} \right]\tilde{c} - \bar{\rho }\tilde{u}_{j} \widetilde{{c^{2} }}} \right] = \bar{\rho }\frac{{\Delta^{2} }}{12}\frac{{\partial \tilde{u}_{j} }}{{\partial x_{k} }}\frac{{\partial \sigma_{v}^{2} }}{{\partial x_{k} }} $$

(8)

The predictions of the model given by Eq. 8 are also shown in Fig. 5, which shows that this model is indeed capable of capturing both gradient and counter-gradient behaviours of SGS flux of variance. The quantitative agreement of Eq. 8 with DNS data is better than the model given by Eq. 7 but the performance of Eq. 8 deteriorates with increasing $\Delta$. In order to avoid this problem, a model expression, which was proposed by Chakraborty and Swaminathan (Chakraborty and Swaminathan 2011) for the Reynolds flux of variance in the context of RANS, has been extended for LES in the following manner:

$$ F_{j} = \left[ {\overline{{\rho u_{j} c}} - \bar{\rho }\tilde{u}_{j} \tilde{c}} \right]\left[ {1 - 2\tilde{c}\left\{ {\frac{{\sigma_{v}^{2} }}{{\tilde{c}\left( {1 - \tilde{c}} \right)}}} \right\}^{0.3} } \right]\frac{{2\sigma_{v}^{2} }}{{\sigma_{v}^{2} + \tilde{c}\left( {1 - \tilde{c}} \right)}} $$

(9)

It can be seen from Fig. 5 that Eq. 9 satisfactorily predicts both qualitative and quantitative predictions of $F_{1}$ and $F_{2}$ for both (and all) filter widths in all cases considered here. It is worth noting that the satisfactory performance of Eq. 9 depends on appropriate modelling of the SGS scalar flux $\left[ {\overline{{\rho u_{j} c}} - \bar{\rho }\tilde{u}_{j} \tilde{c}} \right]$, which has been addressed by the authors elsewhere (Gao et al. 2015a, b; Klein et al. 2016, 2018) and thus is not repeated here. Interested readers are referred to Refs. (Gao et al. 2015a, b; Klein et al. 2016, 2018) for suitable models for the SGS scalar flux $\left[ {\overline{{\rho u_{j} c}} - \bar{\rho }\tilde{u}_{j} \tilde{c}} \right]$ for turbulent premixed flames.

The modelling of $T_{3} = 2\left[ {\overline{{\dot{w}c}} - \bar{\dot{w}}\tilde{c}} \right]$ depends on the closures of $\bar{\dot{w}}$ and $\overline{{\dot{w}c}}$. In the context of flamelet closure for unity Lewis number flames, one gets:

$$ \bar{\dot{w}} = \mathop \int \limits_{0}^{1} \left[ {\dot{w}} \right]_{L} P_{\Delta } \left( c \right)dc\;\;{\text{and}}\;\;\overline{{\dot{w}c}} = \mathop \int \limits_{0}^{1} \left[ {\dot{w}c} \right]_{L} P_{\Delta } \left( c \right)dc $$

(10)

where $\left[ q \right]_{L}$ is the value of a general variable $q$ obtained from a flamelet table and $P_{\Delta } \left( c \right)$ is the presumed sub-grid PDF of $c$, which is often parameterized with the help of $\tilde{c}$ and $\sigma_{v}^{2}$ (Langella and Swaminathan 2016; Langella et al. 2016). Alternatively, it is possible to express $\overline{{\dot{w}c}}$ as: $\overline{{\dot{w}c}} = \bar{\dot{w}}c_{m}$ where $c_{m} = \mathop \int \limits_{0}^{1} \left[ {\dot{w}c} \right]_{L} f\left( c \right)dc/\mathop \int \limits_{0}^{1} \left[ {\dot{w}} \right]_{L} f\left( c \right)dc$ (Bray 1980), where $f\left( c \right)$ is the burning mode PDF, which can be taken as any appropriate continuous function according to Bray (Bray 1980). This thermo-chemical parameter $c_{m}$ varies between 0.7 and 0.9 and for the present thermo-chemistry $c_{m}$ has been found to be 0.84. Accordingly, $T_{3}$ can be expressed as: $T_{3} = 2\bar{\dot{w}}\left[ {c_{m} - \tilde{c}} \right]$. The variations of mean values of $T_{3}$ conditional upon bins of $\tilde{c}$ are shown in Fig. 6 for cases A–D for $\Delta /\delta_{th} = 0.8$ and 4.6. Figure 6 shows that $T_{3} = 2\bar{\dot{w}}\left[ {c_{m} - \tilde{c}} \right]$ with $\bar{\dot{w}}$ extracted from DNS data satisfactorily captures both qualitative and quantitative behaviours of $T_{3}$ for all filter widths for all cases considered here. The filtered reaction rate $\bar{\dot{w}}$ can be expressed using the Favre-filtered SDR $\tilde{N}_{c}$ as (Gao et al. 2014; Ma et al. 2014):

$$ \overline{{\dot{w}}} = \frac{{2\bar{\rho }\widetilde{N}_{c} }}{{2c_{m} - 1}}\left[ {1 - \exp \left( { - \varPhi \frac{\Delta }{{\delta_{th} }}} \right)} \right] + f_{1} \left( {\bar{\rho },\tilde{c}} \right)\exp \left( { - \varPhi \frac{\Delta }{{\delta_{th} }}} \right) $$

(11)

where $f_{1} \left( {\bar{\rho },\tilde{c}} \right)$ is a function such that $\dot{w} = f_{1} \left( {\rho ,c} \right)$ and $\phi = 0.56\delta_{th} S_{L} /\alpha_{T0}$ is a model parameter with $\alpha_{T0}$ being the thermal diffusivity in the unburned gas. It has been shown elsewhere (Gao et al. 2014) that Eq. 11 satisfactorily predicts $\bar{\dot{w}}$ for $\Delta > \delta_{th}$ and the quantitative prediction improves with increasing filter width. Interested readers are referred to Refs. (Dunstan et al. 2013; Gao et al. 2014; Ma et al. 2014) for physical explanations behind this behaviour. Based on Eq. 11, it is possible to express $T_{3}$ as:

Fig. 6

Variations of mean values of $T_{3} \times \delta_{th} /\rho_{0} S_{L}$, $2\bar{\dot{w}}\left[ {c_{m} - \tilde{c}} \right] \times \delta_{th} /\rho_{0} S_{L}$ and the predictions of Eq. 12 (for $\tilde{N}_{c}$ extracted from DNS data and according to the prediction of Eq. 13) conditional upon bins of $\tilde{c}$ for $\Delta /\delta_{th} = 0.8$ and 4.6 in cases A–D

$$ T_{3} = 2\left[ {\frac{{2\bar{\rho }\widetilde{N}_{c} }}{{2c_{m} - 1}}\left\{ {1 - \exp \left( { - \varPhi \frac{\Delta }{{\delta_{th} }}} \right)} \right\} + f_{1} \left( {\bar{\rho },\tilde{c}} \right)\exp \left( { -\Phi \frac{\Delta }{{\delta_{th} }}} \right)} \right]\left( {c_{m} - \tilde{c}} \right) $$

(12)

The predictions of Eq. 12 are also shown in Fig. 6, which shows that $T_{3}$ can be reasonably predicted by this expression and the quantitative agreement improves with increasing $\Delta$. This behaviour arises due to improved prediction of $\bar{\dot{w}}$ by Eq. 11 for large filter widths (Dunstan et al. 2013; Gao et al. 2014; Ma et al. 2014).

Gao et al. (2014) proposed an algebraic model of $\widetilde{N}_{c}$ in the following manner:

$$ \widetilde{N}_{c} = \tilde{D}\nabla \tilde{c} \cdot \nabla \tilde{c} + \left( {1 - f_{b} } \right)\left[ {\frac{{2K_{c}^{*} S_{L} }}{{\delta_{th} }} + \left( {C_{3}^{*} - \tau {\text{Da}}_{\Delta } C_{4}^{*} } \right)\frac{{2u_{\Delta }^{'} }}{3\Delta }} \right]\frac{{\tilde{c}\left( {1 - \tilde{c}} \right)}}{{\beta_{c} }} $$

(13i)

where $K_{c}^{ *} = (\delta_{th} /S_{L} )\mathop \int \limits_{0}^{1} \left[ {\rho N_{c} \nabla \cdot \vec{u}} \right]_{L} f\left( c \right)dc/\mathop \int \limits_{0}^{1} \left[ {\rho N_{c} } \right]_{L} f\left( c \right)dc$ is a thermo-chemical parameter, which is found to be $K_{c}^{*} = 0.77\tau$ for the present thermo-chemistry and $f_{b} = \exp \left[ { - 0.7\left( {\Delta /\delta_{th} } \right)^{1.7} } \right]$ is a bridging function, which ensures that $\widetilde{N}_{c}$ approaches $D\nabla c \cdot \nabla c$, when the flame is fully resolved (i.e. $lim_{\Delta \to 0} \widetilde{N}_{c} = D\nabla c \cdot \nabla c$). The model parameters $C_{3}^{ *} , C_{4}^{*}$ and $\beta_{c}$ are given by the following expressions (Gao et al. 2014):

$$ \begin{aligned} & C_{3}^{ *} = \frac{{2.0\sqrt {{\text{Ka}}_{\Delta } } }}{{1.0 + \sqrt {{\text{Ka}}_{\Delta } } }};\quad C_{4}^{ *} = \frac{{1.2\left( {1.0 - \tilde{c}} \right)^{0.2} }}{{Le^{2.57} \left( {1 + {\text{Ka}}_{\Delta } } \right)^{0.4} }} ; \\ & \beta_{c} = \hbox{max} \left( {\frac{2}{{2c_{m} - 1}},\left[ {1.05\frac{\tau }{\tau + 1} + 0.51} \right]^{4.6} } \right) \\ \end{aligned} $$

(13ii)

where ${\text{Ka}}_{\Delta } = \left( {u_{\Delta }^{'} /S_{L} } \right)^{3/2} \left( {\Delta /\delta_{th} } \right)^{ - 1/2}$ is the sub-grid Karlovitz number. Interested readers are referred to Refs. (Gao et al. 2014a, b; Ma et al. 2014) regarding the derivation of the model given by Eq. 13. The predictions of Eq. 13 are compared to the mean values of $\widetilde{N}_{c}$ extracted from DNS conditional upon the bins of $\tilde{c}$ in Fig. 7 for $\Delta /\delta_{th} = 0.8$ and 4.6, which shows that this model satisfactorily predicts $\widetilde{N}_{c}$ for all cases. It can be seen from Fig. 6 that using Eq. 13 in Eq. 12 provides satisfactory prediction of $T_{3}$ for all cases for the range of filter widths analysed here. Finally, as $\tilde{D}\nabla \tilde{c} \cdot \nabla \tilde{c}$ is a resolved quantity, Eq. 13 can be used to close the molecular dissipation term $D_{v} = - \bar{\rho }\tilde{\epsilon}_{c} = - \bar{\rho }\left[ {\widetilde{N}_{c} - \tilde{D}\nabla \tilde{c} \cdot \nabla \tilde{c}} \right]$.

Fig. 7

Variations of mean values of $\widetilde{N}_{c} \times \delta_{th} /S_{L}$ and the prediction of Eq. 13 conditional upon bins of $\tilde{c}$ for $\Delta /\delta_{th} = 0.8$ and 4.6 in cases A–D

4 Conclusions

The modelling of SGS variance of reaction progress variable has been considered in this paper based on a priori analysis of DNS data of freely propagating statistically planar turbulent premixed flames with different turbulence intensities. It has been found that the sub-grid PDF of reaction progress variable cannot be approximated by a presumed bi-modal distribution with impulses at $c = 0.0$ and 1.0. An algebraic expression of SGS variance has been found to significantly overpredict the corresponding quantity extracted from DNS data. The statistical behaviours of the unclosed terms of the SGS variance transport equation have been analysed based on explicitly filtered DNS data. It has been found that the reaction rate and molecular dissipation contributions act as leading order source and sink terms for all filter widths for all cases considered here, which has been explained based on detailed scaling arguments. The reaction rate and molecular dissipation contributions remain in approximate equilibrium for large filter widths. The modelling of the unclosed terms of the SGS variance transport equation has been addressed based on a priori DNS analysis and suitable models have been identified for the SGS flux of variance, reaction rate contribution and scalar dissipation rate. It is important to note that a recent analysis by Nilsson et al. (2019) reported the budgets of unclosed terms of the SGS variance equation based on a detailed chemistry DNS database of CH₄-air mixture for ${\text{Ka}} = 4 - 4100$. The behaviour of the unclosed terms of the transport equation of SGS variance for H₂O mass fraction based reaction progress variable, especially the approximate balance between $T_{3}$ and $\left( { - D_{v} } \right)$ for large filter widths (i.e. $\Delta \gg \delta_{th}$), are found to be qualitatively similar to the budget of the unclosed terms of the transport equation of the SGS variance shown in this analysis (see Fig. 4). Moreover, the scaling arguments provided here for explaining the behaviours of the unclosed terms of the SGS variance equation in Eq. 6 are also valid for the data reported by Nilsson et al. (2019). The magnitude of SDR $\widetilde{N}_{c}$ is also likely to be affected by the presence of detailed chemistry and transport. However, the SDR $\widetilde{N}_{c}$ closure used in this paper (see Eq. 13i) has previously been assessed for both simple (Gao et al. 2014a, b) and detailed chemistry (Gao et al. 2016; Nilsson et al. 2019) DNS data. It has been demonstrated by Gao et al. (Gao et al. 2016) that the leading order balance of the different unclosed terms of the transport equation of $\widetilde{N}_{c}$ for $\Delta \gg \delta_{th}$ in the context of detailed chemistry DNS remains qualitatively similar to the corresponding findings from simple chemistry DNS data (Gao et al. 2014a, b), and this also validates the assumption based on which the algebraic $\widetilde{N}_{c}$ closure (Eq. 13i) was originally derived (Gao et al. 2014a, b). It is also worth considering the fact that the leading order competition between the reaction rate contribution $T_{3}$ and dissipation contribution $\left( { - D_{v} } \right)$ originating from the approximate balance between reaction rate $\dot{w}$ and molecular diffusion rate $\nabla \cdot \left( {\rho D\nabla c} \right)$, is typical of major reactants and products in premixed combustion in the flamelet combustion regime, which is also valid in the presence of detailed chemistry (Chakraborty et al. 2018). However, this approximate equilibrium may not be maintained for intermediate species depending on their characteristic chemical timescales even within the flamelet regime of combustion. Therefore, the models discussed in this paper need to be used with care while extrapolating them for the purpose of modelling the SGS variance of intermediate species. As a result, further analyses based on detailed chemistry and transport will be necessary for further validation of the model expressions identified in this analysis. Furthermore, the modelling methodology proposed here is only valid for unity Lewis number, globally adiabatic, perfectly premixed turbulent combustion. Although the expressions for the SDR based reaction rate closure (Eq. 11) and algebraic closure of SDR (Eq. 13i) have already been assessed for non-unity Lewis number flames (Gao et al. 2014) and they have been demonstrated to provide satisfactory predictions, further work will be needed in terms of modelling of the unclosed terms of the transport equation of the SGS variance for non-unity Lewis number, non-adiabatic and partially-premixed/stratified mixture combustion. Finally, the proposed model expressions need to be implemented in actual LES for the purpose of a posteriori assessment, which will form the basis of future analyses.