2.1. General evolution equation for surface-averaged quantities

We begin with some convenient definitions. A long-hat over any expression denotes the ensemble-averaged value, i.e.

(2.1)\begin{equation} \hat{\phi} \equiv \lim_{M \to \infty} \frac{1}{M} \sum_{i=1}^{M} \phi_{(i)}, \end{equation}

and an overline denotes ensemble and volume averages taken simultaneously, i.e.

(2.2)\begin{equation} \bar{\phi} \equiv \lim_{M \to \infty} \frac{1}{ M} \sum_{i=1}^{M} \frac{1}{V}\iiint\limits_{V} \phi_{(i)}(t,{\boldsymbol{x}}) \,\mathrm{d} {\boldsymbol{x}}, \end{equation}

where $M$ is the number of realizations in the ensemble, $V$ is the domain volume and $\phi _{(i)}$ pertains to the $i$th realization. An instantaneous surface average of a quantity $\phi$ conditioned on the iso-surface $c({\boldsymbol {x}},t)=\hat {c}$ is defined (Veynante & Vervisch Reference Veynante and Vervisch2002) as

(2.3)\begin{equation} \left \langle \phi \right \rangle_{{s}} |_{\hat{c},t} \equiv \overline{ \phi \lvert \boldsymbol{\nabla} c \rvert \delta(c-\hat{c}) } / \overline{ \lvert \boldsymbol{\nabla} c \rvert \delta(c-\hat{c}) }, \end{equation}

where $\hat {c}$ is a reference value of the reaction progress variable and $\delta (c-\hat {c})$ is the Dirac delta function.

Let $S|_{\hat {c},t}$ denote the iso-surface defined by $c({\boldsymbol {x}},t)=\hat {c}$ whose total area is equal to $A|_{\hat {c},t} \equiv \iint _{S|_{\hat {c},t}} \,\mathrm {d} s$. The Dirac delta function in (2.3) can now be removed by converting volume integrals into surface integrals using the identity of

(2.4)\begin{equation} \iiint\limits_V \phi \lvert \boldsymbol{\nabla} c \rvert \delta( c- \hat{c}) \,\mathrm{d} {\boldsymbol{x}} = \iint\limits_{S|_{\hat{c},t}} \phi \,\mathrm{d} s \end{equation}

(Maz'ja Reference Maz'ja1985; Kollmann & Chen Reference Kollmann and Chen1994; Vervisch et al. Reference Vervisch, Bidaux, Bray and Kollmann1995). We arrive at

(2.5)\begin{equation} \left \langle \phi \right \rangle_{{s}} |_{\hat{c},t} = \widehat{ \displaystyle \iint_{S|_{\hat{c},t}} \phi \,\mathrm{d} s } / \hat{A}|_{\hat{c},t}. \end{equation}

It is well known (Pope Reference Pope1988; Candel & Poinsot Reference Candel and Poinsot1990) that the stretch rate, defined as

(2.6)\begin{equation} \mathcal{K} \equiv a_t - S_d \boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{n}} , \end{equation}

controls the rate of change of the area of an infinitesimal surface element of an iso-scalar surface according to

(2.7)\begin{equation} \mathcal{K} = \frac{1}{ \mathrm{d} s} \frac{\partial^{*} }{\partial^{*}t} \mathrm{d} s. \end{equation}

The change of the total area is thus

(2.8)\begin{equation} \frac{\partial }{\partial t} A|_{\hat{c},t}= \iint\limits_{S|_{\hat{c},t}} \mathcal{K} \,\mathrm{d} s. \end{equation}

Application of the surface average operator to $\mathcal {K}$ yields

(2.9)\begin{equation} \left \langle \mathcal{K} \right \rangle_{{s}} |_{\hat{c},t} =\widehat{\displaystyle \iint_{S|_{\hat{c},t}} \mathcal{K} \,\mathrm{d} s } / \hat{A}|_{\hat{c},t} =\frac{1}{\hat{A}|_{\hat{c},t}} \frac{\partial }{\partial t} \hat{A}|_{\hat{c},t} . \end{equation}

Here, $a_t\equiv \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {u}}- a_n$ is the tangential strain rate and $a_n \equiv {\boldsymbol {n}} {\boldsymbol {n}} \boldsymbol {:} \boldsymbol {\nabla } {\boldsymbol {u}} = n_i n_j \nabla _j u_i$ is the normal strain rate, the nabla operator with a subscript denotes spatial derivatives as $\nabla _j \equiv \partial /\partial x_j$, and the summation convention applies to repeated indexes.

As shown by Yu & Lipatnikov (Reference Yu and Lipatnikov2019c), the time derivative of (2.5) can be expanded as

(2.10)\begin{equation} \frac{\partial \left \langle \phi \right \rangle_{{s}} }{\partial t} + \left \langle \phi \right \rangle_{{s}} \frac{1}{\hat{A}} \frac{\partial \hat{A} }{\partial t} = \frac{1}{\hat{A}} \widehat{\displaystyle \iint_{S} \frac{\partial^{*} \phi }{\partial^{*}t} \mathrm{d} s } +\frac{1}{\hat{A}} \widehat{\displaystyle \iint_{S} \phi \mathcal{K} \,\mathrm{d} s } . \end{equation}

Equation (2.10) can be rewritten to give a general evolution equation for the surface average of any quantity $\phi ({\boldsymbol {x}},t)$,

(2.11)\begin{equation} \frac{\partial \left \langle \phi \right \rangle_{{s}} |_{\hat{c},t} }{\partial t} = \left. \left \langle \frac{\partial^{*} }{\partial^{*}t} \phi \right \rangle_{{s}} \right|_{\hat{c},t} + \left. \left \langle \phi \mathcal{K} \right \rangle_{{s}} \right|_{\hat{c},t} - \left. \left \langle \phi \right \rangle_{{s}} |_{\hat{c},t} \left \langle \mathcal{K} \right \rangle_{{s}} \right|_{\hat{c},t}, \end{equation}

which holds on all iso-surfaces defined by $\hat {c}\in (0,1)$ and at all time instants.

2.2. Surface-following evolution equations for $S_d$ and its constituents

In the following we will derive surface-averaged evolution equations for the displacement speed $S_d$ and its constituents, i.e. the four right-hand side terms in

(2.12)\begin{equation} S_d = S_d^c + S_d^n + S_d^W + S_d^o = \frac{ \mathbb{D} + \mathbb{W} }{ \lvert \boldsymbol{\nabla} c \rvert }. \end{equation}

These terms are: a curvature (tangential diffusion) contribution

(2.13)\begin{equation} S_d^c \equiv \mathcal{D} \boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{n}} =\frac{ \mathcal{D} (\nabla^2 c - \nabla_n \nabla _n c) }{ \lvert \boldsymbol{\nabla} c \rvert } , \end{equation}

