Assessment of an Evolution Equation for the Displacement Speed of a Constant-Density Reactive Scalar Field

The displacement speed that characterises the self-propagation of isosurfaces of a reaction progress variable is of key importance for turbulent premixed reacting flow. The evolution equation for the displacement speed was derived in a recent work of Yu and Lipatnikov (Phys Rev E 100:013107, 2019a) for the case where the flame is described by a transport equation for single reaction progress variable assuming simple transport and one-step chemistry. This equation represents interaction of a number of complex coupled mechanisms related to straining by the velocity field, surface curvature and the scalar gradient. The aim of the current work is to provide detailed physical explanations of the displacement speed equation and its various terms, and to provide a new perspective to understand the mechanisms responsible for observed variations in the displacement speed. The equation is then used to analyze the propagation of a statistically planar reaction wave in homogeneous isotropic constant-density turbulence using direct numerical simulations. Additional emphasis is put on retracting surface segments that have a negative displacement speed, a phenomenon that commonly occurs at high Karlovitz numbers.


3
where t, u , , and denote time, flow velocity vector, diffusion term, and the reaction rate term, respectively. One quantity of particular interest is the displacement speed S d of the scalar field c relative to the local fluid. The displacement speed is controlled by the combined effect of diffusion and reaction and is defined as (Veynante and Vervisch 2002) In a turbulent reacting flow the displacement speed is defined for any point (x, t) for which 0 < c(x, t) < 1 . The displacement speed contains information that characterizes the selfpropagation of the local reaction wave. Statistical information on the displacement speed can therefore benefit premixed combustion modeling, based on for example the level-set (G-equation) formulation (Williams 1985) or the flame surface density (FSD) approach for reaction rate closure (Candel and Poinsot 1990). Furthermore, the variation of displacement speed across different zones inside the reaction wave can provide information regarding the changes in the internal wave structure induced by turbulence. Such information can be particularly useful when studying heavily disturbed flames characterized by a large Karlovitz number (Ka).
Statistical behavior of the displacement speed has been the topic of a number of studies based on direct numerical simulations (DNS) of various turbulent reacting flows, see Gran et al. (1996), Echekki and Chen (1996), Chen and Im (1998), Peters et al. (1998), Echekki and Chen (1999), Chakraborty and Cant (2005), Dopazo et al. (2007), Wang et al. (2017), and Luca et al. (2019) and references therein. In the literature the displacement speed is often examined by correlating it with other quantities characterizing the flame surface or turbulent reacting flow, e.g. with the local flame curvature Chen and Im 1998;Echekki and Chen 1999;Wang et al. 2017;Luca et al. 2019;Chakraborty 2007). It is also common to decompose the displacement speed into three component contributions due to reaction, diffusion in the normal direction, and tangential diffusion induced by curvature (Chen and Im 1998;Peters et al. 1998;Echekki and Chen 1999). Particular focus has been put on the probability of finding locally negative displacement speed (Gran et al. 1996;Echekki and Chen 1996;Peters et al. 1998;Echekki and Chen 1999;Chakraborty and Cant 2005;Wang et al. 2017;Luca et al. 2019;Chakraborty 2007;Sankaran et al. 2015;Cecere et al. 2016;Trisjono et al. 2016;Cecere et al. 2018) as well as the sign of the averaged displacement speed. While previous studies that compare quantities extracted from DNS data contributed significantly to the understanding of turbulence/flame interactions, the mechanism governing the dynamic evolution of displacement speed is still not well understood. This may be attributed to the lack of a quantitative framework describing the evolution of displacement speed that advanced analyses may rely on.
To better understand the self propagation of turbulent reaction waves characterized by the displacement speed, an evolution equation for S d conditionally averaged on the reaction progress variable c was recently derived ) for a simplified reaction wave. Using this equation, shown as Eq. (13) below, the temporal evolution of the averaged displacement speed can be attributed to the terms on the right hand side (rhs) which represent different physical effects. Each term in this equation can be numerically extracted from DNS. It was demonstrated  using several DNS cases that the left hand side (lhs) term match well with the sum of all rhs terms even for a turbulent case with a high Karlovitz number.
The previous paper  focused on theoretical aspects of the equation derivation and its numerical verification. In this work we turn the focus toward the physical interpretation of the different terms that appear in the equation. We will examine their behaviors in a turbulent reacting wave, both before and after it evolves into statistical equilibrium, with particular emphasises on isosurface segments that propagate with negative displacement speed. As shown later, the displacement speed equation involves some high order terms that may lead to numerical complications if multiple factors affecting a realistic flame are simultaneously present. We therefore limit this work to a study of constant-density reaction waves where the thermal expansion effect is inhibited. While such a study is relevant for fundamental topics such as turbulent mixing in passive scalar, this work is also aimed to constitute a basis for future investigation of the displacement speed in more complicated cases where additional realistic effects (e.g. combustion heat release, multi-step chemistry, complex transport, etc.) can be taken into account. In this regard, it is worth remembering that we recently explored the evolution of various measures of local flame thickness (e.g. |∇c| or 1∕|∇c| ) conditionally averaged on the reaction progress variable c by analyzing DNS data obtained from (i) a constant-density single-reaction wave (the same DNS data are considered in the present work), (ii) two single-step-chemistry premixed turbulent flames with heat release and density variations, and (iii) two complexchemistry premixed turbulent flames (Yu et al. 2019). As shown in the cited paper, see p. 240, "the major features of the transient evolutions of all thicknesses are similar in different (constant density and single-step-chemistry, variable density and single-step chemistry, variable density and complex chemistry) cases". These recent results obtained for the conditioned |∇c| and 1∕|∇c| suggest that results obtained for the similarly conditioned S d in this study, i.e. in the simplest case (i), will be relevant not only to constant-density singlereaction waves, but also to premixed turbulent flames provided that the aforementioned limitations of this study are borne in mind.

