A DNS Study of Closure Relations for Convection Flux Term in Transport Equation for Mean Reaction Rate in Turbulent Flow

The present work aims at modeling the entire convection flux \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\overline {\rho \mathbf {u}W}$\end{document}ρuW¯ in the transport equation for a mean reaction rate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\overline {\rho W}$\end{document}ρW¯ in a turbulent flow, which (equation) was recently put forward by the present authors. In order to model the flux, several simple closure relations are developed by introducing flow velocity conditioned to reaction zone and interpolating this velocity between two limit expressions suggested for the leading and trailing edges of the mean flame brush. Subsequently, the proposed simple closure relations for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\overline {\rho \mathbf {u}W}$\end{document}ρuW¯ are assessed by processing two sets of data obtained in earlier 3D Direct Numerical Simulation (DNS) studies of adiabatic, statistically planar, turbulent, premixed, single-step-chemistry flames characterized by unity Lewis number. One dataset consists of three cases characterized by different density ratios and is associated with the flamelet regime of premixed turbulent combustion. Another dataset consists of four cases characterized by different low Damköhler and large Karlovitz numbers. Accordingly, this dataset is associated with the thin reaction zone regime of premixed turbulent combustion. Under conditions of the former DNS, difference in the entire, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\overline {\rho {u}W}$\end{document}ρuW¯, and mean, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\tilde {u}\overline {\rho W}$\end{document}ũρW¯, convection fluxes is well pronounced, with the turbulent flux, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\overline {\rho u^{\prime \prime }W^{\prime \prime }}$\end{document}ρu′′W′′¯, showing countergradient behavior in a large part of the mean flame brush. Accordingly, the gradient diffusion closure of the turbulent flux is not valid under such conditions, but some proposed simple closure relations allow us to predict the entire flux \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\overline {\rho \mathbf {u}W}$\end{document}ρuW¯ reasonably well. Under conditions of the latter DNS, the difference in the entire and mean convection fluxes is less pronounced, with the aforementioned simple closure relations still resulting in sufficiently good agreement with the DNS data.

pronounced, with the turbulent flux, ρu W , showing countergradient behavior in a large part of the mean flame brush. Accordingly, the gradient diffusion closure of the turbulent flux is not valid under such conditions, but some proposed simple closure relations allow us to predict the entire flux ρuW reasonably well. Under conditions of the latter DNS, the difference in the entire and mean convection fluxes is less pronounced, with the aforementioned simple closure relations still resulting in sufficiently good agreement with the DNS data.

Nomenclature b
Model constant c Combustion progress variable c 1 , c 2 Boundaries of the reaction zone D Molecular diffusivity of a species Da th = τ t /τ f Damköhler   Subscripts and superscripts q Reynolds-averaged value of a quantity q q = ρq/ρ Favre-averaged value of a quantity q q Either Reynolds or Favre-averaged value of a quantity q q|c 1 < c < c 2 Value of a quantity q, conditionally averaged at c 1 < c < c 2 q f = q / Value of a quantity q conditioned to flamelets q r = ρqW / Value of a quantity q conditioned to reaction zone ρW q = q −q Fluctuations with respect to Favre-averaged values b Burned f Flamelet r Reaction zone T Turbulent u Unburned