a normal diffusion contribution

(2.14)\begin{equation} S_d^n \equiv \frac{\mathcal{D} \nabla_n \nabla_n c }{ \lvert \boldsymbol{\nabla} c \rvert } = \frac{\mathcal{D} {\boldsymbol{n}} \boldsymbol{\cdot} \boldsymbol{\nabla} ( \lvert \boldsymbol{\nabla} c \rvert ) }{ \lvert \boldsymbol{\nabla} c \rvert }, \end{equation}

a reaction contribution

(2.15)\begin{equation} S_d^W \equiv \frac{ \mathbb{W} } { \lvert \boldsymbol{\nabla} c \rvert } , \end{equation}

as well as a contribution

(2.16)\begin{equation} S_d^o \equiv S_d -S_d^W - S_d^c - S_d^n = \frac{1}{\rho} \nabla_n (\rho \mathcal{D}) \end{equation}

due to variations of $\rho \mathcal {D}$ (often, $S_d^o$ is merged with $S_d^n$). Here, $\nabla _n\equiv n_j \nabla _j$ denotes the spatial derivative along the normal direction. Under the assumption that $\rho \mathcal {D}=\textrm {const}$ the last contribution vanishes, i.e. $S_d^o=0$.

Note that $\lvert \boldsymbol {\nabla } c \rvert$ appears in the denominator in each of the four quantities $\phi \in [ \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}, S_d^n/\mathcal {D}, S_d^W, S_d]$. Consequently, the surface-following derivative of these four quantities has the following general form:

(2.17)\begin{equation} \frac{\partial^{*} }{\partial^{*}t} \phi =\underbrace{-\phi \frac{\partial^{*} }{\partial^{*}t} \ln \lvert \boldsymbol{\nabla} c \rvert }_{ \mathbb{A} (\phi)} +\underbrace{\frac{1}{ \lvert \boldsymbol{\nabla} c \rvert } \frac{\partial^{*} }{\partial^{*}t} (\phi \lvert \boldsymbol{\nabla} c \rvert )}_{ \mathbb{B} (\phi)}. \end{equation}

Here, the term $\mathbb {A}$ comes from the denominator of $\phi$ and represents contribution from the surface-following rate of change in $\lvert \boldsymbol {\nabla } c \rvert$, or in other words, decreasing/increasing separation distance between neighbouring iso-surfaces. The term $\mathbb {B}$ comes from the numerator of $\phi$ and is associated with the surface-following rate of change of the ‘diffusion’ contained inside $(\phi \lvert \boldsymbol {\nabla } c \rvert )$. When $\phi$ designates either curvature or $S_d^n/\mathcal {D}$, the ‘diffusive’ nature of the numerator can be readily seen from the second spatial derivative of $c$ in (2.13) and (2.14), respectively. For $S_d^W$, if it is assumed that the reaction rate depends solely on $c$, the ‘diffusion’ vanishes, i.e.

(2.18)\begin{equation} \frac{\partial^{*} }{\partial^{*}t} \mathbb{W} =\frac{ \partial \mathbb{W} }{ \partial c } \frac{\partial^{*} }{\partial^{*}t} c = 0. \end{equation}

Finally, for the full displacement speed $S_d$, the additional assumption that $\rho \mathcal {D}=\textrm {const}$ with $\mathcal {D}$ depending solely on $c$ gives a standard Laplacian ‘diffusion’ of $\nabla ^2 c$. As a summary, the corresponding $\mathbb {B}$ terms under the aforementioned assumptions become

(2.19)\begin{gather} \mathbb{B} ( S_d ) = \frac{\mathcal{D}}{ \lvert \boldsymbol{\nabla} c \rvert } \frac{\partial^{*} }{\partial^{*}t} (\nabla^2 c ) , \end{gather}

(2.20)\begin{gather}\mathbb{B} \left( \frac{S_d^n}{\mathcal{D}} \right) = \frac{1}{ \lvert \boldsymbol{\nabla} c \rvert } \frac{\partial^{*} }{\partial^{*}t} \left(\nabla_n \nabla_n c \right), \end{gather}

(2.21)\begin{gather}\mathbb{B} \left(\frac{S_d^c}{\mathcal{D}} \right) = \frac{1}{ \lvert \boldsymbol{\nabla} c \rvert } \frac{\partial^{*} }{\partial^{*}t} \left(\nabla^2 c -\nabla_n \nabla_n c \right) \end{gather}

and $\mathbb {B} ( S_d^W )=0$.

In appendix A we further expand the terms $\mathbb {A}$ and $\mathbb {B}$ and derive surface-following evolution equations for the four quantities $\phi \in [ \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}, S_d^n/\mathcal {D}, S_d^W,S_d ]$. All these equations have the same general form

(2.22)\begin{align} \frac{\partial^{*} \phi }{\partial^{*}t} &= \underbrace{ \mathbb{A} _0 (\phi) + \mathbb{A} _{1u}(\phi)+\mathbb{A} _2(\phi) }_{\mathbb{A} (\phi)}\nonumber\\ &\quad +\,\underbrace{ \mathbb{B} _{1u}(\phi)+\mathbb{B} _{2u}(\phi) + \mathbb{B} _{3u}(\phi) +\mathbb{B} _4(\phi)+\mathbb{B} _5(\phi)+\mathbb{B} _6(\phi) }_{\mathbb{B} (\phi)} , \end{align}

with

(2.23)\begin{gather} \mathbb{A} _0 (\phi) \equiv -\phi \mathcal{K} , \end{gather}

(2.24)\begin{gather}\mathbb{A} _{1u}( \phi ) \equiv \phi ( \boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{u}} ), \end{gather}

(2.25)\begin{gather}\mathbb{A} _2( \phi ) \equiv -\phi \boldsymbol{\nabla} \boldsymbol{\cdot} ( S_d {\boldsymbol{n}} ). \end{gather}

The complete equations are given by (A13), (A17), (A18) and (A24), respectively. For convenience, those terms are also summarized in table 1. Here $S_{ij}\equiv (\nabla _i u_j + \nabla _j u_i)/2$ is the flow strain rate tensor and the subscript ‘$u$’ in the terms of $\mathbb {A} _{1u}$, $\mathbb {B} _{1u}$, $\mathbb {B} _{2u}$ and $\mathbb {B} _{3u}$ highlights that these terms involve flow velocity. We also want to emphasize that, for the quantities $\boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}$ and $S_d^n/\mathcal {D}$, (2.22) has been derived without further assumptions and is valid for any reaction wave described by (1.1). On the other hand, for (2.22) to be valid for $\phi =S_d^W$, the assumption of $(\partial ^*/\partial ^*t)\mathbb {W}=0$, which is fulfilled when $\mathbb {W}$ depends solely on $c$, is needed. For (2.22) to be valid for $\phi =S_d$, two assumptions of (i) $\rho \mathcal {D}=\textrm {const}$ and (ii) $\mathcal {D}$ depends solely on $c$ are needed.

