**Rixin Yu <sup>1</sup> and Andrei N. Lipatnikov 2,\***


Received: 28 January 2019; Accepted: 15 February 2019; Published: 19 February 2019

**Abstract:** Propagation of either an infinitely thin interface or a reaction wave of a nonzero thickness in forced, constant-density, statistically stationary, homogeneous, isotropic turbulence is simulated by solving unsteady 3D Navier–Stokes equations and either a level set (G) or a reaction-diffusion equation, respectively, with all other things being equal. In the case of the interface, the fully developed bulk consumption velocity normalized using the laminar-wave speed SL depends linearly on the normalized rms velocity u /SL. In the case of the reaction wave of a nonzero thickness, dependencies of the normalized bulk consumption velocity on u /SL show bending, with the effect being increased by a ratio of the laminar-wave thickness to the turbulence length scale. The obtained bending effect is controlled by a decrease in the rate of an increase δAF in the reaction-zone-surface area with increasing u /SL. In its turn, the bending of the δAF(u /SL)-curves stems from inefficiency of small-scale turbulent eddies in wrinkling the reaction-zone surface, because such small-scale wrinkles characterized by a high local curvature are smoothed out by molecular transport within the reaction wave.

**Keywords:** reaction waves; turbulent reacting flows; turbulent consumption velocity; bending effect; reaction surface area; molecular transport; direct numerical simulations

#### **1. Introduction**

In turbulent reacting flows, the so-called bending effect consists in decreasing the rate (dUT)/du of an increase in turbulent consumption velocity UT by the rms turbulent velocity u with increasing u , i.e., the second derivative of the function UT(u ) is negative. As illustrated in Figure 1, due to this effect, the curve plotted in an orange solid line is bent and, at high u , shows significantly lower consumption velocities when compared to the straight dashed red line. Since this basic phenomenon was documented, e.g., in premixed turbulent flames [1,2], it has been challenging the research community and different approaches to explaining and modeling the bending effect have been put forward.

The most recognized approach consists in highlighting the so-called stretch effect, i.e., variations in the local structure of reaction wave (e.g., flame) and the local consumption velocity uc, caused by turbulent stretching of the wave [3–5]. Here, uc is a properly normalized rate of production of a major reaction product, integrated along the local normal to a thin reaction zone, which is assumed to be inherently laminar within the framework of the discussed concept. The straightforward manifestation of the stretch effect consists in changing the mean value uc of the local consumption velocity with increasing u , followed by eventual local reaction extinction at sufficiently high u . A decrease in uc and the local reaction extinction can also affect the area AF of the reaction-wave surface, but this manifestation of the discussed mechanism is indirect, i.e., it is a consequence of the dependence of uc or uc on the local stretch rate or u , respectively. According to the theory of stretched laminar premixed

flames [6,7], which is well developed in the case of single-step chemistry and the asymptotically high activation energy E of the global combustion reaction, the dependence of uc on the local stretch rate s is controlled by differences in molecular diffusivities of a fuel, DF, oxygen, DO, and the heat diffusivity κ of the mixture. In the equidiffusive (DF = DO = κ) adiabatic case, uc does not depend on s if E tends to infinity, but uc can depend on s and the flame extinction by the stretch rate can occur at a finite activation energy [8]. The reader interested in further discussion of the stretch effect in premixed flames is referred to review paper [9].

**Figure 1.** A sketch of the bending effect for turbulent consumption velocity UT as a function of the rms turbulent velocity u .

Even if the mean uc is close to the speed SL of the unperturbed laminar reaction wave, the bending effect can occur. For instance, a recent Direct Numerical Simulation (DNS) study [10] has shown that the bending effect can be controlled by the mean flame surface area AF, i.e., the second derivative of the function AF(u ) can be negative in spite of uc ≈ SL. Under such conditions, the bending effect might be attributed to various physical mechanisms, e.g., statistical equilibrium between small-scale turbulent eddies and reaction rate [11], collisions of reaction waves [12,13], or smoothing of small-scale wrinkles of the reaction-wave surface due to its propagation [14].

While all the aforementioned approaches [3–5,11–14] were developed by studying premixed flames, they place the focus of consideration on the influence of turbulence on combustion, but disregard the back influence of the combustion on the turbulence. However, phenomena caused by combustion-induced thermal expansion can also contribute to the bending effect. For instance, small-scale turbulent eddies may be inefficient in wrinkling reaction-zone surface, because they disappear due to dilatation and an increase in the kinematic viscosity of the preheated mixture [15,16]. The reader interested in further discussion of the thermal expansion effects is referred to review papers [17,18].

