# A Layered Description of a Premixed Flame Stabilized in Stagnating Turbulence

Arnaud Mura, Vincent Robin, Kim Q N Kha, Michel Champion

► To cite this version:

Arnaud Mura, Vincent Robin, Kim Q N Kha, Michel Champion. A Layered Description of a Premixed Flame Stabilized in Stagnating Turbulence. Combustion Science and Technology, 2016, Recent Advances in Turbulent Premixed Combustion Modeling, 188 (9), pp.1592-1618. 10.1080/00102202.2016.1195822. hal-05502137

HAL Id: hal-05502137

https://hal.science/hal-05502137v1

Submitted on 10 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# A Layered Description of a Premixed Flame Stabilized in Stagnating Turbulence

Arnaud Mura, Vincent Robin, Kim Q. N. Kha & Michel Champion

**To cite this article:** Arnaud Mura, Vincent Robin, Kim Q. N. Kha & Michel Champion (2016) A Layered Description of a Premixed Flame Stabilized in Stagnating Turbulence, *Combustion Science and Technology*, 188:9, 1592-1618, DOI: 10.1080/00102202.2016.1195822

**To link to this article:** https://doi.org/10.1080/00102202.2016.1195822

<figure><img src="image_2.png" /></figure>

Published with license by Taylor & Francis
Group, LLC© Arnaud Mura, Vincent Robin,
Kim Q. N. Kha, and Michel Champion

<figure><img src="image_3.png" /></figure>

Published online: 02 Jun 2016.

<figure><img src="image_4.png" /></figure>

Submit your article to this journal

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

Article views: 777

<figure><img src="image_7.png" /></figure>

View related articles

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

View Crossmark data

<figure><img src="image_10.png" /></figure>

CrossMark
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

OPEN ACCESS

# A Layered Description of a Premixed Flame Stabilized in Stagnating Turbulence

Arnaud Mura, Vincent Robin, Kim Q. N. Kha, and Michel Champion

Institut Pprime - UPR 3346 - CNRS - ISAE-ENSMA - Université de Poitiers, Futuroscope, France

## ABSTRACT

In the thin-flame regime of turbulent premixed combustion, instantaneous progress variable gradients are essentially fixed by propagating flamelets. The laminar flamelet internal structure thus imposes, at least to some extent, the progress variable one-point one-time statistics, i.e., the progress variable PDF $P(c; x, t)$. In the present study a generalized presumed probability density function (PDF) shape is considered to account for possible departures from the thin-flame limit. The parameters that characterize this PDF shape depend on the ratio of the Kolmogorov length scale to laminar flame width. The corresponding framework may be associated to a layered description of the turbulent flame brush (TFB) with the thicknesses of the different sub-layers determined from the knowledge of both the Reynolds and Karlovitz numbers, $Re_T$ and $Ka$. It incorporates boundary layers on both sides of the TFB, which are associated to the finite thickness of the local flamelets, i.e., small but finite values of the Karlovitz number. This description leads to a modeling proposal for premixed turbulent combustion. A tabulation is constructed based on canonical one-dimensional planar unstrained laminar premixed flame computations and the different quantities are tabulated as functions of $Re_T$ and $Ka$. The final closure is applied to the numerical simulations of premixed flames stabilized in stagnating turbulence.

## ARTICLE HISTORY

Received 13 October 2015
Revised 7 January 2016
Accepted 19 February 2016

## KEYWORDS

Flame brush; Flamelet regime; Leading edge; Premixed flame; Stagnating turbulence; Turbulent combustion

## Introduction

The development of computational fluid dynamics (CFD) tools over the last decades has allowed their current applications to the numerical simulations of various industrial turbulent flows featuring complex geometries and realistic conditions. However, the physical phenomena involved in turbulent reactive flows display such a wide range of characteristic scales that the statistical, i.e., Reynolds Averaged Navier-Stokes (RANS), representation or the filtered, i.e., large eddy simulations (LES), description of these flows are still required to proceed with their numerical simulations. The direct numerical simulation (DNS) of such reactive flows in complex geometry remains indeed unlikely in the foreseeable future. The use of LES may provide a satisfactory description of the large-scale turbulent dynamics at the price of important computational efforts and costs so that its practical use for engineering purposes still remains limited. Accordingly, RANS and U-RANS approaches remain the reference tools to proceed with the industrial design

CONTACT Arnaud Mura

<figure><img src="image_13.png" /></figure>

arnaud.mura@ensma.fr

<figure><img src="image_14.png" /></figure>

Institut Pprime, UPR 3346 CNRS, Université de Poitiers, ENSMA,

Téléport 2-1, Avenue Clément Ader, Futuroscope 86961, France.

Color versions of one or more of the figures in the article can be found online at www.tandfonline.com/gcst.
<figure><img src="image_15.png" /></figure>

of large-scale energy conversion devices. In practice, the resort to any of these two frameworks (i.e., RANS or LES) does involve the use of some numerical and modeling parameters so that the corresponding computations cannot be considered as fully predictive thus far. Therefore, the models involved in these numerical simulations and their behaviors must be perfectly understood so as to obtain results that recover relevant trends and can be useful to complement experimental data. The objective of the present study is to provide the analysis and detailed sensitivity investigation of a recent modeling proposal through its application to the numerical simulation of an impinging premixed turbulent flame.

One specific difficulty associated with the numerical simulation of turbulent reactive flows lies in the domains of validity of turbulent combustion models, which remain quite restricted. For instance, the main part of the flame brush can be often described as belonging to the flamelet regime, which is characterized by infinitely small chemical time scales compared to characteristic turbulent time scales. However, on the one hand, the local increase of turbulent kinetic energy from one given location to the other may lead to some thickening of the flame associated to finite-rate chemistry effects. On the other hand, towards the fresh reactants side, the characteristic chemical time becomes infinitely large in comparison with the characteristic turbulent time scales, thus making irrelevant the flamelet assumption to describe this part of the flame brush, i.e., the leading edge. Therefore, the turbulent combustion models proposed for the numerical simulation of premixed reactive flowfields must be as simple as possible but must also be realistic enough to handle such a wide range of variations of the turbulent combustion regime. Here, we choose to analyze a simple model for the mean reaction rate, able to change its functional form inside the computational domain, depending on the turbulent combustion regime encountered. The relevant criterion accounting for these possible changes is based herein on a layered description of the turbulent flame brush structure that depends upon both Reynolds and Karlovitz numbers.

## Finite thickness flamelets statistical description

In this study, fully premixed conditions are studied so that a unique progress variable $c(x, t)$, defined to be $c=0$ in the fresh unburned mixture and $c=1$ in combustion products, is considered to characterize the composition of the mixture. Accordingly, Lewis numbers for all species are close to unity, and heat losses as well as composition inhomogeneities are not addressed. Thus, the Favre averaged progress variable $\tilde{c} = \overline{\rho c}/\overline{\rho}$ is the solution of the following balance equation:

$$ \frac{\partial \bar{\rho} \tilde{c}}{\partial t} + \frac{\partial}{\partial x_k} \left( \bar{\rho} \tilde{v}_k \tilde{c} + \overline{\rho v_k'' c''} - \overline{\rho D \frac{\partial c}{\partial x_k}} \right) = \bar{\Omega}_c \quad (1) $$

where $v_k$ denotes the velocity field, $D$ the molecular diffusion coefficient, and $\Omega_c = \rho\omega_c$ is the chemical reaction rate expressed in kg·m⁻³·s⁻¹. In the above equation, the quantities $c'' = c - \tilde{c}$ and $v_k'' = v_k - \tilde{v}_k$ denote Favre fluctuations of the progress variable and velocity components, respectively. The quantity $\overline{\rho v_k'' c''}$ denotes the turbulent fluxes of the progress variable; the reader is referred to Robin et al. (2011) and references therein for further details concerning its modeling. Under the low Mach number assumption, the
<figure><img src="image_16.png" /></figure>