Table 1. Terms in (2.22) for different choices of $\phi$.

As noted above, under the assumption of $\rho \mathcal {D}=\text {const}$, the decomposition of displacement speed simplifies from (2.12) to $S_d = S_d^n + \mathcal {D} \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}} + S_d^W$. By examining table 1, it can be verified that such a simple decomposition relation also holds for all terms and their parts in (2.22). More specifically, the relations

(2.26)\begin{equation} \left.\begin{gathered} \mathbb{A} _i(S_d) =\mathcal{D}\cdot\mathbb{A} _i(S_d^n/\mathcal{D}) +\mathcal{D}\cdot \mathbb{A} _i( \boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{n}}) + \mathbb{A} _i(S_d^W), \\ \mathbb{B} _i(S_d) =\mathcal{D}\cdot\mathbb{B} _i(S_d^n/\mathcal{D}) + \mathcal{D}\cdot \mathbb{B} _i( \boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{n}}) + \mathbb{B} _i(S_d^W) \end{gathered}\right\}\end{equation}

are satisfied for all terms $\mathbb {A} _i$ and $\mathbb {B} _i$ in table 1.

2.3. Various terms in the transport equations

For each of the four quantities $\phi \in [ \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}, S_d^n/\mathcal {D}, S_d^W,S_d ]$ shown in table 1, each right-hand side term and its sub-parts in (2.22) contribute to the surface-following variation in $\phi$, i.e. $(\partial ^*/\partial ^*t) \phi$. Additional explanations and implications for those sub-terms may be summarized as follows.

1. (a.i) The first part of $\mathbb {A}$, i.e. term $\mathbb {A} _0(\phi )$, defined as a product of the stretch rate and $\phi$, disappears when $(\partial ^*/\partial ^*t)\phi$ is inserted into (2.22) to obtain the surface-averaged equations (2.29).

2. (a.ii) The second term $\mathbb {A} _{1u}(\phi )$ stems from the dilatation $\boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {u}}$. This term with $\phi$ set as $S_d$ may stay non-negative for certain simple reacting flows, e.g. the non-conservative form of continuity equation (3.7) shown later yields $\mathbb {A} _{1u}(S_d)=S_d^2 \lvert \boldsymbol {\nabla } c \rvert / (c+ 1/\theta )\geqslant 0$, with $\theta$ being a positive constant.

3. (a.iii) The third term $\mathbb {A} _2(\phi )$ indicates a correlation between $\phi$ and the divergence of a vector field that contains the displacement speed, i.e. the vector $S_d{\boldsymbol {n}}$.

1. (b.i) The first three parts of $\mathbb {B}$ contain contributions from the flow. The first term $\mathbb {B} _{1u}(S_d) = \mathcal {D} {\boldsymbol {n}}\boldsymbol {\cdot } \nabla ^2 {\boldsymbol {u}}$ represents the normal projection of a flow Laplacian vector and is controlled by the curvature, i.e. $\mathbb {B} _{1u}(S_d)= \mathcal {D} \mathbb {B} _{1u}( \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}})$ and $\mathbb {B} _{1u}(S_d^n/D) )=0$.

2. (b.ii) The second term $\mathbb {B} _{2u}(S_d)=-2\mathcal {D}({\boldsymbol {\nabla }\boldsymbol {\nabla } c}/{ \lvert \boldsymbol {\nabla } c \rvert }):\boldsymbol {\nabla } {\boldsymbol {u}} = -2\mathcal {D}({\nabla _i\nabla _j c}/{ \lvert \boldsymbol {\nabla } c \rvert }) S_{ij}$ is a product of the symmetrical strain rate tensor $S_{ij}= (\nabla _i {\boldsymbol {u}}_j+\nabla _j {\boldsymbol {u}}_i)/2$ and another tensor $(\mathbb {S}_d^D)''_{ij}\equiv \mathcal {D}({\nabla _i\nabla _j c}/{ \lvert \boldsymbol {\nabla } c \rvert })$ extended from the total diffusion displacement speed, i.e. $S_d^D=S_d^c+S_d^n=S_d - S_d^W=\mathcal {D} ({\nabla _i\nabla _i c}/{ \lvert \boldsymbol {\nabla } c \rvert })$, with the latter being the trace of the former, i.e. $S_d^D = (\mathbb {S}_d^D)''_{ii}$. Similarly, term $\mathbb {B} _{2u}( \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}})$ for the curvature part $\nabla _i n_i$ or term $\mathbb {B} _{2u}(S_d^n/\mathcal {D})$ for the normal diffusion part $S_d^n/\mathcal {D}=n_i\nabla _i \ln \lvert \boldsymbol {\nabla } c \rvert$ is a product of $S_{ij}$ and the curvature tensor $\nabla _i n_j$ or the normal diffusion tensor $(\mathbb {S}_d^n)''_{ij} \equiv \mathcal {D} n_i \nabla _j (\ln \lvert \boldsymbol {\nabla } c \rvert )$), respectively.

3. (b.iii) The third term vanishes for the full displacement speed, i.e. $\mathbb {B} _{3u}(S_d)=0$, due to cancellation of two contributions of the opposite signs. For the curvature, $\mathbb {B} _{3u}( \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}} )= - \mathbb {B} _{3u}(S_d^n / \mathcal {D}) = \nabla _n a_n$ represents the normal gradient of the normal strain rate. Note that application of the curvature equation, i.e. (2.22) for $\boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}$, to an one-dimensional (1-D) flame reveals a pair of cancelled, but non-zero terms, i.e. $\mathbb {B} _{3u}( \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}})$ and $\mathbb {B} _{1u}( \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}})$, as shown in figure 5(k) discussed later.

4. (b.iv) The fourth term $\mathbb {B} _4(S_d) = \mathcal {D}\nabla ^2 S_d$ represents diffusion of $S_d$ calculated as a sum of the second derivatives along the three orthogonal directions. Among these three directions, solely the normal direction ${\boldsymbol {n}}$ appears in $\mathbb {B} _4(S_d^n/\mathcal {D})= n_i n_j\nabla _i\nabla _j S_d$ (note that $n_i n_j\nabla _i\nabla _j \psi$ is not equal to $\nabla _n\nabla _n \psi$ for arbitrary $\psi$). The sum of the two directions orthogonal to ${\boldsymbol {n}}$ appears solely in the curvature equation, i.e. $\mathbb {B} _4( \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}) = \nabla ^2 S_d - n_i n_j\nabla _i\nabla _j S_d$.

