Effects of Body Forces on Turbulent Kinetic Energy Transport in Premixed Flames

Abstract

The effects of buoyancy on turbulent premixed flames are expected to be strong due to the large changes in density between the unburned and fully burned gases. The present work utilises three-dimensional direct numerical simulations of statistically planar turbulent premixed flames under decaying turbulence to study the influence of buoyancy on the evolution of turbulent kinetic energy within the flame brush. Four sets of turbulence parameters have been studied, with four different values of Froude number for each set. It is found that for a given set of turbulence parameters, flame wrinkling increases with an increase in body force magnitude in the case of unstable stratification, which is also reflected in the increased values of turbulent burning velocity and flame surface area. An increase in body force magnitude in the case of stable stratification acts to reduce the extent of flame wrinkling. Turbulent kinetic energy and its dissipation rate are found to be affected by both the magnitude and direction of the body force. For low turbulence intensities considered here, turbulent kinetic energy increases from the leading edge of the flame brush before decaying eventually towards the product side of the flame brush. For high turbulence intensities, the turbulent kinetic energy is found to decay across the flame brush, and it is also found that the effect of body force on the evolution of turbulent kinetic energy is marginal in the case of high turbulence intensities. The effects of body force magnitude and direction on the statistical behaviours and closures of the various terms of the turbulent kinetic energy transport equation are analysed, and existing models have been modified to account for Froude number effects, where necessary.

1 Introduction

Flame generated turbulence, first proposed by Karlovitz et al. (1953), is a characteristic of premixed combustion and refers to the influence of the premixed flame on the turbulent flow field within the combustion domain. Bray and Libby (1976) used physical arguments to associate the flame generated turbulence to the mean velocity gradient across the flame brush and this hypothesis was confirmed experimentally by Moreau and Boutier (1977), who also provided evidence of counter-gradient transport of scalars in turbulent premixed combustion. The role of the fluctuating pressure gradient in turbulent kinetic energy transport, in addition to that of the mean pressure gradient, was indicated by Kuznetsov (1979) and Strahle (1983) and this was subsequently confirmed by Direct Numerical Simulation (DNS) studies (Rutland and Cant 1994; Zhang and Rutland 1995; Nishiki et al. 2002; Chakraborty et al. 2011a, b; Lai et al. 2017; Ahmed et al. 2019). Bray et al. (1981) demonstrated the role of buoyancy on flame generated turbulence in premixed turbulent combustion and their models in the context of Reynold Averaged Navier Stokes (RANS) simulations showed satisfactory agreement with experimental data.

In premixed flames, density changes by a significant amount between the unburned and fully burned gases and therefore, body forces are expected to affect the flame wrinkling and flame normal acceleration under turbulent conditions. It has been demonstrated by Veynante and Poinsot (1997) that buoyancy significantly influences statistical behaviours of turbulent scalar flux and flame wrinkling in premixed turbulent combustion using two-dimensional direct numerical simulations (DNS) data. Although an increased level of gradient transport and an augmentation of flame wrinkling have been reported for buoyantly unstable flows (Veynante and Poinsot 1997), its implication on the turbulent kinetic energy transport has not been analysed in detail in the existing literature. This study aims to address this gap in existing literature by using a three-dimensional DNS database involving a number of statistically planar turbulent premixed flames subjected to different values of the body force in the direction of mean flame propagation for different turbulence intensities. Thus, the main objectives of the present study are as follows:

• 1. To investigate the effects of body force on turbulent kinetic energy transport within the turbulent premixed flame brush.

• 2. To evaluate the performance of the existing closure models for the different terms of the turbulent kinetic energy transport equation in the context of RANS simulations and propose suitable modifications, wherever necessary, to incorporate the effects of body forces.

The remainder of this paper is organised as follows. The next section describes the essential mathematical background which is followed by the details of the numerical implementation. The results are subsequently presented and discussed. The final section provides a summary of the main findings and conclusions that are drawn from the present study.

2 Mathematical Background