It is also worth noting that Damköhler [19] has arrived at the following scaling UT∝SL(u L/κ) 1/2 by reducing the influence of very intense turbulence on reaction wave to enhancement of heat and mass transfer within the wave by turbulent eddies. Here, L is an integral turbulent length scale. From the purely mathematical viewpoint, this scaling results in the bending, but the physical mechanism hypothesized by Damköhler [19] is associated an increase in uc when compared to SL.

Finally, the bending effect can be pronounced differently in different reaction waves. For instance, levelling-off of UT(u )-curves, followed by a decrease in UT with increasing u , is well documented in expanding statistically spherical flames [1,2], but the positive dUT/du was obtained from statistically stationary flames at much higher values of u /SL [20,21].

Besides the reaction-wave-configuration effects, which are beyond the scope of the present study, the aforementioned physical mechanisms of the bending may be divided into three groups; (A) mechanisms that highlight differences in the molecular transport coefficients, i.e., the stretch effect [3–5], (B) mechanisms that highlight thermal expansion effects in flames, and (C) mechanisms that address equidiffusive flames, but disregard the thermal expansion effects [11–14,19]. The physical mechanisms of the bending may also be divided into two other groups; (a) mechanisms that highlight variations in local consumption velocity uc due to the stretch effect [3–5] or the transport enhancement [19] and (b) mechanisms that highlight the bending of AF(u )-curves.

To conclude this brief introduction, it is worth noting that the fact that various physical mechanisms of the bending effect are discussed in the literature does not mean that all relevant physical mechanisms have already been revealed.

The present communication does not aim at comparing all the aforementioned physical mechanisms. The goals of the communication are solely restricted to (i) supporting recent finding [10] that the bending of UT(u )-curves can be controlled by the bending of AF(u )-curves, (ii) comparing physical mechanisms from Group (C), and (iii) emphasizing a physical mechanism that controls the bending of AF(u )-curves under conditions of the present study, but has yet been outside of the mainstream discussions, to the best of the present authors' knowledge. It is worth stressing that, under certain conditions in turbulent flames, the emphasized physical mechanism might play a less important role when compared to preferential diffusion or thermal expansion effects, which are not addressed in the present study. This reservation should be borne in mind when applying the reported results to modeling premixed turbulent combustion.

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

For these purposes, a DNS study of propagation of (i) an infinitely thin interface, see fine red line in Figure 1, and (ii) a single-reaction wave of a finite thickness, see thick orange shape, in statistically the same (in both cases) homogeneous, isotropic, forced, constant-density turbulence affected neither by the interface nor by the reaction was performed.

The constant-density turbulence is described by the continuity

$$\nabla \cdot \mathbf{u} = 0 \tag{1}$$

and Navier–Stokes equations

$$
\partial \mathbf{u} / \partial \mathbf{t} + (\mathbf{u} \cdot \nabla) \mathbf{u} = -\rho^{-1} \nabla \mathbf{p} + \mathbf{v} \nabla^2 \mathbf{u} + \mathbf{f},\tag{2}
$$

where t is time, **u** is the flow velocity vector, ρ, ν, and p are the density, kinematic viscosity, and pressure, respectively, and a vector-function **f** is added in order to maintain constant turbulence intensity by using energy forcing at low wavenumbers.

Propagation of the infinitely thin interface is modeled by level set (or G) equation [22]

$$
\partial \mathbf{G} / \partial \mathbf{t} + \mathbf{u} \cdot \nabla \mathbf{G} = \mathbf{S}\_{\mathcal{L}} \lor \nabla \mathbf{G} \mid \tag{3}
$$

where G is a signed distance function to the closest interface associated with G(**x**,t) = 0. Attributes, methodology, and results of the simulations that dealt with Equations (1)–(3) were already discussed by us in details in References [23,24].

Moreover, propagation of a reaction wave of a non-zero thickness is modeled by the following reaction-diffusion equation

$$
\partial \mathbf{c} / \partial \mathbf{t} + \mathbf{u} \cdot \nabla \mathbf{c} = \mathbf{D} \Delta \mathbf{c} + \mathbf{W}\_{\prime} \tag{4}
$$

for a scalar field c, which is equal to zero and unity in fresh reactants and products, respectively. Here,