Evolution Equation for Displacement Speed
A brief description will first be given for the derivation of the displacement speed evolution equation, both in an "averaged" as well as a "local" sense. More details on the derivation can be found in Yu and Lipatnikov (2019a). The explanations that follow will be directed toward physical interpretation of equation terms.
For an arbitrary quantity (x, t) an isosurface-following time derivative operator * * t can be defined by The composite velocity vector is defined as with n ≡ c∕|∇c| being the normal to the local isosurface. Using these definitions, Eq. (1) reduces to * * t c = 0 , reflecting the fact that c remains constant for any point that follows the isosurface.
In a general reacting flow characterized by Eq. (1), an instantaneous surface-average of a quantity conditioned on the isosurface c(x, t) =ĉ can be defined (Veynante and Vervisch 2002) as Here, ĉ is a reference value of the reaction progress variable on an isosurface, (c −ĉ) is the Dirac delta function and over-line denotes ensemble and volume averages taken simultaneously, i.e., ≡ lim M→∞ where M is the number of realizations in the ensemble, V is domain volume, and (i) pertains to the i-th realization.
Note that the Dirac delta function in Eq. (5) can be removed if the identity ∭ V |∇c| (c −ĉ)dx = ∬ S|̂c ,t ds is used (Maz'ja 1985;Kollmann and Chen 1994;Vervisch et al. 1995). This gives Here, S|̂c ,t is the isosurface defined by c(x, t) =ĉ whose total area is equal to A|̂c ,t ≡ ∬ S|̂c ,t ds . The long-hat operator over any expression represents the ensemble-averaged value, i.e.
The rate of change of the total area is related to the stretch rate by t A|̂c ,t = ∬ S|̂c ,t Kds , or written with surface-average since the stretch rate, defined as is well-known (Pope 1988;Candel and Poinsot 1990) to control the rate of change in the area (ds) of an infinitesimal element of an iso-scalar surface, i.e. K = 1 ds * * t ds . Here a t ≡ ∇ ⋅ u − n i n j ( u i ∕ x j ) is the tangential strain rate. The time derivative of Eq. (6) can now be expanded using the product rule to yield [see Ref. Yu and Lipatnikov (2019a) for a more formal derivation] Equation (9) can be rewritten to give a general evolution equation for surface-averaged value of the quantity , which holds on all isosurfaces ĉ ∈ (0, 1) and at time instants and for an arbitrary variable . If is the displacement speed, S d , an equation for surface-averaged S d can then be obtained as long as * * t S d is known. For this purpose, we assume a general diffusion term of the form = 1 ⋅ ( D c) and rewrite Eq.
(2) as with and D being density and diffusivity respectively. Now make the two assumptions that D =constant and that the reaction rate depends solely on c, which is likely to be valid for low Mach number globally adiabatic unity Lewis number flames. Apply the operator * * t to Eq. (11), to get because * * t = c * * t c = 0 . Equation (12) indicates that the isosurface-following rate of change in S d is attributed to two terms T 1 and T 2 : the term T 1 is a contribution from the isosurface-following rate of change in |∇c| , or in other words, decreasing/increasing separation between neighboring isosurfaces; while the term T 2 is due to the isosurface-following rate of change in diffusion. We then substitute Eq. (12) into Eq. (10) with = S d to get where the last averaged term is 1 which is related to a stretch-rate-induced difference between the time derivative of the isosurface-averaged value and the averaged isosurface-following derivative, i.e.
In the appendixes B and C of Ref. Yu and Lipatnikov (2019a) it is shown that the two terms in T 1 and T 2 can be recast into a form containing only spatial derivatives which facilitates their extraction from numerical simulations: Expand u * in Eqs. (15, 16) using Eq. (4), then Eq. (12) for * * t S d can be re-expressed as 1 3 In Eq. (17), the term T u comprises all contributions to * * t S d due to non-stationary flow motion, while four of its constituent parts, T u0 , T u1 , T u2 and T u3 , denote contributions due to flow-dilatation, tangential flow strain, velocity diffusion, and the product of c-diffusion tensor with velocity strain rate tensor, respectively. In fact the link between S d and the fluid velocity field through dilatation rate has been reported in previous works of Chakraborty and Cant (2004), and Chakraborty and Cant (2005). The term T S d in Eq. (17) denotes the remaining contribution to * * t S d due to non-uniform spatial distribution of S d . Comparing Eq. (12) with Eq. (17) yields It might be worth noting, that although the definition of S d by Eq.
(2) does not involve velocity, the S d -change following a local isosurface element, * * t S d , receives a direct contribution from velocity through the term T u . To better understand the physical meaning of two terms T u and T S d in Eq. (17), we consider a scenario in a periodic channel starting with a perfectly-planar, steady one-dimensional (1D), unstretched reaction wave solution imposed on a turbulent flow field with |u| ≠ 0 (for simplicity, we further assume that the flow has constant density, this scenario will be adopted in the later results section). The term T S d may be viewed as a history contribution to S d , originating from the velocity field. In fact, any later evolved non-zero value of T S d | 0<t≤∞ can be traced back in time to an initial nonzero value of T u | t=0 , as compared to a zero valued T S d | t=0 resulting from the initial 1D steady wave solution with a constant displacement speed. Therefore, in terms of contribution to the change of isosurface propagation speed according to Eq. (17), T u may be understood as an "immediate" contribution exerted by the velocity field through interaction with a temporally accumulated non-uniform distribution in S d . On the other hand, we may consider a second evolution scenario starting from a fully-developed turbulent reaction wave with its flow field being forced to stationary yielding zero T u for t ≥ 0 . Such a setup S d still evolves according to Eq. (17), however, with a single non-zero term T S d . A possible solution of such a wave evolution after long enough time is t S d | t ∞ = 0 and T S d | t ∞ = 0 , i.e. the wave being "attracted" toward to a planar steady-1D wave solution with a uniform S d . We also note that the terms T u and T S d defined in (17) do not explicitly involve the density, therefore the density variation can only influence the above two terms indirectly, through (i) density dependence of the flow field therefore affecting T u , e.g. the term T u0 vanishes in a constant-density flow, and (ii) historical contribution through T S d .
It might also be of interest to identify the velocity contribution to changes in the "averaged" displacement speed, t ⟨S d ⟩ s . This can be achieved by taking into account the additional velocity contribution in ⟨T 3 ⟩ s . Inserting Eq. (8) into Eq. (14), Eq. (13) can then be rewritten as where N u summarizes all the direct velocity contribution to t ⟨S d ⟩ s . We note here that the above derivations were performed for the simplified case of constant density, but the derivation is also valid in the more general case of D=const where the density is not required to be constant. When dealing with variable density reacting flow, Eq. (1) is often expressed in an equivalent form as by utilizing the continuity equation t + ∇ ⋅ ( u) = 0 . From the perspective of combustion modeling it appears more appropriate to treat the composite density-weighted displacement speed, ( S d ) , in the rhs of Eq. (20). For instance, in a steady 1D laminar flame with density variation ( S d ) is constant in space while S d is not. The isosurface-following derivative of ( S d ) can be calculated as * * t ( S d ) = S d * * t + * * t S d = * * t S d by assuming that is a function that only depends on c, which is again a reasonable assumption for low Mach number globally adiabatic unity Lewis number flames. This gives * * t = c * * t c = 0 and the time derivative of averaged ( S d ) simply becomes These simple relations suggest that the derived S d -equations can be directly applied to varying-density flames by adding a density scaling. Moreover, as the current analysis considers constantdensity conditions, the transport equation of ⟨ S d ⟩ s can be obtained by multiplying both sides by . It is worth noting that both ⟨ S d ⟩ s and ⟨S d ⟩ s are necessary for the FSD-based closure of the mean reaction rate and the terms of the FSD transport equation (Chakraborty andCant 2007, 2009). Thus, it is fundamentally important to consider the transport equations of ⟨S d ⟩ s and ⟨ S d ⟩ s due to their modelling importance instead of the transport equation of instantaneous S d .
We also want to point out a second route for density variations in flames to affect the S d -evolution; through the thermal-expansion acting on the flow to inflict T u . Consider the aforementioned example 1D steady flame whose velocity field must change across the flame. T u and its three components that involve velocity derivative across flame, T u0 , T u2 and T u 3 , are already non-zero before applying any disturbance. (Note that, in contrast T 1 and T 2 are zero in the same example). Moreover when considering the evolution equation (13) for average S d , its three rhs terms ⟨T 1 ⟩ s , ⟨T 2 ⟩ s and ⟨T 3 ⟩ s contains different velocity con- respectively; the thermal-expansion therefore exert a direct influence to term ⟨T 1 ⟩ s through the flow dilatation T u0 , and as well as play an indirect modulation of terms ⟨T 2 ⟩ s and ⟨T 3 ⟩ s . The additional mechanism for density variations in flames to affect S d suggests that care must be taken when extrapolating results obtained from constant-density cases. At the same time, it is worth noting that while various thermal expansion effects were revealed in premixed turbulent flames, the latest review articles (Lipatnikov and Chomiak 2010;Sabelnikov and Lipatnikov 2017) do not refer to any study that documents a substantial influence of thermal expansion on the density-weighted displacement speed ( S d ) . This fact, as well as the recent DNS study (Yu et al. 2019) discussed in the end of the previous section, suggest that results obtained from a constant-density turbulent reaction wave can be useful for understanding the influence of turbulence on a premixed flame.

3
2019b; Elperin et al. 2016;Sabelnikov et al. 2019a, b). In those DNSs, propagation of a statistically 1D, initially planar, single-reaction wave in a homogeneous, isotropic, statistically stationary forced turbulence was simulated in the case of a dynamically passive wave, i.e. a wave that does not change the fluid density and viscosity and does therefore not affect turbulence. The invoked simplifications allowed us to sample more statistics, as will be discussed later, and to investigate a large number of substantially different cases. Moreover, the simplifications significantly facilitate analysis and interpretation of the DNS data. Numerical results reported in the following are relevant largely to constant-density turbulent reacting flows, however those results also provides a basis for further investigation of more realistic premixed flames.