thermodynamic contribution $p_0(x, t)$ to the pressure field $p(x, t)$ is considered to be homogeneous, i.e., $p_0(x, t) = p_0(t)$. Therefore, density may be expressed as a function of the composition variable only, i.e., $\rho(c) = \rho_u/(1 + \tau c)$, where $\tau = (\rho_u/\rho_b) - 1$ is the expansion factor.

As it is well known, the closure of the mean chemical rate $\overline{\Omega}_c$ appearing in the right-hand side (RHS) of Eq. (1) is one of the most important sub-problems in the field of turbulent combustion modeling and various expressions or functions of the mean progress variable $\tilde{c}$, i.e., $\overline{\Omega}_c(\tilde{c}) = \overline{\Omega}_c(\tilde{c}) = \bar{\rho}(\tilde{c})\tilde{\omega}_c(\tilde{c})$, can be found in the literature. Most of the available closures have been developed under the flamelet hypothesis (Bray, 2016). In the corresponding flamelet regime of turbulent combustion (Borghi, 1985; Peters 1986; Williams, 1976), it is standardly assumed that the molecular processes, i.e., molecular diffusion and chemical reactions, are the dominant processes. However, the laminar flame thickness $l_f$ does not provide any firm information about the real thickness of flamelets present in the turbulent flows. Therefore, we proposed recently to introduce a characteristic length scale $l_x$, defined as the relevant turbulent flamelet length scale at which molecular processes involved in the flame become the dominant processes (Kha et al., 2016). In the first step of the analysis, the value of $l_x$ has been approximated by the Kolmogorov length scale $\eta_K$. The ratio between the laminar flame thickness and the turbulent flamelets length scale $Rf_x$ is thus considered to be of the order of the square root of the Karlovitz number: $Rf_x = l_f/l_x = O(l_f/\eta_K) = O(\sqrt{Ka})$. Accordingly, the extreme values of the progress variable that delineate the flamelet profile are no longer zero and unity but will be denoted $c_{min}$ and $c_{max}$. Beyond these limits, the laminar flame structure is indeed expected to be altered by the turbulent flowfield (Anand and Pope, 1985; Mura et al., 2003). These values will depend on the Karlovitz number. For instance, as established in the recent analysis by Kha et al. (2016), we have $c_{min} = \exp(-1/Rf_x)$ and $c_{max} - 1$ in the limit of infinitely

<figure><img src="image_17.png" /><figcaption>**Figure 1.** Evolutions of $\delta_c = c_{max} - c_{min}$ as a function of the Karlovitz number. These results are obtained in the limit of large values of the Zel'dovich number $\beta$ and for two distinct approximations of $l_x$: $l_x = \eta_K$ (dashdotted line) and $l_x = l_G$ (continuous line).</figcaption></figure>
<figure><img src="image_18.png" /></figure>

large values of the Zel'dovich number, such that $\delta_c = c_{\max} - c_{\min} = 1 - \exp(-1/Rf_x)$. Figure 1 displays the evolution of $\delta_c$ as a function of $Ka$ for two distinct approximations of $l_x$. The first corresponds to the one discussed above, i.e., $l_x = \eta_K$, which leads to $\delta_c = 1 - \exp(-1/\sqrt{Ka})$. The second assumption implies $l_x = l_G$ with $l_G$ being the Gibson scale introduced by Peters (1986, 1992). It may be worth recalling here that the Gibson scale is determined by setting the eddy turnover velocity to the laminar flamelet propagation velocity $S_L$, in such a manner that $l_G = l_T$. $(S_L/v_{RMS})^3$ with $v_{RMS}$ the characteristic turbulent velocity fluctuation such that $Re_T = v_{RMS} \cdot l_T/v$. In this case, it can be readily shown that $\delta_c = 1 - \exp(-1/Ka^2)$.

Considering that the flamelets characteristic scales $l_x$ may differ from the laminar flame thickness, $l_f$, the following expression of the progress variable PDF is introduced:

$$ 
\begin{aligned}
\tilde{P}(c; \mathbf{x}, t) &= \delta(c - c_{\min}), \quad \forall \tilde{c} \in [0, c_{\min}] \\
\tilde{P}(c; \mathbf{x}, t) &= \alpha_c \delta(c - c_{\min}) + \beta_c \delta(c - c_{\max}) + \gamma_c \tilde{f}(c), \quad \forall \tilde{c} \in [c_{\min}, c_{\max}] \\
\tilde{P}(c; \mathbf{x}, t) &= \delta(c - c_{\max}), \quad \forall \tilde{c} \in [c_{\max}, 1]
\end{aligned}
\quad (2) $$

where $\alpha_c$ and $\beta_c$ denote the probability that $c$ equals $c_{\min}$ and the probability that $c$ equals $c_{\max}$, respectively. This PDF can be thought as a generalization of the well-known BML form (Bray and Moss, 1977), which is recovered in the limit associated to infinitely large values of the Damköhler number ($\gamma_c \to 0$) and infinitely small values of the Karlovitz number ($c_{\min} \to 0$) and ($c_{\max} \to 1$). In the above expressions, the function $\tilde{f}(c)$ represents the PDF conditioned inside the internal structure of flamelets. This function can be deduced from the structure of a one-dimensional planar unstrained premixed laminar flame:

$$ \tilde{f}(c) = C_1^* \left| \frac{dc}{dx^*} \right|^{-1} \quad (3) $$

where $C_1^* = C_1 l_f$ is a nondimensionalized normalization constant:

$$ 
\begin{aligned}
C_1^* &= \left( \int_{c_{\min}}^{c_{\max}} \left| \frac{dc}{dx} \right|^{-1} dc \right)^{-1} \cdot l_f \\
&= \left( \int_{x_{\min}}^{x_{\max}} dx \right)^{-1} \cdot l_f = \frac{l_f}{l_x} = Rf_x = O(\sqrt{Ka})
\end{aligned}
\quad (4) $$

As emphasized above, the values of $c_{\min}$ and $c_{\max}$ may differ from zero and unity, i.e., the Karlovitz number is not infinitely small. Thus, the progress variable PDF, as defined by Eq. (2), as well as the internal PDF and its nondimensionalized normalization constant $C_1^*$, take into account the finite thickness of flamelets.

Accordingly, first- and second-order moments of the progress variable are:

$$ \tilde{c} = \alpha_c c_{\min} + \beta_c c_{\max} + \gamma_c I_1, \quad \forall \tilde{c} \in [c_{\min}, c_{\max}] \quad (5) $$

$$ \tilde{c}'' = \alpha_c c_{\min}^2 + \beta_c c_{\max}^2 + \gamma_c I_2 - \tilde{c}^2, \quad \forall \tilde{c} \in [c_{\min}, c_{\max}] \quad (6) $$
with

$$I_1 = \int_{c_{\min}}^{c_{\max}} c \tilde{f}(c) dc \qquad (7)$$

$$I_2 = \int_{c_{\min}}^{c_{\max}} c^2 \tilde{f}(c) dc \qquad (8)$$

and

$$\int_{c_{\min}}^{c_{\max}} \tilde{f}(c) dc = 1 \qquad (9)$$

