**Andrei N. Lipatnikov 1,\*, Shinnosuke Nishiki <sup>2</sup> and Tatsuya Hasegawa <sup>3</sup>**


Received: 11 February 2019; Accepted: 4 March 2019; Published: 7 March 2019

**Abstract:** In this study, closure relations for total and turbulent convection fluxes of flame surface density and scalar dissipation rate were developed (i) by placing the focus of consideration on the flow velocity conditioned to the instantaneous flame within the mean flame brush and (ii) by considering the limiting behavior of this velocity at the leading and trailing edges of the flame brush. The model was tested against direct numerical simulation (DNS) data obtained from three statistically stationary, one-dimensional, planar, premixed turbulent flames associated with the flamelet regime of turbulent burning. While turbulent fluxes of flame surface density and scalar dissipation rate, obtained in the DNSs, showed the countergradient behavior, the model predicted the total fluxes reasonably well without using any tuning parameter. The model predictions were also compared with results computed using an alternative closure relation for the flame-conditioned velocity.

**Keywords:** turbulent flame; premixed turbulent combustion; countergradient transport; flame surface density; scalar dissipation rate; modeling; direct numerical simulations

## **1. Introduction**

Among various approaches to modeling premixed turbulent combustion, concepts that deal with a transport equation for the mean flame surface density (FSD) Σ = |∇c| [1,2] or the Favre-averaged scalar dissipation rate (SDR) <sup>χ</sup><sup>~</sup> <sup>=</sup> <sup>ρ</sup>D∇c·∇c/ <sup>ρ</sup> [3] have become particularly popular over the past decade [4–6]. Here, c is the combustion progress variable, which characterizes the state of a reacting mixture in a flame, and is equal to zero and unity in unburned reactants and equilibrium combustion products, respectively, ρ is the mixture density, D is the molecular diffusivity of c, and <sup>q</sup> and <sup>q</sup><sup>~</sup> <sup>=</sup> <sup>ρ</sup>q/ <sup>ρ</sup> designate the Reynolds and Favre (i.e., mass-weighted), respectively, mean values of an arbitrary quantity q. As reviewed elsewhere [4–6], both transport equations involve a number of terms that should be modeled, but the present communication is solely restricted to total transport terms **u**Σ = **u** Σ + **u** Σ and <sup>ρ</sup>**u**χ <sup>=</sup> <sup>ρ</sup> **u**<sup>~</sup> <sup>χ</sup><sup>~</sup> <sup>+</sup> <sup>ρ</sup>**u**"χ", turbulent contributions to which (i.e., **u** Σ and ρ**u**"χ") require closure relations. Here, u is the flow velocity vector, and <sup>q</sup> = q − q and q" = q − q<sup>~</sup> are fluctuations of q with respect to its Reynolds and Favre-averaged values, respectively.

As reviewed elsewhere [4–6], in applications, the turbulent transport terms are commonly modeled invoking a paradigm of gradient diffusion, for example, **u** Σ = −Dt∇ Σ and ρ**u**"χ" = − ρDt∇ χ~, where Dt > 0 is turbulent diffusivity. However, at least for the turbulent flux <sup>ρ</sup>**u**"c"

of the combustion progress variable c, such a paradigm is challenged by countergradient turbulent transport [7,8], that is, <sup>ρ</sup>**u**"c"·∇ c<sup>~</sup> > 0, which is well documented in various flames, as reviewed elsewhere [9–11]. For the FSD or SDR concept, an issue of eventually countergradient flux of **u** Σ or ρ**u**"χ" is of particular importance, because the corner-stone assumption of the two concepts (i.e., a linear relation between the mean mass rate <sup>W</sup> of creation of c and <sup>ρ</sup><sup>u</sup> <sup>Σ</sup> [1,2] or <sup>ρ</sup> χ<sup>~</sup> [3]) is best justified under conditions of the flamelet combustion regime [12], which are commonly associated with the countergradient transport [9–11]. Accordingly, the goal of the present work is to suggest and validate simple closure relations capable of predicting the total convection fluxes **u**Σ = **u** Σ + **u** Σ and <sup>ρ</sup>**u**χ <sup>=</sup> <sup>ρ</sup> **u**<sup>~</sup> <sup>χ</sup><sup>~</sup> <sup>+</sup> <sup>ρ</sup>**u**"χ" and, in particular, the countergradient turbulent fluxes **u** Σ and ρ**u**"χ" in the flamelet combustion regime.

