DNS Study of the Bending Effect Due to Smoothing Mechanism

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 S_L depends linearly on the normalized rms velocity u′/S_L. In the case of the reaction wave of a nonzero thickness, dependencies of the normalized bulk consumption velocity on u′/S_L 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 δA_F in the reaction-zone-surface area with increasing u′/S_L. In its turn, the bending of the δA_F(u′/S_L)-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.

1. Introduction

In turbulent reacting flows, the so-called bending effect consists in decreasing the rate (dU_T)/du′ of an increase in turbulent consumption velocity U_T by the rms turbulent velocity u′ with increasing u′, i.e., the second derivative of the function U_T(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 u_c, caused by turbulent stretching of the wave [3,4,5]. Here, u_c 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 〈u_c〉 of the local consumption velocity with increasing u′, followed by eventual local reaction extinction at sufficiently high u′. A decrease in 〈u_c〉 and the local reaction extinction can also affect the area A_F of the reaction-wave surface, but this manifestation of the discussed mechanism is indirect, i.e., it is a consequence of the dependence of u_c or 〈u_c〉 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 u_c on the local stretch rate $\dot{\mathrm{s}}$ is controlled by differences in molecular diffusivities of a fuel, D_F, oxygen, D_O, and the heat diffusivity κ of the mixture. In the equidiffusive (D_F = D_O = κ) adiabatic case, u_c does not depend on $\dot{\mathrm{s}}$ if E tends to infinity, but u_c can depend on $\dot{\mathrm{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].

Even if the mean 〈u_c〉 is close to the speed S_L 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 〈A_F〉, i.e., the second derivative of the function 〈A_F〉(u′) can be negative in spite of 〈u_c〉 ≈ S_L. 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,4,5,11,12,13,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 U_T∝S_L(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 〈u_c〉 when compared to S_L.

Finally, the bending effect can be pronounced differently in different reaction waves. For instance, levelling-off of U_T(u′)-curves, followed by a decrease in U_T with increasing u′, is well documented in expanding statistically spherical flames [1,2], but the positive dU_T/du′ was obtained from statistically stationary flames at much higher values of u′/S_L [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,4,5], (B) mechanisms that highlight thermal expansion effects in flames, and (C) mechanisms that address equidiffusive flames, but disregard the thermal expansion effects [11,12,13,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 u_c due to the stretch effect [3,4,5] or the transport enhancement [19] and (b) mechanisms that highlight the bending of 〈A_F〉(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 U_T(u′)-curves can be controlled by the bending of 〈A_F〉(u′)-curves, (ii) comparing physical mechanisms from Group (C), and (iii) emphasizing a physical mechanism that controls the bending of 〈A_F〉(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

∇·u = 0

(1)

and Navier–Stokes equations

∂u/∂t + (u·∇)u = −ρ⁻¹∇p + ν∇²u + f,

(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]

∂G/∂t + u·∇G = S_L|∇G|,

(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

∂c/∂t + u·∇c = DΔc + W,

(4)

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

W = (1 − c)exp[−(Ze(1 + τ)²)/τ(1 + τc)]/[τ_r(1 + τ)]

(5)

is the reaction rate, τ = 6, and Ze = 6. DNS cases are set up (i) by specifying the reaction-wave speed S_L and thickness δ_F = D/S_L and (ii) by finding required constant values of D and reaction time scale τ_r 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 Λ_x × Λ × Λ. It is discretized using a uniform staggered Cartesian grid of N_x × N × N cells with N_x = N(Λ_x/Λ) = 4N. Therefore, spatial resolution Δx = Λ_x/N_x = Λ/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,28,29]. As shown earlier [23,24,25,26], the turbulence achieves statistical stationarity, homogeneity, and isotropy over the entire domain, with correlations R_xx(r) = 〈u(x,y,z,t)u(x + r,y,z,t)〉, R_yy(r) = 〈v(x,y,z,t)v(x,y + r,z,t)〉, and R_zz(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. Studied cases.

Case	Re0	Nx × N × N	η/Δx	L11/η	u′/SL	L11/δF	Da	Ka	〈UT,c〉/u′	〈UT,G〉/u′
A1	50	256 × 642	0.68	12	0.5	2.0	4.1	1.34	2.18	3.66
A2	-	-	-	-	1.0	-	2.0	2.69	1.31	2.68
A3	-	-	-	-	2.0	-	1.0	5.38	0.89	2.19
A4	-	-	-	-	5.0	-	0.4	13.4	0.57	1.83
A5	-	-	-	-	10.0	-	0.2	26.9	0.40	1.67
B1	100	512 × 1282	0.87	17	0.5	3.7	7.5	0.84	2.45	3.66
B2	-	-	-	-	1.0	-	3.8	1.67	1.69	2.71
B3	-	-	-	-	2.0	-	1.9	3.34	1.27	2.23
B4	-	-	-	-	5.0	-	0.8	8.36	0.80	1.87
B5	-	-	-	-	10.0	-	0.4	16.7	0.55	1.77
C1	200	1024 × 2562	1.07	25	0.5	6.7	14.	0.55	2.77	3.70
C2	-	-	-	-	1.0	-	6.7	1.10	2.06	2.80
C3	-	-	-	-	2.0	-	3.4	2.21	1.61	-
C4	-	-	-	-	5.0	-	1.4	5.52	1.11	-
C5	-	-	-	-	10.0	-	0.7	11.0	0.76	1.87
TRW	-	-	-	-	-	2.1	0.2	36.2	0.42	-

. Here, L₁₁ is the longitudinal integral length scale of the turbulence evaluated by the integrating the correlation R_xx(r) over distance r, η = (ν³/〈〈ε〉〉)^1/4 is the Kolmogorov length scale, 〈〈ε〉〉 is the dissipation rate ε = 2νS_ijS_ij averaged over the computational volume and time, and S_ij = (∂u_i/∂x_j + ∂u_j/∂x_i)/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 L₁₁ and the initial Reynolds number Re₀ = u′Λ/(4ν). In other words, L₁₁ and Re₀ 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/τ_f and Ka = τ_f/τ_η are the Damköhler and Karlovitz numbers, respectively, τ_f = δ_F/S_L, τ_t = L₁₁/u′, and τ_η = (ν/〈〈ε〉〉)^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′/S_L are studied by varying S_L. 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 δ_F retains the same value in all 15 cases A, B, and C in spite of variations in S_L. In order to keep δ_F constant, the Schmidt number Sc = ν/D is changed. A ratio of L₁₁/δ_F 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 L₁₁/δ_F is decreased by increasing the thickness δ_F when compared to Case C5. It is worth remembering that the values of L₁₁/δ_F, Da, and Ka, reported in Table 1, characterize the reaction waves described with Equations (1 and 2) and (4 and 5), whereas L₁₁/δ_F = Da = ∞ and Ka = 0 for an infinitely thin interface, which is characterized solely with u′/S_L in the present study.

In order to initiate the studied process, either an interface G(x₀,y,z,t) = 0 or a planar wave c(x,t) = c_L(ξ) is released at x₀ = Λ_x/2 and t = 0 so that the value of the combustion progress variable integrated over the half-space of x < x₀ to be equal to the value of c integrated over the half-space of x > x₀. Here, ξ = x − x₀ and c_L(ξ) 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]

〈U_T〉 = Λ⁻²d/dt〈∭ c(x,t)dx〉,

(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

〈U_T〉 = Λ⁻²〈∭ W(x,t)dx〉

(7)

is also applied. Because Equations (6) and (7) yield very close results in all studied cases, solely values of 〈U_T〉, 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 τ⁰ = Λ/(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 D_F ≠ D_O ≠ κ 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 U_T, e.g., see review papers [30,31], do not invoke the density ratio σ, thus, implying a weak influence of σ on U_T or S_T, (ii) recent target-directed experiments [32], as well as earlier measurements [33], did not reveal a substantial influence of σ on U_T 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 〈U_T〉/S_L 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 〈U_T,G〉 is significantly higher than 〈U_T,c〉, with the difference in 〈U_T,G〉 and 〈U_T,c〉 being increased by u′/S_L. Therefore, the nonzero thickness of the reaction wave reduces its bulk consumption velocity when compared to the infinitely thin interface. Second, if u′/S_L is kept constant, then, the ratio of 〈U_T,c〉/S_L is increased by L₁₁/δ_F, cf. filled blue squares, black circles, and red triangles. Third, while 〈U_T,G〉/S_L depends linearly on u′/S_L, see open symbols and dashed lines, the solid curves, which fit to the DNS data on 〈U_T,c〉/S_L, are clearly bent. Fourth, the bending effect is more pronounced at a lower L₁₁/δ_F, see filled blue squares.

It is worth noting that the computed increase in 〈U_T,c〉/S_L by L₁₁/δ_F or u′/S_L agrees with available experimental data, at least qualitatively. In particular, solid lines in Figure 2 show power-law fits a(u′/S_L)^b, where b = 0.44, 0.51, and 0.58 for L₁₁/δ_F = 2.0, 3.7, and 6.7, respectively. Similar values of the scaling exponent q in U_T∝u′^q 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 U_T by an integral length scale of turbulence.

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

δA_F = Λ⁻²(c₂ − c₁)⁻¹〈∭ |∇c|_f(x,t)dx〉

(8)

in the area A_F of the reaction-zone surface depends on u′/S_L very similarly to 〈U_T,c〉/S_L. Here, |∇c|_f is the value of the Flame Surface Density |∇c| conditioned to the reaction zone c₁ ≤ c ≤ c₂, whose boundaries c₁ and c₂ are given by W(c₁) = W(c₂) = max{W(c)}/2. In fact, the bending of the δA_F(u′/S_L)-curves (dotted-dashed lines) is even more pronounced when compared to 〈U_T,c〉/S_L-curves (solid lines), thus, implying a slight increase in 〈u_c〉/S_L by u′/S_L. This observation could be attributed to (i) an increase in the reaction rate integrated over the mixing zones (c ≤ c₁), 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 〈u_c〉/S_L by u′/S_L 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 〈U_T,c〉/S_L(u′/S_L)-curves is mainly controlled by the bending of δA_F(u′/S_L)-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 〈U_T,c〉/S_L(u′/S_L)-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 〈U_T,c〉/S_L ≈ δA_F, thus, implying that, under conditions of the present study, the influence of turbulence on 〈U_T,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 〈U_T,c〉/S_L∝ (u′L₁₁/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 〈U_T,c〉 on L₁₁ when compared to the Damköhler’s scaling).

Second, eventual reduction in δA_F 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 |ηh_m|>b of the Kolmogorov length scale and the local iso-surface curvature conditioned to c₁ ≤ c ≤ c₂ (i) is weakly affected by L₁₁ if the laminar-wave thickness δ_F is kept constant, cf. dotted-dashed lines which show results obtained in cases A5, B5, and C5, but (ii) is significantly increased when δ_F 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, h_m = ∇·n_G and n_G = −∇G/|∇G| or h_m = ∇·n_c and n_c = −∇c/|∇c| in the simulations that deal with Equation (3) or Equations (4) and (5), respectively. Therefore, when the thickness δ_F is decreased, the probability of appearance of small-scale highly curved wrinkles on the reaction-zone surface is increased, thus, increasing the surface area A_F and, hence, the turbulent consumption velocity 〈U_T,c〉. The highest values of A_F and 〈U_T,G〉 are associated with the infinitely thin interface and 〈U_T,G〉∝u′ in this case, see open symbols and dashed lines in Figure 2. An increase in δ_F results in decreasing 〈U_T,c〉 when compared to 〈U_T,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 |δ_Fh_m| > b at c₁ ≤ c ≤ c₂ (i) is decreased with increasing L₁₁/δ_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′/S_L, L₁₁/δ_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 δ_F in cases C5 and TRW result from variations in the diffusivity D, the significant dependence of the curvature magnitude on δ_F 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 A_F of the surface. Such a hypothesis is further supported by approximately equal values of 〈U_T,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).

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

∂c/∂t + u(∂c/∂r) − 2(D/r)(∂c/∂r) = D(∂²c/∂r²) + W

(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 r_w of the wave curvature is on the order of the Kolmogorov length scale η, the magnitude D/r_w 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 S_L 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 r_w. Moreover, while a turbulent eddy whose length scale is significantly smaller than δ_F 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/r_w = O(S_L) and r_w = 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 η/δ_F is decreased due to an increase in u′/S_L. This physical mechanism acts to reduce increasing the wave surface area and turbulent consumption velocity with increasing u′/S_L, 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 S_d = (D∇²c + W)/|∇c| of the reaction zone is strongly increased by h_m if |δ_Fh_m| < 1. The probability of S_d/S_L < −1 is larger than 60% if δ_Fh_m > 2. Accordingly, if the local curvature of the reaction zone is positive and sufficiently high (δ_Fh_m > 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 (δ_Fh_m < −1), the zone moves to reactants (the probability of finding positive S_d > S_L 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.

On the contrary, violet dashed or blue dotted-dashed line shows that the probability of finding S_d < −S_L or S_d > S_L, respectively, in the reaction zone (c₁ ≤ c ≤c₂) depends weakly on the total strain S²=S_ij S_ij. 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 〈U_T,c〉/S_L(u′/S_L)-curves, computed in the case of δ_F > 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′/S_L (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 δ_F 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 U_T 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 S_L was substituted with S_L−νSc⁻¹ ∇·n_G. Indeed, results obtained at Sc = ∞, i.e., using the unperturbed laminar-wave speed S_L indicate a linear dependence of δA_F on u′ [45], with U_T = S_LδA_F being very close to 〈U_T,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 A_F(u′)-curves reported in [45] is well pronounced due to the physical mechanism highlighted above, as the second term in S_L−νSc⁻¹ ∇·n_G 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,4,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 A_F independently of a ratio of 〈u_c〉/S_L. On the contrary, the stretch effect (i) is mainly controlled by the differences in D_F, D_O, and κ in the adiabatic case and (ii) manifests itself in significant variations in the mean local consumption velocity 〈u_c〉 independently of the area A_F. 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 r_wS_L/D. If the radius r_w of the curvature created by the smallest-scale turbulent eddies is significantly less than δ_F = D/S_L, 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 <U_T,c> (i) is reduced when a ratio L₁₁/δ_F of the longitudinal integral length scale L₁₁ of the turbulence to the laminar wave thickness δ_F is decreased and (ii) is significantly lower than <U_T,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 δ_F of the reaction wave.

First, the computed <U_T,c>/S_L(u′/S_L)-curves show bending. The bending effect is less pronounced at higher L₁₁/δ_F 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 δA_F in the reaction-zone-surface area with increasing u′/S_L. In its turn, the bending of the δA_F(u′/S_L)-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 S_L in statistically the same turbulence.