$$\mathcal{W} = (1 - \mathbf{c}) \exp[- (\mathbf{Z} \mathbf{e} (1 + \tau)^2) / \tau (1 + \tau \mathbf{c})] / [\tau \mathbf{r} (1 + \tau)] \tag{5}$$

is the reaction rate, τ = 6, and Ze = 6. DNS cases are set up (i) by specifying the reaction-wave speed SL and thickness δ<sup>F</sup> = D/SL and (ii) by finding required constant values of D and reaction time scale τ<sup>r</sup> in pre-simulations of a planar 1D laminar-wave modeled by Equations (4) and (5).

Because the attributes, methodology, and results of the simulations that dealt with Equations (1), (2), (4), and (5) were already discussed by us in detail [25,26], we will restrict ourselves to a very brief summary of those DNSs and refer the interested reader to the cited papers.

The computational domain is a rectangular box of size of Λ<sup>x</sup> × Λ × Λ. It is discretized using a uniform staggered Cartesian grid of Nx × N × N cells with Nx = N(Λx/Λ) = 4N. Therefore, spatial resolution Δx = Λx/Nx = Λ/N = Δy = Δz is the same in the axial (x) and transverse (y and z) directions. The boundary conditions are periodic in all three directions, thus, enabling a piece of reaction zone that comes to the left boundary (x = 0) at certain t, y and z to enter the computational domain through the right boundary (x = Λx) at the same t, y and z, respectively.

Turbulence is generated and forced using techniques discussed elsewhere [27–29]. As shown earlier [23–26], the turbulence achieves statistical stationarity, homogeneity, and isotropy over the entire domain, with correlations Rxx(r) = u(x,y,z,t)u(x + r,y,z,t), Ryy(r) = v(x,y,z,t)v(x,y + r,z,t), and Rzz(r) = w(x,y,z,t)w(x,y,z + r,t) being very close to each over and vanishing at r = Λ/2. Here, the mean values · are evaluated by averaging the velocity fields over transverse coordinates and time.