Introduction
The critical point of the turbulent combustion theory stems from averaging reaction rates subject to fluctuations in the local temperature T and concentrations. This is an issue of severe importance, because (i) the rates of reactions that control heat release depend on T in a highly non-linear manner, (ii) the magnitudes of the temperature fluctuations are typically large in turbulent flames, and (iii) these fluctuations exhibit a wide range of length and time scales. As reviewed elsewhere [1][2][3][4][5][6][7], among the most widely used approaches to solving this highly non-linear and multiscale problem, there are methodologies that deal with the Reynolds-averaged or filtered form of (i) the following transport equation for a combustion progress variable c, which is assumed to completely characterize the state of a reacting mixture, and (ii) an extra transport equation for either mean scalar dissipation rateχ = D∇c · ∇c [6,8] or mean flame surface density¯ = |∇c| [3,4,[9][10][11], with the mean rate W being hypothesized to be linearly related to eitherχ or¯ . More specifically,ρ W is often considered to be equal to ρ u S L [3,4,[9][10][11] and the following linear relation W =χ/ (2c m − 1) was derived by Bray [12] by assuming that the probability of finding intermediate (between unburned and fully burned) states of a reacting mixture is much less than unity. Here, t is time, u is the flow velocity vector, ρ is the mixture density, the combustion progress variable may be defined as follows c = (T − T u )/ (T b − T u ) in the simplest case of adiabatic burning of an equidiffusive mixture (i.e. D n = D for all species n = 1, . . . , N) characterized by unity Lewis number (i.e. Le = κ/D = 1) and a low Mach number, κ is the heat diffusivity of the mixture, W is the rate of production of the combustion progress variable in the flame, S L is the laminar flame speed, c m = ρcW /ρW is commonly assumed to be constant [6,12], q andq = ρq/ρ are the Reynolds and Favreaveraged values of a quantity q, respectively, and subscripts u and b designate unburned and burned mixture, respectively. It is worth noting, however, that even the precise knowledge ofχ or¯ does not allow us to precisely evaluate W , because the linear relations between these three quantities are just assumptions, which are best justified if unburned and fully burned gases are separated by a thin zone (flamelet) that retains the structure of the laminar flame. Nevertheless, recent analysis [13,14] of Direct Numerical Simulation (DNS) data shows that, even under such conditions (i.e. in the flamelet regime of premixed turbulent combustion [1][2][3][4][5][6][7][8]), the linear relations do not hold in certain flame zones. For instance, a ratio ofρ W / ρ u S L¯ can be significantly larger than unity at c > 0.8 [14] and c m strongly varies at c 1 [13] (it is worth remembering that the leading edge of a premixed turbulent flame brush may play a crucial role in the flame propagation, as discussed in detail elsewhere [7]). Accordingly, it would be of interest to straightforwardly evaluate the mean rate W by solving an appropriate transport equation without invoking extra assumptions regarding a relation between W and χ or¯ .
This requirement is addressed in recent papers by the present authors [15,16] where the following transport equations for the instantaneous and mean rates of product creation, i.e. W and W have been derived starting from Eq. 1 and assuming that W is solely controlled by c, i.e. W = W (c). This assumption holds, in particular, in the case of adiabatic burning, low Mach number, single-step chemistry, equal diffusivities of fuel and oxidant, and Le = 1.0. The same simplifications are commonly invoked by models that deal with transport equations forχ [6,8] or¯ [3,4,[9][10][11]. Here, q = q −q for any quantity q. The statistical behaviors of the transport equations of W and¯ have been compared elsewhere [16] in the context of RANS simulations and the interested reader is referred to Ref. [16] for further information.
It is worth noting that (i) counterparts of Eq. 2 can be derived in more challenging cases (e.g., Le = 1 or complex combustion chemistry) and (ii) Eq. 2 may also be filtered [16] for Large Eddy Simulation (LES) of turbulent flames. Nevertheless, before addressing more complicated problems, it is worth analyzing the newly introduced approach [15,16] from the simplest case by thoroughly investigating the basic features of Eq. 3.
In order to apply Eq. 3 to simulations of premixed turbulent flames, closure relations for the third term on the Left Hand Side (LHS) and for the three terms on the Right Hand Side (RHS) should be developed. In this regard, the biggest challenge consists of the fact that the magnitudes of the second and third terms on the RHS are much higher than the magnitudes of other terms under typical conditions, whereas the signs of the two dominant terms are opposite [15,16]. Accordingly, if separate closure relations are developed for each of the two dominant terms, even small uncertainties that stem from the closure relations can yield large residuals for Eq. 3. This challenging problem was resolved in Ref. [15], where the following relation was proposed to jointly close the three terms on the RHS of Eq. 3. Then, Eq. 3 reads and Eq. 5 was validated [15,16] with respect to the DNS data, which will be discussed later. Here, ṡ|c 1 < c < c 2 denotes stretch rateṡ = ∇ · u − nn : ∇u + S d ∇ · n conditioned to the reaction zone, which is bounded by c 1 and c 2 such that ρW (c 1 ) = ρW (c 2 ) = max {ρW (c)}/2, the unit vector n = −∇c/ |∇c| is locally normal to the instantaneous flame surface, and S d = [∇ · (ρD∇c) + ρW ] / (ρ |∇c|) is the local displacement speed. Nevertheless, two terms ρu W and ṡ|c 1 < c < c 2 should still be modeled in Eq. 5. As the discussed approach was put forward very recently [15,16], the present communication is restricted to analyzing a single unclosed term in Eq. 5. This contributes to the ultimate goal to develop closure relations for the entire Eq. 3 in order to subsequently use it in RANS and LES studies of flames investigated in experiments. Thus the particular goal of the present work is to assess several simple closure relations for the former (turbulent flux) term by analyzing two sets of DNS data obtained from statistically planar flames associated with the flamelet and thin reaction zone regimes [1] of premixed turbulent combustion.
It is worth noting that assumptions invoked in the following to close the turbulent flux term are more clear if they are applied to the entire convection term ρuW . Accordingly, in the rest of the paper, we will address closure relations for the entire convection term, whereas the counterpart closure relation for the turbulent transport term can be obtained using the following identity ρu W = ρuW −ρũ W and the closure relation for ρuW . Such an approach is fully justified, because the final goal consists in evaluating the sum of the mean convectionρũ W and turbulent transport ρu W terms on the RHS of Eq. 5, rather than each term separately.
It is also worth noting that the three terms on the LHS of Eq. 5 are similar to the unsteady, mean convection, and turbulent transport terms on the LHSs of the transport equations for χ [6,8] or¯ [3,4,[9][10][11], but the RHSs of the latter two equations involve several unclosed terms, contrary to a single unclosed term on the RHS of Eq. 5. This single term (i) differs substantially from any term on the RHS of a transport equation forχ [6,8] or¯ [3,4,[9][10][11] and (ii) offers an opportunity to attain a new insight into flame-turbulence interaction, as discussed in detail elsewhere [15,16].
In the next section, several simple closure relations for the entire convection flux ρuW will be suggested. In the third section, the attributes of DNS data will be reported. Assessment of the aforementioned closure relations using the DNS data will be discussed in the fourth section, followed by conclusions.

Simple Closure Relations
Utilizing analogy with quantities q f ≡ q /¯ conditioned to flamelets within a mean flame brush [3,4], let us consider q r ≡ qρW /ρW to be a value of a quantity q, conditioned to reaction zones. Then, in the statistically planar case addressed in the present paper, a closure relation for the term ρuW = ρW u r is required. Here, u is the x-component of the flow velocity vector u and the x-axis is normal to the mean flame brush.
Let us study whether or not the conditioned velocity u r can be evaluated invoking models proposed to determine the conditioned velocity u f . It is worth noting that modeling of the mean rate ρW is beyond the scope of the present paper and this quantity will be extracted from the DNS data in the following. The reader interested in modeling ρW is referred to [1][2][3][4][5][6][7][8] and references quoted therein.
As a starting point, let us invoke the following simple relation where σ = ρ u /ρ b is the density ratio, u u and u b are velocities conditioned to unburned and burned mixture, respectively, and c designates either the Reynolds-averaged, c, or the Favre-averaged,c, combustion progress variable, as will be discussed later. Equation 6 was earlier proposed to model the conditioned velocity u f [17]. This simple model is based on a paradigm of infinitely thin flamelets and on the following reasoning. First, if flamelets are infinitely thin, events associated with arrival of a piece of a flamelet to the trailing edge of a mean flame brush should be accompanied by arrival of unburned gas to the same volume and vice versa. Therefore, u f → u u at c → 1, at least if the density is constant. Second, based on similar arguments, we arrive at Thus, in such a case, Eq. 6 is nothing more than a linear interpolation between two limiting conditions. Nevertheless, this simple model is well supported by the results of a recent DNS study of self-propagation of an infinitely thin interface in constant-density turbulence [18].
If flamelet thickness does not vanish and σ > 1, we may expect that the aforementioned limit (c → 1 or c → 0) relations are approximate in the best case, while u f > |u u | at c → 1 and u f < |u b | at c → 0 due to an increase in flow velocity from the unburned to the burned sides of flamelets. Equation 6 addresses such effects in part by multiplying u b with σ −1 . This modification offers an opportunity to allow for the influence of thermal expansion on the limit (c → 0) relation between u b and u f,u conditioned to the leading edge of flamelets, but u f,u < u f due to the aforementioned thermal expansion effects. Accordingly, Eq. 6 exhibits worse agreement with DNS data obtained from flames characterized by σ > 1 [19,20] when compared to the DNS data obtained in the case of σ = 1 [18].
It is also worth noting that, first, Eq. 6 is a linear interpolation between two limit points, and thus both c = c and c =c may be used for the interpolation. Second, in Eq. 6, c and (1 − c ) may be substituted with bridging functions f ( c ) and g ( c ) such that f ( c → 0) → 0 and f ( c → 1) → 1 whereas g ( c → 0) → 1 and g ( c → 1) → 0. Third, the conditioned velocities u u and u b may be multiplied with factors b u (σ ) ≥ 1 and b b (σ ) ≤ 1 that allow for u f > |u u | at c → 1 and u f < |u b | at c → 0 due to the influence of thermal expansion on velocity within flamelets. Thus, Eq. 6 may be generalized as follows For instance, because a reaction zone is substantially thinner than a flamelet that contains the zone [7,21], variations in the flow velocity within the zone along the normal to it could be neglected to the leading order. If the variations in the local flamelet structure are neglected, then one obtains ρ u (u · n) u = ρ r (u · n) r = ρ b (u · n) b due to the local mass conservation. Consequently, b u = ρ u /ρ r and b b = ρ b /ρ r . Here, ρ r is the density within the reaction zone, e.g. the density conditioned to the peak value of ρW .
Strictly speaking, even if f ( c ) , g ( c ) , b u , and b b are known, Eq. 7 does not solve the problem of modeling the velocity u r , because the conditioned velocities u u and u b also require closure relations. However, at least under conditions associated with the flamelet regime of premixed turbulent combustion, the two velocities may be evaluated (i) using the following well-known Bray-Moss-Libby (BML) equations [22,23] (ii) solving the Favre-averaged Navier-Stokes equations in order to determineũ, and (iii) invoking some of available models of the flux ρu c , which are reviewed elsewhere [24,25]. Moreover, models for u u have also been developed, e.g. [26][27][28]. Accordingly, the conditioned velocities u u and u b are considered to be known in the present paper and are extracted from the DNS data. Because turbulent scalar fluxes are often modeled invoking an assumption of gradient diffusion (e.g., to the best of the present authors' knowledge, solely such models have yet been applied to terms ρu and ρu χ in transport equations for mean flame surface density¯ and mean scalar dissipation rateχ , respectively, [3,4,6]), the following closure relations were also tested. Here, W designates either W or W , D t = C μk 2 /ε is turbulent diffusivity, C μ is a constant,k = ρu · u /2ρ andε = 2ρνS ij S ij /ρ are the Favre-averaged turbulent kinetic energy and its dissipation rate, respectively, which were extracted from the DNS, ν is the kinematic viscosity of the mixture, S ij = 0.5 ∂u i /∂x j + ∂u j /∂x i is the rate-ofstrain tensor, u i is the i-th component of the velocity vector u, and summation convention applies to repeated indexes i and j .

DNS Attributes
In order to assess closure relations given by Eq. 7 or 9, we analyzed DNS data obtained earlier in two sets of simulations, which were consistent with the framework of the present study (adiabatic burning, low Mach number, single-step chemistry, unity Lewis number). One DNS database (flames H, M, and L) was created by Nishiki et al. [29,30] by simulating weakly turbulent combustion in the flamelet regime and was analyzed in a number of recent papers [13-16, 27, 31-38]. Another DNS database (flames B, C, D, and E) was created by Chakraborty et al. [39,40] by simulating combustion in small-scale intense turbulence (the thin-reaction-zone regime [1] of premixed burning) and was also analyzed in a number of recent papers cited elsewhere [15,16,27,41]. Because the DNS attributes were already discussed in detail in the literature, we will restrict ourselves to a very brief summary of the simulations. In both sets of DNS studies, unsteady 3D balance equations for mass, momentum, energy, and mass fraction of the deficient reactant were numerically solved and the ideal gas state equation was used. Combustion chemical mechanism was simplified by a single Arrhenius type irreversible chemical reaction. The Lewis and Prandtl numbers were taken to be equal to 1.0 and 0.7, respectively. Other basic flame characteristics are reported in Table 1 Table 1 are the turbulence characteristics at t = 0) and averaging was performed over transverse planes at t/τ t = 2 (cases D), 3 (cases C and E), or 4.34 (case B), which amount to one chemical time scale τ f . At those instants, the values of u /S L decayed by 61% (case B), 45% (case C), 24% (case D), and 34% (case E) in comparison to the initial values reported in Table 1.
In cases H, M, and L, homogeneous isotropic turbulence was generated in a separate box, was injected into the computational domain at x = 0, and decayed along the direction x (reported in Table 1 are the turbulence characteristics at the inlet). At the leading edges of the mean flame brushes, which were located about 1 mm downstream of the inlet in flames H, M, and L, turbulent kinetic energy was decreased by a factor of about two or three in case H or L, respectively. Accordingly, the Damköhler numbers evaluated at the leading edges were larger than Da th reported in Table 1 by a factor of about 1.5. In case L, the turbulence decay was more pronounced, because the inlet bulk flow velocity was lower (0.8 m/s) than in case H (1.15 m/s). The point is that, in each case, the bulk velocity was equal to the mean turbulent flame speed in order to retain the flame brush within the computational domain for a sufficiently long time [29,30]. In spite of the turbulence decay, the computed mean turbulent flame speeds were higher than S L by a factor of about two in cases H, M, and L [29,30].
After a transition time, which was longer than 1.5τ t and much longer than τ f , both turbulent burning velocity U T = ρ −1 u 0ρ W dx and mean flame brush thickness δ T = 1/max |∇c|, evaluated in cases H, M, and L, oscillated around statistically steady values [38]. Accordingly, averaging of various quantities q (x, t) was performed over transverse planes and over time after the transition period and during the statistically steady period, which was longer than τ t and much longer than τ f . Subsequently, the obtained dependencies of q (x) were transformed to dependencies of q (c) exploiting the monotonicity of the computed axial profiles c (x) of the Reynolds-averaged combustion progress variable.
All quantities reported in the rest of the paper are normalized using ρ u if applicable, i.e. ρ will designate density normalized using the unburned gas density and ρ u = 1 in the following.

Weakly turbulent flames H, M, and L
Let us begin with discussing data obtained from flames H, M, and L, associated with the flamelet regime of premixed turbulent combustion. It is worth noting that these data are analyzed in the coordinate framework attached to the mean flame brush such that dc/dx ≥ 0. Figure 1 shows that the conditioned velocity u r = ρuW /ρW introduced above is indeed very close to the velocity conditioned straightforwardly to the reaction zone determined with the following constraint ρW > max {ρW }/2, cf. black symbols and red dashed lines.
Second, at c > 0.8, the entire ρuW and meanρũ W fluxes are very close to one another, thus, indicating that the turbulent flux ρu W is negligible when compared toρũ W under conditions of the DNS cases considered here.
Third, the entire flux ρuW =ρũ W + ρu W is significantly larger than the mean fluxρũ W at c < 0.8, cf. black symbols and blue dotted-dashed lines, with the difference being increased by σ , cf. cases H, M, and L. Therefore, ρu W (c < 0.8) > 0, whereas Eq. 9 with D t > 0 yields the opposite sign of ρuW − ρũ W at c < 0.  especially, H is expected, because the flux ρu c shows countergradient behavior in these weakly turbulent flames [30].
As the gradient-diffusion closure given by Eq. 9 is contradicted by the present DNS data at c < 0.5, let us test Eq. 7. For this purpose various bridging functions f ( c ) and g ( c ) and various factors b u and b b were invoked. In particular, In the following, we will restrict ourselves to discussing a few closure relations that yield the best agreement with the DNS data on the entire flux ρuW .
Red double-dotted-dashed lines in Fig. 3 show that even the simplest linear interpolation u r =cu u + (1 −c) u b yields acceptable agreement with the DNS data. At c < 0.5 in cases M and, especially, H, the agreement is substantially improved by multiplying the conditioned velocity u b with a factor of b b = ρ b /ρ r , see blue dotted-dashed lines, as suggested in the second section. However, the latter model yields slightly worse results at high c values. The agreement with the DNS data at large c can be improved by substitutingc with a larger bridging function c ≥c in the first termcu u on the RHS in order to increase the magnitude of this term, but such a modification, see violet dashed lines, yields worse results at lower c.
A fact that substitution of a factor of b b = ρ b /ρ r with unity makes the agreement with the DNS data notably worse at c < 0.5 in cases M and, especially, H, cf. curves 3 and 1, is worth further discussing. On the one hand, in that range of the mean flame brush, u r is mainly controlled by the u b -term in Eq. 7, whereas the u u -term plays a minor role, because (i) u u is significantly lower than u b under conditions of the analyzed DNS, see Fig. 4, and (ii)c < (1 −c) if c < 0.5. Accordingly, a comparison of curves 1 and 3 with the DNS data supports reasoning for introducing a factor of b b = ρ b /ρ r into Eq. 7, as discussed in the second section. On the other hand, for the same reasoning, a factor of b u = ρ u /ρ r could also be introduced into the u u -term, but such a model (not shown) significantly overestimates the DNS data, especially at large c. Therefore, it is worth discussing why u r at c → 1 is more close to u u than to u u ρ u /ρ r .
In order to clarify the issue, the following feature of premixed turbulent flames appears to be of importance. As discussed in detail elsewhere [38,[42][43][44], local flamelet structure is strongly perturbed at large c, where the probability of finding highly curved tips of unburned mixture fingers [38] or cusps [42][43][44] is substantial. Accordingly a simple relation of ρ u (u · n) u = ρ r (u · n) r = ρ b (u · n) b , which is valid in unperturbed flamelets and is invoked to introduce b u = ρ u /ρ r and b b = ρ b /ρ r in the second section, does not seem to hold near the finger tips or cusps, as the local flamelet curvature is highly negative and the difference in u u , u r , and u b is significantly less pronounced at c → 1 when compared to the unperturbed laminar flame. Indeed, Fig. 4 shows that in the trailing part (c > 0.9) of mean flame brush, the axial velocity of unburned gas is strongly increased, see blue dotted-dashed lines. As discussed elsewhere [38], this increase in u u is more pronounced at larger density ratios (case H) and is controlled by the axial pressure gradient induced due to combustion in surrounding flamelets. Consequently, at c → 1, the conditioned velocity u u reaches values sufficiently close to the mean velocity u in products, which is equal to u b (c → 1). For instance, a ratio of u b /u u is less than two at c → 1 in flame H, whereas the magnitude |u · n| of flow velocity is increased by a factor of 6.5 in the counterpart laminar flame in the coordinate framework attached to it. Thus, at c → 1, the local axial acceleration of the flow within flamelets is weakly pronounced and the difference in u r and u u is significantly smaller when compared to the laminar flame. Unless a model or theory capable for predicting a ratio of u r /u u at c → 1 is developed, the use of the simplest assumption of u r (c → 1) → u u (c → 1) appears to be a better choice when compared to invoking an alternative assumption of u r (c → 1) → u u (c → 1) ρ u /ρ r .
Nevertheless, it is worth remembering that the real u r (c → 1) should be between the two limiting expressions (more close to the former one) and, therefore, u r (c → 1) > u u (c → 1). Substitution of the Favre-averagedc with a larger Reynolds-averaged c in the u u -term in Eq. 7, see curves 2 in Fig. 3, appears to be the simplest way to mimic the difference in u r (c → 1) and u u (c → 1), but such a modification does not solve the problem, as shown in Fig. 3.
Alternatively, the u u -term could involve an empirical factor b, which is tuned to improve agreement with the DNS data, see orange solid lines in Fig. 3. Such a method allows us to get very good agreement with the DNS data obtained from all three weakly turbulent flames H, M, and L at various c, with the required tuning of b = 1.4 ± 0.1 being more than modest. Figure 4 shows two more points that are worth noting. First, even at c → 0, a ratio of u b /u u or u r /u u is significantly lower than the counterpart quantity in the unperturbed laminar flame. However, this effect is unlikely to stem from weak local acceleration of the flow in flamelets. Indeed, if we consider a piece of flamelet x f (y, z, t) that (i) is locally normal to the x-axis, i.e. ∂x f /∂y = ∂x f /∂z = 0, (ii) moves towards the leading edge of the mean flame brush in the coordinate framework attached to it, i.e. ∂x f /∂t < 0,and (iii) retain the structure of the unperturbed laminar flame, i.e. the local flow acceleration is as strong as in the laminar flame, then, even in such a case, the axial local flow velocity at the burned side of this element should be less than σ S L (or equal to σ S L if ∂x f /∂t = 0). Accordingly, we could assume that u b (c → 0) ≤ σ S L for the above purely kinematic reasoning. This limit value is equal to 4.5, 2.6, or 1.0 m/s in case H, M, or L respectively, and is sufficiently close to u b (c → 0) plotted in Fig. 4, see red double-dotted-dashed lines.
Second, comparison of black solid and orange dashed lines in Fig. 4 shows that u r does tend to ρ b u b /ρ r at c → 0, in line with simple reasoning discussed in the second section.

Thin reaction zone regime flames B-E
Let us consider DNS data obtained from four turbulent flames B-E characterized by low Damköhler numbers and associated with thin reaction zone regime [1] of premixed turbulent combustion. Again, the data are analyzed in the coordinate framework attached to the mean flame brush such that dc/dx ≥ 0. It is also worth noting that results obtained in those simulations were normalized using the laminar flame speed and thickness.
First, similarly to Fig. 1, Fig. 5 also shows that the conditioned velocity u r = ρuW /ρW is very close to the velocity conditioned straightforwardly to the reaction zone determined with the following constraint ρW > max {ρW }/2, cf. black symbols and red dashed lines, respectively.
Second, similarly to Fig. 1, Fig. 5 also shows that the mean convection fluxρũ W and the entire flux ρuW are close to one another in the trailing parts (0.6 < c) of all four mean flame brushes, cf. blue dotted-dashed lines with black symbols, respectively.
Third, at lower c, the difference inρũ W and ρuW is most pronounced in cases D and E, characterized by the highest ratios of u /S L , see Table 1. On the contrary, in cases B and C, the difference between the two fluxes is low in spite of the facts that (i) these two flames are characterized by significantly larger ratios of u /S L when compared to flames H, M, and L, discussed earlier, but (ii) the difference inρũ W and ρuW is well pronounced in flame M, which is characterized by σ = 5.0 close to σ = 5.5 in cases B and C. Because flame M is also characterized by a significantly larger ratio of L/δ th when compared to flames B and C, we could assume that the DNS data obtained from all seven flames and considered all together indicate that a role played by the turbulent flux ρu W (i.e. relative magnitude of the difference inρũ W and ρuW ) is increased by the density ratio (flames H, M, and L), u /S L (flames B-E), and L/δ th (flames M, B, and C). Definitely, more simulations performed by independently varying u /S L and L/δ th are required to thoroughly validate such an assumption.
Fourth, in the range of c < 0.5, associated with d W /dx = (d W /d c) (d c/dx) > 0 and dρW /dx = dρW /d c (dc/dx) > 0, the mean convection fluxρũ W is slightly lower  Table 1, impedes the countergradient turbulent flux ρu W , as could be expected via analogy with the turbulent flux ρu c , which shows the countergradient (gradient) behavior at lower (larger) values of u / (σ S L ) [45]. Figure 6 shows that the two simplest expressions, i.e. u r =cu u +(1 −c) u b and u r = cu u + (1 −c) (ρ b /ρ r ) u b yield very reasonable agreement with the DNS data obtained from all four flames, with the latter equation, see blue double-dotted-dashed lines, doing a little better job when compared to the former one, see red dotted-dashed lines. It is worth noting that Fig. 6 reports DNS data obtained solely in the range of c where the number of sample points is sufficient in order to evaluate conditioned quantities in a statistically meaningful manner.
Contrary to flames H, M, and L, discussed earlier, multiplication of u u with a factor of b = 1 + 0.04σ = 1.42 results in a worse agreement with the DNS data at 0.5 < c. This difference between results obtained in the two sets of DNS databases could be either attributed to different DNS setups or associated with an eventual decrease in b with decreasing Damköhler number. Further target-directed DNSs performed by varying Da th in a wide range of values are definitely required to find the proper explanation.

Conclusions
The turbulent flux ρu W shows the countergradient behavior in the leading parts (c < 0.5) of turbulent flame brushes characterized by sufficiently low ratios of u /S L (cases H, M, L, and B). An increase in u /S L results in transition from the countergradient to the gradient turbulent flux ρu W (cases D and E) at c < 0.5.
In the trailing parts (0.8 < c) of all seven investigated flame brushes, the magnitude of the turbulent flux ρu W is much smaller than the magnitude of the mean convection flux ρũ W .
The following simple closure relation which does not involve an empirical tuning parameter, offers an opportunity to reasonably well model the entire convection flux ρuW in the transport equation for the mean reaction rate under substantially different conditions (both weakly turbulent combustion associated with the flamelet regime and burning in small-scale intense turbulence, associated with the thin reaction zone regime). Under conditions associated with the flamelet regime of premixed turbulent combustion, even better agreement with the DNS data can be obtained by introducing a single empirical parameter b = 1.2 + 0.04σ = 1.4 ± 0.1 into the following closure relation ρuW =ρ cbu u + (1 −c) (ρ b /ρ r ) u b W , (11) which reduces to Eq. 10 if the tuning constant b is skipped. This parameter is required to mimic an increase in flow velocity from the unburned side of a flamelet to its reaction zone due to combustion-induced thermal expansion. A fact that the value of this model parameter is significantly lower than initially proposed ratio ρ u /ρ r of the gas densities in the unburned reactants and reaction zone, respectively, is attributed to damping the local flow acceleration due to strong negative curvature of flamelets that are highly probable at the trailing edge of a premixed turbulent flame brush.