DNS Database
Since the DNS attributes are discussed in detail elsewhere (Yu et al. , 2015(Yu et al. , 2019Lipatnikov 2017a, b, 2019b;Elperin et al. 2016;Sabelnikov et al. 2019b), we will restrict ourselves to a brief summary of the simulations.
The propagation of a single-reaction wave is governed by the equation where is the reaction rate, R is a constant reaction time scale, = 6 and Ze = 6 in order for the rate W to depend on c in a non-linear manner. Equation (21) is a simplified form of Eq. (1) with = D∇ 2 c and = W. The wave propagates in a forced, homogeneous and isotropic turbulence described by the Navier-Stokes equations where p is the pressure and the density and kinematic viscosity are constants. Therefore, the flow is not affected by the wave, as mentioned earlier. A function is used to maintain the turbulence intensity by applying energy forcing at low wavenumbers (Lamorgese et al. 2005).
The wave evolves in a rectangular box with size of x × × represented using a uniform grid of N x × N × N cubic cells. The boundary conditions are periodic in all directions. In other words, when the reaction wave reaches the left boundary ( x = 0 ) of the computational domain, the reaction wave enters the domain through its right boundary ( x = x ). This method allows greatly improved sampling of statistics by simulating many cycles of wave propagation through the computational domain, but the method may only be used in the case of = const and = const and provided that the mean wave brush thickness is smaller than the length of the computational domain. These constraints are satisfied in all present simulations.
An initial turbulence field is generated by synthesizing prescribed Fourier waves  with an initial rms velocity u 0 and the integral scale 0 = ∕4 . The initial turbulent Reynolds number Re 0 = u 0 0 ∕ can be adjusted by changing the domain width . Subsequently, a non-decaying incompressible turbulent field is obtained by integrating Eq. (23). At the same Re 0 , turbulent fields characterized by the same rms velocity u � ≈ u 0 , but different longitudinal integral length scales L 11 are generated (Yu and Lipatnikov 2017a) by appropriately adjusting following Lamorgese et al. Lamorgese et al. (2005).
The governing equations are solved using an in-house DNS solver (Yu et al. 2012) developed for low Mach number reacting flows.
In the present work, both fully-developed (statistically stationary) and transient reaction waves are studied. Simulations are started from the pre-computed laminar-wave profile of c L ( ) initially ( t = 0 ) released at x 0 = x ∕2 and the subsequent evolution of the field c( , t) is simulated by solving Eq. (21). Sampling of the fully-developed statistics starts when a statistically stationary state has been reached and sampling is done for a duration of at least 50 0 t . In order to study transient turbulent reaction waves, several copies of the same pre- Totally 45 cases characterized by the Damköhler number Da = t ∕ F = 0.01-24.7, the Karlovitz number Ka = F ∕ = 0.36-587, u � ∕S L = 0.5-90, and L 11 ∕ F = 0.39-12.4 were simulated, with a few cases designed to show weak sensitivity of computed results to grid resolution, ∕L 11 , etc. (Yu and Lipatnikov 2017a). Here, F = F ∕S L is the wave time scale, t = L 11 ∕u � and = ( ∕ ) 1∕2 are integral and Kolmogorov time scales of the turbulence, respectively, and = 2 S ij S ij with S ij ≡ ( x j u i + x i u j )∕2 is the dissipation rate averaged over the computational domain and for at least 50 0 t after turbulence has evolved into statistical stationarity.

Case Setup
In the present paper, the focus is on the results obtained in a representative highly turbulent and well resolved case, whose major characteristics are reported in Table 1 (case C), where Pe = u � L 11 ∕D is the turbulent Péclet number and a ratio of F ∕ x characterizes the grid resolution in terms on number of grid points per the laminar wave thickness. Moreover, in case C, L 11 ∕ = 0.11 , 0 t ∕ t = 2.3 , and ∕ x = 1.1 . Here, = ( 3 ∕ ) 1∕4 is the Kolmogorov length scale. Figure 1 shows a snapshot with multiple coexisting c t -fields and one c s field in case C.
One issue of concern when considering a highly turbulent reaction wave is the appearance of zero-gradient points (Gibson 1968) with |∇c|(x, t) = 0 . Here we should point out that zero gradient points have been excluded from the definitions of surface average in Eq. (5), see Yu and Lipatnikov (2019a), and Yu and Lipatnikov (2019c), however certain quantities of interest, e.g. S d and K , may locally grow unboundedly large in the neighborhood of a zero-gradient point, which poses a challenge for numerical calculation of surface-averaged quantities. To explore the eventual influence of such points and their neighborhoods on the accuracy of evaluation of various terms in the evolution equation, two supplementary 2D cases (case A and B) were designed.
Case A is largely identical to case C, but the turbulent field is replaced with a frozen shear flow, i.e. u(x, y, z, t) = −u 0 cos(2 y∕ ) , v = w = 0 , and the momentum Eq. (23) is not solved. In case A the isosurfaces are only curved, there is no zero-gradient points in the computational domain. Characteristics of case A, reported in Table 1, are calculated using L 11 = ∕2 , u � = u 0 , and t = t = L 11 ∕u 0 .
Case B is designed to have a small number of zero-gradient points in the computational domain. Similarly to case A, case B is also based on a frozen velocity field u(x, y, z, t) = u �� (x, y) , v(x, y, z, t) = v �� (x, y) , w = 0 , and ⋅ u �� = 0 . The field u ′′ represents a 2D, zero-mean, spatially fluctuating velocity field generated using a reduced version of the method for synthesizing the initial 3D turbulence. Characteristics of case B, reported in Table 1, are calculated by analyzing the 2D velocity field ( u 0 and L 11 ) and using ≈ t = L 11 ∕u 0 to evaluate Ka. Case B mimics a moderately disturbed reaction wave propagating through vortexes and allowing the appearance of zero-gradient points during collisions of reaction-waves. Comparison of results computed in cases A (no zero-gradient points), B (a few zero-gradient points), and C (no restrictions on appearance of the zerogradient points) offers an opportunity to estimate the influence of such points on the accuracy of evaluation of various terms appearing in the S d -evolution equation.
In the two frozen-velocity cases A and B, simulations of multiple transient waves c t are performed largely similarly to case C (further details on the difference are discussed in Yu and Lipatnikov (2019a)) but the duration of transient sampling is changed from 2 0 t to 2 F . It is worth noting that the transient data not only is of interest in itself because the vast majority of premixed turbulent flames are developing flames, as discussed in detail elsewhere (Lipatnikov and Chomiak 2002;Lipatnikov 2012), but also offer an opportunity to control the following numerical issue. As discussed earlier, in very intense turbulence, a scalar isosurface of c(x, t) =ĉ can become very complicated and can contain multiple zero-gradient points, close to which a quantity of interest can be unboundedly large. As a result, certain surface-averaged might be unbounded. Since the transient waves c t m (x, t) =ĉ begin their evolution from a regular flat initial surface, monitoring evolution of (i) isosurfaces of the transient fields c t m (x, t) =ĉ and (ii) the relevant surface-averaged quantities offer an opportunity to detect any anomaly in the developing surface-averaged

Results and Discussions
DNS data has been used to examine the evolution equation for surface-averaged S d , Eq. (13) or an equivalent Eq. (19). These equations should hold for all isosurfaces of c(x, t) =ĉ such that ĉ ∈ (0, 1) and at all time instants t. Our interest is first placed on how various terms contribute to the change in averaged displacement speed.
In each simulated case (A, B, or C), the following terms were computed: (a) the left hand side (lhs) term, i.e. the time derivative of ⟨S d ⟩ s , (b) the three rhs terms in Eq. (13), and (c) the sum of all the rhs terms, i.e. ∑ rhs . Note that the lhs term t ⟨S d ⟩ s (red dots in all the figures) is evaluated by first obtaining the sequence of transient-evolving ⟨S d ⟩ s calculated at 20 sampled time instants t i = (i 2 ∕200) * F for i = (1, .., 20) and then applying a discrete approximation of the time derivative. Figure 2 shows the temporal evolution of these terms normalized using S L ∕ 2 F and conditioned on three values of ĉ representing the preheat zone c(x, t) = 0.1 , a middle zone c(x, t) = 0.5 and the reaction zone c(x, t) = 0.88 where W(0.88) = max(W) . For comparison of transient statistics between cases, the time t is normalized by * F which is set to F in the two cases A and B and 0 t ( ≈ F ∕22.5 ) in case C. Figure 3 shows the ĉ-isosurfacedependence of all terms at three representative time instants: an early instant at 0.045 * F , a middle instant 1.125 * F , and the fully developed state at t ∞ . Furthermore, the averaged displacement speed ⟨S d ⟩ s is shown in Fig. 4 at similar conditions as those in Figs. 2 and 3. Figures 2 and 3 show that the difference between t ⟨S d ⟩ s and ∑ rhs is small for all cases A-C for almost all isosurfaces and time instants. This fact supports Eq. (13) and also shows that all terms have been computed with high enough numerical accuracy.
Before proceeding to a discussion of Eq. (13), one may notice in Fig. 4 that the averaged displacement speed is negative in the reaction zone in the high-Ka case C during the transient wave evolution. Nevertheless, for the fully developed reaction wave in all three cases, the averaged displacement speed recovers to be positive throughout all wave zones and its value stays slightly above unity, see Fig. 4c. The above observations of opposite signs of displacement speed agree quialitatively with DNS data by Chakraborty (2007) who simulated single-step chemistry turbulent flames associated with the thin reaction zone regime and reported high probability of finding negative local S d but positive averaged displacement speed throughout the entire flame brushes. 1 3

Trends During Transient Wave Evolution
For temporal evolution of the terms in Eq. (13) one may notice that the early evolution of the time derivative term t ⟨S d ⟩ s is different in the preheat zone ( ĉ = 0.1 ) than in the reaction zone ( ĉ = 0.88 ). In the preheat zone the term t ⟨S d ⟩ s rises from an initial zero to a positive peak, seen for 0 < t∕ * F ≈ 0.3 in Fig. 2a, d for cases A and B, as well as for 0 < t∕ * F ≈ 0.1 for case C in Fig. 2g. On the contrary the term t ⟨S d ⟩ s in the reaction zone ĉ = 0.88 drops from an initial zero to a negative peak for the same early duration as shown in Fig. 2c, f, i. These early trends are further reflected by the negative slope in the profile of t ⟨S d ⟩ s along ĉ , which is shown at a representative early instant t∕ * F = 0.045 in Fig. 3a, d, g. The observed early negative slope in t ⟨S d ⟩ s is mainly attributed to the diffusion contribution term ⟨T 2 ⟩ s ; compare ⟨T 2 ⟩ s with t ⟨S d ⟩ s in Fig. 2a, d, g. Among the three time-dependent terms ⟨T 1 ⟩ s , ⟨T 2 ⟩ s and ⟨T 3 ⟩ s shown in Fig. 2, the overall temporal evolution of t ⟨S d ⟩ s is best mimicked by the diffusion contribution term ⟨T 2 ⟩ s , even though significant difference between the other two terms always exists at t ∞ . The two terms ⟨T 1 ⟩ s and ⟨T 3 ⟩ s tend to develop values of significant magnitude with opposite sign during later time of wave evolution, compare the two curves plotted in triangles both in Fig. 2 for t∕ * F > 0.3 and in Fig. 3. Each of the two terms ⟨T 1 ⟩ s and ⟨T 3 ⟩ s conditioned at the reaction zone ( ĉ = 0.88 ) can even change its sign twice during its evolution, as shown in Fig. 2.i for the high Ka case C. Such an oscillating behavior is more clearly observed by comparing Fig. 3g-i, which show that the overall slope of the profile ⟨T 1 ⟩ s along ĉ at three time instants changes from positive to negative and then finally back to positive. Moreover, unlike being One observation related to the above phenomenon is that the maximum magnitude of the variation in all displayed terms during the transient evolution conditioned at the middle zone, shown in Fig. 2h for ĉ = 0.5 , is much smaller than the ones conditioned at the two other zones, shown in Fig. 2g, i for ĉ = 0.1 and ĉ = 0.88 , respectively. Such observation may be related to the fact that in a reference steady laminar wave, the denominator term |∇c| in the S d -definition Eq. (2) attains its largest value at the middle ĉ , which is less prone to perturbation than values at two ends of ĉ . For this reason, further conditioned analysis in Sect. 4.2.1 will be based on the isosurface at ĉ = 0.5.
Regarding the velocity-contribution to t ⟨S d ⟩ s , Fig. 3a, d, g show for three cases A-C that the ĉ-profile of N u at t∕ * F = 0.045 matches well with the one for t ⟨S d ⟩ s . The match is particularly good for ĉ < 0.3 in cases A and B in which the "effective" evolution time t∕ F is significantly shorter than in case C. This reflects a simple fact that it is the nonuniform velocity field which kick-starts the initial variation in ⟨S d ⟩ s � t=0 and continues for a short while to play a major role for the change in ⟨S d ⟩ s during early transient wave evolution. Compared with t ⟨S d ⟩ s , the term N u changes from being larger in the preheat zone ( ĉ < 0.4 ) to being smaller toward the reaction zone ( ĉ > 0.8 ) in all cases A-C for later times ( 1.125 * t ≤ t ≤ t ∞ ). One interpretation of the observation at t ∞ is that, for a hypothetical zero-N u scenario in which the flow field of a fully developed wave is suddenly changed to a uniform flow, the averaged displacement speed ⟨S d ⟩ s must respond by decreasing in  . 3 The ĉ-dependencies of the lhs term, three individual rhs terms, and sum of all rhs terms from Eq. (13), as well as a velocity contribution term N u in Eq. (19) for cases A (top row), B (middle row) and C (bottom row). Results computed at t = 0.045 * F and t = 1.125 * F are plotted in the left and middle columns, respectively, with * F = F in cases A and B or * F = F ∕22.5 in case C. All terms are normalized with S L ∕ 2 F the preheat zone and increasing in the reaction zone respectively, since t ⟨S d ⟩ s = N S d and its rhs term stores a negative copy of the non-zero valued N u prior to the imposed velocity modification. The observed sign-flipping variation of N u along ĉ is (observed in our data, not shown) largely attributed to the component T u 3 which involves a diffusion tensor D∇∇c∕|∇c| , which in a 1D flame becomes a normal-diffusion contribution to displacement speed (i.e D∇ 2 c∕|∇c| ) and must flip its sign along a monotonic 1D solution of c L (x).

Trend in Fully Developed Wave
When a reacting wave has evolved into a statistically stationary state, or is fully developed, the time derivative of any statistic must be zero. Therefore, for Eq. (13) at t ∞ , the lhs term as well as the sum of all rhs terms must be zero for all ĉ ∈ (0, 1) . This is clearly seen for the three cases A, B and C in Fig. 3c, f, i (blue open circles and red dots). However, each of the three individual rhs terms, ⟨T 1 ⟩ s , ⟨T 2 ⟩ s and ⟨T 3 ⟩ s , are not necessarily zero at t ∞ . Note that ⟨K⟩ s �̂c ,t ∞ = 0 according to statistical stationarity and Eq. (7) [confirmed by the data shown in Fig. 7 in Ref. Yu and Lipatnikov (2019a)].
In fact, Fig. 3c, f, i show, in the fully developed cases A-C, that ⟨T 3 ⟩ s stays negative and ⟨T 1 ⟩ s stays positive for all ĉ . The value of ⟨T 2 ⟩ s generally rises with an increase in ĉ , the value of ⟨T 2 ⟩ s in the reaction zone ( ĉ > 0.8 ) is positive, however in the far upstream preheat zone ( ĉ = 0.1 ) its value drops slightly below zero when the reaction wave becomes more disturbed by the flow (compare cases A and C). In the two low Ka cases A and B, the maximum absolute values of all three (normalized) terms at t ∞ across all ĉ stay below 4; in the high Ka case C this value significantly increases above 1000. This can be explained by that the Kolmogorov eddy should be the dominant cause for change in isosurface movement (this will become clear in following discussions when the Ka is used for normalization of relevant terms in case C).

Averages Conditioned at Fast/Slow Isosurfaces
To better understand the causes for changes in the displacement speed we define a S d -conditioned isosurface-average as Following a similar procedure as from Eqs. (5) to (6), the average can be re-expressed as which is related to the average over the entire isosurface by where an isosurface-conditioned probability density function (pdf) is defined for a quantity by with pdf( , |∇c|, c;t, x) being a point-wise joint pdf [see appendix in Ref. Yu and Lipatnikov (2019c)] and n s is a normalization constant to enforce ∫ +∞ −∞ pdf s ( )d = 1. Figure 5 shows the isosurface-conditioned displacement-speed pdf, i.e. pdf s (S d )|̂c ,t ∞ , and the S d -dependent surface-average of curvature , i.e. [∇ ⋅ n] s |Ŝ d ,ĉ,t ∞ , for three fully-developed cases A, B and C on a representative isosurface with ĉ = 0.5 . For reference, raw sampling values of curvature ∇ ⋅ n and its isosurface-average, i.e. ⟨∇ ⋅ n⟩ s �̂c ,t ∞ are also shown in Fig. 5. Similarly, Fig. 6 shows for three fully-developed cases A-C, the S d -dependent isosurfaceaverages [ ] s |Ŝ d ,ĉ,t ∞ for each ∈ [∇ ⋅ n, T 1 + T 2 , T 1 , S d K, T u , T u 1 , T u 2 ] , i.e. for terms relevant to changes in the displacement speed of an isosurface element. The averaged terms ⟨ ⟩ s �̂c ,t ∞ are also plotted in Fig. 6 for reference. Note, that all quantities displayed in the two Figs. 5, 6 for the two low-Ka number cases A and B are appropriately normalized based on a length unit of laminar flame thickness F and a velocity unit of laminar flame speed S L . However, for the high Ka case C, the normalization was performed using the Kolmogorov length scale ∝ F ∕

√
Ka and the Kolmogorov velocity scale ∕ ∝ S L √ Ka to indicate an important role played by the Kolmogorov-scale eddies. Figure 5 shows that, although the curvature averaged over the whole isosurface ⟨∇ ⋅ n⟩ s is close to zero in all three cases, the S d -dependent surface-average of curvature [∇ ⋅ n] s |Ŝ d ,ĉ,t ∞ deviates from zero and its value generally rises with an increase in S d , reflecting the well-known fact that the displacement speed contains a direct contribution from curvature, e.g. S d = S c d + S n d + S W d with S c d = D∇ ⋅ n (assuming constant D , the rest becomes S n d = D∇ n ∇ n c∕|∇c| and S W d = ∕|∇c| where ∇ n = n ⋅ ). For cases A and B with a low Ka, Fig. 5 show pdf s (S d ) peaked at S L , which indicates that a major part of the isosurface has a displacement speed around S L ; there is also a small amount of fast propagating isosurface with S d significantly larger than S L . Those isosurfaces are on average associated with large curvature, (cf. [∇ ⋅ n] s ≫ −1 F for S d ≫ S L indicated by black lines in Fig. 5 for cases A and B). As shown by the cyan circles in Fig. 5 for case B, for faster propagating isosurfaces at larger S d the curvature values develop a wider spread and can even become negative. Nevertheless S d is positive due to other contributions to the displacement speed, i.e. S n d and S W d . In all cases, there exist a probability of finding "retreating" isosurface elements that propagate with negative S d and they are generally associated with negative curvature. The occurrence of negative S d become significant in the high-Ka case C shown by pdf s (S d ) in Fig. 5c. As far as the present simulations of negative displacement speed are concerned, these results obtained in the simplest case of a constant-density, single-reaction wave are consistent with earlier computations of negative S d in 2D DNS of complex-chemistry premixed turbulent flames (Gran et al. 1996;Echekki and Chen 1996;Peters et al. 1998;Echekki and Chen 1999;Cecere et al. 2018), 3D DNS of single-step-chemistry premixed turbulent flames (Chakraborty and Cant 2005;Chakraborty 2007), and 3D DNS of complex-chemistry premixed turbulent flames (Wang et al. 2017;Luca et al. 2019;Sankaran et al. 2015;Cecere et al. 2016;Trisjono et al. 2016).
It can be of particular interest to examine various terms contributing to changes in propagating speed of an isosurface element, see Fig. 6. For instance, according to Eq. (12) the summation term of T 1 + T 2 is equal to * * t S d . Its S d -conditioned isosurface-average, i.e T 1 + T 2 s |Ŝ d means the rate of change in propagating speed averaged for all isosurface elements propagating at a same speed of Ŝ d . Figure 6 shows for three cases A-C a clear trend of T 1 + T 2 s |Ŝ d rise toward larger positive value of Ŝ d , which, together with curvature shown in Fig. 5, suggests an accelerated advancement of positively curved isosurfaces associated with a typical "consumption" scenario of cusps or small fuel pockets surrounded by products.
For retreating isosurface with negative S d , the average T 1 + T 2 s stays negative in cases B, suggesting an acceleration of a "quenching" process for a product pocket surrounded by fuel. Such a process is more evident at large negative Ŝ d for case C in which significant amount of product can be transported over to the fuel side due to intensive turbulent mixing. This behaviour is also in agreement with curvature dependence of displacement speed reported by previous variable density single-step (Chakraborty 2007;Chakraborty and Cant 2004) and detailed chemistry Peters et al. 1998;Echekki and Chen 1999) DNS results.
According to Eq. (26), for any quantity the two averages [ ] s |Ŝ d and ⟨ ⟩ s are related by a weight of pdf s (S d ) . In the two low-Ka cases A and B shown in Fig. 6, the value of ⟨T 1 + T 2 ⟩ s differs significantly from the value of T 1 + T 2 s |Ŝ d at Ŝ d ≈ S L where pdf s (S d ) reaches its peak, seen in Fig. 5. This indicates that the combined T 1 + T 2 term averaged over the entire isosurface receives significant contribution from a small amount of . For the high-Ka number case C, the absolute value of T 1 + T 2 s |Ŝ d is still large at large positive and negative values of Ŝ d , however the whole average ⟨T 1 + T 2 ⟩ s stays small as a result of overall cancellation T 1 + T 2 s weighted by a pdf s (S d ) which is symmetric around zero. Figure 6 shows that at ĉ = 0.5 the term T 1 s |Ŝ d ,ĉ approximately follows the term T 1 + T 2 s |Ŝ d ,ĉ in the three fully developed cases, which can be attributed to a low magnitude of ∇ 2 c at ĉ = 0.5 . The term T 2 becomes more important near reaction zone ( ĉ > 0.8 ), see ⟨T 2 ⟩ s �̂c ,t ∞ in Fig. 3f. For the high-Ka case C, Fig. 6 . This corresponds to a negative value of [K] s , suggesting a net reduction in the summed area of fast-propagating isosurfaces. This observation is not in conflict with the fact that there is no change in the area of fully-developed isosurfaces (due to statistical stationary). The inset in Fig. 6 shows that [K] s , conditioned at the low-speed isosurface elements with −1 <Ŝ d ∕(S L √ Ka) < 1 , changes to become positive for area production of low speed isosurfaces. Now we examine the velocity contribution to the changes of S d . First, retreating isosurfaces with a large negative S d are promoted by the velocity field, as indicated in Fig. 6 by a negative T u s at negative S d . Compared with aforementioned negative term of T 1 + T 2 s at negative S d , the velocity-contribution term T u s has a larger magnitude in the small Ka case B, and has a smaller magnitude in the large Ka case C. This suggests that the importance of the "direct" role played by velocity in creating retreating isosurfaces [compared with its "indirect" historic role played through T S d in Eq. (18)] decreases with an increase in Ka. Regarding the isosurfaces with large positive S d (i.e. S d > 3S L max( √ Ka, 1) ), the direct Fully-developed [ ] s |Ŝ d ,ĉ,t ∞ , ⟨ ⟩ s �̂c ,t ∞ at ĉ = 0.5 for each ∈ [∇ ⋅ n, T 1 + T 2 , T 1 , S d K, T u , T u1 , T u2 ] in cases A (left), B (middle) and C (right). All terms are normalized by S 2 L ∕ F (A, B) and Ka 2 S 2 L ∕ F (C); normalization of S d is same as in Fig. 5 velocity contribution T u s stays positive, therefore, also promotes their existence. However, in all three cases this direct velocity contribution is small compared to the combined terms T 1 + T 2 s , or compared to its indirect contribution through T S d . Furthermore, in the two "turbulent" cases B and C, the aforementioned promotion effect to retreating isosurfaces is largely attributed to the term T u 3 , since T u 3 = T u − T u1 − T u2 in constant-density wave and the magnitudes of T u1 s and T u2 s are relatively small compared to T u s at large negative S d , see Fig. 6.