The simulations are performed using three velocity fields A, B, and C, whose characteristics are reported in Table 1. Here, L11 is the longitudinal integral length scale of the turbulence evaluated by the integrating the correlation Rxx(r) over distance r, <sup>η</sup> = (ν3/ <sup>ε</sup>) 1/4 is the Kolmogorov length scale, ε is the dissipation rate ε = 2νSijSij averaged over the computational volume and time, and Sij = (*∂*ui/*∂*xj + *∂*uj/*∂*xi)/2 is the rate-of-strain tensor. The major difference between the three velocity fields consists of the width Λ of the computational domain, which controls the length scale L11 and the initial Reynolds number Re0 = u Λ/(4ν). In other words, L11 and Re0 are increased by increasing Λ, whereas u or ν remain the same in the simulations.



All cases studied by solving either Equations (1–3) or Equations (1 and 2) and (4 and 5) and using velocity field A, B, or C are reported in Table 1, where Da = τt/τ<sup>f</sup> and Ka = τf/τ<sup>η</sup> are the Damköhler and Karlovitz numbers, respectively, τ<sup>f</sup> = δF/SL, τ<sup>t</sup> = L11/u , and τ<sup>η</sup> = (ν/ ε) 1/2 are the reaction wave, eddy-turn-over, and Kolmogorov time scales, respectively. For each velocity field, five cases characterized by different ratios of u /SL are studied by varying SL. When propagation of a reaction wave of a nonzero thickness is addressed using Equations (1 and 2) and (4 and 5), the laminar-wave thickness δ<sup>F</sup> retains the same value in all 15 cases A, B, and C in spite of variations in SL. In order to keep δ<sup>F</sup> constant, the Schmidt number Sc = ν/D is changed. A ratio of L11/δ<sup>F</sup> is increased by increasing the width Λ of the computational domain from field A to field C. In a single TRW (thick reaction wave) case, the ratio of L11/δ<sup>F</sup> is decreased by increasing the thickness δ<sup>F</sup> when compared to Case C5. It is

worth remembering that the values of L11/δF, Da, and Ka, reported in Table 1, characterize the reaction waves described with Equations (1 and 2) and (4 and 5), whereas L11/δ<sup>F</sup> = Da = ∞ and Ka = 0 for an infinitely thin interface, which is characterized solely with u /SL in the present study.

In order to initiate the studied process, either an interface G(x0,y,z,t) = 0 or a planar wave c(x,t) = cL(ξ) is released at x0 = Λx/2 and t = 0 so that the value of the combustion progress variable integrated over the half-space ofx<x0 to be equal to the value of c integrated over the half-space ofx>x0. Here, ξ = x − x0 and cL(ξ) is the pre-computed laminar-wave profile. Subsequently, evolution of the field G(x,t) or c(x,t) is simulated by solving Equation (3) or Equations (4) and (5), respectively. In the former case, c(x,t) = H[G(x,t)], where H is Heaviside function.

Mean bulk turbulent consumption velocity is evaluated as follows [24]

$$
\langle \mathbf{U}\_{\rm T} \rangle = \Lambda^{-2} \mathbf{d} / \det \langle \iiint \mathbf{c}(\mathbf{x}, \mathbf{t}) d\mathbf{x} \rangle \tag{6}
$$

because this method can be applied to both sets of simulations. When processing the DNS data obtained by solving the reaction-diffusion Equation (3), the following alternative method

$$
\langle \mathsf{U}\_{\mathsf{T}} \rangle = \Lambda^{-2} \langle \iiint \, \mathsf{W}(\mathsf{x}, \mathsf{t}) \mathsf{d}\mathsf{x} \rangle \tag{7}
$$

is also applied. Because Equations (6) and (7) yield very close results in all studied cases, solely values of UT, obtained using Equation (6), will be reported in the following. Here, · designates a mean value averaged over time from t\* till t\* + Δt, the starting instant t\* allows the forced turbulence to reach statistical stationarity, and Δt is longer than 50 time scales τ<sup>0</sup> = Λ/(4u ). Reported in the two right columns in Table 1 are normalized fully developed mean bulk turbulent consumption velocities obtained by solving Equations (4 and 5) and (3), respectively. Cases C3 and C4 were not simulated by solving the G-equation (3) in Reference [24].

It is worth stressing that the major difference between the two sets of the simulations consists in solving either the G-equation (3) or the reaction-diffusion Equation (4), whereas numerical methods, turbulence characteristics, etc. are basically similar in both sets. Accordingly, in the present communication, the focus of consideration is placed on phenomena that stem from the finite thickness of the reaction wave.

As far as the physical mechanisms discussed in Section 1 are concerned, the present simulations allow for all mechanisms from Group (C). Moreover, the stretch effect is also addressed. Indeed, variations in the local consumption velocity in and the local quenching of stretched inherently laminar flames can occur in the studied adiabatic equidiffusive flames, because the Zeldovich number Ze in Equation (5) is finite [8]. However, it is worth remembering that the stretch effect can be much more pronounced if DF = DO = κ or if heat losses play a substantial role. The thermal expansion effects are not taken into account in the present simulations, because the physical mechanisms from Group (C) do not allow for them. We may also note that (i) the vast majority of approximations of experimental data on UT, e.g., see review papers [30,31], do not invoke the density ratio σ, thus, implying a weak influence of σ on UT or ST, (ii) recent target-directed experiments [32], as well as earlier measurements [33], did not reveal a substantial influence of σ on UT either, and (iii) recent DNS studies, e.g., Figures 10 and 11 in [34] or Figure 2a in [35], do not indicate such an influence.

#### **3. Results**

Normalized fully developed turbulent consumption velocities UT/SL obtained by solving either the G-equation (2) or the reaction-diffusion Equation (3), with all other things being equal, are plotted in open or filled symbols, respectively, in Figure 2. Comparison of the two sets of results indicates, first, that UT,G is significantly higher than UT,c, with the difference in UT,G and UT,c being increased by u /SL. Therefore, the nonzero thickness of the reaction wave reduces its bulk consumption velocity when compared to the infinitely thin interface. Second, if u /SL is kept constant, then, the ratio of

 UT,c/SL is increased by L11/δF, cf. filled blue squares, black circles, and red triangles. Third, while UT,G/SL depends linearly on u /SL, see open symbols and dashed lines, the solid curves, which fit to the DNS data on UT,c/SL, are clearly bent. Fourth, the bending effect is more pronounced at a lower L11/δF, see filled blue squares.

**Figure 2.** Normalized fully developed turbulent consumption velocity UT/SL (symbols) and a relative increase δA in the wave surface area (dotted-dashed lines) vs normalized root mean square (rms) velocity u /SL. Open symbols show DNS data obtained in the case of an infinitely thin interface, while dashed lines are linear fits to these data. Filled symbols show results simulated in the case of a reaction wave of a nonzero thickness, while solid lines are power-law fits to these data.

It is worth noting that the computed increase in UT,c/SL by L11/δ<sup>F</sup> or u /SL agrees with available experimental data, at least qualitatively. In particular, solid lines in Figure 2 show power-law fits a(u /SL) b, where b = 0.44, 0.51, and 0.58 for L11/δ<sup>F</sup> = 2.0, 3.7, and 6.7, respectively. Similar values of the scaling exponent q in UT∝u<sup>q</sup> were documented in various experimental studies, e.g., q ≈ 0.5 [30,31,36,37], q = 0.56, see data by Kido et al. [38] fitted in Reference [31], q = 0.63 [39], or q = 0.67, see data by Kobayashi et al. [40,41] fitted in Reference [31]. All these experimental databases also indicate an increase in UT by an integral length scale of turbulence.

Curves plotted in dotted-dashed lines in Figure 2 show that a relative increase δAF

$$
\delta \mathbf{A}\_{\mathbf{F}} = \Lambda^{-2} (\mathbf{c}\_2 - \mathbf{c}\_1)^{-1} \langle \iiint \nabla \mathbf{c} \, \mathbb{I}\_{\mathbf{f}} (\mathbf{x}, \mathbf{t}) d\mathbf{x} \rangle \tag{8}
$$

in the area AF of the reaction-zone surface depends on u /SL very similarly to UT,c/SL. Here, |∇c|f is the value of the Flame Surface Density |∇c| conditioned to the reaction zone c1 ≤ c ≤ c2, whose boundaries c1 and c2 are given by W(c1) = W(c2) = max{W(c)}/2. In fact, the bending of the δAF(u /SL)-curves (dotted-dashed lines) is even more pronounced when compared to UT,c/SL-curves (solid lines), thus, implying a slight increase in uc/SL by u /SL. This observation could be attributed to (i) an increase in the reaction rate integrated over the mixing zones (c ≤ c1), as such zones are expanded by small-scale turbulent eddies [19], or (ii) an increase in the local reaction rate in the vicinity of cusps [42]. Because the increase in uc/SL by u /SL does not contribute to the bending effect, but weakly resists it, further discussion of this trend is beyond the scope of the present communication.

All in all, DNS data reported in Figure 2 indicate that, under conditions of the present study, the bending of UT,c/SL(u /SL)-curves is mainly controlled by the bending of δAF(u /SL)-curves, in line with the recent finding by Nivarti and Cant [10]. The present DNS data imply that the bending effect stems from the finite thickness of the reaction wave. Indeed, because the sole difference between the previous [23,24] and present simulations consists in substituting the G-equation (3) with Equations (4) and (5), the obtained difference between the linear and bent UT,c/SL(u /SL)-curves should not be attributed to the exploited method of turbulence generation, numerical resolution or scheme, insufficiently large width of the computational domain, etc. The difference between the linear and bent curves can solely be controlled by δF. Let us discuss physical mechanisms that could control this difference.

First, the DNSs yield UT,c/SL ≈ δAF, thus, implying that, under conditions of the present study, the influence of turbulence on UT,c is not controlled by an increase in local burning rate due to enhancement of the local heat and mass transfer by turbulent eddies. It is worth noting, however, that DNS data computed by us at significantly higher Ka 1 [26,43] support scaling of UT,c/SL∝ (u L11/D)1/2, as predicted by Damköhler [19] by reducing the influence of turbulence on burning rate to the enhancement of the local heat and mass transfer by the turbulence. In other words, the mechanism by Damköhler [19] can play an important role, but under conditions that differ significantly from the conditions of the present study (e.g., the present DNS data show a substantially weaker dependence of UT,c on L11 when compared to the Damköhler's scaling).

Second, eventual reduction in δAF due to interface-interface or wave-wave collisions [12,13] is taken into account when numerically solving Equation (3) or (4), respectively. Nevertheless, the bending effect is not pronounced in the former case, thus, implying that the collisions do not control the effect under conditions of the present study. If the collisions were of importance, they would yield the bending in the simulations with Equation (3) also.

Third, a physical mechanism that controls the bending effect under conditions of the present study is revealed in Figure 3. In particular, DNS data plotted in lines in Figure 3 indicate that the cumulative probability of finding highly curved reaction zones characterized by a sufficiently large product |ηhm|>b of the Kolmogorov length scale and the local iso-surface curvature conditioned to c1 ≤ c ≤ c2 (i) is weakly affected by L11 if the laminar-wave thickness δ<sup>F</sup> is kept constant, cf. dotted-dashed lines which show results obtained in cases A5, B5, and C5, but (ii) is significantly increased when δ<sup>F</sup> is decreased, cf. (ii.a) violet dashed line (case TRW) and the three dotted-dashed lines or (ii.b) magenta solid and black double-dashed-dotted lines, which show results computed by solving Equations (3)–(5), respectively, in case C5. Here, b is a positive number of unity order, hm = ∇·**n**<sup>G</sup> and **n**<sup>G</sup> = −∇G/|∇G| or hm = ∇·**n**<sup>c</sup> and **n**<sup>c</sup> = −∇c/|∇c| in the simulations that deal with Equation (3) or Equations (4) and (5), respectively. Therefore, when the thickness δ<sup>F</sup> is decreased, the probability of appearance of small-scale highly curved wrinkles on the reaction-zone surface is increased, thus, increasing the surface area AF and, hence, the turbulent consumption velocity UT,c. The highest values of AF and UT,G are associated with the infinitely thin interface and UT,G∝u in this case, see open symbols and dashed lines in Figure 2. An increase in δ<sup>F</sup> results in decreasing UT,c when compared to UT,G∝u , i.e., the bending effect.

DNS data plotted in symbols in Figure 3 support such a scenario by indicating that the cumulative probability of finding highly curved reaction zones characterized by |δFhm| > b at c1 ≤ c ≤ c2 (i) is decreased with increasing L11/δF, cf. filled symbols, but (ii) is weakly affected by variations in the velocity field (A or C), Ka, and Kolmogorov scales, provided that u /SL, L11/δF, and, hence, the Damköhler number are kept constant, cf. red triangles (case A5) and violet crosses (case TRW). Therefore, comparison of lines and symbols in Figure 3 implies that the curvature magnitude is primarily controlled by the thickness δF. Because variations in δ<sup>F</sup> in cases C5 and TRW result from variations in the diffusivity D, the significant dependence of the curvature magnitude on δ<sup>F</sup> implies that it is the molecular transport that impedes wrinkling reaction-zone surface by small-scale eddies in intense turbulence.

In other words, the DNS data plotted in Figure 3 imply that the magnitude of the local curvature of the reaction zone is bounded by the reaction-wave thickness δF, whereas smaller-scale turbulent eddies are inefficient in wrinkling the reaction-zone surface. Accordingly, the small-scale range of the turbulence spectrum weakly contributes to an increase in the area AF of the surface. Such a hypothesis is further supported by approximately equal values of UT,c obtained in cases A5 and TRW, associated with almost the same Da (which characterizes large-scale eddies), but substantially different Ka (which characterizes the smallest-scale turbulent eddies).

**Figure 3.** Cumulative probability that the normalized absolute value of local curvature of reaction zone is larger than a positive threshold number b. Lines and symbols show results obtained by normalizing the curvature using the Kolmogorov length scale η and the laminar-wave thickness δF, respectively. All results were obtained in cases characterized by the same u /SL = 10, but different L11/δ<sup>F</sup> = 2.1 (cases A5 and TRW), L11/δ<sup>F</sup> = 3.7 (case B5), or L11/δ<sup>F</sup> = 6.7 (case C5). In case TRW, the laminar reaction-wave thickness δ<sup>F</sup> is larger by a factor of four when compared to a reference value used in all other cases associated with Equations (4) and (5). Magenta solid line shows DNS data obtained by solving the level set Equation (3).

To illustrate a mechanism that bounds |hm| by δF, let us follow Zel'dovich et al. [44] and rewrite Equation (4) as follows

$$\left(\partial\mathbf{\hat{c}}/\partial\mathbf{t} + \mathbf{u}(\partial\mathbf{c}/\partial\mathbf{r}) - 2(\mathbf{D}/\mathbf{r})(\partial\mathbf{c}/\partial\mathbf{r}) = \mathbf{D}(\partial^2\mathbf{c}/\partial\mathbf{r}^2) + \mathbf{W} \tag{9}$$

in the spherical coordinate framework. For an expanding reaction wave, the last term on the left-hand side reduces the wave speed when compared to the counterpart planar wave. If the local radius rw of the wave curvature is on the order of the Kolmogorov length scale η, the magnitude D/rw of the negative speed resulting from the considered term is on the order of the Kolmogorov velocity (if the Schmidt number Sc = O(1)) and can be much larger than SL if Ka 1. This curvature-induced speed is negative (positive) for a local wrinkle with the curvature center in products (reactants) and, therefore, tends to damp the local wrinkle. This physical mechanism is controlled by the molecular diffusion and acts even if the local consumption velocity does not depend on the curvature radius rw. Moreover, while a turbulent eddy whose length scale is significantly smaller than δ<sup>F</sup> perturbs the local wave surface during a short lifetime of the eddy, the considered mechanism can smooth the perturbation even after disappearance of the eddy until D/rw = O(SL) and rw = O(δF). Accordingly, if η < δF, the small-scale range of the entire turbulence spectrum appears to be inefficient in wrinkling the wave surface due to such a smoothing mechanism and this inefficient range expands to smaller length scales when η/δ<sup>F</sup> is decreased due to an increase in u /SL. This physical mechanism acts to reduce increasing the wave surface area and turbulent consumption velocity with increasing u /SL, thus, causing the bending effect.

Indeed, filled black triangles in Figure 4 show that, in case TRW associated with the most pronounced bending effect, the probability of finding high locally negative displacement speed Sd = (D∇2c + W)/|∇c| of the reaction zone is strongly increased by hm if |δFhm| < 1. The probability of Sd/SL < −1 is larger than 60% if δFhm > 2. Accordingly, if the local curvature of the reaction zone

is positive and sufficiently high (δFhm > 1), the zone statistically tends (i) to move to products, i.e., to the curvature center, and, therefore, (ii) to smooth out the local wrinkle of the zone surface. If the local curvature of the reaction zone is negative and sufficiently high (δFhm < −1), the zone moves to reactants (the probability of finding positive Sd > SL is almost equal to unity, see red stars), i.e., to the curvature center, and, therefore, smooths out the local wrinkle on the zone surface again.

**Figure 4.** Probabilities of Sd/SL < −1 (violet dashed line and black triangles) and Sd/SL > 1 (blue dotted-dashed line and red circles), conditioned to (i) the local curvature (symbols) hm of the reaction zone, normalized using <sup>δ</sup>F, or (ii) the local total strain (lines) S<sup>2</sup> normalized using <sup>ε</sup>/(2ν). Case TRW.

On the contrary, violet dashed or blue dotted-dashed line shows that the probability of finding Sd < −SL or Sd > SL, respectively, in the reaction zone (c1 ≤ c ≤c2) depends weakly on the total strain S2=Sij Sij. This result implies that turbulent strain rates affect propagation of the reaction zone with respect to the local flow weakly from the statistical viewpoint.

#### **4. Discussion**

The above analysis of the present DNS data indicates that the bending of UT,c/SL(u /SL)-curves, computed in the case of δ<sup>F</sup> > 0 (see filled symbols in Figure 2) is controlled by the bending of the mean area of the reaction-zone surface as a function of u /SL (see dotted-dashed lines in Figure 2). The latter bending is controlled by the following physical mechanism. When a reaction front has a negligible thickness, turbulent eddies of various scales can wrinkle the front surface, increase its area, and, hence, increase turbulent consumption velocity. However, if the local thickness δ<sup>F</sup> of the reaction wave is comparable with or larger than the Kolmogorov length scale η, the local molecular transport efficiently smooths out small-scale wrinkles of the reaction-zone surface. Therefore, the local molecular transport impedes increasing the surface area due to the highest local stretch rates created by the smallest turbulent eddies. Consequently, the turbulent consumption velocity UT is reduced when compared to the case of the infinitely thin interface.

It is worth noting that the above analysis is consistent with DNS data by Wenzel and Peters [45], who simulated self-propagation of a passive interface in constant-density turbulence by solving Equation (3), where SL was substituted with SL−νSc−<sup>1</sup> ∇·**n**G. Indeed, results obtained at Sc = <sup>∞</sup>, i.e., using the unperturbed laminar-wave speed SL indicate a linear dependence of δAF on u [45], with UT = SLδAF being very close to UT,G obtained by us by solving Equation (3), cf. stars with other symbols in Figure 4 in Reference [24]. However, at Sc = 0.25 or 0.50, the bending of AF(u )-curves reported

in [45] is well pronounced due to the physical mechanism highlighted above, as the second term in SL−νSc−<sup>1</sup> ∇·**n**<sup>G</sup> models the smoothing effect of the curvature-induced displacement speed.

The smoothing effect emphasized in the present work differs fundamentally from the well-recognized stretch effect [3–5]. In particular, the smoothing effect (i) is controlled by the molecular diffusivity D of the deficient reactant and (ii) manifests itself in the reduction in AF independently of a ratio of uc/SL. On the contrary, the stretch effect (i) is mainly controlled by the differences in DF, DO, and κ in the adiabatic case and (ii) manifests itself in significant variations in the mean local consumption velocity uc independently of the area AF. In the present study, the stretch effect is not pronounced, because the differences in the molecular transport coefficients are not allowed for.

At first glance, the smoothing mechanism emphasized in the present communication appears to be similar to the smoothing effect of flame propagation, which is considered to play an important role by fractal models of premixed turbulent combustion. For instance, the latter mechanism controls the Gibson length scale [46]. However, there are fundamental differences between this well-known propagation-smoothing mechanism and the diffusion-smoothing mechanism emphasized in the present communication. For instance, first, the former and latter mechanisms are associated with the right-hand side and the third term on the left-hand side of Equation (9), respectively. Accordingly, a ratio of magnitudes of effects caused by the two mechanisms may be estimated as follows rwSL/D. If the radius rw of the curvature created by the smallest-scale turbulent eddies is significantly less than δ<sup>F</sup> = D/SL, the diffusion-smoothing mechanism should dominate. Second, the propagation-smoothing mechanism can act not only in the case of a reaction wave of a finite thickness, but also in the case of self-propagation of an infinitely thin interface. The fact that the simulations with the G-equation do not yield the bending effect, see open symbols and dashed lines in Figure 2, implies that the propagation-smoothing mechanism is of minor importance under conditions of the present study. It is also worth noting that the present results are qualitatively consistent with experimental data [47] that indicate that the inner cut-off length scale of wrinkles on a flame surface scales as the laminar flame thickness, rather than the Gibson length scale.

In a variable density case, the efficiency of small-scale eddies in wrinkling the reaction zone can be reduced not only due to the smoothing mechanism emphasized in the present communication, but also due to the disappearance of the small-scale eddies due to dilatation and an increase in the mixture viscosity by the temperature. In order to understand what mechanism (smoothing or thermal expansion) dominates and under which conditions, the present DNS study should be extended to flames characterized by substantial density variations.

#### **5. Conclusions**

A DNS study of propagation of either an infinitely thin interface or a reaction wave of a nonzero thickness in forced, constant-density, statistically stationary, homogeneous, isotropic turbulence was performed by solving Navier–Stokes equations and either a level set or a reaction-diffusion equation, respectively. In the latter case, the computed mean wave speed <UT,c> (i) is reduced when a ratio L11/δ<sup>F</sup> of the longitudinal integral length scale L11 of the turbulence to the laminar wave thickness δ<sup>F</sup> is decreased and (ii) is significantly lower than <UT,G> simulated in the former case, with all other things being equal. Moreover, the following results obtained in the present work stem from the finite thickness δ<sup>F</sup> of the reaction wave.

First, the computed <UT,c>/SL(u /SL)-curves show bending. The bending effect is less pronounced at higher L11/δ<sup>F</sup> and vanishes in the case of the infinitely thin interface.

Second, under conditions of the present study, the bending effect is controlled by a decrease in the rate of an increase δAF in the reaction-zone-surface area with increasing u /SL. In its turn, the bending of the δAF(u /SL)-curves stems from inefficiency of small-scale turbulent eddies in wrinkling the reaction-zone surface, because small-scale wrinkles are smoothed out by molecular transport within the local reaction wave. Such a smoothing effect is not pronounced in the case of self-propagation of the infinitely thin interface at a constant speed SL in statistically the same turbulence.

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

**Funding:** The first author (R.Y.) gratefully acknowledges the financial support by the Swedish Research Council (VR) and the National Centre for Combustion Science and Technology (CeCOST). The second author (A.N.L.) gratefully acknowledges the financial support by Chalmers Area of Advance Transport and the Combustion Engine Research Center (CERC).

**Acknowledgments:** The simulations were performed using the computer facilities provided by the Swedish National Infrastructure for Computing (SNIC) at Beskow-PDC Center. Some results of this work were orally presented at the 10th Mediterranean Combustion Symposium and were reported in the Proceedings of the Symposium, see http://ircserver2.irc.cnr.it/.

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

## **References**


#### *Fluids* **2019**, *4*, 31


© 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/).