In the present analysis, DNS of turbulent premixed combustion has been carried out in three dimensions with simplified chemistry due to the computational requirements of detailed chemistry simulations, which make them unsuitable for a parametric study involving several simulations, such as the present one. Hence, a single-step Arrhenius type irreversible chemical reaction has been used in the present study where 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}$ as $c = \left( {Y_{P} - Y_{P0} } \right)/\left( {Y_{P\infty } - Y_{P0} } \right)$. Here $0$ and $\infty$ refer to the unburned and fully burned gases, respectively and under adiabatic conditions $c$ is identical to the non-dimensional temperature $T = \left( {\hat{T} - T_{0} } \right)/\left( {T_{ad} - T_{0} } \right)$ (Poinsot and Veynante 2021), where $\hat{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 a body force, the momentum conservation equation in the i^th direction takes the following form (Veynante and Poinsot 1997):

$$ \frac{{\partial \left( {\rho u_{i} } \right)}}{\partial t} + \frac{{\partial \left( {\rho u_{i} u_{j} } \right)}}{{\partial x_{j} }} = - \frac{\partial p}{{\partial x_{i} }} + S_{i} + \frac{{\partial \tau_{ij} }}{{\partial x_{j} }} $$

(1)

where $\rho$ is the density, $u_{i}$ is the i^th component of fluid velocity, $p$ is the pressure, $\tau_{ij}$ is the component of viscous stress tensor and $S_{i}$ is the source/sink term in the i^th direction. The source term is assumed to take non-zero values only in the $x_{1}$-direction which is taken to align with the mean direction of flame propagation. Here, $S_{1}$ is expressed as: $S_{1} = \rho {\Gamma } = \rho g^{*} S_{L}^{2} /\delta_{Z}$ where ${\Gamma }$ represents the externally imposed dimensional acceleration/deceleration, $\delta_{Z} = \alpha_{T0} /S_{L}$ is the Zel’dovich flame thickness, $S_{L}$ is the unstrained laminar burning velocity, $\alpha_{T0}$ is the thermal diffusivity in the unburned gas, and $g^{*}$ stands for the inverse of Froude number-squared (i.e., $g^{*} = {\Gamma }\delta_{Z} /S_{L}^{2} = Fr^{ - 2}$). A positive (negative) value of $g^{*}$ indicates an externally imposed flow acceleration (deceleration). Alternatively, a positive value of $g^{*}$ represents an unstable stratification where the heavier cold unburned reactants reside on top of lighter hot burned products, whereas a negative value of $g^{*}$ signifies a situation of stable stratification where the heavier cold unburned reactants reside underneath the lighter hot burned products.

Using Eq. (1) the transport equation for the turbulent kinetic energy $\tilde{k} = 0.5\overline{{\rho u_{i}^{^{\prime\prime}} u_{i}^{^{\prime\prime}} }} /\overline{\rho }$ can be derived as (Nishiki et al. 2002; Chakraborty et al. 2011a, b; Lai et al. 2017; Ahmed et al. 2019; Jones and Launder 1973; Durbin and Reif 2010; Wilcox 2002):

$$ \begin{aligned} \frac{{\partial \left( {\overline{\rho }\tilde{k}} \right)}}{\partial t} + \frac{{\partial \left( {\overline{\rho }\tilde{u}_{j} \tilde{k}} \right)}}{{\partial x_{j} }}{ } & = \underbrace {{ - \overline{{\rho u_{i}^{{^{\prime\prime}}} u_{j}^{{^{\prime\prime}}} }} \frac{{\partial \tilde{u}_{i} }}{{\partial x_{j} }}}}_{{T_{1} }}{ }\underbrace {{ - \overline{{u_{i}^{{^{\prime\prime}}} }} \frac{{\partial \overline{p}}}{{\partial x_{i} }}}}_{{T_{2} }} + \underbrace {{\overline{{p^{\prime}\frac{{\partial u_{k}^{{^{\prime\prime}}} }}{{\partial x_{k} }}}} }}_{{T_{3} }} + \underbrace {{\overline{{u_{i}^{{^{\prime\prime}}} \frac{{\partial \tau_{ij} }}{{\partial x_{j} }}}} }}_{{T_{4} }}\underbrace {{ - \frac{{\partial \left( {p^{\prime}u_{i}^{{^{\prime\prime}}} } \right)}}{{\partial x_{i} }}}}_{{T_{4} }}{ } \\ & \quad \underbrace {{ - \frac{\partial }{{\partial x_{i} }}\left( {\overline{{\frac{1}{2}\rho u_{i}^{{^{\prime\prime}}} u_{k}^{{^{\prime\prime}}} u_{k}^{{^{\prime\prime}}} }} } \right)}}_{{T_{6} }} \\ \end{aligned} $$

(2)

Here, $\overline{q}, \tilde{q} = \overline{\rho q} /\overline{\rho }$ and $q^{\prime\prime} = q - \tilde{q}$ refer to the Reynolds averaged, Favre-averaged and Favre fluctuation of a general quantity $q$, respectively. In Eq. (2), the unclosed terms $T_{1}$ to $T_{6}$ are the mean velocity gradient, mean pressure gradient, pressure dilatation, viscous dissipation, pressure transport and the turbulent transport terms, respectively. Effectively, the term $T_{1}$ is closed within the framework of $k -\epsilon$ model, whereas the terms $T_{2}$ to $T_{6}$ need to be modelled. For modelling purposes, it is useful to express the viscous dissipation term $T_{4}$ as (Nishiki et al. 2002):

$$ {\text{T}}_{4} { } = \overline{{u_{i}^{{^{\prime\prime}}} \frac{{\partial \tau_{ij} }}{{\partial x_{j} }}}} = - \overline{\rho }\tilde{ \epsilon} + \underbrace {{\left[ { \overline{{u_{i}^{^{\prime\prime}} \frac{\partial }{{\partial x_{k} }}\left( {\mu \frac{{\partial u_{k}^{^{\prime\prime}} }}{{\partial x_{i} }}} \right)}} - \frac{2}{3}\overline{{u_{i}^{^{\prime\prime}} \frac{\partial }{{\partial x_{i} }}\left( {\mu \frac{{\partial u_{k}^{^{\prime\prime}} }}{{\partial x_{k} }}} \right)}} } \right]}}_{{T_{V} }} + \frac{\partial }{{\partial x_{j} }}\left( {\mu \frac{{\partial \tilde{k}}}{{\partial x_{j} }}} \right) $$

(3)

The first and second terms on the right-hand side of Eq. (3) are unclosed but the first term involves the dissipation rate of turbulent kinetic energy $\tilde{\epsilon} = \overline{{\mu \left( {\partial u_{i}^{^{\prime\prime}} /\partial x_{j} } \right)\left( {\partial u_{i}^{^{\prime\prime}} /\partial x_{j} } \right)}} /\overline{\rho }$ whereas the second term $T_{V}$ needs to be modelled. The last term is negligible for high values of turbulent Reynolds numbers and is already closed. Equation (3) suggests that the complete closure of $T_{4}$ is dependent on the closures of $\tilde{\epsilon} = \overline{{\mu \left( {\partial u_{i}^{^{\prime\prime}} /\partial x_{j} } \right)\left( {\partial u_{i}^{^{\prime\prime}} /\partial x_{j} } \right)}} /\overline{\rho }$ and $T_{V}$. It is worth noting that an alternative definition of dissipation rate of turbulent kinetic energy $\tilde{\epsilon}_{comp} = \left[ {0.5\overline{{\mu \left( {\partial u_{i}^{^{\prime\prime}} /\partial x_{j} + \partial u_{j}^{^{\prime\prime}} /\partial x_{i} } \right)\left( {\partial u_{i}^{^{\prime\prime}} /\partial x_{j} + \partial u_{j}^{^{\prime\prime}} /\partial x_{i} } \right)}} - 2/3\overline{{\mu \left( {\partial u_{k}^{^{\prime\prime}} /\partial x_{k} } \right)^{2} }} } \right]/\overline{\rho }$ exists for compressible flows (Wilcox 2002) but it does not offer any extra benefit in terms of the closure of $T_{4}$ because $\tilde{\epsilon}$ appears explicitly in the expression for $T_{4}$ (see Eq. (3)). The same approach was adopted in previous analyses by Nishiki et al. (2002) and Chakraborty et al. (2011a, b). Moreover, it has been found that there are no significant qualitative and quantitative differences between $\tilde{\epsilon}$ and $\tilde{\epsilon}_{comp}$ for this database and thus $\tilde{\epsilon}_{comp}$ is not explicitly shown here for the sake of conciseness and the implications of different definitions of the dissipation rate of kinetic energy will not be discussed further in this paper.

It is also important to note that a volume-integral of $T_{5}$ for wall bounded flows yields $\mathop \int \limits_{V} T_{5} dV = - \mathop \int \limits_{V} \overline{{\partial \left( {p^{\prime}u_{i}^{^{\prime\prime}} } \right)}} /\partial x_{i} dV = \mathop {\oint }\limits_{S} \left( {p^{\prime}u_{i}^{^{\prime\prime}} } \right).n_{i} dS$ according to the divergence theorem where $n_{i}$ is the i^th component of the surface normal. This suggests that for wall-bounded flows $\mathop \int \limits_{V} T_{5} dV$ vanishes because $u_{i}^{^{\prime\prime}}$ vanishes on the wall surface. This suggests that $T_{5}$ does not generate or destroy the turbulent kinetic energy but only redistributes it. By contrast, $\mathop \int \limits_{V} T_{3} dV = \mathop \int \limits_{V} \overline{{p^{\prime}\left( {\partial u_{i}^{^{\prime\prime}} /\partial x_{i} } \right)}} dV \ne 0$, which suggests that the pressure dilatation term $T_{3}$ can be a source/sink term. Therefore, the modelling of $T_{3}$ and $T_{5}$ are treated separately, as done in previous DNS studies (Zhang and Rutland 1995; Nishiki et al. 2002; Chakraborty et al. 2011a, b; Lai et al. 2017; Ahmed et al. 2019), which dealt with turbulent kinetic energy transport in premixed turbulent flames.

The statistical behaviour of the unclosed terms, their response to magnitude and direction of body forces and implications of body force on the closure modelling will be discussed further in Sect. 4.

3 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$ (i.e., $Pr = 0.7)$ and Zel’dovich number ${\upbeta } = T_{ac} \left( {T_{ad} - T_{0} } \right)/T_{ad}^{2}$ $\left( {\text{i.e.,}}\;{{ \upbeta }} = 6.0 \right)$, where $T_{ac}$ is the activation temperature. The heat release parameter ${\uptau } = \left( {T_{ad} - T_{0} } \right)/T_{0}$ is considered to be $4.50$. These values of $\beta$ and $\tau$ are representative of methane-air flames preheated to $T_{0}$ = 415 K. The simulations have been carried out on a rectangular domain of size $140.4{\updelta }_{{\text{Z}}} \times 70.2{\updelta }_{{\text{Z}}} \times 70.2{\updelta }_{{\text{Z}}}$, which is discretised by a grid size of 800 × 400 × 400, ensuring that the grid spacing remains smaller than the Kolmogorov length scale and 10 grid points are accommodated within the thermal flame thickness $\delta_{th} = \left( {T_{ad} - T_{0} } \right)/\max \left| {\nabla \hat{T}} \right|_{L}$. 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 the Batchelor-Townsend kinetic energy spectrum (Batchelor and Townsend 1948) is used to initialise the turbulent velocity field and the turbulence decays with time. The flame and reacting scalar fields have been initialised by a steady-state one-dimensional unstretched laminar premixed flame solution. The initial values of normalised root-mean-square turbulent velocity fluctuation $u^{\prime}/S_{L}$, integral length scale to flame thickness ratio $l_{T} /\delta_{th}$, Damköhler number $Da = l_{T} S_{L} /u^{\prime}\delta_{th}$ ($Da_{Z} = l_{T} S_{L} /u^{\prime}\delta_{Z}$) and Karlovitz number $Ka = \left( {u^{\prime}/S_{L} } \right)^{1.5} \left( {l_{T} /\delta_{th} } \right)^{ - 0.5}$ ($Ka_{Z} = \left( {u^{\prime}/S_{L} } \right)^{1.5} \left( {l_{T} /\delta_{Z} } \right)^{ - 0.5}$) are presented in Table 1. For the values of $u^{\prime}/S_{L}$ and $l_{T} /\delta_{th}$ 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, where $\left| {g^{*} } \right| = 3.12$ corresponds to a pressure gradient of $2000 Pa/m$ (Shepherd et al. 1982). This value is much greater than the typical pressure gradient observed in open premixed flames such as those in Bunsen burners, where the pressure gradient is of the order of 200 Pa/m (one-tenth of the value used in the present study). To provide more perspective, the vertical pressure gradient in the earth’s atmosphere is of the order of 10 Pa/m.

Table 1 Initial values of numerical parameters considered for the simulations

The simulation time is taken to be ~ 2.3 times the chemical time scale $\tau_{chem}$ (where $\tau_{chem} = \delta_{th} /S_{L}$), and also greater than 3.0 times the eddy turnover time $\tau_{eddy}$. This equates to 3.0–8.0 times $\tau_{eddy}$ for Set-A to Set-D, respectively. 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^{\prime}$ decayed by about 58% and $l_{T}$ increased by about 26% in comparison to their initial values in the unburned gas. The Reynolds/Favre-averaged values of a general quantity $q$ is evaluated by ensemble averaging the relevant quantities in $x_{2} - x_{3}$ planes at a given $x_{1}$ location following several previous analyses (Rutland and Cant 1994; Zhang and Rutland 1995; Nishiki et al. 2002; Chakraborty et al. 2011a, b; Lai et al. 2017; Ahmed et al. 2019; Bray et al. 1981; Veynante and Poinsot 1997).

4 Results and Discussion

The contour plots of reaction progress variable $c$ at the $x_{1} - x_{2}$ midplane of the simulation domain for $g^{*} = - 3.12, 0.0$ and $3.12$, for all sets of turbulence parameters considered, at the time when statistics were extracted are shown in Fig. 1. It can be seen from Fig. 1 that the value of $g^{*}$ clearly affects the flame morphology and Rayleigh–Taylor type of instability possibly contributes to the augmentation of flame wrinkling with increasing $g^{*}$ values (i.e., from $- 3.12$ to $3.12$). The increase in flame wrinkling with increasing $g^{*}$ and $u^{\prime}/S_{L}$ 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. 2. The turbulent flame surface area and flame speed have been evaluated using the volume integrals $A_{T} = \mathop \int \limits_{V} \left| {\nabla c} \right|dV$ and $S_{T} = \mathop \int \limits_{V} \dot{w}dV/{\uprho }_{0} A_{P}$, respectively, where ${\uprho }_{0}$ is the unburned gas density, $A_{P}$ is the projected flame surface area in the direction of mean flame propagation, and the subscripts $T$ and $L$ refer to the turbulent and laminar flame conditions, respectively. These observations are consistent with the increased flame wrinkling observed in the two-dimensional DNS study by Veynante and Poinsot (1997) and with the theory proposed by Libby (1989).

Fig. 1

Contours of reaction progress variable c for $g^{*} = - 3.12, 0.0$ and $3.12$ for all sets of turbulence parameters considered

Fig. 2

Variation of $S_{T} /S_{L}$ and $A_{T} /A_{L}$ with $g^{*}$ for all sets of turbulence parameters considered

4.1 Effects of ${{g}}^{{*}}$ on the Variations of  $\tilde{{k}}$ and ${\tilde{\epsilon}}$

The variations of normalised turbulent kinetic energy $\tilde{k}/\tilde{k}_{0}$ and its dissipation rate $\tilde{\epsilon}/\tilde{\epsilon}_{0}$ with Favre averaged reaction progress variable $\tilde{c}$ for different values of $g^{*}$, for Set-A to Set-D are exemplarily shown in Fig. 3, where $\tilde{k}_{0}$ and ${\tilde{\upepsilon}}_{0}$ are values of turbulent kinetic energy and dissipation rate at $\tilde{c} = 0.005$. It can be seen from Fig. 3 that the variations of turbulent kinetic energy and its dissipation rate within the flame brush are affected by both the magnitude and direction of $g^{*}$ (i.e., from $- 3.12$ to $3.12$). It is evident from Fig. 3 that $\tilde{k}/\tilde{k}_{0}$ increases significantly from the leading edge of the flame brush before decaying towards the product side of the flame brush for small $g^{*}$ values in Set-A (i.e., $g^{*}$ = $- 3.12$, $- 1.56$ and $0.0$), whereas a slight decay in $\tilde{k}/\tilde{k}_{0}$ across the flame brush is observed in Set-D for larger $g^{*}$ values. Moreover, in Set-D, the influence of $g^{*}$ on the variation of kinetic energy across the flame brush has been found to be much weaker than in Set-A. This may be attributed to the weakening of the body force effects with an increase in $u^{\prime}/S_{L}$, which can be explained by considering the ratio of body forces to the inertial forces due to turbulent velocity fluctuations as: $Fr_{T}^{ - 2} = {\Gamma }l_{T} /u^{^{\prime}2} = \left( {g^{*} S_{L}^{2} /{\updelta }_{Z} } \right)/\left( {u^{^{\prime}2} /l_{T} } \right)\sim g^{*} \left( {l_{T} /\delta_{Z} } \right)/\left( {u^{^{\prime}2} /S_{L}^{2} } \right)$, where $Fr_{T} = u^{\prime}/\sqrt {{\Gamma }l_{T} }$ is the turbulent Froude number. Accordingly, an increase in $u^{\prime}/S_{L}$ for a given value of $l_{T} /\delta_{Z}$ and $g^{*}$ leads to a decrease in the influence of body forces. It can also be observed from Fig. 3 that ${\tilde{\upepsilon}}/{\tilde{\upepsilon}}_{0}$ exhibits significant augmentation within the flame brush in Set-A, whereas this tendency in Set-D is much weaker than in Set-A.

Fig. 3

Variations of $\tilde{k}/\tilde{k}_{0}$ and $\tilde{\epsilon}/\tilde{\epsilon}_0$ with $\tilde{\user2{c}}$ for different values of $g^{*}$ for all sets of turbulence parameters considered

The case with $g^{*} = - 3.12$ for Set-A remains weakly turbulent with a small extent of flame wrinkling (see Figs. 1, 2). The stabilising effects of buoyancy acts to weaken turbulent fluctuations with an increase in magnitude of negative $g^{*}$, and these effects are the strongest in Set-A because of the smallest initial turbulence intensity. This leads to a smaller peak $\tilde{k}/\tilde{k}_{0}$ value for $g^{*} = - 3.12$ within the flame brush than in the $g^{*} = - 1.56$ case (see Fig. 3). This non-monotonic trend between $g^{*} = - 3.12$ and 1.56 is also observed in the variation of turbulent scalar flux $\overline{{\rho u_{1}^{^{\prime\prime}} c^{\prime\prime}}}$ for Set-A [reported elsewhere (Varma et al. 2021)]. However, this behaviour is not observed for higher values of turbulence intensity as the inertial effects tend to dominate over buoyancy effects where $\tilde{k}/\tilde{k}_{0}$ shows a monotonic trend with the variations of $g^{*}$.

4.2 Statistical Behaviour of Turbulent Kinetic Energy Transport

In order to understand the $g^{*}$ dependence of turbulent kinetic energy distribution within the flame brush, it is useful to consider the statistical behaviour of the terms $T_{1}$ to $T_{6}$ of the turbulent kinetic energy transport equation. The variation of these terms, normalised by $\delta_{Z} /\rho_{0} S_{L}^{3}$, for all the cases considered in the present study are shown in Fig. 4. The following observations can be made from Fig. 4:

• The pressure dilatation term $T_{3}$ is the major source term across all the cases. The magnitude of $T_{3}$ remains largely unaffected by $g^{*}$ but its relative importance decreases with increasing $u^{\prime}/S_{L}$.

• The viscous dissipation term $T_{4}$ is the major sink term in all the cases and its magnitude decreases with increasing $g^{*}$ (i.e., from $- 3.12$ to $3.12$), especially in cases with small initial values of $u^{\prime}/S_{L}$ (Set-A and Set-B).

• The mean velocity gradient term $T_{1}$ also acts as a sink term across all cases and its magnitude decreases with increasing $g^{*}$, thereby becoming negligible at high $g^{*}$ values (i.e., for $g^{*} = 1.56$ and $3.12$). For negative $g^{*}$ values, the magnitude of $T_{1}$ remains comparable to that of $T_{4}$.

• The mean pressure gradient term $T_{2}$ acts as a source term in all the cases and its magnitude decreases with increasing $g^{*}$. Therefore, $T_{2}$ acts as a significant leading order term for negative $g^{*}$ values (i.e., for $g^{*} = - 3.12$ and $- 1.56$) but becomes negligible at high $g^{*}$ values (i.e., for $g^{*} = 3.12$). It can also be seen that the magnitude of $T_{2}$ increases with increasing initial $u^{\prime}/S_{L}$ values.

• The pressure transport term $T_{5}$ and the turbulent transport term $T_{6}$ remain significant only for negative $g^{*}$ values (i.e., $g^{*} = { } - 3.12$ and $- 1.56$). It can be seen that $T_{5}$ predominantly assumes positive values towards the reactant side of the flame brush before becoming negative towards the product side. The term $T_{6}$ shows an opposite qualitative behaviour to that of $T_{5}$, and assumes negative values towards the unburned gas side and positive values towards the burned gas side.

Fig. 4

Variations of $\left\{ {T_{1} - T_{6} } \right\} \times \delta_{Z} /\rho_{0} S_{L}^{3}$ with $\tilde{c}$ for all cases considered

It can be seen from Figs. 3 and 4 that the combined contribution of the source terms is higher than that of the sink terms towards the unburned gas side of the flame brush, leading to an augmentation of turbulent kinetic energy magnitude, especially for lower $g^{*}$ values (compare $g^{*} = - 3.12, - 1.56$ and $0.0$). For higher $g^{*}$ values, the contributions of source and sink terms balance each other, leading to no significant production or decay of turbulent kinetic energy within the flame brush. This is more evident for cases with low initial values of $u^{\prime}/S_{L}$ (i.e., Set-A).

4.3 Modelling of the Mean Velocity Gradient term ${{\varvec{T}}}_{1}$

For modelling the mean velocity gradient term $T_{1}$, it is necessary to model the Reynolds stress $\overline{{\rho u_{i}^{^{\prime\prime}} u_{j}^{^{\prime\prime}} }}$, which is commonly modelled as the Favre averaged Reynolds stress for non-reacting turbulent flows and this is given as (Jones and Launder 1973; Durbin and Reif 2010; Wilcox 2002):

$$ \frac{{\overline{{\rho u_{i}^{^{\prime\prime}} u_{j}^{^{\prime\prime}} }} }}{{\overline{\rho }}} = - \nu_{t} \left( {\frac{{\partial \tilde{u}_{i} }}{{\partial x_{j} }} + \frac{{\partial \tilde{u}_{j} }}{{\partial x_{i} }}} \right) + \frac{2}{3}\delta_{ij} \left( {\nu_{t} \frac{{\partial \tilde{u}_{k} }}{{\partial x_{k} }} + \tilde{k}} \right) $$

(4)

where $\nu_{t}$ is the kinematic eddy viscosity given by $\nu_{t} = C_{\mu } \tilde{k}^{2} /\tilde{\epsilon}$ with the standard model constant $C_{\mu } = 0.09$. Using Eq. (4), the Reynolds stress component $\overline{{\rho u_{1}^{^{\prime\prime}} u_{1}^{^{\prime\prime}} }}$ is given as:

$$ \frac{{\overline{{\rho u_{1}^{^{\prime\prime}} u_{1}^{^{\prime\prime}} }} }}{{\overline{\rho }}} = - \frac{4}{3}\nu_{t} \frac{{\partial \tilde{u}_{1} }}{{\partial x_{1} }} + \frac{2}{3}\tilde{k} $$

(5)

The predictions of the model given by Eq. (5) are shown in Fig. 5 alongside the $\overline{{\rho u_{1}^{^{\prime\prime}} u_{1}^{^{\prime\prime}} }} /\overline{\rho }$ extracted from DNS data. It can be seen from Fig. 5 that this model underpredicts the magnitude of $\overline{{\rho u_{1}^{^{\prime\prime}} u_{1}^{^{\prime\prime}} }} /\overline{\rho }$ for most of the cases, especially for lower $g^{*}$ values (i.e., $g^{*} = - 3.12, - 1.56,$ and $0.0$), whereas the predictions of Eq. (5) improve significantly for high $g^{*}$ values (i.e., $g^{*} = 1.56$ and $3.12$). The probability density function (pdf) of reaction progress variable $c$ and the joint pdf between the velocity vector $\vec{u}$ and $c$ may be represented according to the Bray–Moss–Libby (BML) analysis (Bray et al. 1985) as:

$$ \begin{aligned} & P\left( c \right) = \alpha_{c} \delta \left( c \right) + \beta_{c} \delta \left( {1 - c} \right) + O\left( {\gamma_{c} } \right)\;{\text{and}} \\ & P\left( {\vec{u},c} \right) = \alpha_{c} P_{R} \left( {\vec{u}} \right)\delta \left( c \right) + \beta_{c} P_{P} \left( {\vec{u}} \right)\delta \left( {1 - c} \right) + O\left( {\gamma_{c} } \right) \\ \end{aligned} $$

(6)

where $\alpha_{c} , \beta_{c} ,$ and $\gamma_{c}$ are the coefficient functions corresponding to the contributions of reactants, products and reacting gas, respectively, whereas $P_{R} \left( {\vec{u};\vec{x}} \right)$ and $P_{P} \left( {\vec{u};\vec{x}} \right)$ refer to the pdfs of  $\vec{u}$ in the unburned and fully burned gases, respectively. The contribution of the term $O\left( {\gamma_{c} } \right)$ is considered to be negligible under the strict flamelet assumption (Bray et al. 1985). Using Eq. (6), the mean velocity in the $x_{1}$-direction may be written as (Bray et al. 1985):

$$ \tilde{u}_{1} = \left( {1 - \tilde{c}} \right) \overline{u}_{{1}_{R}} + \tilde{c}\overline{u}_{{1}_{P}} + O\left( {\gamma_{c} } \right) $$

(7)

where $\overline{u}_{{1}_{R}}$ and $\overline{u}_{{1}_{P}}$ are the conditional mean velocities in the $x_{1}$-direction in the reactants and products, respectively. Using Eqs. (6) and (7), $\overline{{\rho u_{1}^{^{\prime\prime}} u_{1}^{^{\prime\prime}} }} /\overline{\rho }$ can be expressed as:

$$ \frac{{\overline{{\rho u_{1}^{^{\prime\prime}} u_{1}^{^{\prime\prime}} }} }}{{\overline{\rho }}} = \left( {\overline{u}_{{1}_{P}} - \overline{u}_{{1}_{R}} } \right)^{2} \tilde{c}\left( {1 - \tilde{c}} \right) + \left( {1 - \tilde{c}} \right)\overline{{\left( {u_{1}^{^{\prime}2} } \right)}}_{R} + \tilde{c}\overline{{\left( {u_{1}^{^{\prime}2} } \right)}}_{P} + O\left( \gamma \right) $$

(8)

where $\overline{{\left( {u_{1}^{^{\prime}2} } \right)}}_{R}$ and $\overline{{\left( {u_{1}^{^{\prime}2} } \right)}}_{P}$ are the conditionally averaged mean-squared velocity fluctuations in the $x_{1}$-direction. The second and third terms on the right-hand side of Eq. (8) represent the non-reacting contribution and this can be modelled for unity Lewis number cases as (Chakraborty et al. 2011a, b):

$$ \left( {1 - \tilde{c}} \right)\overline{{\left( {u_{1}^{^{\prime}2} } \right)}}_{R} + \tilde{c}\overline{{\left( {u_{1}^{^{\prime}2} } \right)}}_{P} = \frac{2}{3}\tilde{k} - \frac{4}{3}\frac{{C_{\mu } \rho_{0} \tilde{k}^{2} }}{{\overline{\rho }\tilde{\epsilon}}}\frac{{\partial \tilde{u}_{1} }}{{\partial x_{1} }} $$

(9)

where $\overline{\rho }\tilde{\epsilon}/\rho_{0}$ is the density-weighted dissipation rate which is used to account for the changes in kinematic viscosity with temperature. The contribution of the first term on the right-hand side of Eq. (8) can be modelled using expressions for the turbulent scalar flux and the variance of the reaction progress variable, both of which may be obtained using the BML framework (Bray et al. 1985) and are given as:

$$ \frac{{\overline{{\rho u_{1}^{^{\prime\prime}} c^{\prime\prime}}} }}{{\overline{\rho }}} = \tilde{c}\left( {1 - \tilde{c}} \right)\left( {\overline{u}_{{1}_{P}} - \overline{u}_{{1}_{R}} } \right) + O\left( {\gamma_{c} } \right) $$

(10)

$$ \frac{{\overline{{\rho c^{\prime\prime}c^{\prime\prime}}} }}{{\overline{\rho }}} = \tilde{c}\left( {1 - \tilde{c}} \right) + O\left( {\gamma_{c} } \right) $$

(11)

Using Eqs. (9)–(11), a model for $\overline{{\rho u_{1}^{^{\prime\prime}} u_{1}^{^{\prime\prime}} }} /\overline{\rho }$ may be given as (Chakraborty et al. 2011a, b):

$$ \frac{{\overline{{\rho u_{1}^{^{\prime\prime}} u_{1}^{^{\prime\prime}} }} }}{{\overline{\rho }}} = - \frac{4}{3}\frac{{C_{\mu } \rho_{0} \tilde{k}^{2} }}{{\overline{\rho }\tilde{\epsilon}}}\frac{{\partial \tilde{u}_{1} }}{{\partial x_{1} }} + \frac{2}{3}\tilde{k} + \frac{{\overline{{\rho u_{1}^{^{\prime\prime}} c^{\prime\prime}}} \overline{{\rho u_{1}^{^{\prime\prime}} c^{\prime\prime}}} }}{{\overline{\rho } \overline{{\rho c^{\prime\prime}c^{\prime\prime}}} }} $$

(12)

The predictions of the model given by Eq. (12) are shown in Fig. 5 where $\tilde{\epsilon}$, $\overline{{\rho u_{1}^{^{\prime\prime}} c^{\prime\prime}}}$ and $\overline{{\rho c^{\prime\prime}c^{\prime\prime}}}$ are extracted from the DNS data in the spirit of a priori analysis, and it can be seen that this model satisfactorily predicts the variation of $\overline{{\rho u_{1}^{^{\prime\prime}} u_{1}^{^{\prime\prime}} }} /\overline{\rho }$ across the flame brush for all the cases considered in the present study. However, it is worthwhile to note that $\tilde{\epsilon}$, $\overline{{\rho u_{1}^{^{\prime\prime}} c^{\prime\prime}}}$ and $\overline{{\rho c^{\prime\prime}c^{\prime\prime}}}$ are unclosed quantities and therefore the accuracy of the prediction of Eq. (12) in an actual RANS simulation depends on the modelling accuracies of $\tilde{\epsilon}$, $\overline{{\rho u_{1}^{^{\prime\prime}} c^{\prime\prime}}}$ and $\overline{{\rho c^{\prime\prime}c^{\prime\prime}}}$.

Fig. 5

Variations of $\overline{{\rho u_{1}^{^{\prime\prime}} u_{1}^{^{\prime\prime}} }} /\overline{\rho }S_{L}^{2}$ with $\tilde{c}$ along with the predictions of Eq. (5) and Eq. (12) for all cases considered

4.4 Modelling of the Mean Pressure Gradient term ${{T}}_{2}$

For unity Lewis number flames, the density variation across the flame brush can be expressed as (Bray et al. 1985):

$$ \frac{1}{\rho } = \frac{{\left( {1 + \tau c} \right)}}{{\rho_{0} }} = \frac{c}{{\rho_{P} }} + \frac{1 - c}{{\rho_{0} }} $$

(13)

where $\rho_{P} = \rho_{0} /\left( {1 + \tau c} \right)$ is the density of the products. Equation (13) can be used to obtain a relation for $\overline{{u_{1}^{^{\prime\prime}} }}$ in the following manner (Nishiki et al. 2002):

$$ \overline{{u_{1}^{^{\prime\prime}} }} = \mathop \int \limits_{0}^{1} u_{1}^{^{\prime\prime}} P\left( c \right)dc = \mathop \int \limits_{0}^{1} \rho u_{1}^{^{\prime\prime}} \frac{P\left( c \right)}{\rho }dc = \left( {\frac{1}{{\rho_{P} }} - \frac{1}{{\rho_{0} }}} \right)\overline{{\rho u_{1}^{^{\prime\prime}} c^{\prime\prime}}} $$

(14)

Equation 14 gives rise to the model for the mean pressure gradient term $T_{2}$ by Nishiki et al. (2002), which can be written as:

$$ T_{2} = - \overline{{u_{1}^{^{\prime\prime}} }} \frac{{\partial \overline{p}}}{{\partial x_{1} }} = - \frac{\tau }{{\rho_{0} }}\overline{{\rho u_{1}^{^{\prime\prime}} c^{\prime\prime}}} \frac{{\partial \overline{p}}}{{\partial x_{1} }} $$

(15)

The predictions of the model given by Eq. (15) with $\overline{{\rho u_{1}^{^{\prime\prime}} c^{\prime\prime}}}$ extracted from the DNS data are shown in Fig. 6, alongside $T_{2}$ extracted from DNS data. Since Eq. (15) is an exact form of $T_{2}$ subject to the BML assumption (Bray et al. 1985), the predictions of the model show excellent agreement with the DNS data. However, $\overline{{\rho u_{1}^{^{\prime\prime}} c^{\prime\prime}}}$ is an unclosed quantity and therefore the successful prediction of $T_{2}$ will depend on the closure of $\overline{{\rho u_{1}^{^{\prime\prime}} c^{\prime\prime}}}$, which is beyond the scope of the current analysis but interested readers are referred to Veynante and Poinsot (1997) for further information related to the buoyancy effects on the modelling of turbulent scalar flux in the context of RANS.

Fig. 6

Variations of $T_{2} \times \delta_{Z} /\rho_{0} S_{L}^{3}$ with $\tilde{c}$ along with the predictions of Eq. (15) for all cases considered

4.5 Modelling of the Pressure Dilatation Term $T_{3}$

According to Zhang and Rutland (1995), the pressure dilatation term may be expressed in the following manner:

$$ T_{3} = \overline{{p^{\prime}\frac{{\partial u_{i}^{^{\prime\prime}} }}{{\partial x_{i} }}}} = \frac{1}{V}\mathop \int \limits_{V}^{ } p^{\prime}\frac{{\partial u_{i}^{^{\prime\prime}} }}{{\partial x_{i} }}dV = \frac{1}{V}\mathop \int \limits_{S}^{ } \mathop \int \limits_{\delta }^{ } p^{\prime}\frac{{\partial u_{i}^{^{\prime\prime}} }}{{\partial x_{i} }}d\xi .dS $$

(16)

where $dS$ and $dV$ refer to the elemental surface and volume elements. Considering the definition of generalised flame surface density which is given by ${\Sigma } = S/V = \overline{{\left| {\nabla c} \right|}}$ (Boger et al. 1998; Candel and Poinsot 1990) and based on Eq. (16) an expression for $T_{3}$ may be given as $T_{3} = I_{s} {\Sigma }$, where $I$ is given as (Lai et al. 2017):

$$ \begin{aligned} & I = \mathop \int \limits_{\delta }^{ } \left( {p - \overline{p}} \right)\left( {\frac{{du_{n} }}{d\xi } - \frac{{\partial \tilde{u}_{i} }}{{\partial x_{i} }}} \right)d\xi \\ & {\text{or}}\;I = \mathop \int \limits_{\delta }^{ } \left( {p - \overline{p}} \right)\left( {\frac{{d\left( {u_{n} - u_{0} } \right)}}{d\xi } - \frac{{\partial \tilde{u}_{i} }}{{\partial x_{i} }}} \right)d\xi \\ \end{aligned} $$

(17)

where $\xi$ is the local flame normal direction, $u_{n}$ is the velocity component in the local flame normal direction and $u_{0}$ is the component of flame normal velocity in the reactants. A one-dimensional evaluation along the flame normal direction gives (Zhang and Rutland 1995):

$$ p - p_{0} = \rho_{0} u_{0} \left( {u_{n} - u_{0} } \right) = \rho_{0} u_{0} {\Delta }u $$

(18)

Substituting Eq. (18) into Eq. (17) followed by surface averaging gives:

$${\langle I \rangle}_{s} = {\langle\left( {p_{0} - \overline{p}} \right){\Delta u}\rangle}_{s} - \frac{1}{2}{\langle \rho_{0} u_{0} \left( {{\Delta }u} \right)^{2}\rangle}_{s} -{ \Bigg\langle \left( {p - \overline{p}} \right)\left( \frac{\partial \tilde{u}_{i} }{\partial x_{i} } \right)\delta \Bigg\rangle}_{s} $$

(19)

where $\delta$ is a length scale that characterises the flame thickness. For unity Lewis number flames, the first and second terms on the right-hand side of Eq. (19) may be scaled as follows (Lai et al. 2017):

$${\langle \left(p_{0} - {\overline{p}} \right){\Delta}u \rangle}_{s} \sim \tau^{2} {\rho_{0}} {S_{L}^{3}} \, ;\,{ \langle {\rho_{0}} u_{0} \left({\Delta u} \right)^{2} \rangle}_{s} \sim \tau^{2} {\rho_{0}} {S_{L}^{3}}$$

(20)

whereas the last term is scaled according to:

$${\langle \left( {p - \overline{p}} \right)\left( \frac{\partial \tilde{u}_{i} }{\partial x_{i} } \right)\delta \rangle_s} \sim \tau \rho_{0} S_{L}^{3} \frac{U_{mean} }{u^{\prime}} \frac{1}{Da} $$

(21)

where $U_{mean}$ is the characteristic mean velocity. Therefore, this last term remains negligible for high Damköhler number flames but may remain significant for low Damköhler number flames. Zhang and Rutland (1995) proposed the following model for the pressure dilatation term $T_{3}$:

$$ T_{3} = \frac{1}{2}C_{z} \rho_{0} S_{L}^{3} \tau^{2} \tilde{c}{\Sigma } $$

(22)

where $C_{z}$ is a model parameter which takes a value of $1.35$ and $\tilde{c}$ on the right-hand side accounts for the assumed dependence of the pressure fluctuation on $\tilde{c}$. The model given by Eq. (22) will henceforth be referred to as the PDZ model. Nishiki et al. (2002) proposed another model for the pressure dilatation term $T_{3}$ using a FSD-based closure for the mean reaction rate (i.e., $\overline{{\dot{w}}} = \rho_{0} S_{L} {\Sigma }$) and this is given as:

$$ T_{3} = C_{3} S_{L}^{2} \tau^{2} \overline{{\dot{w}}} $$

(23)

where $C_{3}$ is a model parameter and takes a value of 0.35. This model will henceforth be referred to as the PDN model. The predictions of the PDZ and PDN models have been shown in Fig. 7 alongside the pressure dilatation term $T_{3}$ extracted from the DNS data. It can be seen from Fig. 7 that the PDN model satisfactorily predicts the variations of $T_{3}$ across all the cases, whereas the PDZ model tends to underpredict the magnitude of $T_{3}$ towards the unburned gas side of the flame brush. Both the PDZ and PDN models tend to overpredict the magnitude of $T_{3}$ for cases with negative $g^{*}$ values, especially for lower initial $u^{\prime}/S_{L}$ values (i.e., Set-A and Set-B). It is worth noting that in Fig. 7 both ${\Sigma }$ and $\overline{{\dot{w}}}$ are extracted from DNS data for the PDZ and PDN models, respectively, in the spirit of a priori analysis. However, both ${\Sigma }$ and $\overline{{\dot{w}}}$ need closure in the context of RANS simulations and it is possible to model the mean reaction rate $\overline{{\dot{w}}}$ in terms of ${\Sigma }$. The effects of $g^{*}$ on the closure of the FSD ${\Sigma }$ and the FSD based mean reaction rate $\overline{{\dot{w}}}$ have been addressed elsewhere (Varma et al. 2021) by the present authors and the interested readers are referred to Varma et al. (2021) for further information in this regard.

Fig. 7

Variations of ${\varvec{T}}_{3} \times \varvec{ \delta }_{{\varvec{Z}}} /{\varvec{\rho}}_{0} {\varvec{S}}_{{\varvec{L}}}^{3}$ with $\tilde{c}$ along with the PDZ and PDN model predictions for all cases considered

It is worth noting that the contribution of the term shown in Eq. (21) may play a significant role for low Damköhler number (i.e., $Da \ll 1$) flames, but it can be seen from Fig. 7 that for the case with the lowest Damköhler number in the present study, i.e., Set-D with $u^{\prime}/S_{L} = 10.0$, the PDN model satisfactorily predicts the magnitude and variation of the pressure dilatation term across the flame brush and therefore the contribution of the term given by Eq. (21) remains weak for the range of the parameters considered here. Therefore, a modification of the PDN model was not attempted in this analysis.

4.6 Modelling of the Viscous Dissipation and Molecular Diffusion term $T_{4}$

The variations of $\nabla \left( {\mu \nabla \tilde{k}} \right)$, $\left( { - \overline{\rho }\tilde{\epsilon}} \right)$, $T_{V}$ and ${\text{T}}_{4}$ for all cases considered in the present study are shown in Fig. 8. It can be observed from Fig. 8 that the $T_{V}$ term acts as a sink term for all the cases, which is consistent with the findings of Nishiki et al. (2002). Figure 8 also shows that the contributions of $\nabla \left( {\mu \nabla \tilde{k}} \right)$ and $T_{V}$ are negligible in comparison to that of $\left( { - \overline{\rho }\tilde{\epsilon}} \right)$ for most of the cases, except for the cases with negative $g^{*}$ values with small initial $u^{\prime}/S_{L}$ (e.g., for $g^{*} = - 3.12$ and $- 1.56$ in Set-A) where the effects of turbulence are rendered weak due to the stable density stratification. According to Nishiki et al. (2002), the order of magnitude of $T_{V}$ may be given as:

$$ T_{V} = \overline{{u_{i}^{^{\prime\prime}} \frac{\partial }{{\partial x_{k} }}\left( {\mu \frac{{\partial u_{k}^{^{\prime\prime}} }}{{\partial x_{i} }}} \right)}} - \frac{2}{3}\overline{{u_{i}^{^{\prime\prime}} \frac{\partial }{{\partial x_{i} }}\left( {\mu \frac{{\partial u_{k}^{^{\prime\prime}} }}{{\partial x_{k} }}} \right)}} \sim - O\left( {\tilde{k}^{1/2} \rho_{0} \tau S_{L}^{2} {\Sigma }} \right) $$

(24)

where $u_{i}^{^{\prime\prime}}$, $\mu$ and $\partial u_{k}^{^{\prime\prime}} /\partial x_{k}$ are scaled using $\tilde{k}^{1/2}$, $\rho_{0} S_{L} \delta_{L}$ and $\tau S_{L} {\Sigma }\sim \tau S_{L} /\delta_{th}$, respectively (Nishiki et al. 2002). This gives rise to the following model for $T_{V}$ given by Nishiki et al. (2002):

$$ T_{V} = - C_{add} \rho_{0} \tau \tilde{k}^{1/2} S_{L}^{2} {\Sigma } $$

(25)

where $C_{add}$ = 0.25 (Nishiki et al. 2002) is a model constant. The model given by Eq. (25) will henceforth be denoted as the VN model and the predictions of this model when ${\Sigma }$ is extracted from DNS data are shown in Fig. 9 alongside the $T_{V}$ term extracted from DNS data. It can be seen from Fig. 9 that the VN model gives satisfactory predictions of the $T_{V}$ term for the cases with low initial $u^{\prime}/S_{L}$ (Set-A and Set-B), whereas for cases with high initial $u^{\prime}/S_{L}$ (Set-C and Set-D) the VN model overpredicts the magnitude of $T_{V}$ towards the unburned gas side of the flame brush. It is worthwhile to note that in the VN model, the spatial derivative $\partial u_{k}^{^{\prime\prime}} /\partial x_{k}$ is scaled with respect to the length scale given by the thermal flame thickness $\delta_{th}$, which may not be valid in the thin reaction zones regime of premixed combustion where the enhanced turbulent mixing in the preheat zone causes the effects of dilatation to spread over a larger distance (Chakraborty et al. 2011a). Chakraborty et al. (2011a) used the integral length scale of turbulence $l_{T} \sim \tilde{k}^{3/2}/\tilde{\epsilon}$ to scale the spatial derivative $\partial u_{k}^{^{\prime\prime}} /\partial x_{k}$ and obtained an order of magnitude estimate for $T_{V}$ which is given as:

$$ T_{V} \sim O\left( { - \tilde{\epsilon}\mu_{0} {\Delta }u_{s} /\tilde{k}\delta_{th} } \right) $$

(26)

where ${\Delta }u_{s}$ is an appropriate slip velocity which can be expressed as ${\Delta }u_{s} = \tau S_{L}$ for unity Lewis number flames. Equation (26) may then be re-written as (Chakraborty et al. 2011a):

$$ T_{V} \sim O\left( { - \tilde{k}^{1/2} \rho_{0} \tau S_{L}^{2} {\Sigma } \times Re_{t}^{ - 1/2} Da^{ - 1/2} } \right) $$

(27)

Equation (27) suggests that the VN model is expected to overpredict the magnitude of $T_{V}$ when the turbulent Reynolds number $Re_{t}$ is higher, which is observed in the case of Set-D in Fig. 9. Using Eq. (27), Chakraborty et al. (2011a) proposed another model for $T_{V}$ and for unity Lewis numbers, this may be given as:

$$ T_{V} = - C_{\nu } \frac{{\tilde{\epsilon}}}{{\tilde{k}}}\mu_{0} \tau S_{L} \frac{{\tilde{c}^{p} \left( {1 - \tilde{c}} \right)^{q} }}{{\delta_{th} }} $$

(28)

where $C_{\nu }$ is a model parameter and $\tilde{c}^{p} \left( {1 - \tilde{c}} \right)^{q}$ is used to capture the variation of $T_{V}$ across the flame brush. This model will henceforth be referred to as the VCKC model. In the present study, it has been found that for $p = 1.8$, $q = 2.1$ and $C_{\nu } = 0.03 + 0.0145$ erfc $\left( {0.75g^{*} + 1.223} \right)$, the VCKC model predicts the variation of $T_{V}$ satisfactorily for all cases considered, which can be observed from Fig. 9. It should be noted that the effects of the body forces have been incorporated into the model parameter $C_{\nu }$ to capture the variation of $T_{V}$ with $g^{*}$.

Fig. 8

Variations of $\nabla \left( {{\varvec{\mu}}\nabla \tilde{\user2{k}}} \right)$, $\left( { - \overline{\rho }{\tilde{\epsilon}}} \right)$, $T_{V}$ and $T_{4}$ with $\tilde{c}$ for all the cases considered. All terms have been normalized with respect to $\delta_{Z} /\rho_{0} S_{L}^{3}$

Fig. 9

Variations of $T_{V} \times \delta_{Z} /\rho_{0} S_{L}^{3}$ with $\tilde{c}$ along with the predictions of the VN model and the VCKC model for all cases considered

It is worth noting that the closures of Eqs. (25) and (28) depend also on the closures of ${\Sigma }$ and $\tilde{\epsilon}$, which are also unclosed terms. In the context of $k -\epsilon$ model (Jones and Launder 1973) a separate modelled transport equation is solved for $\tilde{\epsilon}$ (Jones and Launder 1973; Durbin and Reif 2010; Wilcox 2002), whereas the effects of buoyancy on the modelling of ${\Sigma }$ for this database have been addressed elsewhere (Varma et al. 2021) by the present authors.

4.7 Modelling of the Pressure Transport Term $T_{5}$

The pressure transport term $T_{5}$ may be written as (Rogallo 1981):

$$ T_{5} = - \frac{{\overline{{\partial \left( {p^{\prime}u_{i}^{^{\prime\prime}} } \right)}} }}{{\partial x_{i} }} = - \overline{{u_{i}^{^{\prime\prime}} \frac{{\partial p^{\prime}}}{{\partial x_{i} }} }} - \overline{{p^{\prime}\frac{{\partial u_{i}^{^{\prime\prime}} }}{{\partial x_{i} }}}} = - \overline{{u_{i}^{^{\prime\prime}} \frac{{\partial p^{\prime}}}{{\partial x_{i} }} }} - T_{3} $$

(29)

where $- \overline{{u_{i}^{^{\prime\prime}} \partial p^{\prime}/\partial x_{i} }}$ is the fluctuating pressure gradient term. The modelling of $T_{3}$ has been discussed previously and hence for modelling of $T_{5}$ it is necessary to model the fluctuating pressure gradient term. Launder et al. (1975) proposed a model for this term, which will henceforth be referred to as the PTL model, and this may be written as (Bray et al. 2000):

$$ - \overline{{u_{i}^{^{\prime\prime}} \frac{{\partial p^{\prime}}}{{\partial x_{i} }} }} = C_{1} \overline{{\rho u_{i}^{^{\prime\prime}} u_{j}^{^{\prime\prime}} }} \frac{{\partial \tilde{u}_{i} }}{{\partial x_{j} }} - C_{2} \overline{\rho }\tilde{\epsilon} $$

(30)

where $C_{1}$ and $C_{2}$ are model parameters that take values of 1.5 and 0.2, respectively. Strahle (1983) modelled the fluctuating pressure gradient term as:

$$ - \overline{{u_{i}^{^{\prime\prime}} \frac{{\partial p^{\prime}}}{{\partial x_{i} }} }} = \frac{1}{2}C_{st} \overline{{\rho u_{i}^{^{\prime\prime}} u_{j}^{^{\prime\prime}} }} \frac{{\partial \tilde{u}_{i} }}{{\partial x_{j} }} $$

(31)

where $C_{st}$ is a model constant of the order of unity ($C_{st} = 2.5$ in the present study) and this model will hereafter be called the PTS model. Another model, proposed by Zhang and Rutland (1995) (and referred to as the PTZR model), can be written as:

$$ - \overline{{u_{i}^{^{\prime\prime}} \frac{{\partial p^{\prime}}}{{\partial x_{i} }} }} = \frac{1}{2}\rho_{0} S_{L}^{3} {\Sigma }\tau^{3} \frac{{\tilde{c}\left( {1 - \tilde{c}} \right)}}{{\left( {1 + \tau } \right)}} $$

(32)

Domingo and Bray (2000) used an alternate approach and modelled the quantity $- \overline{{u_{i}^{^{\prime\prime}} \partial p/\partial x_{j} }}$ in the following manner:

$$ \begin{aligned} - \overline{{u_{i}^{^{\prime\prime}} \frac{\partial p}{{\partial x_{j} }}}} & = - \frac{{\overline{{\rho u_{i}^{^{\prime\prime}} c^{\prime\prime}}} }}{{\overline{\rho }\left( {1 + \tau \tilde{c}} \right)}}\left[ { - \frac{{\partial \overline{p}_{R} }}{{\partial x_{i} }} + \left( {1 + \tau } \right)\frac{{\partial \overline{p}_{P} }}{{\partial x_{i} }}} \right] - \frac{{\left( {1 - \tilde{c}} \right)}}{{\left( {1 + \tau \tilde{c}} \right)}}\overline{{u_{R}^{^{\prime}} \frac{{\partial p_{R}^{^{\prime}} }}{{\partial x_{i} }}}} \\ & \quad - \frac{{\left( {1 + \tau } \right)\tilde{c}}}{{\left( {1 + \tau \tilde{c}} \right)}}\overline{{u_{P,i}^{^{\prime}} \frac{{\partial p_{P}^{^{\prime}} }}{{\partial x_{i} }}}} + \frac{1}{2}\rho_{0} \tau S_{L}^{2} \Sigma \langle \vec{N} \cdot \vec{M}_{i} \rangle \frac{{\overline{{\rho u_{i}^{^{\prime\prime}} c^{\prime\prime}}} }}{{\overline{\rho }\left( {1 - \tilde{c}} \right)}} \\ & \quad + 0.35\rho_{0} \tau^{2} S_{L}^{3} {\Sigma }\langle(\vec{N} \cdot \vec{M}_{i} )\rangle\langle(\vec{N} \cdot \vec{M}_{i} )\rangle \\ \end{aligned} $$

(33)

where $\vec{N} = - \nabla c/\left| {\nabla c} \right|$ is the local flame normal vector and $\vec{M}_{i}$ is the component of the unit normal vector that describes the direction of mean flame propagation. The quantity $- \overline{{u_{i}^{^{\prime\prime}} \partial p/\partial x_{j} }}$ may be written as:

$$ - \overline{{u_{i}^{^{\prime\prime}} \frac{\partial p}{{\partial x_{j} }}}} = - \overline{{u_{i}^{^{\prime\prime}} }} \frac{{\partial \overline{p}}}{{\partial x_{i} }} - \overline{{u_{i}^{^{\prime\prime}} \frac{{\partial p^{\prime}}}{{\partial x_{i} }}}} = T_{2} - \overline{{u_{i}^{^{\prime\prime}} \frac{{\partial p^{\prime}}}{{\partial x_{i} }} }} $$

(34)

The modelling of the $T_{2}$ term has been discussed previously and the model given by Eqs. (33) and (34) will be referred to as the PTDB model. Since the PTDB model was derived for flames in the corrugated flamelets regime of premixed turbulent combustion, Chakraborty et al. (2011a) proposed a new model for flames in the thin reaction zones regime. This model was developed as a combination of the PTL and PTDB models, using a linear bridging based on the segregation factor $g = \overline{{\rho c^{\prime\prime}c^{\prime\prime}}} /\overline{\rho }\tilde{c}\left( {1 - \tilde{c}} \right)$. This model may be given as (Chakraborty et al. 2011a):

$$ \begin{aligned} - \overline{{u_{i}^{^{\prime\prime}} \frac{{\partial p^{\prime}}}{{\partial x_{i} }}}} & = - g\frac{{\overline{{\rho u_{i}^{^{\prime\prime}} c^{{\prime \prime }} }} }}{{\overline{\rho }\left( {1 + \tau \tilde{c}} \right)}}\left[ { - \frac{{\partial \overline{p}_{R} }}{{\partial x_{i} }} + \left( {1 + \tau } \right)\frac{{\partial \overline{p}_{P} }}{{\partial x_{i} }}} \right] - g\frac{{\left( {1 - \tilde{c}} \right)}}{{\left( {1 + \tau \tilde{c}} \right)}}\overline{{u_{R}^{^{\prime}} \frac{{\partial p_{R}^{^{\prime}} }}{{\partial x_{i} }}}} \\ & \quad - g\frac{{\left( {1 + \tau } \right)\tilde{c}}}{{\left( {1 + \tau \tilde{c}} \right)}}\overline{{u_{P,i}^{^{\prime}} \frac{{\partial p_{P}^{^{\prime}} }}{{\partial x_{i} }}}} + \frac{1}{2}g\rho_{0} \tau S_{L}^{2} \Sigma \langle \vec{N}.\vec{M}_{i} \rangle \frac{{\overline{{\rho u_{i}^{^{\prime\prime}} c^{{\prime \prime }} }} }}{{\overline{\rho }\left( {1 - \tilde{c}} \right)}} \\ & \quad + 0.35\rho_{0} g\tau^{2} S_{L}^{3} \Sigma\langle (\vec{N}.\vec{M}_{i} )\rangle\langle(\vec{N}.\vec{M}_{i} )\rangle + \left( {1 - g} \right)C_{1N} \overline{{\rho u_{i}^{^{\prime\prime}} u_{j}^{^{\prime\prime}} }} \frac{{\partial \tilde{u}_{i} }}{{\partial x_{j} }} \\ & \quad - \left( {1 - g} \right)C_{2N} \overline{\rho }\tilde{\epsilon} + (g\tau /\rho_{0} )\overline{{\rho u_{i}^{^{\prime\prime}} c^{\prime\prime}}} \partial \overline{p}/\partial x_{i} \\ \end{aligned} $$

(35)

where the model parameters $C_{1N} = 12.64 + 8.42$ erf $\left( {0.67g^{*} - 1.0} \right)$ and $C_{2N} = 1.07 - 0.931$ erf $\left( {1.12g^{*} + 1.31} \right)$ have been modified compared to the original model (Chakraborty et al. 2011a) to incorporate the dependence on the value of $g^{*}$. This new model will hereafter be called the PCKC model.

The predictions of the PTL, PTS, PTZR, PTDB and PCKC models are shown in Fig. 10, alongside the fluctuating pressure gradient term extracted from the DNS data. It can be seen from Fig. 10 that only the PCKC model satisfactorily predicts the variation of the fluctuating pressure gradient across the flame brush for all cases, although a slight over-prediction is observed for high $g^{*}$ values ( $g^{*} = 3.12$).

Fig. 10

Variations of $- \overline{{u_{i}^{^{\prime\prime}} \partial p^{\prime}/\partial x_{i} }} \times \delta_{Z} /\rho_{0} S_{L}^{3}$ with $\tilde{c}$ along with the predictions of the PTL, PTS, PTZR, PTDB and PCKC models for all cases considered

It should further be noted that the predictions of the PTDB and PCKC models depend on the accuracy of the sub-models for $\overline{p}_{R}$, $\overline{p}_{P}$ and the term $- \left[ {\left( {1 - \tilde{c}} \right) \times \overline{{u_{R,i}^{^{\prime}} \partial p_{R}^{^{\prime}} /\partial x_{i} }} + \left( {1 + \tau } \right)\tilde{c} \times \overline{{u_{{P,{\text{i}}}}^{^{\prime}} \partial p_{P}^{^{\prime}} /\partial x_{{\text{i}}} }} } \right]/\left( {1 + \tau \tilde{c}} \right)$, which is beyond the scope of the present work and is therefore not covered here and these quantities are extracted from DNS data for the predictions of these models in Fig. 10. Moreover, the PTL and PTZR model predictions are shown in Fig. 10 for $\tilde{\epsilon}$ and ${\Sigma }$ extracted from DNS data but their closures in RANS simulations are likely to affect the predictions of these pressure transport models.

4.8 Modelling of the Turbulent Transport Term $T_{6}$

For statistically planar flames, it is necessary to model the turbulent flux of kinetic energy $\overline{{\rho u_{1}^{^{\prime\prime}} k}} = \overline{{\rho u_{1}^{^{\prime\prime}} u_{j}^{^{\prime\prime}} u_{j}^{^{\prime\prime}} }} /2$ to model the turbulent transport term $T_{6}$. According to the BML analysis (Bray et al. 1985), this term is given as:

$$ \overline{{\rho u_{1}^{^{\prime\prime}} k}} = \frac{1}{2}\overline{{\rho u_{1}^{^{\prime\prime}} u_{j}^{^{\prime\prime}} u_{j}^{^{\prime\prime}} }} = \overline{\rho }\left( {1 - \tilde{c}} \right)\overline{{u_{1R}^{^{\prime}} k_{R} }} + \overline{\rho }\tilde{c}\overline{{u_{1P}^{^{\prime}} k_{P} }} + \overline{\rho }\tilde{c}\left( {1 - \tilde{c}} \right)\left( {\overline{k}_{P} - \overline{k}_{R} } \right) + \frac{1}{2}\overline{\rho }\tilde{c}\left( {1 - \tilde{c}} \right)\left( {\overline{u}_{1P} - \overline{u}_{1R} } \right)^{3} \left( {1 - 2\tilde{c}} \right) $$

(36)

where $\overline{k}_{R}$ and $\overline{k}_{P}$ are the conditional values of the turbulent kinetic energy in the reactants and products, respectively. The turbulent flux of kinetic energy is usually modelled using a gradient hypothesis in non-reacting turbulent flows (Jones and Launder 1973; Durbin and Reif 2010; Wilcox 2002) which may be expressed as:

$$ \overline{{\rho u_{1}^{^{\prime\prime}} k}} = - C_{T1} \overline{\rho }\frac{{\tilde{k}^{2} }}{{\tilde{\epsilon}}}\frac{{\partial \tilde{k}}}{{\partial x_{1} }} $$

(37)

where $C_{T1} = C_{\mu } = 0.09.$ This approach may be used to model the non-reacting contribution, represented by the first and second terms on the right-hand side of Eq. (36), which can then be written as:

$$ \overline{\rho }\left( {1 - \tilde{c}} \right)\overline{{u_{1R}^{^{\prime}} k_{R} }} + \overline{\rho }\tilde{c}\overline{{u_{1P}^{^{\prime}} k_{P} }} = - C_{T2} \overline{\rho }\frac{{\tilde{k}^{2} }}{{\tilde{\epsilon}}}\frac{{\partial \tilde{k}}}{{\partial x_{1} }} $$

(38)

The slip velocity in Eq. (36) may be modelled using the turbulent scalar flux $\overline{{\rho u_{1}^{^{\prime\prime}} c^{\prime\prime}}}$ as given by the BML analysis (Bray et al. 1985) and this may be expressed as:

$$ \overline{u}_{1P} - \overline{u}_{1R} \approx \frac{{\overline{{\rho u_{1}^{^{\prime\prime}} c^{\prime\prime}}} }}{{\overline{{\rho c^{^{\prime\prime}2} }} }} $$

(39)

whereas the difference in mean turbulent kinetic energy between the products and reactants may be modelled as (Chakraborty et al. 2011a):

$$ \overline{k}_{P} - \overline{k}_{R} = - C_{T3} m_{i} \frac{{\partial \tilde{k}}}{{\partial x_{i} }}\frac{1}{{\Sigma }} $$

(40)

where $m_{i} = - \left( {\partial \tilde{c}/\partial x_{i} } \right)/\left| {\nabla \tilde{c}} \right|$ is the resolved flame normal vector component in the i^th direction and $C_{T3}$ is a model parameter. Combining Eqs. (36) and (38)–(40) gives a model for the turbulent flux of kinetic energy, which can be written as (Chakraborty et al. 2011a):

$$ \overline{{\rho u_{1}^{^{\prime\prime}} k}} = - \overline{\rho }C_{T2} \frac{{\tilde{k}^{2} }}{{\tilde{\epsilon}}}\frac{{\partial \tilde{k}}}{{\partial x_{1} }} - C_{T3} \overline{{\rho u_{1}^{^{\prime\prime}} c^{\prime\prime}}} m_{i} \frac{{\partial \tilde{k}}}{{\partial x_{i} }}\frac{1}{{\Sigma }} + \frac{1}{{2\left( {\overline{{\rho c^{^{\prime\prime}2} }} } \right)^{2} }}\left( {\overline{{\rho u_{1}^{^{\prime\prime}} c^{\prime\prime}}} } \right)^{3} \left( {1 - 2\tilde{c}} \right) $$

(41)

It should be noted that when the first two terms on the right-hand side of Eq. (41) become negligible with respect to the last term, the transition of $\overline{{\rho u_{1}^{^{\prime\prime}} k}} /\rho_{0} S_{L}^{3}$ from positive to negative occurs at $\tilde{c} = 0.5$, which may not be true for flames in the thin reaction zones regime since the pdf of $c$ cannot be approximated by a bimodal distribution as used in the derivation of Eq. (36). Chakraborty et al. (2011a) modified Eq. (41) for applicability within the thin reaction zones regime and this modified model may be given as:

$$ \overline{{\rho u_{1}^{^{\prime\prime}} k}} = - \overline{\rho }C_{T2} \frac{{\tilde{k}^{2} }}{{\tilde{\epsilon}}}\frac{{\partial \tilde{k}}}{{\partial x_{1} }} - C_{T3} \overline{{\rho u_{1}^{^{\prime\prime}} c^{\prime\prime}}} m_{i} \frac{{\partial \tilde{k}}}{{\partial x_{i} }}\frac{1}{{\Sigma }} + \frac{1}{{2\left( {\overline{{\rho c^{^{\prime\prime}2} }} } \right)^{2} }}\left( {\overline{{\rho u_{1}^{^{\prime\prime}} c^{\prime\prime}}} } \right)^{3} \left( {1 - 2\tilde{c}g^{s} } \right) $$

(42)

where $g = \overline{{\rho c^{\prime\prime}c^{\prime\prime}}} /\overline{\rho }\tilde{c}\left( {1 - \tilde{c}} \right)$ is the segregation factor and $s$ is a model constant. In the present study, the effects of body forces have been incorporated into the model parameters $C_{T2}$ and $C_{T3}$, which are given as $C_{T2} = 0.45 + 0.15$ erf $\left( {1.624g^{*} - 0.31} \right)$ and $C_{T3} = 0.305 + 0.295$ erf $\left( {1.95g^{*} - 0.684} \right)$, and the value of the model constant $s$ has been taken to be $0.2$. The predictions given by Eqs. (37) and (42) are shown in Fig. 11 alongside the $\overline{{\rho u_{1}^{^{\prime\prime}} k}} /\rho_{0} S_{L}^{3}$ term extracted from the DNS data. It can be seen that Eq. (37) fails to accurately predict the qualitative or quantitative behaviour of $\overline{{\rho u_{1}^{^{\prime\prime}} k}} ,$ and moreover predicts the opposite sign of $\overline{{\rho u_{1}^{^{\prime\prime}} k}}$, which is indicative of the occurrence of counter-gradient transport of turbulent kinetic energy. The model given by Eq. (42) satisfactorily predicts the variation of $\overline{{\rho u_{1}^{^{\prime\prime}} k}}$ across the flame brush, although underprediction of the magnitude of $\overline{{\rho u_{1}^{^{\prime\prime}} k}}$ is observed for negative $g^{*}$ values at high initial $u^{\prime}/S_{L}$ cases (i.e., Set-C and Set-D).

Fig. 11

Variations of $\overline{{\rho u_{1}^{^{\prime\prime}} k}} /\rho_{0} S_{L}^{3}$ with $\tilde{c}$ along with the predictions of the models given by Eqs. (37) and (42) for all cases considered

It is worth noting the model predictions are shown in Fig. 11 for $\tilde{\epsilon}$, $\overline{{\rho c^{^{\prime\prime}2} }}$ and ${\Sigma }$ extracted from DNS data in the spirit of a priori analysis but in a RANS simulation the accuracy of the closures of these unclosed terms will also affect the fidelity of the modelling of $\overline{{\rho u_{1}^{^{\prime\prime}} k}}$.

4.9 Model Parameterisation and Error Assessment

It can be seen from Figs. 9, 10 and 11 that the magnitudes of $T_{V}$ (in modelling of $T_{4}$), the fluctuating pressure gradient term $- \overline{{u_{i}^{^{\prime\prime}} \partial p^{\prime}/\partial x_{i} }}$ (in modelling of $T_{5}$) and turbulent flux of kinetic energy $\overline{{\rho u_{1}^{^{\prime\prime}} k}}$ varies with $g^{*}$, which necessitates the incorporation of the effects of $g^{*}$ into the models for closure of these terms. Accordingly, the model constants $C_{v}$, $C_{1N}$, $C_{2N}$, $C_{T2}$ and $C_{T3}$ are expressed as functions of $g^{*}$. Although a certain degree of empiricism is involved in the expressions for these model parameters, they satisfy all asymptotic requirements in terms of $g^{*} .$ Moreover, most turbulent combustion models involve a degree of empiricism and even if the expressions are considered to be purely empirical parameterisations of the DNS data, they remain as valuable as any new experimental parameterisation which describes a physical phenomenon that is yet to be analysed in detail. For the purpose of quantification of the error level of the models considered here the average $R^{2}$ values and the average error quantity $E$ (both $R^{2}$ and $E$ are averaged over different $u^{\prime}/S_{L}$ for each $g^{*}$) between the model predictions and the corresponding DNS results for the unclosed terms are listed in “Appendix”, where $E$ is defined as:

$$ E = \sqrt {\frac{{\mathop \sum \nolimits_{i = 1}^{N} \left( {\langle Q_{i}\rangle_{{\tilde{c}}}^{DNS} - \langle Q_{i}\rangle_{{\tilde{c}}}^{Model} } \right)^{2} }}{N}} \times \frac{1}{{max\left| {\langle Q_{i}\rangle_{{\tilde{c}}}^{DNS} } \right|}} \times 100 \% $$

(43)

where $N = 50$ is the number of uniformly spaced bins in the range $0.005 \le \tilde{c} \le 0.995$ for the cases considered in the present study, and $\langle Q_{i}\rangle_{{\tilde{c}}}^{DNS}$ and $\langle Q_{i}\rangle_{{\tilde{c}}}^{Model}$ are the mean values of a general quantity $Q$ from DNS and model expressions condition upon i^th bin of $\tilde{c}$. It can be seen from the “Appendix” that all the recommended models have sufficiently high $R^{2}$ and low $E$ values providing the confidence that the identified model expressions have the potential to be successful in an actual RANS simulation. Moreover, in combination with the figures showing the comparison between the model predictions and the corresponding DNS data conditioned upon $\tilde{c}$ shows the extent to which the model captures the qualitative behaviours of these quantities since it ensures that equal importance is given to all $\tilde{c}$ values.

In order to demonstrate that the performance of the newly developed modelled $\tilde{k}$ equation, RANS simulations of statistically planar flames corresponding to Set-A and Set-D are exemplarily conducted for $g^{*} = - 3.12$ and 3.12 using a well-known finite-volume open-source code Code_Saturne (Archambeau et al. 2004) (see: http://www.code-saturne.org). This code has previously been used in several RANS investigations of premixed turbulent combustion (Ahmed and Prosser 2016, 2018; Dong et al. 2012). A one-dimensional domain of $280.8{\updelta }_{{\text{Z}}}$, discretised by equidistant 1600 grid points (which ensures grid convergence for turbulence as well as the flame), is used for these RANS simulations where a second order central difference scheme is used for the purpose of discretisation. For the RANS simulations, the mean reaction rate is modelled as (Boger et al. 1998):

$$ \overline{{\dot{\omega }_{c} }} \approx \rho_{0} S_{L} {\Sigma } $$

(44)

where ${\Sigma }$ is modelled using an existing algebraic closure (Poinsot and Veynante 2021; Bray et al. 1984; Bray and Libby 1986):

$$ {\Sigma } = \frac{{g_{{\Sigma }} \tilde{c}\left( {1 - \tilde{c}} \right)}}{{\sigma_{y} L_{y} }}\frac{{\left( {1 + \tau } \right)}}{{\left( {1 + \tau \tilde{c}} \right)^{2} }} $$

(45)

where $g_{{\Sigma }} = 1.50$, and the parameters $\sigma_{y}$ and $L_{y}$ are taken to be 0.5 and $0.37[\left( {2\tilde{k}/3} \right)^{3/2} /\tilde{\epsilon}](S_{L} /\sqrt {2\tilde{k}/3} )^{n}$ with $n = 1.0$, following previous suggestions (Poinsot and Veynante 2021; Bray et al. 1984; Bray and Libby 1986). It is worth noting that the applicability of Eqs. (44) and (45) along with its model parameters has not been assessed for the premixed turbulent flames under the influence body force and they are merely used in this part of the analysis to enable the RANS simulations to demonstrate the influence of $\tilde{k}$ modelling. The simulations of statistically planar flames corresponding to Set-A and Set-D for the extreme cases with $g^{*} = - 3.12$ and 3.12 have been conducted for the standard $k -\epsilon$ model where the modelled $\tilde{k}$ and $\tilde{\epsilon}$ equations take the following form (Durbin and Reif 2010; Wilcox 2002):

$$ \frac{{\partial \left( {\overline{\rho }\tilde{k}} \right)}}{\partial t} + \frac{{\partial \left( {\overline{\rho }\tilde{u}_{j} \tilde{k}} \right)}}{{\partial x_{j} }}{ } = + \frac{\partial }{{\partial x_{j} }}\left[ {\left( {\mu + \frac{{\mu_{t} }}{{\sigma_{k} }}} \right)\frac{{\partial \tilde{k}}}{{\partial x_{j} }}} \right]\underbrace {{ - \overline{{\rho u_{i}^{{^{\prime\prime}}} u_{j}^{{^{\prime\prime}}} }} \frac{{\partial \tilde{u}_{i} }}{{\partial x_{j} }}}}_{{P_{\varepsilon } }} - \overline{\rho }\tilde{\epsilon} $$

(46)

$$ \frac{{\partial \left( {\overline{\rho }\tilde{\epsilon}} \right)}}{\partial t} + \frac{{\partial \left( {\overline{\rho }\tilde{u}_{j} \tilde{\epsilon}} \right)}}{{\partial x_{j} }}{ } = \frac{\partial }{{\partial x_{j} }}\left[ {\left( {\mu + \frac{{\mu_{t} }}{{\sigma_{\varepsilon } }}} \right)\frac{{\partial \tilde{\epsilon}}}{{\partial x_{j} }}} \right] + C_{\varepsilon 1} \frac{{\tilde{\epsilon}}}{{\tilde{k}}}P_{\varepsilon } - \overline{\rho }{ }C_{\varepsilon 2} \frac{{\tilde{\epsilon}^{2} }}{{\tilde{k}}} $$

(47)

Here, $\sigma_{k} = 1.0, \sigma_{\varepsilon } = 1.3, C_{\varepsilon 1} = 1.44$ and $C_{\varepsilon 2} = 1.92$ are the model parameters and in the context of the standard $k - \epsilon$ model the eddy viscosity $\mu_{t} = \overline{\rho }C_{\mu } \tilde{k}^{2} /\tilde{\epsilon}$ is considered with $C_{\mu } = 0.09$. Moreover, $\overline{{\rho u_{i}^{^{\prime\prime}} u_{j}^{^{\prime\prime}} }}$ in $P_{\varepsilon }$ is modelled according to Eq. (4). The turbulent scalar flux is modelled according to Veynante et al. (1997) as:

$$ \overline{{\rho u_{{\text{i}}} ^{\prime\prime}c^{\prime\prime}}} = - \frac{{\mu_{t} }}{{\sigma_{t} }}\frac{{\partial \tilde{c}}}{{\partial x_{i} }} - M_{i} b_{s} \overline{\rho }\tau S_{L} \tilde{c}\left( {1 - \tilde{c}} \right) $$

(48)

Here, the turbulent Schmidt number $\sigma_{t}$ is taken to be 1.0 (i.e., $\sigma_{t} = 1.0$), $M_{i} = - \left( {\partial \tilde{c}/\partial x_{i} } \right)/\left| {\nabla \tilde{c}} \right|$ is the i^th component of the resolved flame normal and $b_{s}$ is an appropriate constant which is estimated from DNS data and considered for the RANS simulation. It is worth noting that optimising Eq. (48) for the buoyancy effects is beyond the scope of the current analysis and thus no parameterisation of $b_{s}$ in terms of $g^{*}$ has been attempted here.

The same simulations have been conducted for the modelled $\tilde{k}$ equation according to the current a priori analysis, which takes the following form:

$$ \begin{aligned} & \frac{{\partial \left( {\overline{\rho }\tilde{k}} \right)}}{\partial t} + \frac{{\partial \left( {\overline{\rho }\tilde{u}_{j} \tilde{k}} \right)}}{{\partial x_{j} }} \\ { } & \quad = - \overline{{\rho u_{i}^{{^{\prime\prime}}} u_{j}^{{^{\prime\prime}}} }} \frac{{\partial \tilde{u}_{i} }}{{\partial x_{j} }} - \frac{\tau }{{\rho_{0} }}\overline{{\rho u_{i}^{^{\prime\prime}} c^{\prime\prime}}} \frac{{\partial \overline{p}}}{{\partial x_{i} }} + \frac{\partial }{{\partial x_{j} }}\left( {\mu \frac{{\partial \tilde{k}}}{{\partial x_{j} }}} \right) - \overline{\rho }\tilde{\epsilon} \\ & \quad \quad - C_{\nu } \frac{{\tilde{\epsilon}}}{{\tilde{k}}}\mu_{0} \tau S_{L} \frac{{\tilde{c}^{p} \left( {1 - \tilde{c}} \right)^{q} }}{{\delta_{th} }} + \frac{1}{2}\rho_{0} S_{L}^{3} {\Sigma }\tau^{3} \frac{{\tilde{c}\left( {1 - \tilde{c}} \right)}}{{\left( {1 + \tau } \right)}} \\ & \quad \quad + \frac{\partial }{{\partial x_{i} }}\left( {\overline{\rho }C_{T2} \frac{{\tilde{k}^{2} }}{{\tilde{\epsilon}}}\frac{{\partial \tilde{k}}}{{\partial x_{i} }}} \right) + \frac{\partial }{{\partial x_{i} }}\left( {C_{T3} \overline{{\rho u_{i}^{^{\prime\prime}} c^{\prime\prime}}} m_{j} \frac{{\partial \tilde{k}}}{{\partial x_{j} }}\frac{1}{{\Sigma }}} \right) \\ & \quad \quad - \frac{\partial }{{\partial x_{i} }}\left\{ {\frac{1}{{2\left( {g\overline{\rho }\tilde{c}\left( {1 - \tilde{c}} \right)} \right)^{2} }}\left( {\overline{{\rho u_{i}^{^{\prime\prime}} c^{\prime\prime}}} } \right)\left( {\overline{{\rho u_{j}^{^{\prime\prime}} c^{\prime\prime}}} } \right)\left( {\overline{{\rho u_{j}^{^{\prime\prime}} c^{\prime\prime}}} } \right)\left( {1 - 2\tilde{c}g^{s} } \right)} \right\} \\ \end{aligned} $$

(49)

In Eq. (49), $T_{2} ,T_{3} ,T_{V} , - \overline{{u_{i}^{^{\prime\prime}} \partial p^{\prime}/\partial x_{i} }}$ and $\overline{{\rho u_{i}^{^{\prime\prime}} k}}$ are modelled by Eqs. (15), (23), (28), (32) and (42), respectively. Equation (35) for $- \overline{{u_{i}^{^{\prime\prime}} \partial p^{\prime}/\partial x_{i} }}$ has not been attempted here because it needs additional closures of $\overline{p}_{R}$ and $\overline{p}_{P}$. For solving Eq. (49), $\overline{{\rho u_{i}^{{^{\prime\prime}}} u_{j}^{{^{\prime\prime}}} }}$ is expressed using Eq. (12), and $\overline{{\dot{\omega }_{c} }} ,{{ \Sigma }},$ and $\overline{{\rho u_{i} ^{\prime\prime}c^{\prime\prime}}}$ are modelled by Eqs. (44), (45) and (48) respectively. A modified form of the $\tilde{\epsilon}$ transport equation is used with the new $\tilde{k}$ transport equation, where the mean pressure gradient effects are accounted for by adding an extra term on the right-hand side of Eq. (47) as $- C_{1} \left( {\tau /\rho_{0} } \right)\overline{{\rho u_{i}^{^{\prime\prime}} c^{\prime\prime}}} \left( {\partial \overline{p}/\partial x_{i} } \right) \left( {\tilde{\epsilon}/\tilde{k}} \right)$ following Kolla and Swaminathan (2010). Moreover, the representative segregation factor $g = \overline{{\rho c^{\prime\prime}c^{\prime\prime}}} /\overline{\rho }\tilde{c}\left( {1 - \tilde{c}} \right)$ value is extracted from DNS data instead of solving a transport equation of $\overline{{\rho c^{\prime\prime}c^{\prime\prime}}} /\overline{\rho }$ because the effects of body force on the scalar variance is a separate study in its own merit and is beyond the scope of the current analysis. The formulation with the newly modelled $\tilde{k}$ equation (i.e., Eq. (49)) will henceforth be referred to as the NEW formulation in this paper.

The values of the normalised turbulent kinetic energy $\tilde{k}/\tilde{k}_{0}$ for $g^{*} = - 3.12$ and 3.12 in the case of Set-A and Set-D, as obtained from DNS data, standard $k -\epsilon$ model and NEW formulation are shown in Fig. 12, where $\tilde{k}_{0}$ is the value of turbulent kinetic energy at $\tilde{c} = 0.005$. It can be seen from Fig. 12 that the standard $k -\epsilon$ model does not capture the qualitative behaviour of $\tilde{k}/\tilde{k}_{0}$ within the flame brush, as extracted from DNS data. By contrast, the NEW formulation shows local augmentation of $\tilde{k}/\tilde{k}_{0}$ within the flame brush, which is qualitatively similar to that observed from DNS data. Although there are quantitative differences between the $\tilde{k}/\tilde{k}_{0}$ values obtained from DNS data and the predictions of the NEW formulation, this difference is perhaps not surprising because the performance of Eq. (49) is also dependent on the accuracy of the modelling of $\overline{{\dot{\omega }_{c} }} ,{{ \Sigma }}, \overline{{\rho u_{i} ^{\prime\prime}c^{\prime\prime}}}$ and $\tilde{\epsilon}$, which are not optimised for body force effects and the associated inaccuracies lead to the quantitative differences between the $\tilde{k}/\tilde{k}_{0}$ values between DNS data and the predictions of the NEW formulation. In this context it is important to note that implementation of the proposed model expressions into a RANS code requires particular care since the models proposed in this work are in turn dependent on models for dissipation rate of turbulent kinetic energy, turbulent scalar flux, scalar variance, FSD and conditional mean pressures. Modelling of buoyancy effects on the closures of the dissipation rate of turbulent kinetic energy, scalar variance and conditional mean pressures are yet to be reported in the literature and these aspects are beyond the scope of the present analysis. Therefore, implementation of the proposed models in a RANS code is not a straightforward task. Moreover, it is also worth noting that numerical schemes used in the RANS code are also expected to influence the simulation predictions and hence it is not straightforward to assess the performance of the models for the terms in the turbulent kinetic energy transport equation in isolation. It should also be noted that any model that does not perform well in a priori analyses is unlikely to perform well in actual RANS simulations, whereas models that have been validated based on a priori analyses need to be implemented in actual RANS to assess their performance while interacting with models for other unclosed terms, turbulence models and numerical schemes. It is also worth noting that a priori modelling approach was adopted by several authors in the past (Rutland and Cant 1994; Zhang and Rutland 1995; Nishiki et al. 2002; Chakraborty et al. 2011a, b; Lai et al. 2017; Ahmed et al. 2019; Domingo and Bray 2000), and the same approach has been undertaken in this paper. Moreover, the validation of the models after implementation into RANS code would require an experimental database that provides information on the different terms in the turbulent kinetic energy transport equation.

Fig. 12

Values of the normalised turbulent kinetic energy $\tilde{k}/\tilde{k}_{0}$ for $g^{*} = - 3.12$ and $3.12$ in the case of Set-A and Set-D, as obtained from DNS data, standard $k - \epsilon$ model and NEW formulation

5 Conclusions

The effects of body force on the transport of turbulent kinetic energy have been evaluated using a DNS database of three-dimensional statistically planar flames with unity Lewis number for different sets of initial $u^{\prime}/S_{L}$ values under decaying turbulence. 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. It has been found that the turbulent kinetic energy and its dissipation rate are affected by both the magnitude and direction of the body force, especially for cases with low initial $u^{\prime}/S_{L}$ such as Set-A and Set-B. For higher initial $u^{\prime}/S_{L}$ values, such as Set-C and Set-D, the effects of body force tend to be weak due the relative strengthening of the inertial effects. The effects of body forces on the statistical behaviour of the various terms of the turbulent kinetic energy transport equation have been analysed and it is found that the pressure dilatation term $T_{3}$ and the mean pressure gradient term $T_{2}$ are the major source terms, whereas the viscous dissipation term $T_{4}$ and the mean velocity gradient term $T_{1}$ act as sink terms. The contributions of the pressure transport term $T_{5}$ and turbulent transport term $T_{6}$ are found to be significant only for cases with stable density stratification, i.e., for negative $g^{*}$ values.

The performance of existing closure models for the various terms in the transport equation have also been evaluated based on a priori analysis and appropriate models for the mean velocity gradient, mean pressure gradient, and pressure dilatation terms have been identified. The existing models for the viscous dissipation and molecular diffusion term $T_{4}$, the pressure transport term $T_{5}$ and the turbulent transport term $T_{6}$ have been modified to account for the effects of body forces and it has been shown that the modified models capture the variations of the respective terms satisfactorily for all cases considered. It is worthwhile to note that for the models considered in this work, the terms such as turbulent scalar flux $\overline{{\rho u_{1}^{^{\prime\prime}} c^{\prime\prime}}}$, variance of reaction progress variable $\overline{{\rho c^{^{\prime\prime}2} }}$, flame surface density ${\Sigma }$, and mean chemical reaction rate $\overline{{\dot{w}}}$ have been taken from the DNS data and thus the performance of the respective models in the turbulent kinetic energy transport equation in an actual RANS simulation will therefore also depend on the closures of these unclosed terms. Moreover, it is important to note that this study has been conducted with simple chemistry in the interest of a detailed parametric analysis and future analyses need to include the effects of detailed chemistry and transport in order to validate the modified models proposed in this analysis. A preliminary implementation of the modified $\tilde{k}$ transport equation in RANS simulations shows a significant improvement in the $\tilde{k}$ predictions within the mean flame brush when compared with the DNS data. A more detailed a posteriori assessment of the models based on RANS simulations will also be required to validate the practical applicability of these models.

Appendix: Average \({R}^{2}\) and \(E\) Values Between DNS Data and Model Predictions

Mean velocity gradient $T_{1}$	$g^{*}$	Model—Eq. (5)		Model—Eq. (12)+
		$R^{2}$	$E$	$R^{2}$	$E$
	− 3.12	0.97	40.3	0.97	7.9
	− 1.56	0.97	36.5	0.99	9.1
	0	0.95	29.2	0.99	13.8
	1.56	0.89	22.1	0.97	16.1
	3.12	0.86	25.1	0.88	19.1

Pressure dilatation term $T_{3}$	$g^{*}$	PDZ		PDN+
		$R^{2}$	$E$	$R^{2}$	$E$
	− 3.12	0.4	59.5	0.64	55.9
	− 1.56	0.35	43.3	0.74	35.0
	0	0.33	30.3	0.74	16.7
	1.56	0.31	29.5	0.52	21.1
	3.12	0.44	26.4	0.45	23.4

Viscous dissipation and molecular diffusion term $T_{4}$	$g^{*}$	VN		VCKC+
		$R^{2}$	$E$	$R^{2}$	$E$
	− 3.12	0.72	31.3	0.88	24.5
	− 1.56	0.8	34.7	0.89	18.9
	0	0.5	31.1	0.77	19.2
	1.56	0.44	29.4	0.82	23.5
	3.12	0.4	27.8	0.75	25.7

Pressure transport term $T_{5}$	$g^{*}$	PTS#		PCKC+
		$R^{2}$	$E$	$R^{2}$	$E$
	− 3.12	0.33	53.8	0.59	39.2
	− 1.56	0.36	36.3	0.69	23.2
	0	0.45	32.0	0.84	17.7
	1.56	0.62	42.4	0.84	38.3
	3.12	0.57	46.6	0.73	61.4

Turbulent transport term $T_{6}$	$g^{*}$	Model—Eq. (37)		Model—Eq. (42)+
		$R^{2}$	$E$	$R^{2}$	$E$
	− 3.12	0.92	56.9	0.90	27.4
	− 1.56	0.91	56.2	0.91	24.8
	0	0.79	57.2	0.84	26.5
	1.56	0.31	47.2	0.33	42.1
	3.12	0.13	27.6	0.27	26.1

1. ⁺Recommended models
2. ^#The comparison is shown with the PTS model since it gives the best predictions among the other models