DNS study of the bending effect due to smoothing mechanism

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.


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 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][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, smallscale 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 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.

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 and Navier-Stokes equations 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|, 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 for a scalar field c, which is equal to zero and unity in fresh reactants and products, respectively. Here, 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. Here, L 11 is the longitudinal integral length scale of the turbulence evaluated by the integrating the correlation R xx (r) over distance r, η = (ν 3 / ε ) 1/4 is the Kolmogorov length scale, ε is the dissipation rate ε = 2νS ij S 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 11 and the initial Reynolds number Re 0 = u Λ/(4ν). In other words, L 11 and Re 0 are increased by increasing Λ, whereas u or ν remain the same in the simulations.  (1)(2)(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 11 /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 11 /δ 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 11 /δ F is decreased by increasing the thickness δ F when compared to Case C5. It is worth remembering that the values of L 11 /δ F , Da, and Ka, reported in Table 1, characterize the reaction waves described with Equations (1 and 2) and (4 and 5), whereas L 11 /δ 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 0 ,y,z,t) = 0 or a planar wave c(x,t) = c L (ξ) is released at x 0 = Λ x /2 and t = 0 so that the value of the combustion progress variable integrated over the half-space of x < x 0 to be equal to the value of c integrated over the half-space of x > x 0 . Here, ξ = x − x 0 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 = Λ −2 d/dt Fluids 2019, 4, x 5 of 13 Fluids 2019, 4, x; doi: www.mdpi.com/journal/fluids waves described with Equations (1 and 2) and (4 and 5), whereas L11⁄δF = 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 of x < x0 to be equal to the value of c integrated over the half-space of x > 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] 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 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 τ 0 = Λ⁄(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.

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.
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 Fluids 2019, 4, x 5 of 13 waves described with Equations (1 and 2) and (4 and 5), whereas L11⁄δF = 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 of x < x0 to be equal to the value of c integrated over the half-space of x > 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] 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 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 τ 0 = Λ⁄(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.

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.
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 τ 0 = Λ/(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.

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 11 /δ 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 11 /δ F , see filled blue squares. It is worth noting that the computed increase in 〈UT,c〉/SL by L11⁄δF 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⁄δF = 2.0, 3.7, and 6.7, respectively. Similar values of the scaling exponent q in UT∝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 UT by an integral length scale of turbulence.
Curves plotted in dotted-dashed lines in Figure 2 show that a relative increase δAF 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 It is worth noting that the computed increase in U T,c /S L by L 11 /δ 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 11 /δ 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 Fluids 2019, 4, x 5 of 13 waves described with Equations (1 and 2) and (4 and 5), whereas L11⁄δF = 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 of x < x0 to be equal to the value of c integrated over the half-space of x > 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] 〈UT〉 = Λ -2 d/dt〈∭ c(x,t)dx 〉, 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 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 τ 0 = Λ⁄(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 |∇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 1 ≤ c ≤ c 2 , whose boundaries c 1 and c 2 are given by W(c 1 ) = W(c 2 ) = 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 1 ), 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 11 /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 11 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 1 ≤ c ≤ c 2 (i) is weakly affected by L 11 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.  (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 |δ F h m | > b at c 1 ≤ c ≤ c 2 (i) is decreased with increasing L 11 /δ 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 11 /δ 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  To illustrate a mechanism that bounds |hm| 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(∂ 2 c/∂r 2 ) + 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 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 δ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⁄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 η⁄δF 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∇ 2 c + 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 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 /S L = 10, but different L 11 /δ F = 2.1 (cases A5 and TRW), L 11 /δ F = 3.7 (case B5), or L 11 /δ F = 6.7 (case C5). In case TRW, the laminar reaction-wave thickness δ F 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 |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(∂ 2 c/∂r 2 ) + 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∇ 2 c + W)/|∇c| of the reaction zone is strongly increased by h m if |δ F h m | < 1. The probability of S d /S L < −1 is larger than 60% if δ F h m > 2. Accordingly, if the local curvature of the reaction zone is positive and sufficiently high (δ F h 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 (δ F h 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. Fluids 2019, 4, x; doi: www.mdpi.com/journal/fluids 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. 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 S 2 =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.

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 δ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'⁄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 δ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 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 -1 ∇•nG. Indeed, results obtained at Sc = ∞, 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 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 1 ≤ c ≤c 2 ) depends weakly on the total strain S 2 =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.

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 −1 ∇·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 −1 ∇·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 w S 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.

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 11 /δ F of the longitudinal integral length scale L 11 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 11 /δ 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.