In the next section, the developed model is presented and the direct numerical simulation (DNS) data used to access it are summarized. In the third section, results of the model validation using the DNS data are discussed and the model is compared with alternative approaches. Conclusions are drawn in the fourth section.

#### **2. Method of Research**

#### *2.1. Model*

Using velocities

$$<\langle \mathbf{u} \rangle\_{\mathbf{f}, \Sigma} = \langle \mathbf{u} \Sigma \rangle / \langle \Sigma \rangle \tag{1}$$

and

$$
\langle \mathbf{u} \rangle\_{\mathbf{f}, \chi} = \langle \mathfrak{p} \mathbf{u} \chi \rangle / \langle \mathfrak{p} \chi \rangle \tag{2}
$$

conditioned to flamelets, modeling of the total fluxes **<sup>u</sup>**Σ <sup>=</sup> **<sup>u</sup>**f,<sup>Σ</sup> <sup>Σ</sup> and <sup>ρ</sup>**u**χ <sup>=</sup> **<sup>u</sup>**f,<sup>χ</sup> <sup>ρ</sup> χ<sup>~</sup> in the transport equations for <sup>Σ</sup> and <sup>χ</sup>~, respectively, is reduced to modeling the conditioned velocities, whereas closure relations for the turbulent fluxes **u** Σ and ρ**u**"χ" are not required.

In order to model the conditioned velocities, let us start with a discussion of a constant-density flow. In such a case, the following simple closure relation,

$$
\langle \mathbf{u} \rangle\_{\mathbf{f}} = (\mathbf{1} - \langle \mathbf{c} \rangle) \langle \mathbf{u} \rangle\_{\mathbf{f}} + \langle \mathbf{c} \rangle \langle \mathbf{u} \rangle\_{\mathbf{u}},\tag{3}
$$

was proposed [13] and validated in a recent DNS study [14] of self-propagation of an infinitely thin interface in a constant-density turbulent flow. Here, **u**<sup>f</sup> is velocity conditioned to the interface, subscripts u and b designate quantities conditioned to reactants and products, respectively, which are separated by the interface. Note that the conditioned velocities **u**<sup>u</sup> and **u**<sup>b</sup> vary along the normal to the mean reaction wave. For instance, in the statistically 1D, planar, and stationary case, **u**<sup>u</sup> = **u**u(x) and **u**<sup>b</sup> = **u**b(x) if the x-axis is normal to the mean wave, with averaging being performed over transverse planes.

Equation (3) is based on the following simple physical reasoning. For an infinitely thin interface, its arrival to the trailing edge of a mean reaction-wave brush is always accompanied by the arrival of reactants to the same trailing point. Therefore, **u**f→ **u**<sup>u</sup> at c→1. Due to similar arguments, **u**f→ **u**<sup>b</sup> at c→0. Thus, Equation (3) is nothing more than a linear interpolation between two limiting relations.

Let us assume that the same interpolation holds both for **u**f,<sup>Σ</sup> and **u**f,χ, that is, **u**f,<sup>Σ</sup> = **u**f,<sup>χ</sup> = **u**f, in the case of a constant-density "flamelet" of a finite thickness. Such an assumption is based on a small thickness of a flamelet when compared to a mean flame brush thickness in a typical case.

If we allow for combustion-induced thermal expansion, effects due to the finite flamelet thickness appear to be of more importance due to significant velocity variations within thin flamelets. In particular, **u**f,b·**n** < **u**f·**n** < **u**f,u·**n**, because local flow acceleration within a flamelet occurs in the direction opposite to the direction of a unit vector **n** = −∇c/|∇c| that is locally normal to the flamelet (i.e., ( **u**f,b − **u**f,u)·**n** < 0). Here, **u**f,u and **u**f,b are evaluated on the reactant and product