It seems worth emphasizing that, according to Eq. (2), we have $\tilde{c} = c_{\min}$ for $\tilde{c} \in [0, c_{\min}]$ and $\tilde{c} = c_{\max}$ for $\tilde{c} \in [c_{\max}, 1]$, while the second-order centered moment $\tilde{c}''^2$ is zero for mean progress variable values $\tilde{c}$ that lie in $[0, c_{\min}] \cup [c_{\max}, 1]$.

The probability to find flamelets is proportional to the intermittency factor between the fresh reactants and fully burned gases $\tilde{c}(1 - \tilde{c})$ and to the departure from the maximum segregation rate value $S_{\max} - S$ [see, for instance, Bray et al. (2006)]:

$$\gamma_c = (S_{\max} - S) \frac{\tilde{c}(1 - \tilde{c})}{I_1 - I_2} \qquad (10)$$

where the segregation rate S reaches its maximum value $S_{\max}$, which, in the limit of large values of the Reynolds number, can be expressed as follows:

$$S_{\max} = \frac{(c_{\max} - \tilde{c})(\tilde{c} - c_{\min})}{\tilde{c}(1 - \tilde{c})} \qquad (11)$$

The above equality indeed provides a generalized expression of the maximum segregation rate consistent with Eq. (2).

Thus, the expression of the mean chemical rate can be written as:

$$\tilde{\omega}_c = \frac{\overline{\Omega}_c}{\overline{\rho}} = \gamma_c I_{\omega_c} = (S_{\max} - S) \frac{\tilde{c}(1 - \tilde{c})I_{\omega_c}}{I_1 - I_2} \qquad (12)$$

where $I_{\omega_c}$ is the mean chemical rate conditioned inside flamelet structures, a quantity obtained from laminar flame computations:

$$I_{\omega_c} = \int_{c_{\min}}^{c_{\max}} \omega_c(c) \tilde{f}(c) dc \qquad (13)$$

Finally, the modified form of the PDF proposed above evidences the existence of boundary layers associated to the flame edges.

## Properties of finite thickness flamelets

The properties of laminar flamelets are now studied in detail. Different levels of description can be retained to conduct the analysis. The limit of infinitely large values of the Zel'dovich number $\beta$ provides a theoretical framework to get analytical results. For instance, in this limit, the progress variable evolution is given by $c(x^*) = \exp(x^*)$, where the normalized coordinate $x^* = x/l_f$ has been introduced.
<figure><img src="image_19.png" /></figure>

For the purpose of analyzing the influence of finite values of the Zel'dovich number, laminar flame calculations based on simplified chemistry and transport properties may be first considered. The numerical calculations are performed by using the Mathematica™ software, and the conservation equations are solved in a frame of reference attached to the flame so that only the steady-state solution is evaluated. Thus, the transport equation of the progress variable writes:

$$ \rho_u v_u \frac{dc}{dx} = \frac{d}{dx} \left( \rho D \frac{dc}{dx} \right) + \Omega_c \quad (14) $$

where $v_u$ denotes the velocity in the fresh reactants, and the reaction rate is given by:

$$ \Omega_c = \mathcal{B} \frac{\rho_u(1-c)}{(1+\tau c)} \exp\left(-\frac{\beta(1-c)}{1-\alpha(1-c)}\right) \quad (15) $$

where $\mathcal{B}$ is the reduced pre-exponential factor and $\alpha = (T_b - T_u)/T_b$ is a temperature factor (Williams, 1985). The chemical reaction rate can be recast in $\Omega_c = \rho_u \mathcal{B} \mathcal{F}(c; \tau, \beta)$ with the function $\mathcal{F}(c; \tau, \beta)$ defined as follows:

$$ \mathcal{F}(c; \tau, \beta) = \frac{(1-c)}{1+\tau c} \exp\left[\frac{-\beta(1-c)}{1-\tau(1-c)/(1+\tau)}\right] $$

Then, the dimensionless form of the transport equation (14) is:

$$ \frac{dc}{dx^*} = \frac{d}{dx^*} \left( \frac{\rho D}{\rho_u D_u} \frac{dc}{dx^*} \right) + \omega_c \quad (16) $$

where $\omega_c = \Omega_c D_u / \rho_u v_u^2 = E \mathcal{F}(c; \tau, \beta)$ is the dimensionless chemical reaction rate. The quantity $D_u$ denotes the molecular diffusion coefficient at the cold edge of the flame and the nondimensional number $E$ is defined by $E = \mathcal{B} D_u / v_u^2$. The above equation is solved in the domain $\epsilon \le c \le 1 - \epsilon$ by following the procedure of Bray et al. (2006). More realistic laminar flame structures may be evaluated from a detailed chemical mechanism using the Cantera library. A methane air flame is considered here for a value of the equivalence ratio $\Phi = 1.0$. Chemistry is described with the detailed GRI3.0 mechanism (Bowman et al., 1997). The resulting profiles of the progress variable $c(x)$ and the corresponding scalar gradient $dc/dx$ are presented in Figure 2 with the progress variable evolution deduced from $\rho = \rho_u / (1 + \tau c)$. The burning velocity value is 38.6 cm/s. The detailed chemical mechanism leads to a slow evolution of the progress variable from $c = 0.9$ to $c = 1.0$ on the burned gases side. This is clearly visible on the scalar gradient profiles (see Figure 2). It is expected that turbulence may alter the laminar flamelet profile in this specific zone where the scalar gradients are rather weak and some experimental data bring support to this picture (Chen and Bilger 2002; Chen et al., 1996).

The conditional scalar moments $I_1$ and $I_2$, see Eqs. (7) and (8), can be calculated from the knowledge of the laminar flame structure and associated internal PDF shape $\tilde{f}(c)$, and the flamelet boundaries $c_{min}$ and $c_{max}$. As emphasized above, these flamelet boundaries $c_{min}$ and $c_{max}$ are related to the value of the Karlovitz number. Therefore, the quantities defined by Eqs. (7), (8), and (13) are related to both the flamelet inner structure and the value of the Karlovitz number. In the case of infinitely large values of the Zel'dovich number, i.e., $\beta \to \infty$, the upper boundary is such that $c_{max} \to 1$ and the lower limit $c_{min}$ is
<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 2.** Laminar flame structure: progress variable $c(x)$ plotted vs. the spatial coordinate $x$ (top), and its gradient $dc/dx$ plotted vs. $c$ (bottom) obtained for a mixture of methane and air (equivalence ratio $\Phi = 1.0$) from a computation based on a detailed description of chemistry.</figcaption></figure>

related to the Karlovitz number and may be expressed from $c_{min} = \exp(-1/\sqrt{Ka})$. Moreover, in this case, the scalar gradient through the laminar flame structure is given by $dc/dx^* = c$, except in the reaction zone ($c \to 1$). Thus, the conditional scalar moments $I_1$ and $I_2$ can be obtained analytically:

$$I_1^{\beta \to \infty} = \int_{c_{min}}^{1} c \tilde{f}(c) dc = C_1^* \int_{c_{min}}^{1} \left| \frac{dc}{dx^*} \right|^{-1} cdc = \sqrt{Ka} \int_{c_{min}}^{1} dc = \sqrt{Ka}(1 - c_{min}) \quad (17)$$

$$I_2^{\beta \to \infty} = \int_{c_{min}}^{1} c^2 \tilde{f}(c) dc = C_1^* \int_{c_{min}}^{1} \left| \frac{dc}{dx^*} \right|^{-1} c^2 dc = \sqrt{Ka} \int_{c_{min}}^{1} cdc = \frac{\sqrt{Ka}}{2} (1 - c_{min}^2) \quad (18)$$

where the expression of the PDF $\tilde{f}(c)$ has been used and, according to Eq. (4), the normalization constant $C_1^*$ has been approximated by $\sqrt{Ka}$.