5. (b.v) The fifth term $\mathbb {B} _5(S_d)$ represents a diffusion-driven ‘convective’ transport of $S_d$ at a generalized velocity based on the normal diffusion contribution to $S_d$. In fact, $\mathbb {B} _5(S_d)=2({\mathbb {S}_d^n})'_i \nabla _i S_d$, where $({\mathbb {S}}_d^n)'_i \equiv \mathcal {D} \nabla _i (\ln \lvert \boldsymbol {\nabla } c \rvert )$. Not surprisingly, $\mathbb {B} _5(S_d^n)=\mathbb {B} _5(S_d)$ and $\mathbb {B} _5( \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}})=\mathbb {B} _5(S_d^W)=0$.

6. (b.vi) The sixth term $\mathbb {B} _6(S_d)$ stems from variations in the normal vector ${\boldsymbol {n}}$ along the normal direction. Its normal diffusion part $\mathbb {B} _6(S_d^n)\equiv S_d(\nabla _i n_j-\nabla _j n_i)\nabla _i n_j$ involves a curl of the normal vector ${\boldsymbol {n}}$, which may not necessarily be rotation free. Nevertheless, as shown later, $\mathbb {B} _6(S_d^n)$ tends to vanish after being averaged over a wave surface.

Another theoretical perspective for understanding the present set of equations is discussed in appendix B, where evolution of an 1-D diffusion-reaction wave in a spherical, cylindrical or planar configuration is studied under assumptions of constant density and zero velocity. In this case, the original wave transport equation (1.1) simplifies to (B1), and the set of equations (2.22) for the four components $S_d$, $S_d^n$, $S_d^c$ and $S_d^W$ reduces to two non-trivial and much simpler equations for $S_d$ and $S_d^n$, see (B12) and (B13), respectively. As later shown in figure 6, the pair of equations for displacement speed and its normal component forms a closed system and can replace the transport equation (B1) for $c$, at least in certain spatial/temporal regions. Interestingly, while (B1) permits different choices of reaction rate function $\mathbb {W}(c)$, this function will solely be translated from the initial conditions when adopting the above replacement, i.e. (B12) and (B13). Furthermore, the simplified 1-D curvature equation (B8) suggests that (i) term $\mathbb {B} _4( \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}})$ is balanced by term $\mathbb {A} ( \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}})$ and (ii) term $\mathbb {B} _6( \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}})$ is balanced by a difference in $(\partial ^*/\partial ^*t) \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}$ and $(\partial /\partial t) \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}$.

2.4. Extension of evolution equations

The decomposition given by (2.22) for the four quantities $\phi \in [ \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}},S_d^n/\mathcal {D},S_d^W, S_d]$ can be extended to functions on $\phi$ and $c$. Consider the surface-following derivative of some function $F(\phi ,c)$,

(2.27)\begin{equation} \frac{\partial^{*} }{\partial^{*}t} F(\phi,c) = \frac{ \partial F}{\partial \phi} \frac{\partial^{*} }{\partial^{*}t} \phi + \frac{ \partial F}{\partial c} \frac{\partial^{*} }{\partial^{*}t} c = \frac{ \partial F}{\partial \phi} \frac{\partial^{*} }{\partial^{*}t} \phi. \end{equation}

Using (2.22), it reads as

(2.28)\begin{equation} \frac{\partial^{*} }{\partial^{*}t}F(\phi,c) = \underbrace{ \frac{ \partial F}{\partial \phi} \mathbb{A} (\phi) }_{\mathbb{A} ( F )}+ \underbrace{ \frac{ \partial F}{\partial \phi} \mathbb{B} (\phi) }_{\mathbb{B} ( F)} \end{equation}

and has a form the same as (2.22). Therefore, the decomposition for $F(\phi ,c)$ is obtained directly from the decomposition for $\phi$ as long as the derivative $\partial F / \partial \phi$ can be evaluated.

Let us consider some functions $F(\phi ,c)$ of physical interest. The ‘absolute’ mean curvature $F (\phi ,c) = | \phi |$ with $\phi = \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}$ is an example of such a quantity. In a commonly adopted configuration of a statistically planar reaction wave, the mean curvature averaged over an entire iso-surface is sufficiently small due to cancellation of contributions from positively and negatively curved surface elements. Conversely, the absolute mean curvature is non-zero as long as some curved elements appear anywhere on the surface. Taking a derivative of the module function gives ${ \partial F}/{\partial \phi }=\sigma _{\pm }(\phi )$ in regions where $|\phi |\neq 0$, and $\sigma _{\pm }(\phi )\equiv \phi /|\phi |$ is the sign of $\phi$. The expanded terms $\mathbb {A} (| \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}|)$ and $\mathbb {B} (| \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}|)$ in (2.22) are also included in table 1, where $\sigma =\sigma _{\pm }( \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}})$ refers to the sign of $\boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}$.

Another function of interest is $F(\phi ,c) = \mathcal {D} \phi$, with the diffusivity $\mathcal {D}(c)$ depending only on $c$. Application of this simple function to $\phi \in [ \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}, | \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}|, S_d^n/\mathcal {D}]$ yields (2.28) with $\mathbb {A} (\mathcal {D} \phi ) = \mathcal {D} \mathbb {A} (\phi )$ and $\mathbb {B} (\mathcal {D} \phi ) = \mathcal {D} \mathbb {B} (\phi )$. Thus, (2.22) holds for the eight quantities $\phi \in [S_d, S_d^W, S_d^n/\mathcal {D}, \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}, | \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}|, S_d^n , \mathcal {D} \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}, \mathcal {D} | \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}| ]$.

2.5. Surface-averaged evolution equations and statistically steady relations

Substitution of (2.22) to (2.11) yields the following evolution equation for surface-averaged quantities:

(2.29)\begin{equation} \frac{\partial }{\partial t} { \left \langle \phi \right \rangle_{{s}} }= \left \langle \mathbb{A} (\phi) \right \rangle_{{s}} + \left \langle \mathbb{B} (\phi) \right \rangle_{{s}} + \underbrace{ \left \langle \phi \mathcal{K} \right \rangle_{{s}} - \left \langle \phi \right \rangle_{{s}} \left \langle \mathcal{K} \right \rangle_{{s}} }_{ \left \langle \mathbb{C}(\phi) \right \rangle_{{s}} }. \end{equation}