sides of the flamelet along the local normal to it. Because the peak value of Σ(c) or χ(c) is shifted to the product side of a typical laminar premixed flame, the simplest way to allow for the discussed difference between the conditioned velocities **u**f,u, **u**f, and **u**f,b appears to consist in (i) neglecting the difference in **u**f,<sup>Σ</sup> or **u**f,<sup>χ</sup> and **u**f,b→ **u**<sup>b</sup> at c→0, but (ii) assuming that both | **u**f,b − **u**f,Σ| | **u**f,b − **u**f,u| and | **u**f,b − **u**f,χ| | **u**f,b − **u**f,u| at c→1. For example,

$$(\langle \mathbf{u} \rangle\_{\mathbf{f}, \mathbf{b}} - \langle \mathbf{u} \rangle\_{\mathbf{f}, \Sigma}) = \sigma^{-1}(\langle \mathbf{u} \rangle\_{\mathbf{f}, \mathbf{b}} - \langle \mathbf{u} \rangle\_{\mathbf{f}, \mathbf{u}}) \text{ or } (\langle \mathbf{u} \rangle\_{\mathbf{f}, \mathbf{b}} - \langle \mathbf{u} \rangle\_{\mathbf{f}, \chi}) = \sigma^{-1}(\langle \mathbf{u} \rangle\_{\mathbf{f}, \mathbf{b}} - \langle \mathbf{u} \rangle\_{\mathbf{f}, \mathbf{u}})$$

in order to avoid invoking tuning parameters. Here, σ = ρu/ρ<sup>b</sup> is the density ratio. Consequently,

$$
\langle \mathbf{u} \rangle\_{\mathbf{f}, \Sigma} \to \sigma^{-1} \langle \mathbf{u} \rangle\_{\mathbf{f}, \mathbf{u}} + (1 - \sigma^{-1}) \langle \mathbf{u} \rangle\_{\mathbf{f}, \mathbf{b}} \text{ and } \langle \mathbf{u} \rangle\_{\mathbf{f}, \chi} \to \sigma^{-1} \langle \mathbf{u} \rangle\_{\mathbf{f}, \mathbf{u}} + (1 - \sigma^{-1}) \langle \mathbf{u} \rangle\_{\mathbf{f}, \mathbf{b}} \text{ at } \langle \mathbf{c} \rangle \to 1.
$$

If we further assume that **u**f,b = **u**<sup>b</sup> at c→1, then,

$$
\langle \mathfrak{u} \rangle\_{\mathfrak{f}, \Sigma} \to \sigma^{-1} \langle \mathfrak{u} \rangle\_{\mathfrak{u}} + (1 - \sigma^{-1}) \langle \mathfrak{u} \rangle\_{\mathfrak{b}} \text{ and } \langle \mathfrak{u} \rangle\_{\mathfrak{f}, \chi} \to \sigma^{-1} \langle \mathfrak{u} \rangle\_{\mathfrak{u}} + (1 - \sigma^{-1}) \langle \mathfrak{u} \rangle\_{\mathfrak{b}} \text{ at } \langle \mathfrak{c} \rangle \to 1.
$$

Note that **u**f,u→ **u**<sup>u</sup> at c→1, because the unburned gas can arrive at the trailing edge of the mean flame brush only together with a flamelet. Finally, since **<sup>u</sup>**<sup>b</sup> <sup>=</sup> **<sup>u</sup>** <sup>=</sup> **<sup>u</sup>**<sup>~</sup> at <sup>c</sup> = 1, we arrive at

$$
\langle \mathfrak{u} \rangle\_{\mathfrak{f}, \Sigma} \to \sigma^{-1} \langle \mathfrak{u} \rangle\_{\mathfrak{u}} + (1 - \sigma^{-1}) \langle \mathfrak{u} \rangle^{\sim} \text{ and } \langle \mathfrak{u} \rangle\_{\mathfrak{f}, \chi} \to \sigma^{-1} \langle \mathfrak{u} \rangle\_{\mathfrak{u}} + (1 - \sigma^{-1}) \langle \mathfrak{u} \rangle^{\sim} \text{ at } \langle \mathfrak{c} \rangle \to 1
$$

