A Direct Numerical Simulation Assessment of Turbulent Burning Velocity Parametrizations for Non-Unity Lewis Numbers

Abstract

The predictions of turbulent burning velocity parameterizations for non-unity Lewis number flames have been assessed based on a single-step chemistry Direct Numerical Simulation (DNS) database of premixed Bunsen flames for different values of characteristic Lewis numbers ranging from 0.34 to 1.2. It has been found that the definition of the turbulent burning velocity is strongly dependent on the choice of projected flame brush area in the Bunsen burner configuration. The highest values of normalized turbulent burning velocity are obtained when the projected flame brush area is evaluated using the area of the isosurface of the Reynolds averaged reaction progress variable of 0.1 out of different options, namely the Favre averaged and Reynolds averaged isosurfaces of reaction progress variable of 0.5 and integral of the gradient of Favre and Reynolds averaged reaction progress variable. Because of the axisymmetric nature of the mean flame brush, the normalized turbulent burning velocity has been found to decrease as the burned gas side is approached, due to an increase in flame brush area with increasing radius. Most models for turbulent burning velocity provide comparable, reasonably accurate predictions for the unity Lewis number case when the projected flame brush area is evaluated using the isosurface of the Reynolds averaged reaction progress variable of 0.1. However, most of these parameterizations underpredict turbulent burning velocity values for Lewis numbers smaller than unity. A scaling relation has been utilized to extend these parameterizations for non-unity Lewis numbers. These revised parameterizations have been shown to be more successful than the original model expressions. These modified expressions also exhibit small values of L₂-norm of the relative error with respect to experimental data from literature for different Lewis numbers, higher turbulence intensity and thermodynamic pressure levels.

1. Introduction

Due to threats of climate change, countries and industries are searching for cleaner and efficient ways of generating energy. Although renewable forms of energy exist, these have their own drawbacks, such as the intermittent generation of energy for wind and solar. Unless and until more efficient methods to store such energy are devised, it is very likely that combustion will remain the principal mode of energy production. Premixed combustion, in which reactants are homogenously mixed, is a good way to reduce the emission of pollutants as it limits the chances of incomplete combustion. Moreover, it is easier to control NOx emissions by optimizing between the peak temperature achieved and power produced. This can be done by either controlling the temperature or the composition of the reactants. The emission of greenhouse gases such as CO₂ can also be reduced by using a premixed combustion of fuels such as hydrogen, ammonia or syngas. Net zero targets by governments can be achieved by the application of low/zero-carbon fuels such as biofuel, hydrogen and ammonia, as proposed by contemporary scenario plans [1].

The application of pure hydrogen enables combustion with the complete elimination of carbon-based greenhouse gas emission [1]. However, the thermochemistry of hydrogen combustion is significantly different from that of hydrocarbon fuels [2], and the presence of lighter chemical species induces significant effects of differential diffusion. This, in turn, impacts upon practical operational burner characteristics such as (i) variations in flame stability maps and (ii) susceptibility to premature combustion (knocking) and flashback. As most existing combustion devices are designed to operate based on hydrocarbon fuels, the switch from hydrocarbon to hydrogen fuel will happen gradually and, in the interim period, hydrogen is expected to be blended in fixed proportions with hydrocarbon fuel to enable a transitional reduction in greenhouse gas emissions. The presence of H₂ in the fuel induces a significant amount of differential diffusion of heat and species due to the non-unity Lewis number. Thus, the non-unity Lewis number effects cannot be ignored in the premixed combustion of High-Hydrogen-Content (HHC) fuels.

It has been shown in several previous analytical [3,4,5,6], experimental [7,8,9,10] and computational [11,12,13,14,15,16,17,18] studies that the turbulent burning velocity in premixed turbulent flames increases with decreasing characteristic Lewis number. As turbulent burning velocity $S_T$ is one of the most important quantities in the modelling of turbulent premixed combustion [7,8,9,19,20,21,22], it is important to consider the modelling of $S_T$ in response to the variations of the characteristic Lewis number. However, a universal scaling law for the turbulent flame speed is yet to be reported [23] and is still a subject of current research. Bray and Cant [24] pointed out that the $S_T/S_L$ is proportional to the product of the flame brush thickness and the maximum of the flame surface density multiplied by the stretch factor $I_0$. These quantities in turn depend on the turbulence structure and the thermophysical properties of the fuel. Often a quadratic scaling of turbulent burning velocity in terms of turbulence intensity is assumed [23]. However, upon onset of the hydro-dynamic instability, a sub-unity scaling exponent for turbulence intensity is reported [23]. Similarly, Kobayashi et al. [25] suggested a correlation of turbulent burning velocity $S_T$ in terms of pressure and turbulence intensity. Driscoll [26] argued that for an ideal geometry independent turbulent flame, $S_T/S_L$ would be expected to be a function of the normalized turbulence intensity and integral length scale as well as the turbulent Markstein number. Further, Driscoll [26] pointed out that for real flames, the wrinkling has a memory effect which will make the functional relationship more complicated. A variety of experimental and numerical data and empirical formulations are discussed in a review article [26]. A rather complex formulation based on measured turbulent burning velocities has been suggested by Filatyev et al. [27] which also provides a fit to the nonlinear bending observed for high turbulence intensities. An important ingredient for this formulation was to include the mean flow velocity and the burner width. The turbulent bending effect has been discussed by several authors, and the reader is referred to [28,29] for further information in this regard. Both the Markstein length and the thermophysical properties, related to the Lewis number of the flame, play an important role. In this respect, it is worth noting that most existing parameterizations of turbulent burning velocity [7,8,9,19,20,21,22] have been proposed for unity Lewis number conditions where the differential diffusion of heat and mass was ignored.

As reactants are composed of different gases with different Lewis numbers, it is difficult to estimate a characteristic mixture Lewis number. However, different methods of estimation of the characteristic Lewis number have been proposed, such as equating the Lewis number of the reactant with the least concentration [30,31], calculating it from measurements made of the heat release rate [8], by numerically estimating it from binary diffusion theory for mixtures [32] or by individually calculating the thermal conductivity and mass diffusivity based on the mole fraction of the constituent species [33]. Bechtold and Matalon [31] have suggested that the mixture Lewis number can be calculated as: $Le=1.0+\left[\left(L{e}_F-1\right)+\left(L{e}_O-1\right){\mathsf{\Lambda }}_{Le}\right]/\left(1+{\mathsf{\Lambda }}_{Le}\right)$ where ${\mathrm{\Lambda }}_{Le}=1+\beta (\mathrm{\Phi }-1)$ with $\mathsf{\Phi }=\varphi H\left(\varphi -1\right)+\left(1/\varphi \right)\left(1-H\left(\varphi -1\right)\right)$ for $\varphi , H\left(x\right)$ and $\beta$ being the equivalence ratio, Heaviside function and Zel’dovich number, respectively. Subscripts F and O are used for fuel and oxidizer, respectively. The present analysis employs the concept of characteristic Lewis number to analyse the effects of $Le$ on the turbulent burning velocity using simple chemistry DNS data of Bunsen burner flames.

To date, most of the turbulent burning velocity parameterizations [7,8,9,19,20,21,22] have been proposed for unity Lewis number conditions and therefore disregard the effects of differential diffusion of heat and species. The performance and applicability of these parameterizations in the flames with the non-unity characteristic Lewis number are yet to be assessed in detail. The present analysis addresses this gap in the literature by considering a DNS database of turbulent premixed Bunsen burner flames with a non-unity characteristic Lewis number to assess the performance of existing parameterizations [7,8,9,19,20,21,22]. The turbulent burning velocity $S_T$ can be defined in terms of the volume-integral in the following manner [34]:

$$S_T=\frac{1}{{\rho }_u{A}_L}{{\displaystyle \int }}_V{\dot{\omega }}_{c}dV$$

(1)

where ${\dot{\omega }}_c$ is the reaction rate of reaction progress variable $c$, ${\rho }_u$ is the unburned gas density and $A_L$ is the projected flame brush area in the direction of mean flame propagation. The evaluation and definition of $A_L$ gives rise to uncertainties to the evaluation of the turbulent burning velocity $S_T$, and the impact of $A_L$ evaluation on $S_T$ will be assessed in this analysis by employing different methodologies to extract the projected flame surface area $A_L$. In this respect, the main objectives of this analysis are:

(a)

To assess the performances of the existing parameterizations of turbulent burning velocity $S_T$ for turbulent premixed flames with characteristic Lewis numbers significantly different from unity.

(b)

To illustrate the impact of the projected flame brush surface area $A_L$ evaluation on turbulent burning velocity $S_T$ for turbulent Bunsen burner flames with different characteristic Lewis numbers.

The rest of the paper will be organized as follows. The mathematical background and numerical implementation pertaining to the current analysis are presented in Section 2 and Section 3 of this paper. The results will be presented in Section 4 of this paper, and finally, the main findings are summarized and conclusions are drawn in Section 5.

2. Mathematical Background

The present analysis deals with the integral quantities such as volume-integrated burning rate and flame surface area, which can be captured reasonably accurately by single-step chemistry. It was shown elsewhere [35,36] that the flame propagation statistics extracted from single-step chemistry DNS is qualitatively similar to that obtained for detailed chemistry DNS and the quantitative differences are of the same order of uncertainty associated with different definitions of reaction progress variable $c$. Thus, a single-step chemistry has been considered in this analysis in the interest of computational economy and to identify the effects of the characteristic Lewis number in isolation. In the context of simple chemistry, the reactive scalar field can be described with the help of reaction progress variable $c$ and non-dimensional temperature $\theta$, which are defined as:

$$c=\frac{Y_{\alpha }-Y_{\alpha u}}{Y_{\alpha b}-Y_{\alpha u}};\theta =\frac{T-T_u}{T_{\alpha d}-T_u}$$

(2)

where $Y_{\alpha }$ is the mass fraction of species $\alpha$ based on which reaction progress variable is defined, the subscripts u and b are used to refer to the values in the unburned gas and burned gas, respectively, and $T_{ad}$ is the adiabatic flame temperature. Accordingly, the reaction rate of the reaction progress variable ${\dot{\omega }}_c$ is defined as:

$${\dot{\omega }}_c=\frac{{\dot{\omega }}_{\alpha }}{Y_{\alpha b}-Y_{\alpha u}}$$

(3)

In the context of Reynolds Averaged Navier–Stokes (RANS) simulations, the mean reaction rate can be modelled in the following manner for high values of Damköhler number (i.e., $Da\gg 1$) [37]:

$${\overline{\dot{\omega }}}_c=2\overline{\rho }{\tilde{\epsilon }}_c/\left(2c_m-1\right)$$

(4)

Here, $\overline{q}$, $\tilde{q}=\overline{\rho q}/\overline{\rho }$ and $q^{″}=q-\tilde{q}$ are the Reynolds averaged, Favre averaged and Favre fluctuation of a general quantity $q$. In Equation (4), ${\tilde{\epsilon }}_c=\overline{\rho D\nabla c^{″}·\nabla c^{″}}/\overline{\rho }$ is the scalar dissipation rate based on Favre fluctuations of reaction progress variable, $D$ is the reaction progress variable diffusivity and $c_m={{\displaystyle \int }}_0^1{[{\dot{\omega }}_{c}c]}_{Lam}f\left(c\right)dc/{{\displaystyle \int }}_0^1{[{\dot{\omega }}_c]}_{Lam}f\left(c\right)dc$ is a thermochemical parameter with $f\left(c\right)$ being the burning mode probability density function, and the subscript ‘Lam’ refers to the unstretched laminar flame quantities. It was discussed by Bray [37] that the assumption of any continuous function for $f\left(c\right)$ is sufficient for the evaluation of $c_m$. Although Equation (4) was originally proposed for $Da\gg 1$ flames, it was demonstrated subsequently by Chakraborty and Cant [17] that this relation holds in an order of magnitude sense for the flames with $Da<1$. However, the scalar dissipation rate ${\tilde{\epsilon }}_c$ is an unclosed quantity and based on the leading order balance of the terms of the transport equation of ${\tilde{\epsilon }}_c$ under $Da\gg 1$ for unity Lewis number conditions, Kolla et al. [22] proposed the following algebraic closure of scalar dissipation rate:

$$\widetilde{{\epsilon }_c}=\frac{1}{{\beta }^{\prime }}\left(2K_c^{*}\frac{S_L}{{\delta }_{th}}+\frac{\tilde{\epsilon }}{\tilde{k}}\left[C_3-\tau C_{4}D{a}_L\right]\right)\tilde{c} \left(1-\tilde{c}\right)$$

(5)

where ${\beta }^{\prime }=6.7$, $C_3=1.5\sqrt{K{a}_L}/\left(1+\sqrt{K{a}_L}\right)$ and $C_4=1.1/{\left(1+K{a}_L\right)}^{0.4}$ are the model parameters, $D{a}_L=S_L\tilde{k}/\tilde{\epsilon }{\delta }_{th}$ is the local Damköhler number, $K{a}_L={\left({\delta }_{th}\tilde{\epsilon }/S_L^3\right)}^{1/2}$ is the local Karlovitz number with ${\delta }_{th},\tilde{k}=\overline{\rho u_i^{″}{u}_i^{″}}/2\overline{\rho }$ and $\tilde{\epsilon }=\overline{\mu \left(\partial u_i^{″}/\partial x_j\right)\left(\partial u_i^{″}/\partial x_j\right)}/\overline{\rho }$ being the thermal flame thickness, turbulent kinetic energy and its dissipation rate, respectively. In Equation (5), $K_c^{*}=\left({\delta }_{th}/S_L\right){{\displaystyle \int }}_0^1{[\left(\nabla ·\overrightarrow{u}\right)\rho D\nabla c·\nabla c]}_{Lam}f\left(c\right)dc/{{\displaystyle \int }}_0^1{[\rho D\nabla c·\nabla c]}_{Lam}f\left(c\right)dc$ is a thermochemical parameter [22] with $\overrightarrow{u}$ being the velocity vector.

Equations (4) and (5) can be utilized to obtain the turbulent burning velocity $S_T$ using the Kolmogorov–Petrovski–Piskunov (KPP) theorem [38,39,40]:

$$S_T=2\sqrt{\frac{{\nu }_t}{S{c}_c{\rho }_u}{\left(\frac{\partial {\overline{\dot{\omega }}}_c}{\partial \tilde{c}}\right)}_{\tilde{c}=0}}$$

(6)

where ${\nu }_t=C_{\mu }{\tilde{k}}^2/\tilde{\epsilon }$ is the eddy kinematic viscosity with $C_{\mu }=0.09$, and $S{c}_c$ is the turbulent Schmidt number, which is of the order of unity. Equation (6) upon using Equations (4) and (5) yields the following expression for $S_T/S_L$ [22]:

$$\frac{S_T}{S_L}=\sqrt{\frac{18C_{\mu }}{\left(2c_m-1\right){\beta }^{\prime }}\left\{\left(2K_c^{*}-\tau C_4\right)\left\{\frac{u^{\prime }l}{S_L{\delta }_{th}}\right\}+\frac{2C_3}{3}\frac{{u^{\prime }}^2}{S_L^2}\right\}}$$