Then, considering the relation between $c_{min}$ and the Karlovitz number, the conditional scalar moments can be expressed as functions of $c_{min}$:
<figure><img src="image_22.png" /></figure>

$$I_1^{\beta \to \infty} = \frac{1 - c_{\min}}{\ln(1/c_{\min})} \qquad (19)$$

$$I_2^{\beta \to \infty} = \frac{1 - c_{\min}^2}{2 \ln(1/c_{\min})} \qquad (20)$$

or as functions of the Karlovitz number:

$$I_1^{\beta \to \infty} = \sqrt{Ka} \left( 1 - \exp \left( - \frac{1}{\sqrt{Ka}} \right) \right) \qquad (21)$$

$$I_2^{\beta \to \infty} = \frac{\sqrt{Ka}}{2} \left( 1 - \exp \left( - \frac{2}{\sqrt{Ka}} \right) \right) \qquad (22)$$

Equations (21) and (22) provide the first- and second-order scalar moments conditioned inside flamelets for large values of the Zel'dovich number as functions of the Karlovitz number only. Flamelet structures obtained from the numerical calculations discussed at the beginning of this section can now be used to evaluate these conditional quantities. Figure 3 displays the integral quantities as functions of the Karlovitz number. It also shows the difference between $I_1$ and $I_2$, which is an important quantity for the evaluation of the mean chemical rate [see Eq. (12)]. It is worth noticing that the ratio $I_2/I_1$ is also an important quantity since, in the limit of infinitely small values of $Ka$, one can show that it imposes the segregation rate at the leading edge<fnref n="1" />:

$$S_{(\tilde{c} \to 0)} = \left[ \frac{c''_{2}}{\tilde{c}(1 - \tilde{c})} \right]_{\tilde{c} \to 0} = \left[ \frac{(I_2/I_1) - \tilde{c}}{(1 - \tilde{c})} \right]_{\tilde{c} \to 0} = \frac{I_2}{I_1} \qquad (23)$$

where the following two identities have been used: $[c''_{2}]_{\tilde{c} \to 0} = \gamma_c I_2 - \tilde{c}^2$ and $[\tilde{c}]_{\tilde{c} \to 0} = \gamma_c I_1$ [see Kha et al. (2016) for further details].

The first part of Figure 3 (on the top) corresponds to the case of large values of the Zel'dovich number, i.e., eqs. (21) and (22), the middle one to the case of finite Zel'dovich number ($\beta = 10$) but simplified chemistry and the third figure (on the bottom) corresponds to detailed chemistry calculations ($\Phi = 1$). These results show that, for vanishingly small values of the Karlovitz number, $I_1$ and $I_2$ vanish for either large values of the Zel'dovich number or finite Zel'dovich number provided that simplified chemistry is considered. However, when detailed chemistry is considered, $I_1$ and $I_2$ tends to unity for vanishingly small values of the Karlovitz number, which is ascribed to the specific behavior of the progress variable gradient observed at the burned gas side when detailed chemistry is considered. It is worth noting that, even in the small Karlovitz number limit, the behaviors of $I_1$ and $I_2$ do depend on the level of description retained for the chemistry, i.e., single step or detailed chemistry. We will see below that this dependence, however, disappears as soon as we consider the difference between $I_1$ and $I_2$, and their ratio $I_2/I_1$.

The quantities $I_1$ and $I_2$ indeed do not appear separately in the mean chemical rate expression but through the difference $I_1 - I_2$ as well as the ratio $I_2/I_1$. Whatever the description of the chemistry, i.e., single step or detailed chemistry, it is indeed

<footnote n="1">¹Note that $\beta_c$ may be neglected compared to $a_c$ and $\gamma_c$ at the leading edge.</footnote>
<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /><figcaption>**Figure 3.** Evolutions of laminar flamelet properties $l_1$, $l_2$, $l_1 - l_2$, and $l_2/l_1$ as functions of the Karlovitz number. The results reported on the top correspond to the infinitely large Zel'dovich number limit ($\beta \to \infty$), those in the middle to one-step chemistry calculations conducted with $\beta = 10$, and those presented at the bottom have been obtained from detailed chemistry computations ($\Phi = 1.0$).</figcaption></figure>
<figure><img src="image_27.png" /></figure>

remarkable that, for small values of the Karlovitz number $Ka$, the difference $I_1 - I_2$ vanishes and the ratio $I_2/I_1$ tends to unity (see second and third graphs of Figure 3). The behavior observed for $I_2/I_1$ in the infinitely large Zel'dovich limit appears to be more specific since, as previously underlined in the recent study by Kha et al. (2016), it tends to 0.5 when $Ka \to 0$. Figure 3 also shows that taking into account the reaction zone, i.e., finite values of the Zel'dovich number, slightly modifies the value of the ratio $I_2/I_1$ only for small values of the Karlovitz number. Another important result is that the difference $I_1 - I_2$ is a constant approximately equal to 0.2 for a large range of Karlovitz number ($0.2 < Ka < 1$). Figure 4 shows the evolution of the flamelet properties $I_1 - c_{min}$ and $c_{max} - I_1$ with the Karlovitz number. These two quantities behave in a way very similar to the one observed for $I_1$ and $I_2$, as shown in Figure 3. Moreover, for large values of the Karlovitz number ($0.5 < Ka < 1$), they are almost constant whatever the laminar flame structure considered. We shall see later on that these two differences $I_1 - c_{min}$ and $c_{max} - I_1$ dictate the thickness of the leading and trailing edge layers of the turbulent flame brush.

<figure><img src="image_28.png" /><figcaption>**Figure 4.** Evolution of two flamelet properties $I_1 - c_{min}$ (top) and $c_{max} - I_1$ (bottom) as functions of the Karlovitz number for different laminar flame structure.</figcaption></figure>
<figure><img src="image_29.png" /></figure>

## Layered description of the turbulent flame and tabulation of finite thickness flamelets characteristics

The turbulent flame brush structure as a whole must be now analyzed in the light of the previous discussion. Figure 5 shows the typical sketch of a turbulent flame brush consistent with the finite thickness flamelets presumed PDF of Eq. (2). Accordingly, the flame brush is bounded by $c_{min} < \tilde{c} < c_{max}$ and we use the following layered description of the turbulent flame brush (Kha et al., 2016):

* Upstream of the turbulent flame brush ($\tilde{c} < c_{min}$), a turbulent preheat layer must be considered. The chemical reactions are exponentially small in this zone and can therefore be discarded from the analysis; the PDF support reduces to $\tilde{P}(c, x) = \delta(c - \tilde{c})$.

* The appearance of the first flamelet structure corresponds to $c_{min} < \tilde{c} < c_{le}$; this is the leading edge of the turbulent flame brush. In this region, the PDF is approximated by $\tilde{P}(c) = (1 - \gamma_c)\delta(c - c_{min}) + \gamma_c\tilde{f}(c)$.

* In the region $c_{le} < \tilde{c} < c_{te}$, the composition changes abruptly from $c = c_{min}$ to $c = c_{max}$. The molecular processes involved in flames occur at the smallest scales of the turbulent flow and the probability $\gamma_c$ to lie inside flamelets is small enough to consider that second-order averaged quantities (variance, scalar turbulent fluxes, etc.) are proportional to $(c_{max} - \tilde{c})(\tilde{c} - c_{min})$. Then the mean chemical rate may be evaluated from Eq. (12).

* The trailing edge of the flame brush ($c_{te} < \tilde{c} < c_{max}$) is associated to the disappearance of the last flamelet structure. In this region, the PDF is approximated by $\tilde{P}(c) = (1 - \gamma_c)\delta(c - c_{max}) + \gamma_c\tilde{f}(c)$.