or

$$(\langle \mathbf{u} \rangle\_{\mathbf{f}, \Sigma} - \langle \mathbf{u} \rangle^{\circ}) \to \sigma^{-1}(\langle \mathbf{u} \rangle\_{\mathbf{u}} - \langle \mathbf{u} \rangle^{\circ}) \text{ and } (\langle \mathbf{u} \rangle\_{\mathbf{f}, \chi} - \langle \mathbf{u} \rangle^{\circ}) \to \sigma^{-1}(\langle \mathbf{u} \rangle\_{\mathbf{u}} - \langle \mathbf{u} \rangle^{\circ}) \text{ at } \langle \mathbf{c} \rangle \to 1.$$

Thus, in the case of σ > 1, we have the following two approximate limiting relations: **u**f,Σ→ **u**<sup>b</sup> or **<sup>u</sup>**f,χ→ **u**<sup>b</sup> at <sup>c</sup>→0 and ( **<sup>u</sup>**f,<sup>Σ</sup> − **u**~)→σ−1( **<sup>u</sup>**<sup>u</sup> − **u**~) or ( **<sup>u</sup>**f,<sup>χ</sup> − **u**~)→σ−1( **<sup>u</sup>**<sup>u</sup> − **u**~) at c→1. Then, linear interpolation results in

$$
\langle \mathbf{u} \rangle\_{\mathbf{f}, \Sigma} - \langle \mathbf{u} \rangle^{\sim} = \langle \mathbf{c} \rangle \sigma^{-1} (\langle \mathbf{u} \rangle\_{\mathbf{u}} - \langle \mathbf{u} \rangle^{\sim}) + (1 - \langle \mathbf{c} \rangle)(\langle \mathbf{u} \rangle\_{\mathbf{b}} - \langle \mathbf{u} \rangle^{\sim}) \tag{4}
$$

and

$$
\langle \mathbf{u} \rangle\_{\mathbf{f}, \chi} - \langle \mathbf{u} \rangle^{\sim} = \langle \mathbf{c} \rangle \sigma^{-1} (\langle \mathbf{u} \rangle\_{\mathbf{u}} - \langle \mathbf{u} \rangle^{\sim}) + (1 - \langle \mathbf{c} \rangle)(\langle \mathbf{u} \rangle\_{\mathbf{b}} - \langle \mathbf{u} \rangle^{\sim}). \tag{5}
$$

In the case of σ = 1 (constant density), Equations (4) and (5) reduce to Equation (3) and yield correct values of **u**f,<sup>Σ</sup> and **u**f,<sup>χ</sup> at least at the boundaries c→0 and c→1 of the mean "flame" brush.

It is worth noting that Equations (4) and (5) extend Equation (3) by allowing for effects due to density variations and finite flamelet thickness. Since such an extension was performed by invoking a few simple assumptions, the developed model required validation and the results of such a validation performed by analyzing DNS data will be reported in Section 3. As we will see in that section, the model works reasonably well, but some differences between the model predictions and the DNS data were observed, probably due to the invoked simplifications. It is also worth remembering that the conditioned velocities **u**<sup>u</sup> and **u**<sup>b</sup> in Equations (4) and (5) are not averaged over the flame brush volume, but vary along the normal to the mean flame brush.

Finally, using the following well-known Bray–Moss–Libby (BML) expressions [8,9,12],

<sup>ρ</sup>u(1 <sup>−</sup>c) = <sup>ρ</sup>(1 <sup>−</sup>c~), <sup>ρ</sup><sup>b</sup> <sup>c</sup> <sup>=</sup> <sup>ρ</sup> c~, (6)

$$
\langle \mathbf{u} \rangle^{\rightharpoonup} = (1 - \langle \mathbf{c} \rangle^{\rightharpoonup}) \langle \mathbf{u} \rangle^{\rightharpoonup} + \langle \mathbf{c} \rangle^{\rightharpoonup} \langle \mathbf{u} \rangle\_{\mathfrak{b} \rhd} \tag{7}
$$

and