(7)

where $S{c}_c=1.0,$ $\tilde{k}=3u^{’2}/2$ and $\tilde{\epsilon }~u^{’3}/l$ are used and these values are understood to be taken at the leading edge of the flame brush. It has been shown elsewhere (Ref. [22]) that the predictions of Equation (7) compare well with experimental data [7,41,42,43] for $S_T$ for premixed flames with a characteristic Lewis number close to unity.

For the sake of completeness, it is useful to consider the other well-known parameterizations of $S_T/S_L$ [7,8,9,19,20,21,22]. Peters [19] proposed the following expression for $S_T/S_L$ based on the leading order balance of the strain rate, kinetic restoration and molecular dissipation of flame surface ratio transport equation in the context of the level-set method:

$$\frac{S_T}{S_L}=1-0.195\frac{l}{{\delta }_z}+\sqrt{{\left(0.195\frac{l}{{\delta }_z}\right)}^2+0.78\left\{\frac{u^{\prime }l}{S_L{\delta }_z}\right\}}$$

(8)

where ${\delta }_z={\alpha }_{Tu}/S_L$ is the Zel’dovich flame thickness with ${\alpha }_{Tu}$ being the unburned gas thermal diffusivity. It is worth noting that Equation (8) is obtained as a positive meaningful root of a quadratic equation, and thus, it is not physically meaningful to modify any model parameters and length scales in isolation in this model expression.

Gülder [20] proposed a parameterization of $S_T/S_L$ based on a large volume of experimental findings in the following manner:

$$\frac{S_T}{S_L}=1+0.62 {\left(\frac{u^{\prime }}{S_L}\right)}^{0.75}{\left(\frac{l}{{\delta }_z}\right)}^{0.25}$$

(9)

Zimont [21] also proposed a similar expression given by:

$$\frac{S_T}{S_L}=1+0.5 {\left(\frac{u^{\prime }}{S_L}\right)}^{0.75}{\left(\frac{l}{{\delta }_z}\right)}^{0.25}$$

(10)

It is important to note that the parameterizations in Equations (7)–(10) are proposed for unity Lewis number conditions, and they do not explicitly account for non-unity Lewis number effects. Moreover, $S_T/S_L=A_T/A_L$ is implicitly assumed in these parameterizations, but this relation is only approximately valid for statistically planar unity Lewis number flames and is rendered invalid for non-unity Lewis number conditions.

Bradley [44] proposed a parameterization of $S_T/S_L$ including the non-unity Lewis number and stretch rate effects in the following manner based on experimental findings:

$$\frac{S_T}{S_L}=1.53 {\left(\frac{u^{\prime }}{S_L}\right)}^{0.55}{\left(\frac{l}{{\delta }_z}\right)}^{0.15}L{e}^{-0.3}$$

(11)

The model expressions provided above are summarised in

Table 1. Summary of normalized turbulent burning velocity $S_T/S_L$ model expressions.

Model	Model Expression
SK model Kolla et al. [22]	$\frac{S_T}{S_L}=\sqrt{\frac{18C_{\mu }}{\left(2c_m-1\right){\beta }^{\prime }}\left\{\left(2K_c^{*}-\tau C_4\right)\left\{\frac{u^{\prime }l}{S_L{\delta }_{th}}\right\}+\frac{2C_3}{3}\frac{{u^{\prime }}^2}{S_L^2}\right\}}$
SP model Peters [19]	$\frac{S_T}{S_L}=1-0.195\frac{l}{{\delta }_z}+\sqrt{{\left(0.195\frac{l}{{\delta }_z}\right)}^2+0.78\left\{\frac{u^{\prime }l}{S_L{\delta }_z}\right\}}$
SG model Gülder [20]	$\frac{S_T}{S_L}=1+0.62 {\left(\frac{u^{\prime }}{S_L}\right)}^{0.75}{\left(\frac{l}{{\delta }_z}\right)}^{0.25}$
SZ model Zimont [21]	$\frac{S_T}{S_L}=1+0.5 {\left(\frac{u^{\prime }}{S_L}\right)}^{0.75}{\left(\frac{l}{{\delta }_z}\right)}^{0.25}$
SB model Bradley [44]	$\frac{S_T}{S_L}=1.53 {\left(\frac{u^{\prime }}{S_L}\right)}^{0.55}{\left(\frac{l}{{\delta }_z}\right)}^{0.15}L{e}^{-0.3}$

for convenience. For statistically planar flames in a canonical configuration, there is no ambiguity in terms of the evaluation of $A_L$ because it is the cross-section of the simulation domain [45]. It is important to appreciate that the experimental evaluations of $A_T$ and $A_L$ are not always straightforward [45]. It has been discussed recently that $A_T$ can be estimated accurately from DNS data using the following expression [45]:

$$A_T={{\displaystyle \int }}_V\left|\nabla c\right|dV$$

(12)

However, the evaluation of $A_L$ gives rise to several possibilities. It is possible to evaluate $A_L$ in the following manner:

$$A_L={{\displaystyle \int }}_V\left|\nabla \overline{c}\right|dV$$

(13)

Another alternative expression for $A_L$ can be obtained as:

$$A_L={{\displaystyle \int }}_V\left|\nabla \tilde{c}\right|dV$$

(14)

Experimental investigations often considered the area of $\overline{c}=0.1$ and 0.5 isosurfaces for the evaluation of $A_L$ [25,41,46,47,48]:

$$A_L=A_{\overline{c}=0.1}\mathrm{or}{A}_L=A_{\overline{c}=0.5}$$

(15)

Borrowing the same concept, $A_L$ could alternatively be calculated based on the area of $\tilde{c}=0.1$ and 0.5 isosurfaces as:

$$A_L=A_{\tilde{c}=0.1}\mathrm{or}{A}_L=A_{\tilde{c}=0.5}$$

(16)

The implications of the assumptions of the evaluation of $A_L$ using Equations (12)–(16) and the performances of $S_T/S_L$ parameterizations based on Equations (7)–(11) will be assessed based on DNS data in Section 4 of this paper. It is worth noting that the turbulent burning velocity in Bunsen flames is often evaluated by using a conventional flame angle method (see Ref. [9] and references therein), which attempts to identify $A_L$ based on the contour of $\overline{c}$ or $\tilde{c}$ in the evaluation of $S_T$. This approach is equivalent to $A_L$ definitions given by Equations (15) and (16). However, the assumption of conical shape due to curvature at the flame tip and estimation of the flame angle based on the tangent to the $\overline{c}$ or $\tilde{c}$ contour introduces additional uncertainties in experimental measurements.

3. Numerical Implementation