* Finally, downstream of the flame brush for $c > c_{max}$, the strength of chemical reactions becomes weaker and weaker, i.e., the characteristic chemical time scale becomes large in comparison with the characteristic turbulent time scale. Therefore, the PDF can be approximated from $\tilde{P}(c, x) = \delta(c - \tilde{c})$.

<figure><img src="image_30.png" /><figcaption>Figure 5. Sketch of the layered description of the turbulent flame brush.</figcaption></figure>
<figure><img src="image_31.png" /></figure>

According to this description of the structure of the turbulent flame brush, the mean chemical rate is expressed in five different sub-layers, respectively: the turbulent preheat layer, the leading edge, the main (intermittent) layer, the trailing edge, and the weak reactive layer. It is necessary to evaluate the boundaries of each of these layers, namely $c_{min}$, $c_{le}$, $c_{te}$, and $c_{max}$. It has been shown earlier that the value of the Karlovitz number determines the flame brush boundaries, i.e., $c_{min}$ and $c_{max}$. Equation (5), once applied to the leading and trailing edges, yields:

$$c_{le} = (1 - \gamma_{le})c_{min} + \gamma_{le}I_1 = c_{min} + \gamma_{le}(I_1 - c_{min}) \quad (24)$$

and

$$c_{te} = (1 - \gamma_{te})c_{max} + \gamma_{te}I_1 = c_{max} - \gamma_{te}(c_{max} - I_1) \quad (25)$$

In these expressions the approximation

$$\gamma_{le} = \gamma_{te} = (\kappa Re_T)^{-3/4} = \kappa^* Re_T^{-3/4} \quad (26)$$

introduced previously by Kha et al. (2016) is used to evaluate $\gamma_{le}$ and $\gamma_{te}$; we assume that the probability to be inside these smallest structures is proportional to the ratio of their thickness $l_x = l_{\eta\kappa}$ to the largest scale of the turbulence, i.e., the integral length scale $l_T$, which is responsible for the largest wrinkles of the flame front. The sensitivity of the computational results to the parameter $\kappa$ will be evaluated afterwards. From these expressions, one can deduce the thickness in the composition space of the leading and trailing edge layers as $\delta_{le} = c_{le} - c_{min} = \gamma_{le}(I_1 - c_{min})$ and $\delta_{te} = c_{max} - c_{te} = \gamma_{te}(c_{max} - I_1)$, respectively.

Considering the dependency of the different boundaries on the Reynolds and Karlovitz numbers, the mean quantities can be tabulated as functions of these two parameters. From a general point of view, tabulations may either be completed before the numerical simulations or generated on-the-fly during the simulation. The tabulation methodology we use here follows the first of these two strategies and is based on the canonical 1D unstrained laminar premixed flame computations that have been presented in the previous section. However, in the present study it is parameterized by the values of the Reynolds and Karlovitz number, an original feature. From the couple $(Ka, Re_T)$, the values $c_{min}(Ka)$, $c_{le}(Ka, Re_T)$, $c_{te}(Ka, Re_T)$, and $c_{max}(Ka)$ are determined as well as the different sub-layers of the turbulent flame brush. These evaluations are based on the Cantera computations of premixed laminar flamelets incorporating detailed chemistry. The corresponding profiles of chemical species and temperature are themselves parameterized as functions of a progress variable $c$, defined as a marker of the heat released by chemical reaction through $\rho = \rho_u/(1 + \tau c)$ and weighted by the progress variable PDF $\tilde{P}(c)$ introduced above. The tabulations includes not only the boundaries of the sub-layers but incorporates also mean species mass fractions as well as any other quantity that is required to evaluate the mean chemical rate, especially $I_1$, $I_2$, and $I_{\omega_c}$.
<figure><img src="image_32.png" /></figure>

From the above layered description, the mean chemical rate expressions are given by:

$$ \tilde{\omega}_c(\tilde{c}) = \begin{cases} \omega_c(\tilde{c}) \approx 0 & \tilde{c} \in [0; c_{\min}] \\ (1 - \gamma_c)\omega_c(c_{\min}) + \gamma_c I_{\omega_c}^* & \tilde{c} \in [c_{\min}; c_{le}] \\ (S_{\max} - S)\tilde{c}(1 - \tilde{c})I_{\omega_c}/(I_1 - I_2) & \tilde{c} \in [c_{le}; c_{te}] \\ (1 - \gamma_c)\omega_c(c_{\max}) + \gamma_c I_{\omega_c}^{**} & \tilde{c} \in [c_{te}; c_{\max}] \\ \omega_c(\tilde{c}) & \tilde{c} \in [c_{\max}; 1] \end{cases} \quad (27) $$

where $I_{\omega_c}^*$ and $I_{\omega_c}^{**}$ are defined by $I_{\omega_c}^* = \int_{c_{\min}}^{c_{\max}} \omega_c(c)\tilde{f}(c)dc$ and $I_{\omega_c}^{**} = \int_{c_{\min}}^{c_{\max}} \omega_c(c)\tilde{f}(c)dc$. However, it has been experienced in the numerical simulations that the leading and trailing edges sub-layers are quite difficult to discriminate from the turbulent preheat layer ($c < c_{\min}$) and weak reactive layer ($c > c_{\max}$). Due to their small thicknesses it is indeed difficult to capture these thin layers on conventional computational meshes. Therefore, a reduced three-layer model of mean chemical rate is used to conduct the numerical simulations:

$$ \tilde{\omega}_c(\tilde{c}) = \begin{cases} \omega_c(\tilde{c}) \approx 0 & \tilde{c} \in [0; c_{le}] \\ \gamma_c I_{\omega_c} & \tilde{c} \in [c_{le}; c_{te}] \\ \omega_c(\tilde{c}) & \tilde{c} \in [c_{te}; 1] \end{cases} \quad (28) $$

Thus, in an attempt to generalize available fast chemistry closures, the present proposal consists in modifying the closure of the mean chemical at the leading and trailing edges of the turbulent flame brush. As emphasized above, in the present application the leading edge sub-layer (LEL) is merged with the turbulent preheat layer (TPL) and it is acknowledged that in the latter the mean progress variable evolution is driven by mean convection and turbulent transport. It is noteworthy that some specific features of the turbulent flame exhibit some sensitivity to the quenching criterion $c_{le}$. For instance, previous one-dimensional traveling wave (TW) analyses have established that, due to the quenching influence, it is quite difficult to recover the classical Kolmogorov-Petrovskii-Piskunov (KPP) solution (Kolmogorov et al., 1937) from a computational estimation of the turbulent flame speed, see Catlin and Lindstedt (1991) and Sabelnikov et al. (1998). The detailed investigation conducted by Corvellec (1998) establishes that the convergence to the theoretical KPP solution, as the quenching value is reduced, is quite slow. Different propagation characteristics of the mean reacting flow may be obtained depending on the accuracy of the numerical integration and on the treatment of the boundary conditions. This study also emphasized that, in practical CFD computations of turbulent premixed flames, i.e., 2D and 3D computations performed with standard discretization schemes, it is almost impossible to get rid of the computational procedure influence and the numerical solution will remain altered by the mesh resolution and boundary conditions (Sabelnikov et al., 1998). The present procedure is quite a contrast to the above discussion. Indeed, the present study is aimed at avoiding the systematic resort to the infinitely fast chemistry approximation (i.e., the telegraphic functional form $\tilde{c}(1 - \tilde{c})$ of the mean chemical rate) in the whole turbulent flame brush thickness (TFB). As emphasized above, such an hypothesis does not hold on the leading and trailing edges of the TFB but is implied in the KPP TW-analyses. A quenching value is finally introduced but it is based on physical grounds and consideration of small but finite values of the Karlovitz number. It involves the thicknesses in the progress variable space of both the turbulent preheat layer, from 0 to
<figure><img src="image_33.png" /></figure>