$$
\langle \rho \mathbf{u}'' \mathbf{c}'' \rangle = \langle \rho \rangle \langle \mathbf{c} \rangle^{\text{\textquotedblleft}} (1 - \langle \mathbf{c} \rangle^{\text{\textquotedblright}}) (\langle \mathbf{u} \rangle\_{\text{\textquotedblleft}} - \langle \mathbf{u} \rangle\_{\text{\textquotedblleft}}),
\tag{8}
$$

*Fluids* **2019**, *4*, 43

we arrive at

$$
\begin{split}
\langle \mathbf{u} \rangle\_{\mathbf{f}, \Sigma} &= \langle \mathbf{u} \rangle\_{\mathbf{f}, \mathbf{\mathcal{X}}} = \langle \mathbf{u} \rangle^{\circ} - \langle \mathbf{c} \rangle \mathbf{\mathcal{F}}^{-1} \langle \rho \mathbf{u}'' \mathbf{c}'' \rangle / \left[ \{\rho\} (1 - \langle \mathbf{c} \rangle^{\circ}) \right] + (1 - \langle \mathbf{c} \rangle) \langle \rho \mathbf{u}'' \mathbf{c}'' \rangle / \left[ \{\rho\} \langle \mathbf{c} \rangle^{\circ} \right] \\ &= \langle \mathbf{u} \rangle^{\circ} - \langle \mathbf{c} \rangle^{\circ} \langle \rho \mathbf{u}'' \mathbf{c} \rangle^{\circ} / \left[ \{\rho\_{\mathbf{u}} (1 - \langle \mathbf{c} \rangle^{\circ}) \} + \langle \rho \mathbf{u}'' \mathbf{c}'' \rangle (1 - \langle \mathbf{c} \rangle^{\circ}) / \left[ \rho\_{\mathbf{u}} \langle \mathbf{c} \rangle^{\circ} \right] \right] \\ &= \langle \mathbf{u} \rangle^{\circ} + (1 - 2 \langle \mathbf{c} \rangle^{\circ}) \langle \rho \mathbf{u}'' \mathbf{c}'' \rangle / \left[ \{\rho\_{\mathbf{u}} \langle \mathbf{c} \rangle^{\circ} (1 - \langle \mathbf{c} \rangle^{\circ}) \} \right].
\end{split} \tag{9}
$$

The problem of modeling the turbulent flux ρ**u**"c" in Equation (9) is not specific to the FSD or SDR concept and should be resolved by any approach that deals with the transport equation for the Favre-averaged <sup>c</sup>~. Accordingly, there are various closure relations for <sup>ρ</sup>**u**"c", as reviewed elsewhere [9–11]. Therefore, modeling of this flux is beyond the scope of the present communication, that is, ρ**u**"c" is considered to be known here. Thus, Equation (9) is the final result of the present model. It is worth stressing that this closure relation does not involve any tuning parameter.

#### *2.2. Direct Numerical Simulations*

In order to test Equation (9), we analyzed DNS data obtained earlier by Nishiki et al. [15,16]. Because these data were used by different research groups in a number of papers [17–36], we will restrict ourselves to a brief summary of these compressible simulations.

The well-known unsteady 3D balance equations for mass, momentum, energy, and mass fraction Y of the deficient reactant were numerically solved. The ideal gas state equation was used. Combustion chemistry was reduced to a single reaction. The molecular transport coefficients were increased by the temperature T (e.g., the kinematic viscosity ν = νu(T/Tu) 0.7). The Lewis and Prandtl numbers were equal to 1.0 and 0.7, respectively. Accordingly, the mixture state was completely characterized with a single combustion progress variablec=1 − Y/Yu = (T − Tu)/(Tb − Tu).

The computational domain was a rectangular box Λ<sup>x</sup> × Λ<sup>y</sup> × Λz, with Λ<sup>x</sup> = 8 mm and Λ<sup>y</sup> = Λ<sup>z</sup> = 4 mm, and was resolved using a uniform rectangular (2Δx = Δy = Δz) mesh of 512 × 128 × 128 points. The x-axis was normal to the mean flame brush and was parallel to the direction of its propagation.