(It is possible to define $\mathbb {C}(\phi ) \equiv \phi ' \mathcal {K}'$ with the prime symbol denoting the fluctuation respective to surface average, $\psi '\equiv \psi - \left \langle \psi \right \rangle _{{s}} |_{\hat {c},t}$.) Henceforth, conditioned subscripts in the surface average notion $\left \langle \cdot \right \rangle _{{s}} |_{\hat {c},t}$ are omitted for brevity.

Following Yu et al. (Reference Yu, Nilsson, Brethouwer, Chakraborty and Lipatnikov2020), (2.29) can be reorganized by grouping all sub-terms involving the velocity into a single term $\mathcal {U}$ in order to highlight the net contribution affected directly by the flow. By separating ${\boldsymbol {u}}$ contained in $\mathcal {K}$ from $\left \langle \mathbb {C} \right \rangle _{{s}}$, the separated ${\boldsymbol {u}}$ terms can be cancelled with $\mathbb {A} _0$. Then, (2.29) reads as

(2.30)\begin{align} \frac{\partial }{\partial t} { \left \langle \phi \right \rangle_{{s}} } &= \underbrace{ \left \langle \mathbb{A} _{1u}(\phi) + \mathbb{B} _{1u}(\phi) +\mathbb{B} _{2u}(\phi)+ \mathbb{B} _{3u}(\phi) \right \rangle_{{s}} - \left \langle \phi \right \rangle_{{s}} \left \langle a_t \right \rangle_{{s}} }_{ \mathcal{U} (\phi) } \nonumber\\ &\quad + \underbrace{ \left \langle \mathbb{A} _2(\phi) + \mathbb{B} _4(\phi)+ \mathbb{B} _5(\phi)+ \mathbb{B} _6(\phi) \right \rangle_{{s}} - \left \langle \phi S_d \boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{n}} \right \rangle_{{s}} + \left \langle \phi \right \rangle_{{s}} \left \langle S_d \boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{n}} \right \rangle_{{s}} }_{ \mathcal{R}(\phi) }, \end{align}

where $\mathcal {R}\equiv \left \langle \mathbb {A} \right \rangle _{{s}} + \left \langle \mathbb {B} \right \rangle _{{s}} + \left \langle \mathbb {C} \right \rangle _{{s}} -\mathcal {U}$ denotes the remaining terms that do not contain the velocity ${\boldsymbol {u}}$ implicitly. Under the assumptions of (i) $\rho \mathcal {D}=\textrm {const}$ and (ii) $\mathcal {D}$ and $\mathbb {W}$ depend solely on $c$, the decomposition rules in (2.26) hold for $\left \langle \mathbb {C} \right \rangle _{{s}}$, $\mathcal {U}$ and $\mathcal {R}$ also, i.e.

(2.31)\begin{equation} \left.\begin{gathered} \left \langle \mathbb{C}(S_d) \right \rangle_{{s}} =\mathcal{D}\cdot \left \langle \mathbb{C}(S_d^n/\mathcal{D}) \right \rangle_{{s}} + \mathcal{D}\cdot \left \langle \mathbb{C}( \boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{n}}) \right \rangle_{{s}} + \big \langle \mathbb{C}(S_d^W) \big \rangle_{{s}} , \\ \mathcal{U}(S_d) =\mathcal{D}\cdot\mathcal{U}(S_d^n/\mathcal{D}) + \mathcal{D}\cdot\mathcal{U}( \boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{n}}) + \mathcal{U}(S_d^W), \\ \mathcal{R}(S_d) =\mathcal{D}\cdot\mathcal{R}(S_d^n/\mathcal{D}) +\mathcal{D}\cdot\mathcal{R}( \boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{n}}) + \mathcal{R}(S_d^W). \end{gathered}\right\}\end{equation}

Finally, by considering a fully developed, statistically stationary wave ($t=\infty$), (2.29) and (2.30) yield the constraints

(2.32)\begin{equation} \left.\begin{gathered} \left \langle \mathbb{A} (\phi) \right \rangle_{{s}} |_{\hat{c},t_\infty}+ \left \langle \mathbb{B} (\phi) \right \rangle_{{s}} |_{\hat{c},t_\infty} + \left \langle \mathbb{C}(\phi) \right \rangle_{{s}} |_{\hat{c},t_\infty} =0, \\ \mathcal{U}(\phi)|_{\hat{c},t_\infty} + \mathcal{R}(\phi)|_{\hat{c},t_\infty} =0, \end{gathered}\right\}\end{equation}

for the fully developed reaction waves. These constraints hold for all $\hat {c}\in (0,1)$ and for $\phi \in [S_d,S_d^W, S_d^n/\mathcal {D}, \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}, | \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}|, S_d^n, S_d^c, |S_d^c|]$. Moreover, since the fully developed surface-averaged stretch rate vanishes due to (2.9), one more constraint of

(2.33)\begin{equation} \left \langle \mathbb{C}(\phi) \right \rangle_{{s}} |_{\hat{c},t_\infty}= - \left \langle \mathbb{A} _0(\phi) \right \rangle_{{s}} |_{\hat{c},t_\infty} \end{equation}

holds.

3.1. Direct numerical simulation database

Since the DNS attributes are discussed in detail elsewhere (Yu et al. Reference Yu, Lipatnikov and Bai2014, Reference Yu, Bai and Lipatnikov2015b; Elperin et al. Reference Elperin, Kleeorin, Liberman, Lipatnikov, Rogachevskii and Yu2016; Yu & Lipatnikov Reference Yu and Lipatnikov2017a,Reference Yu and Lipatnikovb, Reference Yu and Lipatnikov2019a), we will restrict ourselves to a brief summary of the simulations. A unity Lewis number premixed flame is modelled by (1.1) with $\rho \mathcal {D}=\textrm {const}$ and a single-step reaction rate term given by

(3.1)\begin{equation} \mathbb{W} =\frac{1}{(1+\theta)^s}\frac{1-c}{\tau_R} \exp \left[-\frac{Ze (1+\theta)^2}{\theta (1+\theta c)} \right], \end{equation}

where $\tau _R$ is a constant reaction time scale and $Ze=6$ in order for the rate $\mathbb {W}$ to depend on $c$ in a highly nonlinear manner. The density depends on $c$ as

(3.2)\begin{equation} \rho = \frac{ \rho_u }{1+ (1-s) \theta c }, \end{equation}

where two constants $\rho _u$ and $\rho _b$ denote density of unburned and burned gas, respectively. A constant-density reaction wave (i.e. $\rho =\rho _u=\rho _b$) is modelled by setting $s=1$ and $\theta =6$. A flame with heat release and, consequently, $\rho _u/\rho _b > 1$ is modelled by setting $s=0$ and $\theta =\rho _u/\rho _b-1$.

The single-reaction wave propagates in a forced turbulence described by the continuity

(3.3)\begin{equation} \frac{\partial \rho }{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} (\rho {\boldsymbol{u}}) =0, \end{equation}

and the Navier–Stokes equations