$c_{min}$, and leading edge layer, from $c_{min}$ to $c_{le}$. The closure of the mean reaction rate provided by Eq. (28) is now assessed through its application to a practical geometry.

# Computational results

## Description of the numerical setup

The experimental setup that is considered to proceed with the numerical simulations corresponds to a premixed flame stabilized in stagnating turbulence. This geometry has been studied inter alia by Cho et al. (1988, 1989). In these experiments a turbulent stream of reactants at room temperature impinges on a refractive surface so as to stabilize a planar premixed turbulent flame. For many years it has been recognized that such impinging flames provide a reference test case for the validation of premixed turbulent combustion models (Bray et al., 1991, 1994, 2000; Dong et al., 2013; Lee and Huh, 2004; Sabelnikov and Lipatnikov, 2011). In the reference papers by Cho et al. (1988, 1989) data on the velocity and temperature characteristics of these flames are provided. The corresponding set of experimental data is considered as a reference database to assess turbulent combustion closures and it has been used to this purpose in many other studies [see, for instance, the recent analyses conducted by Sabelnikov and Lipatnikov (2011) or Dong et al. (2013)]. A uniform axisymmetric flow of methane and air is provided by a 50-mm diameter nozzle. This jet of reactants is surrounded by a co-flow of air with an outer diameter of 100 mm. The turbulent flames are stabilized in the stagnating flow thanks to a refractive plate placed at a distance of 75 mm above the jet exit plane. The bulk flow velocity of the jet is held constant at 5.0 m/s and the incoming turbulence is generated by a grid. The characteristic integral scale of turbulence $l_T$ can be estimated<fnref n="2" /> to lie between one and two millimeters (Dong et al., 2013). The fluctuating velocity RMS $v_{RMS}$ is approximately equal to 0.3 m/s. The equivalence ratio of the methane-air mixture is unity ($\Phi = 1$).

From the computational point of view, a 2D mesh, made of approximately 15,900 computational cells, is used to represent a small part of the physical domain (see Figure 6). The inlet condition at the bottom of the computational domain is imposed by using the experimental data while an outlet boundary condition is applied on the top right side. Boundary conditions on the left side of the computational domain are given by the symmetry of the flow and on the right side it is approximated by a special condition, which can be either an inlet or an outlet. The numerical simulations are performed with the Code Saturne, a parallel general-purpose 3D low-Mach-number computational fluid dynamics (CFD) solver developed by EDF (Archambeau et al., 2004). It relies on a finite-volume approximation of the Navier-Stokes equations. The set of equations to be solved consists of the averaged Navier-Stokes equations completed with equations for turbulence modeling and for additional scalars (enthalpy, concentration of species, mixture fraction, as well as any other user-defined scalars). The time marching scheme is based on a velocity field prediction step followed by a pressure-correction step. Equations for turbulence and scalars are solved separately afterwards. The discretization in space is based on the fully

<footnote n="2">²Since there are some uncertainties regarding this value, we will evaluate the sensitivity of the computational results to its specification.</footnote>
<figure><img src="image_34.png" /></figure>

<figure><img src="image_35.png" /><figcaption>**Figure 6.** Sketch of the computational setup retained for the numerical simulation of the impinging flame. Top: Computational mesh and boundary conditions description. Bottom: Flow streamline deflection across the flame brush delineated by the averaged heat release rate (HRR).</figcaption></figure>

conservative, unstructured finite-volume framework, with a fully colocated arrangement for all variables. For the present application balance equations for the mean momentum are solved in conjunction with a classical $k-\epsilon$ closure and a transport equation for the mean progress variable [see Eq. (1)]. The boundary conditions for turbulent quantities $k$ and $\epsilon$ at the inlet jet of methane and air are, respectively, $k_0 = 3(u_{\text{RMS}})^2/2$ and $\epsilon_0 = C_{\mu}k_0^{3/2}/l_T$ with $C_{\mu} = 0.09$. A passive scalar, i.e., mixture fraction $\xi$, is also considered to follow the mixing between the jet of reactants and the coflowing jet of air.

## Reference case

To proceed with the numerical simulation, the mean chemical rate in Eq. (1) is represented by the closure discussed in the previous section, see Eq. (28). It is worth noting that, in a first step of the analysis, the progress variable scalar flux $\overline{\rho v_k'' c''}$ is closed by using the standard, i.e., gradient-diffusion law: $\overline{\rho v_k'' c''} = -\bar{\rho} D_T (\partial \tilde{c} / \partial x_k)$ with $D_T = v_T / Sc_T$. The value of the turbulent Schmidt number is set to its reference value $Sc_T = 0.7$ (Launder,
<figure><img src="image_36.png" /></figure>

1976). The present set of results has been obtained with $l_T = 1$ mm and a value of the constant $\kappa$ equal to 0.05, see Eq. (26). The sensitivity of the computational result to the choice of this value will be evaluated subsequently.

It is noteworthy that the model used to perform the present set of numerical simulations has been set forth for fully premixed conditions only. Further work is clearly required to extend it to partially premixed conditions and we therefore did not attempt to describe turbulent combustion in regions where the flow of premixed reactants is subject to dilution by ambient air. The mean chemical rate is thus set to zero as soon as the mixture fraction significantly differs from its value at the nozzle exit. This can be seen in Figures 7 and 8 reporting the field of the mixture fraction $\tilde{\xi}$ superimposed with progress variable isolines and the mean progress variable field $\tilde{c}$ superimposed with mean-chemical rate isolines, respectively. In Figure 8 we also compare the mean field of the progress variable obtained from the transport equation (1) to another expression of the progress variable that can be evaluated from the tabulated detailed chemistry:

$$ \tilde{c}_{Y_{CO}+Y_{CO2}} = \frac{\tilde{Y}_{CO} + \tilde{Y}_{CO2}}{Y_{CO}^{eq} + Y_{CO2}^{eq}} \qquad (29) $$

where $Y_{CO}^{eq}$ and $Y_{CO2}^{eq}$ denote the CO and CO₂ mass fractions at equilibrium.

The profiles of $c_{min}$, $c_{le}$, $c_{te}$, and $c_{max}$, evaluated along the symmetry axis, are reported in Figure 9. This figure shows that the values of $c_{min}$ are too small to resolve the turbulent preheat sublayer [0, $c_{min}$]. It also firmly confirms that it would have been quite difficult to discriminate between $c_{te}$ and $c_{max}$, except at the price of prohibitively high computational costs. The profiles of the mean axial velocity component and mean

<figure><img src="image_37.png" /><figcaption>Figure 7. Averaged mixture fraction field superimposed with mean progress variable iso-lines.</figcaption></figure>

<figure><img src="image_38.png" /><figcaption>Figure 8. Mean progress variable field as obtained from its transport equation (right side) or deduced from the tabulation by using Eq. (29) (left side) superimposed with mean chemical rate iso-lines.</figcaption></figure>
<figure><img src="image_39.png" /></figure>

<figure><img src="image_40.png" /><figcaption>Figure 9. Evolutions of $c_{min}$, $c_{le}$, $c_{te}$, and $c_{max}$ obtained along the symmetry axis of the experimental setup.</figcaption></figure>