Homogeneous isotropic turbulence (the rms turbulent velocity u = 0.53 m/s, an integral length scale of turbulence L = 3.5 mm, the turbulent Reynolds number Ret = u L/ν<sup>u</sup> = 96 [15]) was generated in a separate box and was injected into the computational domain through the left boundary x = 0. In the computational domain, the turbulence decayed along the direction x of the mean flow. The flow was periodic in y and z directions.

At t = 0, a planar laminar flame (c = 0 at x = 0 and c = 1 at x = Λx) was embedded into statistically the same turbulence assigned for the velocity field in the entire computational domain. Subsequently, in order to keep the flame in the computational domain until the end t3 of the simulations, the mean inflow velocity, which was parallel to the x-axis, was increased twice, that is, U(0 ≤ t<t1)=SL < U(t1 ≤ t<t2) < U(t2 ≤ t), with U(t2 ≤ t) being approximately equal to the mean turbulent flame speed ST. Here, SL is the laminar flame speed.

Three cases signifying high (H), medium (M), and low (L) respectively, density ratios were simulated. The basic flame characteristics are reported in Table 1, where Kath = (u /SL) 3/2(L/δth) −1/2 and Dath = (SL/u )(L/δth) are the Karlovitz and Damköhler numbers, respectively, δth = (Tb − Tu)/max|∇T| is the thermal thickness of the laminar flame, and the mean turbulent burning velocity is equal to

$$\rho\_{\mathbf{u}} \langle \mathbf{U}\_{\mathrm{T}} \rangle = \left[ \Lambda\_{\mathrm{y}} \Lambda\_{\mathrm{x}} (\mathbf{t}\_{3} - \mathbf{t}\_{2}) \right]^{-1} \int \int \int \int \mathcal{W}(\mathbf{x}, \mathbf{t}) \mathbf{dx} d\mathbf{t}\_{\mathrm{y}} \tag{10}$$

with the mean turbulent flame speed ST being equal to UT in the considered statistically one-dimensional, planar case.


**Table 1.** Studied cases. H–high, M–medium, L–low.

As discussed in detail elsewhere [25], the present DNS conditions and results are fully consistent with the paradigm of the flamelet regime [12] of premixed turbulent combustion. Accordingly, the DNS data are particularly useful for the goal of the present communication. Indeed, (i) the DNS conditions and data are consistent with the BML, FSD, and SDR concepts, see References [25,29], and [25], respectively, (ii) the axial turbulent flux ρu"c" shows the countergradient behavior [16] due to combustion-induced thermal expansion effects [9–11], and (iii) the use of different density ratios offers an opportunity to vary the magnitude of these effects.

Results presented in the next section were averaged over transverse planes and over time t2 ≤ t<t3 (approximately 200 snapshots). During that time interval, the computed turbulent burning velocity and flame brush thickness oscillated around statistically steady values [26].

#### **3. Results and Discussion**

Results of the developed simple model validation are reported in Figures 1–6. Dotted lines in Figure 1a, Figure 2a, and Figure 3a show the total axial FSD flux uΣ extracted directly from the DNS data, whereas solid lines show the flux uf,<sup>Σ</sup> Σ computed using the model, that is, Equations (1) and (9), with <sup>c</sup>~, <sup>u</sup>~, <sup>ρ</sup>u"c", and <sup>Σ</sup> being extracted from the same DNS data. In all three cases H, M, and L, the agreement between the DNS and model results is very good. Similarly, Figure 1b, Figure 2b, and Figure 3b show the total axial SDR flux ρuχ extracted directly from the DNS data, whereas solid lines show the flux <sup>u</sup>f,<sup>χ</sup> ρχ computed using the model, that is, Equations (2) and (9), with <sup>c</sup>~, <sup>u</sup>~, ρu"c", and ρχ being extracted from the same DNS data. Again, the model is well validated in all three cases.

**Figure 1.** Total axial fluxes of (**a**) flame surface density (FSD) Σ and (**b**) scalar dissipation rate (SDR) χ. Case H.

**Figure 2.** Total axial fluxes of (**a**) FSD Σ and (**b**) SDR χ. Case M.