Summary and Conclusion
For turbulent premixed reacting flow represented by a simplified transport equation for a reaction progress variable, the evolution equation for the displacement speed S d characterizing self-propagation of a reaction progress isosurface was derived in a recent work. In the present work this equation is analyzed with an emphasis on physical interpretation of different factors that affect the displacement speed. Our analysis reveals multiple effects: (i) T 1 , the isosurface-following rate of change in reaction surface density, (ii) T 2 , the isosurface-following rate of change due to diffusion, and (iii) T 3 , a stretch-rate-induced difference between averaged isosurface-following derivative and time derivative of the isosurface averaged value, and (iv) T u or N u , a direct velocity contribution which can be further split into multiple components. The analysis was then applied to study averaged statistics of displacement speed for a planar reaction wave propagating in homogeneous isotropic constant-density turbulence, simulated using a DNS approach. We examine both the transient process when the initial planar wave is disturbed by turbulence, as well as the fully developed wave after all statistics have evolved to a stationary state. During the transient phase, the surface-averaged displacement speed ⟨S d ⟩ s in the reaction zone (preheat zone) experiences a significant decrease (increase). Under highly turbulent conditions, the transient value of ⟨S d ⟩ s in the reaction zone becomes negative. Among the three terms ⟨T 1 ⟩ s , ⟨T 2 ⟩ s and ⟨T 3 ⟩ s , the transient trend of t ⟨S d ⟩ s is largely resembled by the diffusive term ⟨T 2 ⟩ s . The initial variation in ⟨S d ⟩ s is promoted by the velocity contribution N u whose sign changes from positive in the preheat zone to negative in the reaction zone.
In the fully developed wave the averaged displacement speed does not change and its value stays positive and approximates the laminar wave speed. The fully developed wave does yield a zero sum of ⟨T 1 ⟩ s + ⟨T 2 ⟩ s + ⟨T 3 ⟩ s , with ⟨T 1 ⟩ s being positive, ⟨T 3 ⟩ s negative and a ⟨T 2 ⟩ s having the smallest magnitude in the preheat zone. The trend holds no matter the wave evolves in a highly turbulent flow or a simple shear flow.
In a fully developed wave, local analysis of various terms conditioned on a representative isosurface reveals that there exists significant amount of isosurface elements with negative displacement speed in a highly turbulent reaction wave. In such a wave, the magnitude of the terms affecting S d can be taken to scale with the Kolmogorov velocity and length scales. Furthermore, these terms, when averaged over the entire isosurface, can deviate significantly from the corresponding averaged values further conditioned on a fast or slow speed of S d . The isosurfaces propagating with a large positive displacement speed are associated with small fuel pockets surrounded by products. The propagation of those isosurfaces undergoes net acceleration; similarly isosurfaces with net negative displacement speed experience a net negative acceleration. The acceleration process is promoted by the direct velocity contribution T u , mainly through the part of T u3 .