progress variable are presented together with experimental data in Figure 10. A quite satisfactory quantitative agreement is obtained between the numerical simulation results and the experimental data. The sensitivity of the corresponding results to changes in the turbulent Schmidt number value $Sc_T$ is illustrated in Figure 11. The sensitivity of the computational results to a modification of the equivalence ratio $\Phi$ of the incoming stream of reactants has been also examined and, in Figure 12, results obtained for $\Phi = 0.8$ are compared to those obtained for stoichiometric conditions, i.e., $\Phi = 1$, which correspond to experimental operating conditions. It is interesting to see that, for a value of $\kappa$ held constant, the results confirm that the model is sensitive to mixture composition variations. Finally, Figure 13 provides a comparison with another set of results issued from a computation performed with $c_{le} = 0$ and $c_{te} = 1.0$. It evidences the influence of the consideration of the leading and trailing edges sublayers onto the computational results.

### Impact of the CGD modeling on the first moment quantities

Since it is well known that the experimental geometry under consideration is subject to counter-gradient diffusion (Bray et al., 1998), we examine the possible influence of these phenomena on the computational results. Therefore the next set of computational results has been obtained by using an algebraic closure that accounts for the counter-gradient diffusion (CGD) influence. The corresponding closure may be written as follows,

$$ \overline{\rho v_k'' c''} = (\tau S_L \overline{\rho c''^2} \lambda - \bar{\rho} D_T ||\nabla \tilde{c}||) \tilde{M}_k, \quad (30) $$
<figure><img src="image_41.png" /></figure>

<figure><img src="image_42.png" /><figcaption>**Figure 10.** Top: Comparison between the axial profile of the mean progress variable and experimental data documented by Cho et al. (1988). Bottom: Comparison between the axial profile of the mean axial velocity component and experimental data documented by Cho et al. (1988).</figcaption></figure>

where $\tilde{M}_k$ denotes the component of the mean orientation $\tilde{M}$ of the flux, which is approximated by $\tilde{M} = \nabla \tilde{c} / ||\nabla \tilde{c}||$ and $\lambda$ is a parameter related to the fluctuations of the local flame front. Such a closure was previously used by Dong et al. (2013) to conduct the numerical simulation of the same experimental test case.

Figure 14 reports the scalar flux as obtained (i) from the classical gradient-diffusion law (GD), (ii) from the above generalized closure (CGD), and (iii) experimental data. From this figure representing the scalar flux as a function of the mean progress variable, it cannot be disputed that, in the present conditions, CGD processes are
<figure><img src="image_43.png" /></figure>

<figure><img src="image_44.png" /><figcaption>**Figure 11.** Top: Comparison between the axial profile of the mean progress variable and experimental data documented by Cho et al. (1988). Bottom: Comparison between the axial profile of the mean axial velocity component and experimental data documented by Cho et al. (1988). Different values of the turbulent Schmidt number $Sc_T$ are considered.</figcaption></figure>

essential to represent satisfactorily the scalar fluxes. A comparison between the two sets of computational results obtained with a closure that incorporates or not the CGD effects is reported in Figures 15. These figures report the profiles of the mean progress variable and the mean axial component of the velocity respectively. It is quite interesting to see that, compared to the results obtained from the generalized CGD closure, the standard GD closure leads to a satisfactory level of agreement with experimental
<figure><img src="image_45.png" /></figure>

<figure><img src="image_46.png" /><figcaption>**Figure 12.** Top: Comparison between the axial profile of the mean progress variable and experimental data documented by Cho et al. (1988). Bottom: Comparison between the axial profile of the mean axial velocity component and experimental data documented by Cho et al. (1988). Calculations were performed for two distinct values of the equivalence ratio: Φ = 0.8 and Φ = 1.0; the latter corresponds to experimental conditions.</figcaption></figure>

data. Therefore, for the purpose of the present investigation, which is not focused on the representation of counter-gradient diffusion (CGD) and associated flame-generated turbulence (FGT) effects, the standard GD approximation will be kept throughout the rest of this study.
<figure><img src="image_47.png" /></figure>

<figure><img src="image_48.png" /><figcaption>**Figure 13.** Top: Comparison between the axial profile of the mean progress variable and experimental data documented by Cho et al. (1988). Bottom: Comparison between the axial profile of the mean axial velocity component and experimental data documented by Cho et al. (1988). Calculations were performed by considering either the multilayered flame brush representation (continuous line) or by forcing $c_{le}$ and $c_{te}$ to zero and unity, respectively (dotted line).</figcaption></figure>

**Leading edge layer thickness: A sensitivity analysis**

In the model for the mean chemical rate introduced above, an important and critical parameter of the final closure is the value of the constant $\kappa$, which relates through Eq. (26) the probability to lie in the flamelet structure to the turbulent Reynolds number. Equation (25) evidences the explicit dependency of the leading edge layer thickness on
<figure><img src="image_49.png" /></figure>

<figure><img src="image_50.png" /><figcaption>**Figure 14.** Scalar flux plotted vs. the mean progress variable. GD denotes computational results obtained with the standard diffusion-like closure of the progress variable turbulent fluxes. CGD denotes computational results obtained with Eq. (30).</figcaption></figure>

the value of this parameter. It is therefore decided to analyze the dependency of the computational results on the specification of the value of $\kappa$.

Figure 16 reports the profiles of the progress variable fluxes and mean axial component of the velocity respectively. These results have been obtained by using three distinct values of $\kappa$, which are $\kappa = 0.05$, 0.5, and 5. This offers quite a large range of variations for this parameter since the maximal value is 100 times larger than the smallest one. As expected, the obtained results display some dependency on the value used to perform the numerical simulation but it must be acknowledged that the level of sensitivity does not seem to be critical and we will see below that the sensitivity of the computational results to the specification of the inlet boundary conditions is much larger.

**Sensitivity of the computational results to the specification of the inlet flow integral length scale**

Figure 17 reports the profiles of the progress variable fluxes and mean axial component of the velocity respectively. These results have been obtained for three distinct values of the integral length scale imposed at the inlet boundary (jet nozzle exit). The corresponding values are $l_T = 0.5$ mm, $l_T = 1$ mm, and finally $l_T = 2$ mm. Unsurprisingly, the obtained results display an important sensitivity to the specification of $l_T$, which further emphasizes the crucial need for well-documented experimental databases featuring also well-defined boundary conditions.
<figure><img src="image_51.png" /></figure>

<figure><img src="image_52.png" /><figcaption>**Figure 15.** Top: Comparison between the axial profile of the mean progress variable and experimental data documented by Cho et al. (1988). Bottom: Comparison between the axial profile of the mean axial velocity component and experimental data documented by Cho et al. (1988). GD denotes computational results obtained with the standard diffusion-like closure of the progress variable turbulent fluxes. CGD denotes computational results obtained with Eq. (30).</figcaption></figure>

## Conclusions

In this work, a presumed PDF of the progress variable has been used to account for the influence of finite thickness flamelets and associated progress variable boundaries. This leads to the introduction of a layered description of the turbulent flame brush (TFB) based on the consideration of both Reynolds and Karlovitz numbers effects. This description is subsequently used to propose a first-order closure of the mean
<figure><img src="image_53.png" /></figure>

<figure><img src="image_54.png" /><figcaption>**Figure 16.** Sensitivity of the computational results to variations of the modeling parameter κ. Top: Comparison between the axial profile of the mean progress variable and experimental data documented by Cho et al. (1988). Bottom: Comparison between the axial profile of the mean axial velocity component and experimental data documented by Cho et al. (1988).</figcaption></figure>