**Figure 3.** Total axial fluxes of (**a**) FSD Σ and (**b**) SDR χ. Case L.

**Figure 4.** Turbulent axial fluxes of (**a**) FSD Σ and (**b**) SDR χ. Case H.

**Figure 5.** Turbulent axial fluxes of (**a**) FSD Σ and (**b**) SDR χ. Case M.

**Figure 6.** Turbulent axial fluxes of (**a**) FSD Σ and (**b**) SDR χ. Case L.

It is worth noting, however, that the magnitudes of the total fluxes uΣ and ρuχ can mainly be controlled by the mean advection fluxes <sup>u</sup> Σ and <sup>ρ</sup> u<sup>~</sup> <sup>χ</sup>~, respectively, which do not require modeling. Accordingly, to better access Equation (9), turbulent fluxes were also evaluated as follows:

$$
\langle \langle \mathbf{u}' \Sigma' \rangle = \langle \mathbf{u} \Sigma \rangle - \langle \mathbf{u} \rangle \langle \Sigma \rangle,\tag{11}
$$

$$
\langle \rho \mathbf{u}'' \chi'' \rangle = \langle \rho \mathbf{u} \chi \rangle - \langle \rho \rangle \langle \mathbf{u} \rangle^{-} \langle \chi \rangle^{-}. \tag{12}
$$

Results plotted in Figures 4–6 show that Equation (9) predicts the axial turbulent fluxes u Σ and ρu"χ" very well in the largest parts of mean flame brushes, that is, at c < 0.8, in cases H and M, respectively, cf. dotted and solid lines in Figures 4a and 5b, respectively. Moreover, the flux ρu"χ" is well predicted in case L (see Figure 6b), whereas the peak value of the flux is slightly overestimated. As far as flux u Σ in cases M and L or flux ρu"χ" in case H is concerned, Equation (9) yields qualitatively correct results, but some quantitative differences between them and the DNS data are observed, probably due to a simplified treatment of the effects of density variations and finite flamelet thickness in Equations (4) and (5). Nevertheless, these differences do not impede quantitatively predicting the total fluxes <sup>u</sup> Σ and <sup>ρ</sup> u<sup>~</sup> <sup>χ</sup>~, in particular, because <sup>u</sup>Σ is close to the mean flux u Σ in cases M and L, cf. scales of ordinate axes in Figure 2a or Figure 3a and in Figure 5a or Figure 6a, respectively.

It is worth stressing that the turbulent axial fluxes u Σ and ρu"χ" show the countergradient behavior in all studied cases (i.e., u Σ d Σ/dx > 0 and ρu"χ"d ρχ/dx > 0), because d Σ/dx and d ρχ/dx are positive (negative) in the leading (trailing) part of the mean flame brushes (reactants on the left side). Consequently, the widely used concept of gradient diffusion yields wrong direction of the fluxes u Σ and ρu"χ" under conditions of the present DNSs.

It is also worth noting that an alternative model for evaluating the flamelet-conditioned velocity **u**<sup>f</sup> was suggested by Bidaux and Bray in an unpublished work (1994) and was subsequently used by Veynante et al. [37]. The model consists in the following linear interpolation:

$$
\langle \mathbf{u} \rangle\_{\rm f} = (1 - \mathbf{K}) \langle \mathbf{u} \rangle\_{\rm u} + \mathbf{K} \langle \mathbf{u} \rangle\_{\rm b}, \tag{13}
$$

where the constant K could "be related to the iso-c line used to define the flame location" [37]. Equations (7), (8), and (13) yield

$$
\langle \mathbf{u} \rangle\_{\mathbf{f}} = \langle \mathbf{u} \rangle^{\smile} + (\mathbf{K} - \langle \mathbf{c} \rangle^{\smile}) \langle \rho \mathbf{u}'' \mathbf{c}'' \rangle / \lg \langle \mathbf{c} \rangle \langle \mathbf{c} \rangle^{\smile} (1 - \langle \mathbf{c} \rangle^{\smile}) \text{l.} \tag{14}
$$

In the cited paper, the constant K was set equal to 0.5 and we are not aware of any alternative suggestion for K. If K = 0.5, Equation (14) can be rewritten in the following form:

$$
\langle \mathbf{u} \rangle\_{\mathbf{f}} = \langle \mathbf{u} \rangle^{\rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \rm \$$

which is similar to Equation (9), but the last term on the right-hand side (RHS) of Equation (15) is multiplied by a factor of R = ρu/2 ρ. This factor tends to 1/2 and σ/2 at the leading and trailing edges of the mean flame brush. It is equal to unity somewhere in the middle of the mean flame brush.

While Equations (13)–(15) were never applied to model turbulent fluxes of FSD and SDR, to the best of our knowledge, this gap can easily be filled by assuming that **u**f,<sup>Σ</sup> = **u**f,<sup>χ</sup> = **u**<sup>f</sup> and using Equations (1) or (2), (11) or (12), and (13)–(15). The outcomes of such tests of the model by Bidaux and Bray with K = 0.5 are shown in dashed and dotted-dashed lines in Figures 1–6. There are some differences between the results computed using Equations (13) and (15), because Equation (15) is obtained from Equation (13) using Equations (7) and (8), which are exact only in the case of an infinitely thin flame front. In the present DNSs, flamelets have finite thicknesses and, consequently, there are some (sufficiently small) differences between the left-hand side and right-hand side of Equation (7) or (8), extracted from the DNS data (e.g., see Figure 3b in Reference [25]).

The comparison of Equation (13) with Equation (9) shows that Equation (13) with K = 0.5 yields better results for flux ρu"χ" in case H (see Figure 4b), but worse results for flux u Σ in cases H and M (see Figures 4a and 5a). Similarly, when compared to Equation (9), Equation (15) yields better results for flux u Σ in cases M and L (see Figures 5a and 6a), but worse results for both u Σ and ρu"χ" in case H, see Figure 4. Moreover, for both fluxes, Equation (9) is superior to Equation (13) or (15) at the leading edge of the mean flame brush (i.e., at c→0). For instance, Equation (15) yields wrong direction of the fluxes at very low c in cases H and M, whereas Equation (13) underestimates the magnitudes | u Σ | and | ρu"χ"| at very low c in all cases. This difference between Equation (9) and Equation (13) or (15) appears to be of paramount importance, because the leading edge of premixed turbulent flame brush can play the crucial role in its propagation, as discussed in detail elsewhere [38–45]. Furthermore, Equation (13) or (15) yields the same wrong limit behavior of **u**f→ **u**u/2 + **u**b/2, both at the leading and trailing edges of the mean flame brush in the simplest case of the propagation of an infinitely thin, dynamically passive front in the constant-density turbulent flow. On the contrary, in such a case, Equation (9) obtained in the present work (i) reduces to Equation (3), which was well validated in Reference [14], and (ii) yields **u**f→ **u**<sup>b</sup> and **u**f→ **u**<sup>u</sup> at c→0 and c→1, respectively, in line with the simple physical reasoning discussed at the beginning of Section 2.1.

#### **4. Conclusions**

In this study, closure relations for total and turbulent fluxes of flame surface density Σ and scalar dissipation rate χ were developed. The relations did not involve a tuning parameter, very well predicted the total axial fluxes uΣ and ρuχ, and reasonably well predicted turbulent axial fluxes u Σ and ρu"χ" in three flames characterized by significantly different density ratios and

associated with (i) the flamelet regime of premixed turbulent burning, (ii) the validity of the FSD and SDR concepts, and (iii) the countergradient turbulent fluxes of Σ and χ. In the case of a constant density, the developed closure relation reduced to Equation (3), which was recently validated under conditions of gradient turbulent transport [14].

**Author Contributions:** Data curation, S.N. and T.H.; Formal analysis, A.N.L.; Methodology, A.N.L.; Investigation, A.N.L.; Software, S.N.; Writing—original draft preparation, A.N.L.; Writing—review and editing, A.N.L., S.N., and T.H.

**Acknowledgments:** The first author (Andrei N. Lipatnikov) gratefully acknowledges the financial support provided by the Combustion Engine Research Center (CERC).

**Conflicts of Interest:** The authors declare no conflict of interest.

## **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