(3.4)\begin{equation} \frac{\partial \rho {\boldsymbol{u}} }{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} (\rho {\boldsymbol{u}} {\boldsymbol{u}}) = -\boldsymbol{\nabla} p + \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{ \tau} + \boldsymbol{f}, \end{equation}

where $p$ is the pressure, $\tau _{ij} \equiv \mu (\nabla _j u_i +\nabla _i u_j - \frac {2}{3} \delta _{ij} \nabla _k u_k )$ is the viscous stress tensor with a constant dynamic viscosity $\mu$. A function $\boldsymbol {f}\equiv ( \rho / \rho _u ) \boldsymbol {f}_u$ is used to maintain turbulence intensity by applying energy forcing at low wavenumbers, with $\boldsymbol {f}_u$ being simulated using the method developed by Lamorgese, Caughey & Pope (Reference Lamorgese, Caughey and Pope2005) in the case of a constant density. The reader interested in details regarding the characteristics of the forced turbulence and the calculation of the forcing term in wavenumber domain is referred to section 3 in our recent paper (Yu & Lipatnikov Reference Yu and Lipatnikov2017b). When simulating a variable-density reaction wave, the factor $(\rho /\rho _u)$ in $\boldsymbol {f}$ reduces the energy injection on the burned side, thus mimicking the turbulence decay due to dilatation and an increase in the kinematic viscosity ($\nu \equiv \mu /\rho$) in the wave.

A reaction wave evolves in a rectangular box with size $\varLambda _x \times \varLambda \times \varLambda$, represented using a uniform grid of $N_x \times N \times N$ cubic cells. The boundary conditions are periodic in the transverse directions $y$ and $z$ for all simulations. For simulation of a constant-density wave, periodic conditions are also applied in direction $x$ normal to the mean wave surface, because the flow is not affected by the wave propagation in this case. In other words, when a constant-density reaction wave reaches the left boundary ($x=0$) of the computational domain, an identical reaction wave enters the domain through its right boundary ($x=\varLambda _x$). Such a method greatly improves sampling of statistics by simulating many cycles of the wave propagation through the computational domain, but this method may only be used in the case of $\rho =\textrm {const}$ and $\nu =\textrm {const}$, provided that the mean wave brush thickness is smaller than the length of the computational domain. These constraints are satisfied in our simulation of a constant-density wave.

For simulation of a variable-density reaction wave (flame), the inlet and outlet conditions are prescribed at $x=0$ and $x=\varLambda _x$, respectively. Unburned gas is fed into the inlet at $x=0$ with an ‘almost’ periodic copy of the downstream velocity at $x=\varLambda$. More specifically, $u(0,y,z,t) = u(\varLambda ,y,z,t)- ({1}/{\varLambda ^2})\iint _{0,0}^{\varLambda ,\varLambda } u(\varLambda ,y,z,t) \,\mathrm {d} y \,\mathrm {d} z$, $v(0,y,z,t) = v(\varLambda ,y,z,t)$, $w(0,y,z,t) =w(\varLambda ,y,z,t)$ and $c(0,y,z,t) =0$. During the entire simulation, the values of $\rho$ and $c$ remain close to $\rho _u$ and zero, respectively, at $0<x<\varLambda$. Therefore, the near-inlet cubic-box region (i.e. $\forall x,y,z \in [0,\varLambda ]$) resembles an isolated simulation of periodic constant-density turbulence.

During the simulation of the flame C, a control strategy is used to keep the flame inside the computational domain. The strategy, which was adopted in a recent DNS study of intrinsic flame instability (Yu, Bai & Bychkov Reference Yu, Bai and Bychkov2015a), aims at keeping the $x$-position of the mean flame brush, i.e. $X_m(t) \equiv ({1}/{\varLambda ^2})\iiint _{0,0,0}^{\varLambda _x,\varLambda ,\varLambda } (1-c(x,y,z,t) )\,\mathrm {d} x\,\mathrm {d} y \,\mathrm {d} z$, close to a prescribed target position of $X_{T}=0.6\varLambda _x$. For this purpose, (1.1), (3.3) and (3.4) are solved in a coordinate framework that moves at an adjustable axial velocity $U_{m}(t)$. This velocity is controlled in a so-called proportional-integral-derivative fashion, i.e. $U_{m}(t+\Delta t) = U_{m}(t ) + (\mathrm {d}_P + \mathrm {d}_I + \mathrm {d}_D)$, where $\mathrm {d}_P= e/(2\Delta t)$, $\mathrm {d}_I= \int _0^t e(t) \,\textrm {d}t/(10 \Delta t)$, $\mathrm {d}_D = [e(t) - e(t-\Delta t)]/(10 \Delta t)$, $e(t) = X_m(t)- X_{T}$ is the target error and $\Delta t$ is the computational time step. When changing the coordinate framework from $(t,x,y,z)$ to $(t', x', y', z')$, where $t'=t$, $x'=x + \int _0^t U_m \,\mathrm {d} t$, $y'=y$ and $z'=z$, the unsteady terms in (1.1), (3.3) and (3.4) are appropriately changed, i.e. $\partial \phi /\partial t = \partial \phi /\partial t' + U_m \partial \phi /\partial x'$ for any relevant $\phi$. Other terms in these governing equations remain unchanged with the exception of replacement of the coordinates $(t, x, y , z)$ by $(t',x',y',z')$. Accordingly, the same method was adopted to numerically solve the governing equations in both coordinate frameworks.