A DNS database of turbulent premixed Bunsen flames with characteristic Lewis numbers $Le=0.34, 0.6, 0.8, 1.0$ and 1.2 has been considered in this analysis. The simulations have been conducted using a well-known 3D compressible DNS code SENGA+ [49] where all the spatial derivatives for internal grid points are evaluated using a 10th-order central difference scheme, but the order of accuracy gradually drops to a one-sided 2nd-order scheme at the non-periodic boundaries [49]. An explicit 3rd-order low-storage Runge-Kutta scheme [50] is employed for the time-advancement using a generic single-step Arrhenius type irreversible reaction. For the purpose of an extensive parametric analysis in SENGA+, governing equations of mass, momentum, energy and reaction progress variable $c$ are solved in non-dimensional form which are provided elsewhere [51]. The simulation domain is taken to be a cube with each side of $2d_n$ where $d_n$ is the diameter of the nozzle, and a schematic diagram of the computational domain is presented in Figure 1a. The normalized mean inflow velocity $U_B/S_L$, normalized root-mean-square inlet velocity $u^{\prime }/S_L$ and normalized integral length scale of turbulence (i.e., $l/d_n$ and $l/{\delta }_{th}$) and the grid size along with inlet values of bulk Reynolds number $Re={\rho }_0{U}_B{d}_n/{\mu }_0$, Damköhler number $Da=l{S}_L/u^{\prime }{\delta }_{th}$ and Karlovitz number $Ka={\left(u^{\text{'}}/S_L\right)}^{1.5}{\left(l/{\delta }_{th}\right)}^{-0.5}$ are listed in

Table 2. The turbulence inlet flow parameters for the considered cases.

Case	Le	Re	Grid Size	${\mathit{U}}_{\mathit{B}}/{\mathit{S}}_{\mathit{L}}$	${\mathit{u}}_{\mathit{i}\mathit{n}\mathit{l}\mathit{e}\mathit{t}}^{\text{'}}/{\mathit{S}}_{\mathit{L}}$	$\mathit{l}/{\mathit{d}}_{\mathit{n}}$	$\mathit{l}/{\mathit{\delta }}_{\mathit{t}\mathit{h}}$	Ka	Da
A	0.34	1197	250 × 250 × 250	18.0	1.0	1/5	5.20	0.45	5.00
B	0.6	399	250 × 250 × 250	6.0	1.0	1/5	5.20	0.45	5.00
C	0.8	399	250 × 250 × 250	6.0	1.0	1/5	5.20	0.45	5.00
D	1.0	399	250 × 250 × 250	6.0	1.0	1/5	5.20	0.45	5.00
E	1.2	399	250 × 250 × 250	6.0	1.0	1/5	5.20	0.45	5.00

where $l$ is the integral length scale of turbulence, and ${\delta }_{th}=\left(T_{ad}-T_0\right)/\max |{\nabla T|}_L$ is the thermal flame thickness. The grid spacing is sufficient to resolve the thermal flame thickness ${\delta }_{th}$ and the Kolmogorov length scale $\eta$ for all cases considered here. The simulations have been carried out for different values of characteristic Lewis number of the mixture (i.e., $Le=0.34$, 0.6, 0.8, 1.0 and 1.2). For the purpose of isolating the effects of $Le$, the Zel’dovich number $\beta =T_{ac}\left(T_{ad}-T_0\right)/T_{ad}^2$ (=6.0) and heat release parameter $\tau =\left(T_{ad}-T_0\right)/T_0$ (=4.5) are kept unaltered for all cases where $T_{ac}, T_0$ and $T_{ad}$ are the activation temperature, unburned gas temperature and adiabatic flame temperature, respectively. The Prandtl number and the ratio of specific heats assume standard values (i.e., $Pr=0.7$ and $\gamma =1.4)$. All the cases in Table 2 nominally represent the flamelets regime combustion [52], and their position on Borghi–Peters diagram is shown in Figure 1b. The Lewis number 0.34 case is representative of a lean hydrogen–air mixture of equivalence ratio of 0.40. The Lewis number 0.6 and 0.8 cases are representative of hydrogen-blended methane–air mixtures (e.g., 20% and 10% (by volume) hydrogen blended methane–air flames with overall equivalence ratio of 0.6), and the Lewis number 1.2 case is representative of a hydrocarbon–air mixture involving a hydrocarbon fuel which is heavier than methane (e.g., ethylene–air mixture with equivalence ratio of 0.7) [9,10,33,53]. In the $Le=0.34$ case, the inlet velocity was increased from $6S_L$ to $18S_L$ to avoid flashback, which led to an increase in $Re$ [51].

A mean velocity distribution with a hyperbolic tangent-like profile is used for the inlet boundary condition. The digital filter-based method uses filtered random data in order to obtain realistic pseudo-turbulent velocity correlations [54]. In order to overcome the efficiency problems related to the generation of synthetic turbulent inflow data on large-sized distributed grids, several modifications have been implemented: (a) the Gaussian filter in temporal space in this process was substituted by an autoregressive AR1 process requiring only two time levels; (b) the two-dimensional filter kernel after (a) is replaced by the tensor product of two one dimensional filters, which reduces the cost of the filtering operation for a single grid point from $O\left(N^2\right)$ to $O\left(N\right)$, where $N$ is the number of grid points related to the length of the filter in one direction; (c) identical random seeds for generating inflow data in buffer regions that overlap with neighbouring local domains are used, which avoids the necessity of message passing communications; (d) instead of filtering the inflow data for each local domain located in the inflow plane with its allocated CPU, the filtering is done by all available processors. These measures provide an efficient generation of inflow data which takes the order of 1% of the time required for advancing one time step.

All boundaries apart from the one containing the inlet are taken to be partially non-reflecting outflow and are specified using the Navier–Stokes Characteristic Boundary Conditions (NSCBC) technique [55]. The reacting scalars are initialized using an unstrained premixed laminar flame solution, which is specified as a function of radial distance from the centre of the inlet. The statistics in this analysis are recorded after at least two flowthrough times and two initial eddy turnover times. For the purpose of the evaluation of Reynolds/Favre averaged values (i.e., $\overline{q}$ and $\tilde{q}$), the primitive variable $q$ is averaged in time and also in the azimuthal direction using at least 20 statistically independent snapshots for every $Le$ case considered here. Further information on this database can be obtained from Refs. [51,56,57,58,59,60,61].

4. Results and Discussion

The instantaneous views of $c=0.8$ isourfaces (from the product side) coloured by local values of non-dimensional temperature $\theta$ for the cases considered here are shown in Figure 2. It can be seen from Figure 2 that the extent of flame wrinkling increases with decreasing $Le$, although the inlet turbulence intensity $u^{\prime }/S_L$ remains the same. Moreover, the flame wrinkles which are associated with a convex (concave) shape towards the reactants in the $Le<1$ cases are associated with high (low) temperature values, and this tendency is particularly prevalent for the $Le=0.34$ case. By contrast, high (low) temperature values are obtained where flame wrinkles are concave (convex) towards the reactants in the $Le=1.2$ case. This behaviour is well-known and is consistent with several previous analyses [13,15,18]. The focusing of diffusion of fresh reactants into the reaction zone is stronger than the defocusing of heat in the flame surface elements which are convex towards the reactants in the cases with $Le<1$. This gives rise to the simultaneous presence of high values of reactant concentration and temperature at the convexly curved zones towards the reactants, which leads to further increase in reaction rate magnitude and temperature in these zones. Thus, these regions in the $Le<1$ flames propagate faster than the corresponding unstretched planar flames. By contrast, defocusing of diffusion of fresh reactants into the reaction zone is stronger than the focusing of heat in the flame surface elements which are concave towards the reactants in the cases with $Le<1$. As a result, these zones in the $Le<1$ cases are subjected to the simultaneous presence of low reactant concentration and low temperature, which reduces the propagation rate in these regions. This tendency strengthens with decreasing $Le$, and the combination of high flame propagation rates into the reactants for the convexly curved regions and low propagation rates at the concavely curved zones gives rise to increased flame wrinkling with decreasing $Le$ for the cases with $Le<1$. In the $Le=1.2$ case, the focusing of diffusion of fresh reactants into the reaction zone is weaker than the defocusing of heat in the flame surface elements which are convex towards the reactants. This acts to reduce both the reaction rate magnitude and temperature in the regions which are convex towards the reactants in the $Le=1.2$ case. Just the opposite mechanism leads to high reaction rate magnitude and high temperature in the regions which are concave towards the reactants in the $Le=1.2$ case. The combination of low flame propagation rates into the reactants for the convexly curved regions and high propagation rates at the concavely curved zones leads to reduced flame wrinkling in the $Le=1.2$ case in comparison to the corresponding $Le=1.0$ case.

The increase in the extent of flame wrinkling with decreasing $Le$ is reflected in the increases in $A_T/A_L$ with a decrease in $Le$ when all flow parameters remained unchanged, which can be substantiated from Figure 3 where $A_T/A_L$ is evaluated for different choices of $A_L$ (i.e., $A_L={{\displaystyle \int }}_V\left|\nabla \overline{c}\right|dV$, $A_L={{\displaystyle \int }}_V\left|\nabla \tilde{c}\right|dV$, $A_L=A_{\overline{c}=0.1}$, $A_L=A_{\tilde{c}=0.1}$, $A_L=A_{\overline{c}=0.5}$ and $A_L=A_{\tilde{c}=0.5}$). The $A_T/A_L$ values in Figure 3 show that ${{\displaystyle \int }}_V\left|\nabla \overline{c}\right|dV>A_{\overline{c}=0.5}>A_{\overline{c}=0.1}$ and ${{\displaystyle \int }}_V\left|\nabla \tilde{c}\right|dV>A_{\tilde{c}=0.5}>A_{\tilde{c}=0.1}$ for all cases considered here irrespective of the characteristic Lewis number because the value of $A_T$ remains unchanged for a given case. It can further be seen from Figure 3 that the values of $A_T/A_L$ based on $\overline{c}$ provided higher values than the corresponding definition based on $\overline{c}$ irrespective of the value of $Le$.

The higher values of $A_T/A_L$ for smaller values of the Lewis number also leads to increases in $S_T/S_L$ with a decrease in $Le$. This can indeed be substantiated from Figure 4, which shows that $S_T/S_L$ for a given definition increases with decreasing Lewis number $Le$. The increasing trends of $A_T/A_L$ and $S_T/S_L$ are consistent with several previous analytical [3,4,5,6], experimental [7,8,9,10] and computational [11,12,13,14,15,16,17,18] analyses. Consistent with Figure 3, it can further be seen from Figure 4 that $S_T/S_L$ for $A_L={{\displaystyle \int }}_V\left|\nabla \overline{c}\right|dV$ is the smallest and $S_T/S_L$ for $A_L=A_{\overline{c}=0.1}$ is the highest for a given case among the different choices of $A_L$ definitions. The higher values of $S_T/S_L$ for $A_L=A_{\overline{c}=0.1}$ than for $A_L=A_{\overline{c}=0.5}$ are consistent with previous experimental findings [25,46,47,48] for turbulent premixed Bunsen burner flames. The value of $S_T/S_L$ for $A_L=A_{\overline{c}=0.1}$ was found to be 1.74 times that of the corresponding value for $A_L=A_{\overline{c}=0.5}$ for Kobayashi’s experiments [25], whereas Tamadonfar and Gülder [46] reported a ratio of 2–3.6 for the $S_T/S_L$ values for $A_L=A_{\overline{c}=0.05}$ and $A_L=A_{\overline{c}=0.5}$, and Smallwood et al. [47] reported a ratio of 1.2 to 1.5. The DNS data presented here closely correspond with the ratio reported by Smallwood et al. [47] based on their experimental data. Once again, Figure 4 suggests that $S_T/S_L$ based on $\overline{c}$ yields higher values than the corresponding values based on $\tilde{c}$ irrespective of the value of $Le$.

The values of $R=(S_T/S_L)/\left(A_T/A_L\right)$ for the cases considered here are shown in Figure 5, which shows that the ratio $R$ increases with decreasing $Le$. It is worth noting that the ratio $R=(S_T/S_L)/\left(A_T/A_L\right)={{\displaystyle \int }}_V{\dot{\omega }}_{c}dV/{\rho }_u{A}_T$ is independent of $A_L$, and thus, this quantity is independent of the method of evaluation of $A_L$. A unity value of $R=(S_T/S_L)/\left(A_T/A_L\right)$ indicates the perfect validity of Damköhler’s first hypothesis [62]. It can be seen from Figure 5 that $R$ remains greater than unity (i.e., $R>1$) in the unity Lewis number Bunsen flame case (characterized by a global negative mean curvature), and the detailed physical explanations for this behaviour can be found elsewhere [59] and thus will not be repeated here. Figure 5 further suggests that $R$ increases significantly with decreasing $Le$ for the sub-unity Lewis number (i.e., $Le<1$) cases. The ratio $R$ signifies the ratio of consumption rate per unit area between turbulent and laminar flow conditions. Thus, the increasing trends of $R$ with decreasing $Le$ are the manifestations of thermo-diffusive effects induced by the non-unity Lewis number. Under positively strained $Le<1$ flames, the rate of reactants diffusion into the reaction zone dominates over the diffusion of heat out of this zone, which gives rise to an increase in $R$ with decreasing $Le$. This tendency strengthens further with decreasing $Le$ and was reported in previous analyses [15,17,18]. By comparing Figure 3, Figure 4 and Figure 5, it becomes evident that the increase in $S_T/S_L$ is caused to a larger extent by the increase in ${{\displaystyle \int }}_V{\dot{\omega }}_{c}dV$ than the increase in $A_T/A_L$, and both effects together amplify each other.

The flames become thermo-diffusively unstable for $Le<1$ where high fuel concentration and high temperature values are found at the flame surface elements, which are convex towards the reactants (see Figure 2), and this combination locally augments the reaction and propagation rates much greater than the corresponding value for the same reaction progress variable in an unstretched planar flame. By contrast, the flame elements that are concave towards the reactants propagate slower than the corresponding unstretched planar flame due to the combination of low reactant concentration and temperature for $Le<1$ (see Figure 2). The combination of high flame propagation rates at the convexly curved regions and relatively low flame propagation rates at negatively curved regions acts to increase the extent of flame wrinkling and makes the convexly curved regions stable. These mechanisms act to increase the flame surface area (see Figure 3) and the overall burning rate per unit flame surface area (as can be seen from the large value of $R$ in Figure 5) with a decrease in $Le$. As a consequence, turbulent burning velocity exhibits an increasing trend with decreasing $Le$ (see Figure 4). The aforementioned mechanisms become more prominent when $Le<L{e}_c$ where $L{e}_c$ is the threshold Lewis number under which cellular instability becomes triggered [51,63]. For the thermochemistry used here, the cellular instability is obtained for $Le<L{e}_c\approx 0.6$ [51]. Thus, the $Le=0.34$ case shows much higher values of $A_T/A_L$, $S_T/S_L$ and $R=S_T{A}_L/S_L{A}_T$ than the rest of the cases considered here.

By considering Equation (4) and scaling $\overline{\rho }{\tilde{\epsilon }}_c$ using $\overline{\rho }{\tilde{\epsilon }}_c~{\rho }_0{S}_L{\mathsf{\Sigma }}_{gen}/Le$ and ${\overline{\left({\delta }_z\left|\nabla c\right|\right)}}_s~O\left(1\right)$, one gets: ${\overline{\dot{\omega }}}_c~{\rho }_0{S}_L{\mathsf{\Sigma }}_{gen}/Le$ (see Ref. [17]) where ${\mathsf{\Sigma }}_{gen}=\overline{\left|\nabla c\right|}$ is the generalized FSD [64]. Volume integrating both sides of ${\overline{\dot{\omega }}}_c~{\rho }_0{S}_L{\mathsf{\Sigma }}_{gen}/Le$ yields: $S_T/S_L~L{e}^{-1}\left(A_T/A_L\right)$. It can be seen from Figure 5 that $R_{mod}=Le(S_T/S_L)/\left(A_T/A_L\right)$ remains close to unity for all cases considered here, which is consistent with previous findings (i.e., $S_T/S_L~L{e}^{-1}\left(A_T/A_L\right)$) [15,17,18] for the non-unity Lewis number flames. It is worth noting that based on leading edge theory, Lipatnikov and Chomiak [39] proposed a relation which suggests that $S_T/S_L\approx \left[1/\exp \left(Le-1\right)\right]\left(A_T/A_L\right)$. However, the numerical values of $\left[1/\exp \left(Le-1\right)\right]$ and $L{e}^{-1}$ remain comparable for $Le=0.6,0.8,1.0$ and 1.2, but the values are significantly different for $Le=0.34$. Therefore, $S_T/S_L\approx \left[1/\exp \left(Le-1\right)\right]\left(A_T/A_L\right)$ (or $\exp \left(Le-1\right)\left(S_T{A}_L/S_L{A}_T\right)\approx 1.0$) is not maintained for the $Le=0.34$ case and thus will not be considered hereafter in this paper.

The predictions of different models for $S_T/S_L$ are compared to the normalized turbulent burning velocities obtained for different definitions of $A_L$ in Figure 4. The values used for $K_c^{*}/\tau$ and $c_m$ in the case of the SK model are listed in

Table 3. Thermochemical parameters used in the model by Kolla et al. [22].

Le	${\mathit{K}}_{\mathit{c}}^{*}/\mathit{\tau }$	cm
0.34	0.52	0.92
0.60	0.67	0.87
0.80	0.71	0.867
1.00	0.78	0.825
1.20	0.79	0.816

. It can be seen from Figure 4 that the models proposed by Kolla et al. [22], Peters [19], Gülder [20], Zimont [21] and Bradley [44] (Equations (7)–(11), alternatively SK, SP, SG, SZ and SB models in Table 1) yield comparable values to that of $S_T/S_L$ when $A_L=A_{\overline{c}=0.1}$ is used in the unity Lewis number case, for which these models were originally designed and benchmarked. The quantitative agreement is the best for the model proposed by Gülder [20], Bradley [44] and Kolla et al. [22] (i.e., Equations (7), (9) and (11), alternatively SG, SB and SK models), and the models by Peters [19] and Zimont [21] (i.e., SP and SZ) overpredict and underpredict the magnitude of $S_T/S_L$, respectively, when $A_L=A_{\overline{c}=0.1}$ is used. All of these models overpredict $S_T/S_L$ when $A_L={{\displaystyle \int }}_V\left|\nabla \overline{c}\right|dV$, $A_L={{\displaystyle \int }}_V\left|\nabla \tilde{c}\right|dV$, $A_L=A_{\tilde{c}=0.1}$, $A_L=A_{\overline{c}=0.5}$ and $A_L=A_{\tilde{c}=0.5}$ are used for the unity Lewis number case. It is worth noting that the predictions of the models proposed by Peters [19], Gülder [20] and Zimont [21] (Equations (8)–(11), alternatively SP, SG, SZ and SB models in Table 1) do not change for the non-unity Lewis number cases because $u^{\prime }/S_L$ and $l/{\delta }_{th}$ remain identical for all the cases considered here, and these model expressions do not account for the effects of $Le$. The prediction by the model proposed by Kolla et al. [22] varies with $Le$, and the magnitude of $S_T/S_L$ prediction by this model decreases with decreasing $Le$, due to the variation in $c_m$ and $K_c^{*}/\tau$ (see Table 3). The prediction of the model by Bradley [44] (i.e., Equation (11) or SB model) predicts an increase in the magnitude of $S_T/S_L$ with a decrease in $Le$ for a given set of values of $u^{\prime }/S_L$ and $l/{\delta }_{th}$, and the prediction of Equation (11) remains comparable to $S_T/S_L$ when $A_L=A_{\tilde{c}=0.1}$ is used for the cases with $Le=0.6,0.8$ and 1.2 but overpredicts the values of $S_T/S_L$ when $A_L=\underset{V}{{\displaystyle \int }}\left|\nabla \overline{c}\right|dV$, $A_L=\underset{V}{{\displaystyle \int }}\left|\nabla \tilde{c}\right|dV$, $A_L=A_{\overline{c}=0.1}$, $A_L=A_{\overline{c}=0.5}$ and $A_L=A_{\tilde{c}=0.5}$ are used. However, the model by Bradley [44] (i.e., Equation (11) or SB model) underpredicts all possible definitions of $S_T/S_L$ considered here for the $Le=0.34$ case. It is important to note that Equations (7) and (11) predict the unphysical zero value of $S_T$ for laminar conditions (i.e., $u^{\prime }/S_L=1.0$), and adding 1.0 on the right hand side of Equations (7) and (11) (i.e., SK and SB models) gives rise to overpredictions of $S_T/S_L$ when $A_L=A_{\overline{c}=0.1}$ (or $A_L=A_{\tilde{c}=0.1}$) is used for the cases with $Le=0.6, 0.8$ and 1.2. This also does not affect the underpredictions of $S_T/S_L$ for the $Le=0.34$ case.

As the models of $S_T/S_L$ by Kolla et al. [22], Peters [19], Gülder [20] and Zimont [21] (Equations (7)–(10), alternatively SK, SP, SG and SZ models in Table 1) are originally proposed for unity Lewis number flames with an implicit assumption of $S_T/S_L=A_T/A_L$, it is worth considering the performances of these models when the multiplier $L{e}^{-1}$ is used to account for non-unity Lewis number effects, as shown in

Table 4. Summary of modified normalized turbulent burning velocity $S_T/S_L$ model expressions.

Model	Model Expression
SKL model	$\frac{S_T}{S_L}=\frac{1}{Le}\sqrt{\frac{18C_{\mu }}{\left(2c_m-1\right){\beta }^{\prime }}\left\{\left(2K_c^{*}-\tau C_4\right)\left\{\frac{u^{\prime }l}{S_L{\delta }_{th}}\right\}+\frac{2C_3}{3}\frac{{u^{\prime }}^2}{S_L^2}\right\}}$
SPL model	$\frac{S_T}{S_L}=1-0.195L{e}^{-1}\frac{l}{{\delta }_z}+L{e}^{-1}\sqrt{{\left(0.195\frac{l}{{\delta }_z}\right)}^2+0.78\left\{\frac{u^{\prime }l}{S_L{\delta }_z}\right\}}+\left(\frac{1-Le}{Le}\right)\frac{u^{\prime }/S_L}{u^{\prime }/S_L+1}$
SGL model	$\frac{S_T}{S_L}=1+\left[0.62 L{e}^{-1}{\left(\frac{u^{\prime }}{S_L}\right)}^{0.75}{\left(\frac{l}{{\delta }_z}\right)}^{0.25}+\left(\frac{1-Le}{Le}\right)\frac{u^{\prime }/S_L}{u^{\prime }/S_L+1}\right]$
SZL model	$\frac{S_T}{S_L}=1+\left[0.5L{e}^{-1} {\left(\frac{u^{\prime }}{S_L}\right)}^{0.75}{\left(\frac{l}{{\delta }_z}\right)}^{0.25}+\left(\frac{1-Le}{Le}\right)\frac{u^{\prime }/S_L}{u^{\prime }/S_L+1}\right]$
MSKL model	$\frac{S_T}{S_L}=\frac{1}{Le}\sqrt{\frac{18C_{\mu }}{\left(2c_m-1\right){\beta }^{\prime }}\left\{\left(2K_c^{*}-\tau C_4\right)\left\{\frac{u^{\prime }l}{S_L{\delta }_{th}}\right\}+\frac{2C_3}{3}\frac{{u^{\prime }}^2}{S_L^2}\right\}+\frac{L{e}^2}{[u^{\text{'}}l/S_L{\delta }_{th}+1]}}$
MSB model	$\frac{S_T}{S_L}=1.53 {\left(\frac{u^{\prime }}{S_L}\right)}^{0.55}{\left(\frac{l}{{\delta }_z}\right)}^{0.15}L{e}^{-0.3}+\frac{1}{[u^{\text{'}}l/S_L{\delta }_{th}+1]}$

, because $R_{mod}=Le(S_T/S_L)/\left(A_T/A_L\right)$ assumes a value close to unity (see Figure 5). The expressions for the models by Peters [19], Gülder [20] and Zimont [21] in Table 4 (i.e., SPL, SGL and SZL models in Table 4) are modified in such a manner that for large values of $u^{\prime }/S_L$ (i.e., $u^{\prime }/S_L\gg 1$), the expressions in Table 1 multiplied by $L{e}^{-1}$ are obtained, and the original expressions in Table 1 are recovered for $Le=1$. Moreover, these modifications ensure that $S_T=S_L$ is recovered for the laminar condition (i.e., $u^{\prime }/S_L=0$).

The predictions of the modified expressions for the models by Kolla et al. [22], Peters [19], Gülder [20] and Zimont [21] (SKL, SPL, SGL and SZL models in Table 4) are compared for non-unity Lewis number flames considered here in Figure 6. It can be seen from Figure 6 that the SGL models show a good agreement with $S_T/S_L$ when $A_L=A_{\overline{c}=0.1}$ is used for most of the non-unity Lewis number flames but underpredicts for the $Le=1.2$ case. The SPL model slightly overpredicts $S_T/S_L$ when $A_L=A_{\overline{c}=0.1}$ is used for the $Le=0.34, 0.6$ and 0.8 cases but shows reasonable quantitative agreement with $Le=1.0$ and $1.2$ cases. The $S_T/S_L$ values for $A_L=A_{\overline{c}=0.1}$ are underpredicted for $Le=0.34$ cases by the SKL model (see Table 4), although it shows good agreement for $Le=0.6, 0.8, 1.0$ and $1.2$. The SZL model expression shows good agreement with $S_T/S_L$ obtained from DNS data when $A_L=A_{\overline{c}=0.1}$ is used for the $Le=0.6, 0.8, 1.0$ and 1.2 cases, but $S_T/S_L$ is underpredicted for the $Le=0.34$ case.

A modified version of Kolla et al. [22] (i.e., SKL model) can be proposed in the following manner so that the expression for $S_T/S_L$ for $A_L=A_{\overline{c}=0.1}$ can be reasonably predicted for non-unity Lewis number flames while satisfying $S_T/S_L=1.0$ for $u^{\prime }=0$:

$$\frac{S_T}{S_L}=\frac{1}{Le}\sqrt{\frac{18C_{\mu }}{\left(2c_m-1\right){\beta }^{\prime }}\left\{\left(2K_c^{*}-\tau C_4\right)\left\{\frac{u^{\prime }l}{S_L{\delta }_{th}}\right\}+\frac{2C_3}{3}\frac{{u^{\prime }}^2}{S_L^2}\right\}+\frac{L{e}^2}{[u^{\text{'}}l/S_L{\delta }_{th}+1]}}\mathrm{for}{A}_L=A_{\overline{c}=0.1}$$

(17)

Equation (17) will henceforth be referred to as the MSKL model in this paper. The predictions of Equation (17) are also shown in Figure 6, which shows that the MSKL model offers comparable performance to that of the SKL model. Similar to Equation (17), the SB can be modified (henceforth referred to MSB model) in the following manner so that $S_T/S_L$ assumes a value of 1.0 for $u^{\prime }=0$:

$$\frac{S_T}{S_L}=1.53 {\left(\frac{u^{\prime }}{S_L}\right)}^{0.55}{\left(\frac{l}{{\delta }_z}\right)}^{0.15}L{e}^{-0.3}+\frac{1}{[u^{\text{'}}l/S_L{\delta }_{th}+1]}$$

(18)

The prediction of the MSB model remains comparable to that of the SB model. The MSB predicts $S_T/S_L$ reasonably accurately when $A_L=A_{\tilde{c}=0.1}$ is used for the cases with $Le=0.6, 0.8$ and 1.2 but overpredicts the values of $S_T/S_L$ for other definitions of $A_L$. Similar to the SB model, the MSB model also underpredicts all possible definitions of $S_T/S_L$ considered for the $Le=0.34$ case considered here.

The $L_2$-norm of the relative error values $E_{DNS}=\sqrt{{\displaystyle \sum }_{k=1}^N{\left|\left(x_k^{DNS}-x_k^{Model}\right)/x_k^{DNS}\right|}^2}$ (where $x_k$ refers to the value of $S_T/S_L$ for the $k$^th case and ‘DNS’ and ‘Model’ superscripts are used for DNS and model expression values, respectively and $N$ is the total number of different cases considered here) of the model expressions of $S_T/S_L$ listed in Table 1 are shown in Figure 7, which shows that SK, SP, SG, SZ and SB models exhibit comparable $E_{DNS}$, but $E_{DNS}$ values are high for definitions $A_L={{\displaystyle \int }}_V\left|\nabla \overline{c}\right|dV$, $A_L={{\displaystyle \int }}_V\left|\nabla \tilde{c}\right|dV$, $A_L=A_{\overline{c}=0.5}$ and $A_L=A_{\tilde{c}=0.5}$. The corresponding $E_{DNS}$ values of the model expressions for $S_T/S_L$ listed in Table 4 are shown in Figure 8. A comparison between Figure 7 and Figure 8 reveals that the modified expressions in Table 4 significantly decrease the $E_{DNS}$ values when $A_L=A_{\overline{c}=0.1}$ is used for the evaluation of $S_T/S_L$.

It is worth noting that the DNS database used for this analysis considers only a single modest value of turbulence intensity and a single pressure level. Therefore, it is worth considering the predictions of the model expressions summarized in Table 4 for experimental conditions for non-unity Lewis number conditions for different pressures and higher turbulence intensity, which were recently analysed by Lipatnikov et al. [10]. These experimental conditions by Lipatnikov et al. [10] are summarized in

Table 5. Summary of the experimental conditions from Lipatnikov et al. [10].

Conditions	Mixture	u′/SL	l/δth	τ	p	Le
C1	H2-air, $\varphi =0.45$	2.7	63	3.7	1.0 bar	0.35
C2	H2/O2/He, $\varphi =0.45$	2.6	21	3.0	1.0 bar	0.91
C3	CH4-air, $\varphi =1.0$	2.7	69	6.5	1.0 bar	1.00
C4	H2-air, $\varphi =0.45$	4.6	152	3.7	3.0 bar	0.35
C5	H2/O2/He, $\varphi =0.45$	6.3	33	3.0	3.0 bar	0.91
C6	CH4-air, $\varphi =1.0$	4.2	145	6.5	3.0 bar	1.00
C7	H2-air, $\varphi =0.45$	7.0	179	3.7	5.0 bar	0.35
C8	H2/O2/He, $\varphi =0.45$	9.0	19	3.0	5.0 bar	0.91
C9	CH4-air, $\varphi =1.0$	5.4	203	6.6	5.0 bar	1.00

. The predictions of the model expressions in Table 4 are compared to $S_T/S_L$ values reported by Lipatnikov et al. [10] in Figure 9 and the corresponding $L_2$-norm of relative error values $E_{Expt}=\sqrt{{\displaystyle \sum }_{k=1}^N{\left|\left(x_k^{Expt}-x_k^{Model}\right)/x_k^{Expt}\right|}^2}$ (where $x_k$ refers to the value of $S_T/S_L$ for the $k-$th case and ‘Expt’ superscript is used for experimental values) are shown in Figure 10.

It can be seen from Figure 9 that the MSKL, SKL and SPL models predict $S_T/S_L$ values obtained from experimental conditions [10] reasonably well for conditions C1, C2, C5 and C8. The performances of the SGL and SZL models remain comparable to those of the MSKL, SKL and SPL models for C1 and C2 conditions, and the SGL model predicts the experimental value of $S_T/S_L$ satisfactorily for the C5 condition, but both SGL and SZL models underpredict $S_T/S_L$ values significantly for other cases. The MSKL and SKL model expressions overpredict $S_T/S_L$ values for C3, C4, C6, C7 and C9 conditions. Among these conditions, only under the C4 condition, the predictions of the SPL and SZL models remain comparable to those of the MSKL and SKL models and experimental value of $S_T/S_L$, but in other conditions, the SPL, SGL and SZL models underpredict $S_T/S_L$ values. It can indeed be seen from Figure 10 that the SKL, SPL, SGL and MSKL model expressions yield relatively small values of $E_{Expt}$. It is worth noting that the original model expressions by Kolla et al. [22], Peters [19], Gülder [20] and Zimont [21] (i.e., SK, SP, SG and SZ models) underpredict significantly for conditions C1, C2, C4, C5, C7 and C8 and return higher values of $E_{Expt}$ than the SPL, SGL and SZL models and thus are not explicitly shown in Figure 9 and Figure 10. The model MSB underpredicts for all the conditions considered in the experiment.

As is explained in Section 2, most models for turbulent burning velocity were developed for unity Lewis number conditions. As can be seen in Figure 4, while all these models do quite well for $Le=1$, most models used in the present work fail for non-unity Lewis numbers. The models are made to work for non-unity Lewis numbers by taking advantage of the scaling (i.e., $S_T/S_L=L{e}^{-1}{A}_T/A_L$) that is observed in Figure 5, which shows that $R_{mod}=Le{S}_T{A}_L/S_L{A}_T\approx 1$. However, this along with the assumptions in the closure equations is not a precise expression and could have led to the overshoot and undershoot of the predictions. It is worth noting that the experimental methodologies of the measurements of $S_T$ and the projected flame brush area $A_L$ are different from case to case and from the DNS-based methodologies. It has also been noted by Driscoll [26] that the parameterization of turbulent burning velocity might also be geometry dependent to some extent (note experiments in Ref. [10] were carried out for spherical flame kernels and not for Bunsen burner flames). These uncertainties may contribute to the discrepancies between model predictions and $S_T/S_L$ values obtained from DNS/experimental data. However, a careful comparison between Table 5 and Figure 9 reveals that several models provide reasonable predictions when the characteristic Lewis number remains close to unity. Moreover, it can be seen from Figure 10 that the newly proposed modifications to the existing parameterizations of turbulent burning velocity do reasonably well in an average sense, even for non-unity Lewis number flames.

Based on the results shown in Figure 4, Figure 7, Figure 8, Figure 9 and Figure 10, it seems that the SGL, SPL and SZL model expressions provide reasonable estimates of $S_T/S_L$ obtained from both DNS and experiments for a wide range of values of $Le$ and different pressure levels, even though there are quantitative discrepancies for some cases. It is worth noting that the experimental methodologies of the measurements of $S_T$ and the projected flame brush area $A_L$ are different to the DNS-based methodologies. Furthermore, ${Le(S}_T/S_L)/(A_T/A_L)\approx 1.0$ is approximately valid in a scaling sense, but it is not a precise expression. These uncertainties may contribute to the discrepancies between model predictions and $S_T/S_L$ values obtained from DNS/experimental data.

5. Conclusions

The statistical behaviour of the turbulent burning velocity $S_T$ has been analysed based on a Direct Numerical Simulation database of turbulent premixed Bunsen flames for different values of characteristic Lewis number $Le$ ranging from 0.34 to 1.2. It has been found that the turbulent burning velocity normalized by the laminar burning velocity $S_T/S_L$ and the turbulent flame surface area normalized by the projected flame brush area $A_T/A_L$ for a given set of flow parameters increase with decreasing $Le$. Further, the value of $S_T$ depends strongly on the definition of the projected flame brush area $A_L$, and six different alternatives (i.e., $A_L={\int }_V|\nabla \stackrel{-}{c}|dV$, $A_L={\int }_V|\nabla \tilde{c}|dV$, $A_L=A_{\tilde{c}=0.1}$, $A_L=A_{\stackrel{-}{c}=0.5}$ and $A_L=A_{\tilde{c}=0.5}$) have been considered here. The variations of the numerical value of $S_T$ can be up to 50% depending on the definition of $A_L$. It has been found that the highest values of $S_T/S_L$ and $A_T/A_L$ for a given Bunsen premixed flame case are obtained when the projected flame brush area $A_L$ is based on the area of $\stackrel{-}{c}=0.1$ isosurface, and the smallest value is obtained when $A_L$ is evaluated based on the volume-integration of $|\nabla \tilde{c}|$. The choice of higher value of $\stackrel{-}{c}$ or $\tilde{c}$ for the evaluation of $A_L$ has been found to give rise to smaller values of $S_T/S_L$ and $A_T/A_L$ values in this configuration, and the quantitative variation has been found to be as large as 50% in the cases analysed here. This is found to be consistent with previous experimental findings [25,46,47,48] for turbulent premixed Bunsen flames. Moreover, it has been found that $S_T/S_L$ and $A_T/A_L$ values based on $\stackrel{-}{c}$ assume greater values than the corresponding values based on $\tilde{c}$. The predictions of different models for $S_T/S_L$ [19,20,21,22,44] have been compared to the corresponding values obtained from DNS data. It has been found that most available parameterizations of $S_T/S_L$ provide comparable predictions when $A_L$ is evaluated based on the area of the $\stackrel{-}{c}=0.1$ isosurface. However, most of these models were originally proposed for the unity Lewis number [19,20,21,22], where Damköhler’s first hypothesis (i.e., $S_T/S_L=A_T/A_L$) [62] is implicitly assumed. However, Damköhler’s first hypothesis (i.e., $S_T/S_L=A_T/A_L$) is rendered invalid for the non-unity Lewis number flames, and the existing parameterizations tend to underpredict $S_T/S_L$ values for flames with $Le<1$, and this tendency strengthens with a decreasing Lewis number. It has been found that ${Le(S}_T/S_L)/(A_T/A_L)$ remains of the order of unity which has been justified using scaling arguments. This relation has been utilized to extend parameterizations of $S_T/S_L$ (which were originally proposed for the unity Lewis number) for non-unity Lewis number conditions while satisfying the limiting condition $S_T=S_L$ for $u^{\text{'}}=0$. These revised parameterizations based on the models proposed by Peters [19], Gülder [20] and Zimont [21] have been found to be more successful than the other models in terms of $S_T/S_L$ predictions across the datasets considered when $A_L$ is evaluated based on the area of $\stackrel{-}{c}=0.1$ isosurface especially for $Le\ll 1$ flames compared to the original model expressions. The performances of these revised parameterizations have been assessed for a recent experimental dataset [10] for different values of $Le$. It has been found that the expressions, which provide satisfactory quantitative agreement with DNS data, also exhibit small values of $L_2$-norm of the percentage error with respect to the experimental data [10]. However, the performance of these model parameterizations needs to be assessed for different flame configurations under higher turbulence intensity and also in the presence of detailed chemistry and transport, which will form the platform for further investigations.