progress variable. The closure is derived under the assumption of sufficiently large values of the turbulence Reynolds number, with the molecular fluxes neglected with respect to its turbulent counterpart, but accounts for finite values of the Karlovitz number. The corresponding formalism has been used in conjunction with tabulated chemistry to compute Laboratory flames and the preliminary results obtained are rather encouraging.

A sensitivity analysis has also been conducted with respect to the main parameter of the modeling proposal and, as expected, the obtained results display some dependency
<figure><img src="image_55.png" /></figure>

<figure><img src="image_56.png" /><figcaption>**Figure 17.** Sensitivity of the computational results to the specification of the integral length scale $l_T$ at the nozzle exit. Top: Comparison between the axial profile of the mean progress variable and experimental data documented by Cho et al. (1988). Bottom: Comparison between the axial profile of the mean axial velocity component and experimental data documented by Cho et al. (1988).</figcaption></figure>

on the value used to perform the numerical simulation. However, it has been undoubtedly shown that there exists an equivalent sensitivity to the specification of the inlet boundary condition, i.e., the turbulence integral length scale at the nozzle exit. To conclude, it is quite clear that it would be useful to conduct additional validations of the proposed closure against other well-documented experimental data. It would be especially interesting to analyze the response of the computational model in conditions featuring stronger departures from the thin laminar flamelet combustion regime. In
this respect, the highly turbulent premixed V-shaped flames studied by Renou and co-workers (Mazellier et al., 2010) provide an interesting test case so as to proceed with this additional assessment of the proposed model.

## References

Anand, M.S., and Pope, S.B. 1985. Calculations of premixed turbulent flames by PDF methods. *Combust. Flame*, **67**, 127-142.

Archambeau, F., Mehitoua, N., and Sakiz, M. 2004. Code Saturne: A finite volume code for turbulent flows. *Int. J. Finite Volumes*, **1**. Available at: https://hal.archives-ouvertes.fr/hal01115371

Borghi, R. 1985. On the structure and morphology of turbulent premixed flames. In C. Casci (Ed.), *Recent Advances in the Aerospace Science*, Plenum Publishing, New York, pp. 117-138.

Bowman, C.T., Hanson, R.K., Gardiner, W.C., Frenklach, V., Lissianski, M., Goldenberg M., Smith, G.P., Crosley, D.R., and Golden, D.M. 1997. An optimized detailed chemical reaction mechanism for methane combustion and NO formation and reburning, Technical Report No. GRI97/0020. Gas Research Institute, Chicago, IL.

Bray, K.N.C. 2016. Laminar flamelets in turbulent combustion modeling. *Combust. Sci. Technol.*, **188**(9), 1372-1375.

Bray, K.N.C., Champion, M., and Libby, P.A. 1991. Premixed flames in stagnating turbulence. Part I: The general formulation for counterflowing streams and gradient models for turbulent transport. *Combust. Flame*, **84**, 391-410.

Bray, K.N.C., Champion, M., and Libby, P.A. 1994. Flames in stagnating turbulence. In P.A. Libby and F.A. Williams (Eds.), *Turbulent Reacting Flows*, Academic Press, New York.

Bray, K.N.C., Champion, M., and Libby, P.A. 2000. Premixed flames in stagnating turbulence. Part IV: A new theory for Reynolds stresses and Reynolds fluxes applied to impinging flows. *Combust. Flame*, **120**, 1-18.

Bray, K.N.C., Champion, M., Libby, P.A., and Swaminathan, N. 2006. Finite rate chemistry and presumed PDF models for premixed turbulent combustion. *Combust. Flame*, **146**, 665-673.

Bray, K.N.C., and Moss, J.B. 1977. A unified statistical model of the premixed turbulent flame. *Acta Astronaut.*, **4**, 291-319.

Catlin, C.A., and Lindstedt, R.P. 1991. Premixed turbulent burning velocities derived from mixing controlled reaction models with cold front quenching. *Combust. Flame*, **85**, 427-439.

Chen, Y.-C., and Bilger, R. 2002. Experimental investigation of three-dimensional flame-front structure in premixed turbulent combustion. *Combust. Flame*, **131**, 400-435.

Chen, Y.-C., Peters, N., Schneemann, G.A., Wruck, N., Renz, U., and Mansour, M.S. 1996. The detailed flame structure of highly stretched turbulent premixed methane-air flames. *Combust. Flame*, **107**, 223-244.

Cho, P., Law, C., Cheng, R.K., and Shepherd, I. 1989. Velocity and scalar fields of turbulent premixed flames in a stagnation flow. *Proc. Combust. Inst.*, **22**, 739-745.

Cho, P., Law, C., Hertzberg, J.R., and Cheng, R.K. 1988. Structure and propagation of turbulent premixed flames stabilized in a stagnation flow. *Proc. Combust. Inst.*, **21**, 1493-1499.

Corvellec, C. 1998. Etude numérique et analytique des caractéristiques propagatives d'une zone de combustion se développant au sein d'un écoulement turbulent prémélangé modélisé par une approche de type flamelette. PhD thesis. University of Poitiers, Poitiers, France.

Dong, H.Q., Robin, V., Mura, A., and Champion, M. 2013. Analysis of algebraic closures of the mean scalar dissipation rate of the progress variable applied to stagnating turbulent flames. *Flow Turbul. Combust.*, **90**, 301-323.

Kha, K.N.Q., Robin, V., Mura, A., and Champion, M. 2016. Implications of laminar flame finite thickness on the structure of turbulent premixed flames. *J. Fluid Mech.*, **787**, 116-147.

Kolmogorov, A.N., Petrovkii, I., and Piskunov, N. 1937. A study of the diffusion equation with increase in the amount of substance and its application to abiology problem. *Bull. Univ. Moscow Ser. Int. A*, **1**, 116.
<figure><img src="image_57.png" /></figure>

Launder, B.E. 1976. Heat and mass transport. In P. Bradshaw (Ed.), *Turbulence, Topics in Applied Physics*, vol. 12, Springer Verlag, Heidelberg, Germany, pp. 231–287.

Lee, E., and Huh, K.Y. 2004. Zone conditional modeling of premixed turbulent flames at a high Damköhler number. *Combust. Flame*, **138**, 211–224.

Mazellier, N., Danaila, L., and Renou, B. 2010. Multi-scale grid injection: A new tool to generate intense homogeneous and isotropic turbulence for combustion. *J. Turbul.*, **11**, N43.

Mura, A., Galzin, F., and Borghi, R. 2003. A unified PDF-flamelet model for turbulent premixed combustion. *Combust. Sci. Technol.*, **175**, 1573–1609.

Peters, N. 1986. Laminar flamelets concepts in turbulent combustion. *Proc. Combust. Inst.*, **21**, 1231–1256.

Peters, N. 1992. A spectral closure for premixed turbulent combustion in the flamelet regime. *J. Fluid Mech.*, **242**, 611–629.

Robin, V., Mura, A., and Champion, M. 2011. Direct and indirect thermal expansion effects in turbulent premixed flames. *J. Fluid Mech.*, **689**, 149–182.

Sabelnikov, V.A., Corvellec, C., and Bruel, P. 1998. Analysis of the influence of cold front quenching on the turbulent burning velocity associated with an eddy-break-up model. *Combust. Flame*, **113**, 492–497.

Sabelnikov, V.A., and Lipatnikov, A.N. 2011. A simple model for evaluating conditioned velocities in premixed turbulent flames. *Combust. Sci. Technol.*, **183**, 588–613.

Williams, F.A. 1976. Criteria for existence of wrinkled laminar flame structure of turbulent premixed flames. *Combust. Flame*, **26**, 269–270.

Williams, F.A. 1985. *Combustion Theory*, 2nd ed. Benjamin Cummings, San Francisco, CA.