An initial divergence-free turbulence field is generated by synthesizing Fourier waves (Yu & Bai Reference Yu and Bai2014) with an initial r.m.s. velocity $u_0$ and the integral length scale $\ell _0=\varLambda /4$. The initial turbulent Reynolds number $Re_0 = u_0 \ell _0 \rho _u/\mu$ can be changed by changing the domain width $\varLambda$. In the simulations of both constant-density reaction wave B and flame C, a non-decaying turbulent field is obtained by integrating (3.4) to maintain a r.m.s. velocity $u'(t) \approx u_0$. Here, $u'(t)^2 \equiv ({2}/{3 \varLambda ' \varLambda ^2})\iiint _{0,0,0}^{\varLambda ',\varLambda ,\varLambda } k(x,y,z,t) \,\mathrm {d} x \,\mathrm {d} y\,\mathrm {d} z$, $k\equiv (u^2+v^2+w^2)/2$ is the turbulent kinetic energy, $\varLambda ' = \varLambda _x$ if $\rho _u/\rho _b=1$ and $\varLambda '= \varLambda$ when $\rho _u/\rho _b >1$. Different longitudinal integral length scales $L_{11}$ can be generated (Yu & Lipatnikov Reference Yu and Lipatnikov2017b) by appropriately adjusting $\boldsymbol {f}_u$ (and, hence, $\boldsymbol {f}$) using the method by Lamorgese et al. (Reference Lamorgese, Caughey and Pope2005).

3.2. Numerical methods

The governing equations are numerically integrated based on an in-house DNS solver (Yu, Yu & Bai Reference Yu, Yu and Bai2012) developed for low Mach number reacting flows. To treat the stiffness in the calculation of reaction rates in the case of detailed chemistry, the temporal integration of the CDR equation in the original solver is based on a second-order symmetrical Strang splitting method (Strang Reference Strang1968) by placing a full-time step stiff integration of a chemistry calculation in between two half-time step integrations of convection and diffusion terms, while the full-time step $\Delta t$ is determined by the CFL condition. In the present work, which deals with single-step chemistry, the operator splitting strategy is not adopted, because the reaction rate term ($\mathbb {W}$) is non-stiff. Accordingly, temporal integration of the whole CDR equation is performed explicitly over a full-time step. An explicit advancement is proceeded in multiple ($K>1$) sub-time steps of smaller size ($\Delta t/K$) determined by the diffusion-stability-limit, starting from a Runge–Kutta step and followed by Adam–Bashforth steps. More specifically, considering advancement of (1.1) from $t^n=n \Delta t$ to $t^{n+1}=(n+1) \Delta t$, let $\mathbb {E} =-{\boldsymbol {u}}\boldsymbol {\cdot } \boldsymbol {\nabla } c + \mathbb {W}$ denote the combined convection and reaction term. A two-stage Runge–Kutta step is first performed for ($k'=0$, $\alpha =\frac {1}{2}$, $\beta =0$) and, then, for ($k'=1$, $\alpha =1$, $\beta =-\frac {1}{2}$), in order to advance the solution to a first time instant of $t^n+{ \Delta t}/{K}$,

(3.5)\begin{align} \frac{ c^{(({k'+1})/{2K})+n}-c^{({k'}/{2K})+n}}{\Delta t/K} &= \frac{1}{2}\left(\frac{3}{2} \mathbb{E}(c^n,u^{n}) - \frac{1}{2} \mathbb{E}(c^{n-1},u^{n-1}) \right) \nonumber\\ &\quad + (\alpha \mathbb{D}(c^{ ({k'}/{2K})+n} ) + \beta \mathbb{D}(c^{ (({k'-1})/{2K})+n} ) ) , \end{align}

where any superscript multiplied with $\Delta t$ refers to the time instant, followed by $K-1$ sequential Adam–Bashforth steps for $k\in [1, 2,\ldots , K-1]$ to reach the final step at $t^{n+1}$,

(3.6)\begin{align} \frac{ c^{(({k+1})/{K})+n}-c^{({k}/{K})+n} }{\Delta t/K} &= \left(\frac{3}{2} \mathbb{E}(c^n,u^n) - \frac{1}{2} \mathbb{E}(c^{n-1},u^{n-1}) \right) \nonumber\\ &\quad + \left( \frac{3}{2} \mathbb{D} (c^{ ({k}/{K})+n} ) - \frac{1}{2} \mathbb{D}(c^{( ({k-1})/{K})+n} ) \right) . \end{align}

As described by Yu et al. (Reference Yu, Yu and Bai2012), to avoid numerical instability associated with simulation of variable-density flows, the continuity equation (3.3) is rewritten using (1.1) and (3.2) into a non-conservative form as

(3.7)\begin{align} \boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{u}} = - \frac{1}{\rho} \left( \frac{\partial \rho}{\partial t} + {\boldsymbol{u}} \boldsymbol{\cdot} \boldsymbol{\nabla} \rho \right) = - \frac{1}{\rho} \frac{ \partial \rho}{\partial c} \left( \frac{\partial c}{\partial t} + {\boldsymbol{u}} \boldsymbol{\cdot} \boldsymbol{\nabla} c \right)= \begin{cases} 0 & \text{if $\rho_u/\rho_b=$1},\\ \dfrac{ \mathbb{D}+\mathbb{ W} }{ c + 1/\theta } & \text{if $\rho_u/\rho_b >$ 1.} \end{cases} \end{align}

Discretization of (3.7) yields a constraint for the flow velocity at $t^{n+1}$. Thereafter, a standard fractional-step approach for pressure-velocity decoupling results in a variable-coefficient Poisson equation for the pressure, which is solved using Gauss-Seidel iterations with multi-grid acceleration (Yu & Bai Reference Yu and Bai2013). For spatial discretization, a sixth-order central difference scheme is used for all terms containing a spatial derivative with the exception of the convection term in (1.1), which is discretized with a fifth-order WENO scheme (Jiang & Shu Reference Jiang and Shu1996) to avoid numerical overshooting.

3.3. Simulation set-up

In the constant-density case, both developing and fully developed reaction waves are simulated starting from the pre-computed laminar wave profile of $c_L(\xi )$ with $dc_L/d\xi >0$. In order to study the fully developed turbulent reaction wave, a planar wave $c^s(x,0)=c_L(\xi )$ is initially ($t=0$) released at $x_0=\varLambda _x/2$ such that $\int _{-\infty }^{0} c_L(\xi ) \,\textrm {d}\xi = \int _{0}^{\infty }[1-c_L(\xi )]\,\textrm {d}\xi$ and $\xi =x-x_0$. Subsequent evolution of the field $c^s(\boldsymbol {x},t)$ is simulated by solving (1.1). Computation of fully developed statistics with sampling every 100 time steps $\Delta t$ is started after the forced turbulence has reached a statistically stationary state ($t=t^*>3.5\tau _t^0$) and is performed over a time interval longer than 50 $\tau _t^0$. In order to study transient turbulent reaction waves, several copies of the same pre-computed laminar wave profile $c_L (\xi )$ are simultaneously embedded into the turbulent flow in $\mathcal {M}$ equidistantly separated planar zones centred around $x_m/\varLambda _x =(m-0.5)/\mathcal {M}$, i.e. $c^t_m(x,t^*)=c_L(\xi _m)$, where coordinate $\xi _m=x-x_m$ is set using $\int _{-\infty }^{0}c_L(\xi _m)\,\textrm {d}\xi _m = \int _{0}^{\infty }[1-c_L(\xi _m)]\,\textrm {d}\xi _m$ and $m$ is an integer number (1$\leqslant m \leqslant \mathcal {M}=$11). Subsequently, evolutions of $\mathcal {M}$ non-interfering transient fields $c_m^t (\boldsymbol {x},t)$ are simulated by solving $\mathcal {M}$ independent (1.1). The transient simulations are run over 2$\tau _t^0$ before being reset. Subsequently, at $t=t^*+2j\tau _t^0$ with $1\leqslant j \leqslant J$, the flow is again populated by $\mathcal {M}$ new profiles of $c_L (\xi _m )$, the transient simulations are repeated $J$ times, giving an ensemble of $J\times \mathcal {M}$ independent transient waves. For comparison, the fully developed statistics computed using a single $c^s$-field are associated with $2 J\tau ^0_t/(100\Delta t)$ realizations.

Contrary to the constant-density case, the present flame simulations are restricted to a single fully developed turbulent wave $c^s$. The simulations are started by embedding pre-computed laminar flame profiles of $u_L(\xi )$ and $c_L(\xi )$, with $c_L(-\infty )=u_L(-\infty )=0$, into a synthesized turbulence field at $x_0=X_T$. The initial velocity $U_m(0)$ of the moving coordinate framework is equal to the laminar flame speed $S_L$. Subsequently, (1.1), (3.3) and (3.4) are numerically solved using the inlet and outlet boundary conditions, as well as the control strategy for adjusting $U_m(t)$. Similar to the constant-density case, the fully developed statistics are sampled after reaching a statistically stationary numerical solution ($t>t^*$), with the sampling being performed over at least 50 $\tau _t^0$.

The DNS cases are set up combining one of the forced turbulence fields with a reaction wave characterized by a laminar wave speed $S_L$ and thickness $\delta _F=\mathcal {D}/S_L$. The required reaction time scale $\tau _R$ in (3.1) is found through 1-D pre-computations of the laminar wave.

In the present paper the data analysis is mainly focused on results obtained in two representative turbulent and well-resolved cases (constant-density reaction wave B and flame C). The major characteristics are reported in table 2, where $Da=\tau _t/\tau _F$ is the Damköhler number; $Ka=\tau _F/\tau _{\eta }$ is the Karlovitz number; $Pe=u'L_{11}/\mathcal {D}$ is the turbulent Péclet number; $\tau _F=\delta _F/S_L$ is the wave time scale; $\tau _t = L_{11}/u'$ and $\tau _{\eta }=(\nu _u/\bar {\varepsilon }^u)^{1/2}$ are integral and Kolmogorov time scales of the turbulence, respectively; $\bar {\varepsilon }^u = 2 \nu _u \overline {S_{ij}S_{ij} }^u$ is the dissipation rate averaged (i.e. $\overline {\;\cdot \;}^u$) over either the whole domain in the constant-density case or the near-inlet cubic-box region in the variable-density case, and over at least $50\tau _t^0$ after turbulence and the $c^s$-wave has become statistically stationary; and a ratio of $\delta _F/\Delta x$ characterizes the grid resolution in terms of the number of grid points per laminar wave thickness. Moreover, in cases B and C, $L_{11}/\varLambda =0.11$, $\tau _t^0/\tau _t=2.3$ and $\eta /\Delta x=1.1$. Here, $\eta = [\nu _u^3/\bar {\varepsilon }^u]^{1/4}$ is the Kolmogorov length scale. Snapshots for cases B and C are shown in figures 1 and 2, respectively.

Figure 1. (a) Three simultaneously embedded, non-interfering, constant-density developing ($t=0.9\tau _t^0$) reaction waves $c^t_{m}({\boldsymbol {x}},t)$ with $m=1$, 4 and 7 (colour iso-surfaces and fields) and the fully developed reaction wave $c^s({\boldsymbol {x}},t)$ (grey iso-surface and field). Both the peak-reaction-rate iso-surface ($\hat {c}=0.88$) and two-dimensional image of the $c$-field on the $xz$ (the transient waves) or $xy$ (the fully developed wave) boundary surface are shown for each wave. (b) Three-dimensional turbulence field visualized using Lambda2-vortex method using iso-contours of the intermediate eigenvalue of the strain rate tensor (Jeong & Hussain Reference Jeong and Hussain1995) in case B.

Figure 2. Two instantaneous snapshots (two-dimensional cut at the $z$-middle) demonstrate the evolution of a flame in a forced homogeneous turbulence in case C, see table 2. The $z$-vorticity is visualized using different colours. Black iso-lines are shown for nine $c$-levels increasing from left to right. The significant drop in the $z$-vorticity magnitude downstream of the reaction zone reflects turbulence inhibition due to thermal expansion.

Table 2. Three representative DNS cases.

During evolution of a turbulent reaction wave, zero-gradient points characterized by $\lvert \boldsymbol {\nabla } c \rvert ({\boldsymbol {x}},t)=0$ (Gibson Reference Gibson1968) may appear. While the zero-gradient points are excluded from the definition of the surface average by (2.3), see Yu & Lipatnikov (Reference Yu and Lipatnikov2019b,Reference Yu and Lipatnikovc), certain quantities of interest, e.g. $S_d$, $\boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}$ and $\mathcal {K}$, may locally grow unboundedly in the neighbourhood of a zero-gradient point, which pose a challenge for the numerical calculation of surface-averaged quantities. (Note, a surface-averaged value may still remain bounded even when the local value becomes unbounded, see discussion in appendix B of Yu & Lipatnikov (Reference Yu and Lipatnikov2019b).) To explore the eventual influence of such points and their neighbourhoods on the accuracy of evaluation of various terms in the evolution equations, an additional case A is designed. Case A is largely identical to case B, but the turbulent field is replaced with a frozen shear flow, i.e. $u(x,y,z,t)=-u_0 \cos (2{\rm \pi} y/\varLambda )$, $v=w=0$, and the momentum equation (3.4) is not solved. Characteristics of case A, reported in table 2, are calculated using $L_{11}=\varLambda /2$, $u'=u_0$ and $t_\eta = t_\tau =L_{11}/u_0$. Simulation of multiple transient waves $c^t$ in case A is performed largely similarly to case B (some differences between cases A and B are discussed by Yu & Lipatnikov Reference Yu and Lipatnikov2019c) but the duration of transient sampling is changed from $2\tau _t^0$ to $2\tau _F$. As shown in supplemental figure S8, in case A iso-surfaces are only bent, but there is no zero-gradient point in the computational domain. Comparison of results computed in cases A (no zero-gradient points) and B (no restriction on the appearance of zero-gradient points) offers an opportunity to estimate the influence of such points on the evaluation of various surface-averaged terms. Moreover, simulations of developing waves in case B also offer such an opportunity. Indeed, since the transient waves $c_m^t$ begin their evolution from a regular flat initial surface (there is no zero-gradient point at $t=0$), monitoring evolution of (i) iso-surfaces of the transient fields $c_m^t(\boldsymbol {x},t)=\hat {c}$ and (ii) the relevant surface-averaged quantities allows us to detect any anomaly in the developing surface-averaged terms and to see the eventual influence of the zero-gradient points on these terms. Note that transient data can also be of interest in themselves because the vast majority of premixed turbulent flames are developing flames (Lipatnikov & Chomiak Reference Lipatnikov and Chomiak2002; Lipatnikov Reference Lipatnikov2012).