Evolution equations for the decomposed components of displacement speed in a reactive scalar field

Abstract

The study of a turbulent premixed flame often involves analysing quantities conditioned to different iso-surfaces of a reactive scalar field. Under the influence of turbulence, such a surface is deformed and translated. To track the surface motion, the displacement speed ($S_d$) of the scalar field respective to the local flow velocity is widely used and this quantity is currently receiving growing attention. Inspired by the apparent benefits from a simple decomposition of $S_d$ into contributions due to (i) curvature, (ii) normal diffusion and (iii) chemical reaction, this work aims at deriving and exploring new evolution equations for these three contributions averaged over the reaction surface. Together with a previously obtained $S_d$-evolution equation, the three new equations are presented in a form that emphasizes the decomposition of $S_d$ into three terms. This set of equations is also supplemented with a curvature-evolution equation, hence providing a new perspective to link the flame topology and its propagation characteristics. Using two direct numerical simulation databases obtained from constant-density and variable-density reaction waves, all the derived equations and the term-wise decomposition relations are demonstrated to hold numerically. Comparison of the simulated results indicates that the thermal expansion weakly affects the key terms in the considered evolution equations. Thermal expansion can cause variations in the averaged $S_d$ and its decomposed parts through multiple routes more than introducing a dilatation term. The flow plays a major role to influence the key terms in all equations except the curvature one, due to a cancellation between negatively and positively curved surface elements.

1. Introduction

A turbulent premixed flame is often modelled using a reaction progress variable $c({\boldsymbol {x}},t)$. The reaction progress variable is equal to zero and unity on the unburned reactant and burned product sides of a flame, respectively. The evolution of $c$ is described by the convection-diffusion-reaction (CDR) equation

(1.1)\begin{equation} \frac{ \partial c}{ \partial t} + {\boldsymbol{u}} \boldsymbol{\cdot} \boldsymbol{\nabla} c = \underbrace{ \frac{1}{\rho} \boldsymbol{\nabla} \boldsymbol{\cdot}( \rho \mathcal{D} \boldsymbol{\nabla} c ) }_\mathbb{D} + \mathbb{W}, \end{equation}

where $t$, ${\boldsymbol {u}}$, $\rho$, $\mathcal {D}$, $\mathbb {D}$ and $\mathbb {W}$ denote time, velocity, density, diffusivity, the diffusion term and reaction rate, respectively. The term $\mathbb {W}$ often has a highly nonlinear functional dependence on $c$. Consequently, when studying turbulent flames, many quantities of interest exhibit large spatial and temporal fluctuations and tend to possess a strong dependence on $c$. Therefore, both when analysing premixed flames from simulations and experiments, and when formulating models for turbulent premixed flames for use in Reynolds averaged Navier–Stokes or large eddy simulation, it is often useful to consider certain statistics conditioned on $c$. In particular, a quantity $\phi ({\boldsymbol {x}},t)$ averaged over an iso-scalar surface of $c({\boldsymbol {x}},t)=\hat {c}$ at a given $\hat {c}\in (0,1)$ (Veynante & Vervisch 2002) received plenty of attention in the literature.

Consideration of conditional statistics brings forward an analysis perspective of tracing an iso-$c$ surface and studying its evolution. An iso-$c$ surface is formed by the union of points sharing an identical value of $c$, and (1.1) contains all necessary information about how any particular point on an iso-$c$ surface moves with the flow. To see this, we first introduce a displacement speed $S_d$ defined at any point $({\boldsymbol {x}},t)$ for which $0<c({\boldsymbol {x}},t)<1$ (Veynante & Vervisch 2002),

(1.2)\begin{equation} S_d\equiv (\mathbb{D}+\mathbb{W} )/ \lvert \boldsymbol{\nabla} c \rvert . \end{equation}

Then, let us consider a composite ‘total’ velocity ${\boldsymbol {u}}^*\equiv {\boldsymbol {u}}- S_d {\boldsymbol {n}}$ where the local iso-surface normal direction is ${\boldsymbol {n}}\equiv \boldsymbol {\nabla } c/ \lvert \boldsymbol {\nabla } c \rvert$. Finally, a surface-following time derivative operator ${\partial ^{*}}/{\partial ^{*}t}$ is defined as

(1.3)\begin{equation} \frac{\partial^{*} }{\partial^{*}t} \phi \equiv \frac{\partial }{\partial t} \phi+ {\boldsymbol{u}}^*\boldsymbol{\cdot} \boldsymbol{\nabla} \phi \end{equation}

for any quantity $\phi ({\boldsymbol {x}},t)$. Using this operator, (1.1) is reduced to $(\partial ^*/\partial ^*t) c=0$, reflecting the fact that $c$ is constant for any point following the iso-surface. The displacement speed $S_d$ represents the self-propagation speed of the scalar field $c$ relative to the local fluid (i.e. excluding flow convection) due to a combined effect of diffusion and reaction. Thus, $S_d$ is an essential part of the total velocity ${\boldsymbol {u}}^*$ and key information for tracing the movement of the iso-surface points.

One reason to study the displacement speed stems from the fact that it directly affects stretch rate, which in turn controls the rate of change of the local area of an iso-scalar surface within a turbulent reaction wave (Candel & Poinsot 1990). Accordingly, research into the behaviour of $S_d$ in a turbulent reacting flow is of vital importance for understanding the influence of turbulence on the bulk reaction rate, which is proportional to the area of the reaction zone surface. Moreover, studying statistical information about the displacement speed can benefit premixed combustion modelling based on, for example, the level-set (G-equation) formulation (Williams 1985) or the flame surface density approach for reaction rate closure (Candel & Poinsot 1990). Furthermore, variations in the displacement speed across different zones inside a flame can provide additional information about changes in the intrinsic flame structure. Such changes can be induced by turbulence and knowledge about the connection between turbulence and the displacement speed can be particularly useful when studying turbulent flames.

Statistical behaviour of the displacement speed has been the topic of a number of studies based on direct numerical simulations of various turbulent reacting flows, e.g. see Echekki & Chen 1996; Gran, Echekki & Chen 1996; Chen & Im 1998; Peters et al. 1998; Echekki & Chen 1999; Chakraborty & Cant 2005; Dopazo, Martín & Hierro 2007; Wang, Hawkes & Chen 2017a; Luca et al. 2019 and references therein. In particular, there are many studies of correlations between displacement speed and other important quantities characterizing the flame surface or the local turbulent flow. A list of such quantities includes, but is not limited to (i) the local flame curvature (Echekki & Chen 1996; Chakraborty 2007; Sankaran et al. 2015; Wang et al. 2017a; Luca et al. 2019), (ii) strain rate (Chakraborty & Cant 2004; Hawkes & Chen 2006; Kim & Pitsch 2007; Chaudhuri 2015; Cecere et al. 2016), (iii) flame surface topology (Dopazo et al. 2007; Cifuentes et al. 2014), and (iv) alignment characteristics of scalar gradient with principle strain rates (Chakraborty & Swaminathan 2007).

Another research direction consists in developing simple model equations for evaluation of $S_d$ in turbulent reacting flows, e.g. see the review paper by Lipatnikov & Chomiak (2005). Such equations are typically based on the theory of weakly perturbed laminar flames (Matalon & Matkowsky 1982; Pelcé & Clavin 1982; Class, Matkowsky & Klimenko 2003; Kelley, Bechtold & Law 2012), but fail in predicting negative local displacement speeds (Gran et al. 1996; Peters et al. 1998; Chakraborty & Cant 2004; Chakraborty 2007; Wang et al. 2017a; Dave & Chaudhuri 2020) or high local values of $S_d/S_L=O(10)$ documented in recent direct numerical simulation (DNS) studies of premixed turbulent flames (Chaudhuri 2015; Lipatnikov et al. 2015; Im et al. 2016; Uranakara et al. 2016; Lipatnikov et al. 2018). Here, $S_L$ is the laminar flame speed, with $S_d=S_L$ in the unperturbed laminar flame. Recently, significant progress in modelling the high local values of $S_d/S_L$ was made by Dave & Chaudhuri (2020) who explored the evolution of local flame displacement speeds by analysing DNS data and adopting flame tracking techniques developed earlier by Chaudhuri (2015) and Dave, Mohan & Chaudhuri (2018).

It has also been recognized that the displacement speed can naturally be decomposed into three components $S_d^W$, $S_d^n$ and $S_d^c$ (Chen & Im 1998; Peters et al. 1998; Echekki & Chen 1999), which represent reaction, diffusion in the normal direction and tangential diffusion induced due to the surface curvature, respectively. The third component $S_d^c$ links self-propagation characteristics of an iso-surface with its geometry. Such analyses have revealed some physical insights. For instance, Peters (1999) argued that $S_d^W$ and $S_d^n$ dominate in weak turbulence, whereas $S_d^c$ dominates in intense turbulence. Based on this idea, Peters (1999) arrived at an expression for the turbulent flame speed, which is applicable to different combustion regimes. This and other works motivated investigating statistics of $S_d^W$, $S_d^n$ and $S_d^c$, which was addressed in a number of DNS studies (Chakraborty & Cant 2004; Chakraborty 2007; Wang et al. 2017b).

While statistics of the local values of $S_d$, $S_d^W$, $S_d^n$ and $S_d^c$ are addressed in many papers, a transport equation for displacement speed has received little attention yet. Bearing in mind the direct link between $S_d$ and iso-$c$ surface, such an equation should be analysed in the iso-surface-following form, i.e. $(\partial ^*/\partial ^*t)\phi = \sum \text {rhs}$, where rhs abbreviates all right-hand side terms. Evolution equations written in this form were discussed for some quantities such as the absolute gradient $\lvert \boldsymbol {\nabla } c \rvert$ of $c$, which characterizes a separation distance between iso-surfaces (the evolution equation for $\lvert \boldsymbol {\nabla } c \rvert$ is well known and its derivation is given in appendix A, see (A28), for completeness), curvature $\phi = \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}$ (Cifuentes et al. 2018; Dopazo et al. 2018), etc. However, to the best of the authors’ knowledge, such an evolution equation for displacement speed has not been explored until recently. To fill this knowledge gap and to provide a more direct framework for studying displacement speed, we recently derived a surface-following evolution equation for $S_d$ (Yu & Lipatnikov 2019c).

In the same paper we also addressed another problem relevant to previous studies. Given an evolution equation for $\phi$ written in surface-following form, how should a conditioned surface average be applied to ensure that the average of the left-hand side (lhs) terms matches the average of the right-hand side terms? Noticing the key hurdle being the surface average does not commute with the time derivative when the surface itself is changing in time, we derived a general relation, see (2.11) in the next section, to convert the evolution equation for an arbitrary quantity from a non-averaged ‘local’ form to a ‘surface-averaged’ form, thus providing a quantitative framework for numerical study of surface averages in turbulent reacting flows. The above relation was used to generate evolution equations for the surface-averaged $S_d$ and $\lvert \boldsymbol {\nabla } c \rvert$ (Yu et al. 2019). Then, using a DNS database obtained from constant-density turbulent reaction waves, we showed that all terms in the surface-averaged $S_d$-equation could be numerically evaluated with satisfactory accuracy, with a reasonable match between the averaged left-hand and right-hand side terms being achieved. Further numerical investigation (Yu et al. 2020) of the evolution equation for the surface-averaged $S_d$ (for brevity, we will write $S_d$-equation in the following) provided insights into issues such as a retreating iso-surface element moving at a negative displacement speed and relation between rapidly propagating surface elements and curvature. By simply grouping all terms involving ${\boldsymbol {u}}$, the direct effect of the flow velocity on the evolution of the surface-averaged $S_d$ was also explored.

Nevertheless, there is still a number of questions related to the $S_d$-equation that deserve further investigation. First, several terms in the $S_d$-equation contain high-order derivatives and take a relatively complicated form, impeding a better interpretation and understanding of their physical meaning. Second, in a realistic turbulent combustion system, multiple phenomena can affect the self-propagation characteristics of a reaction wave. In order to isolate and study the effect of flow turbulence, our previous numerical assessment on the $S_d$-equation (Yu et al. 2020) was limited to DNS data obtained from a constant-density reaction wave that passively propagated through a background turbulent flow. A simple examination of the derived equation terms (Yu et al. 2020) has identified multiple routes through which the flow dilatation can affect the evolution of $S_d$. However, it is still unknown whether or not such effects are significant enough to change the main trends observed in the constant-density case.

In the present work, inspired by the physical insights Peters (1999) obtained using the simple decomposition of the displacement speed, we derive a set of new surface-averaged evolution equations for each of three decomposed parts of $S_d$. It is found that the corresponding terms in the $S_d$-equation and the equations for its decomposed parts can be organized in a consistent way such that they satisfy the same decomposition rule, thus providing a convenient framework to understand the evolution of $S_d$ and its components. Specifically, and not surprisingly, the equation for the curvature part of $S_d$ is equivalent to the equation for mean curvature. To the best of the authors’ knowledge, the latter equation has not yet been directly linked to the $S_d$-equation. Moreover, the three decomposed parts of $S_d$ provide three new surface-averaged constraint relations that must be satisfied in the statistically fully developed flame. These relations can potentially benefit turbulent combustion modelling. Finally, to answer the questions pertaining to the flow dilatation, we also extend the DNS database by simulating a turbulent flame, with the gas density being decreased due to combustion-induced heat release.

In the rest of this article we first derive evolution equations for the decomposed parts of $S_d$, see § 2. Direct numerical simulation data used for exploring the equations are discussed in § 3. Results are reported in § 4 and conclusions are summarised in § 5.

2. Derivation of equations

In this section we derive new transport and surface-averaged evolution equations for the three components of the displacement speed $S_d$. These components represent reaction, tangential diffusion and normal diffusion contributions to $S_d$. Here, the derivation is described primarily within the ‘surface average’ framework. Derivation of the transport equations within the ‘local’ framework is given in appendix A. For completeness, the appendix also includes the evolution equation for the ordinary displacement speed $S_d$, previously derived by Yu & Lipatnikov (2019c).

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 2002) 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 1985; Kollmann & Chen 1994; Vervisch et al. 1995). 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 1988; Candel & Poinsot 1990) 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 (2019c), 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. (2020), (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. Computational set-up

The above evolution equations derived for different $\phi$, i.e. (2.29) or its equivalent (2.30), were applied to analyse two different DNS data sets obtained from statistically 1-D, initially planar, single-reaction waves propagating in a homogeneous, isotropic, statistically stationary forced turbulence. One data set (case B in the following) was selected from a big DNS database (Yu, Lipatnikov & Bai 2014; Yu, Bai & Lipatnikov 2015b; Elperin et al. 2016; Yu & Lipatnikov 2017a,b; Sabelnikov, Yu & Lipatnikov 2019; Yu & Lipatnikov 2019a) obtained from dynamically passive reaction waves. Such a wave affects neither the fluid density nor the fluid viscosity. Therefore, the wave does not affect the turbulence. These simplifications allowed us to sample more statistics, as will be discussed later, and facilitated analysing and interpreting numerical results.

The second data set (case C) is largely identical to case B, but the thermal expansion was enabled, with the heat diffusivity being equal to the mass diffusivity (i.e. the Lewis number $Le=1$). Accordingly, the wave C affected the fluid density, viscosity and turbulence. In both cases, (i) $\rho \mathcal {D}=\textrm {const}$ and (ii) both $\mathcal {D}$ and $\mathbb {W}$ depend solely on $c$.

The equations derived above can also be applied to realistic premixed flames with complex chemistry and a non-unity Lewis number. As long as we can define a progress variable $c$ based either on the temperature or a species mass fraction, (2.22) and (2.29) hold for both $\boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}$ and $S_d^n/\mathcal {D}$. However, the assumptions invoked to derive (2.22) for $S_d$ or $S_d^W$ do not necessarily hold, e.g. $\rho \mathcal {D}$ may vary within the flame or the quantities $\mathcal {D}$ and $\mathbb {W}$ may depend not only on $c$ but also on other mixture characteristics. Such a complicated problem can be handled using the method developed in the present work, but at a cost of introducing additional terms when extending the derivations leading to (A18) and (A24).

3.1. Direct numerical simulation database

Since the DNS attributes are discussed in detail elsewhere (Yu et al. 2014, 2015b; Elperin et al. 2016; Yu & Lipatnikov 2017a,b, 2019a), 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 (2005) 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 2017b). 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 2015a), 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 2014) 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 2017b) by appropriately adjusting $\boldsymbol {f}_u$ (and, hence, $\boldsymbol {f}$) using the method by Lamorgese et al. (2005).

3.2. Numerical methods

The governing equations are numerically integrated based on an in-house DNS solver (Yu, Yu & Bai 2012) 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 1968) 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. (2012), 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 2013). 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 1996) 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 1995) 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 1968) may appear. While the zero-gradient points are excluded from the definition of the surface average by (2.3), see Yu & Lipatnikov (2019b,c), 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 (2019b).) 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 2019c) 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 2002; Lipatnikov 2012).

4. Results and discussion

In § 4.1 we will use the three DNS cases A–C summarized in table 2 to verify that the numerically computed terms in the averaged (2.29) and its counterpart (2.30) are sufficiently accurate. This will be verified for each of the five quantities $\phi \in [S_d,S_d^n,S_d^W,\mathcal {D} \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}},\mathcal {D}| \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}|]$. In §§ 4.2 and 4.3 the equations and their terms will be analysed from a physical perspective. The analysis mainly consists of three parts.

First, while previous numerical assessment of the $S_d$-equation was limited to constant-density reaction waves (Yu et al. 2020), the present work aims at investigating the effect of thermal expansion on the evolution of $S_d$. For this purpose, results computed in cases B and C will be compared. Second, from a quantifiable geometrical perspective enabled by the new averaged equations for $\boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}$ and $| \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}|$, we will explore the link between the self-propagation characteristics of a wave and its geometrical constraints. Third, we will demonstrate the potential of the extended decomposition rules described by (2.26) and (2.31). These rules offer a deeper explanation for the mechanisms responsible for variation in $S_d$ and its constituents. The rules also allow us to explore common trends between corresponding terms in the different equations.

In the following, results computed in all simulated cases are presented in three figures. First, for each of the five aforementioned quantities $\phi$, the six terms that appear in (2.29) and (2.30) are computed. These terms are: (i) the left-hand side term, i.e. the time derivative of $\left \langle \phi \right \rangle _{{s}}$, (ii) the term $\left \langle \mathbb {A} (\phi ) \right \rangle _{{s}}$, (iii) the term $\left \langle \mathbb {B} (\phi ) \right \rangle _{{s}}$, (iv) the term $\left \langle \mathbb {C}(\phi ) \right \rangle _{{s}}$, (v) the right-hand side of (2.29), i.e. $\sum \text {rhs} = \left \langle \mathbb {A} (\phi ) \right \rangle _{{s}} + \left \langle \mathbb {B} (\phi ) \right \rangle _{{s}} + \left \langle \mathbb {C}(\phi ) \right \rangle _{{s}}$, and (vi) the velocity-related term $\mathcal {U}(\phi )$. All relevant equation terms are listed in table 1. In figure 3 such results are plotted at a representative (transient) time instant of $0.125 \tau _F^*$ in the constant-density cases A and B ($\tau _F^*=\tau _F$ in case A or $\tau _F^*=\tau _t^0 \approx \tau _F/2.25$ in case B). Similar results obtained from the fully developed reaction waves ($t_\infty$) are shown in figure 4 in all three cases. It is emphasized that all cases meet the requirements for the extended decomposition relations described in (2.26) and (2.31) to be valid. Figures 3 and 4 are arranged such that any term selected from a sub-figure in the top row is equal to the sum of the corresponding terms plotted in the same style in the three first sub-figures directly below. For example, the term $\left \langle \mathbb {A} (\phi ) \right \rangle _{{s}}$, shown in green triangles in figure 3(a), is equal to the sum of the curves shown in green triangle in figures 3(c), 3(e) and 3(g).

Figure 3. Various terms in (2.29) and the term $\mathcal {U}$ in (2.30) evaluated for each $\phi \in [ S_d, S_d^n, S_d^W, \mathcal {D} \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}, \mathcal {D}| \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}|]$ (from top to bottom) in two constant-density cases A (a,c,e,g,i) and B (b,d,f,h,j) at a representative time instant of $t=0.125\tau _F^*$. All terms are normalized based on a length unit of $\delta _F$, a speed unit of $S_L$ and a density unit of $\rho _u$.

Figure 4. Various terms in (2.29) and the term $\mathcal {U}$ in (2.30) evaluated for each $\phi \in [ S_d, S_d^n, S_d^W, \mathcal {D} \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}, \mathcal {D}| \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}|]$ (from top to bottom) in two constant-density cases A (a,d,g,j,m) and B (b,e,h,k,n) and in the variable-density case C (c,f,i,l,o) at the fully developed stage of reaction wave evolution. All terms are normalized based on a length unit of $\delta _F$, a speed unit of $S_L$ and a density unit of $\rho _u$.

In § 4.3 we present all sub-terms ($\mathbb {A} _0$, $\mathbb {A} _{1u}$, $\mathbb {A} _{2}$, $\mathbb {B} _{1u}$, $\mathbb {B} _{2u}$, $\mathbb {B} _{3u}$, $\mathbb {B} _{4}$, $\mathbb {B} _{5}$ and $\mathbb {B} _{6}$) for all relevant $\phi$. Plotted in figure 5 are these sub-terms in all three cases at the fully developed state $t_\infty$. For reference, a 1-D steady state version of case C (referred to as case L in the following) is also included in the bottom row of figure 5. Note that all displayed non-zero terms in case L are coupled by three relations specified in the figure caption. Some of the sub-terms are not displayed in figure 5. They either vanish (by definition or due to adopted case assumption, e.g. curvature terms in case L or dilatation terms in cases A and B), or can be trivially deduced using (2.31), e.g. some sub-terms for $\phi =S_d^n$ and $S_d^c$.

Figure 5. Fully developed sub-terms of $\left \langle \mathbb {A} (\phi ) \right \rangle _{{s}}$ (a,d,g,j), first three sub-terms of $\left \langle \mathbb {B} (\phi ) \right \rangle _{{s}}$ (b,e,h,k), and last three sub-terms of $\left \langle \mathbb {B} (\phi ) \right \rangle _{{s}}$ (c,f,i,l) plotted for four $\phi \in [ S_d, S_d^n, \mathcal {D} \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}, S_d^W]$ using solid lines, dashed lines, dot lines with circles and dash–dot lines, respectively. Lines used to show the same sub-term for different $\phi$ are of the same colour. Results computed in cases A, B, C and L are reported in the first, second, third and fourth rows, respectively. All terms are normalized using $\delta _F$, $S_L$ and $\rho _u$.

To evaluate the left-hand side terms $\partial /\partial t \left \langle \phi \right \rangle _{{s}}$ during transient wave evolution, a sequence of transient values of $\left \langle \phi \right \rangle _{{s}}$ for each $\phi$ (see supplemental figures S6 and S7) is calculated at 20 sampling instants $t_i=(i^2/200) \tau _F^*$, where $i = (1,\ldots ,20)$, and a discrete approximation of the time derivative is applied to this sequence. In supplemental material available at https://doi.org/10.1017/jfm.2020.1095, we also include an extra set of figures (S1–S5) showing the temporal evolution of the above terms conditioned to three representative iso-surfaces in cases A and B. Further numerical details regarding evaluation of various terms in (2.29) or (2.30) and calculation of a surface-averaged quantity according to (2.5) or (2.10), respectively, are given in appendices D.3 and D.2 in Yu & Lipatnikov (2019c). Moreover, characteristics of the term fluctuations for the $S_d$-equation have been reported in a recent paper, see figures 5, 6 and § 4.2.1 in Yu et al. (2020).

4.1. Numerical verification

For (2.29) applied to each of the five quantities $\phi \in [S_d, S_d^n, S_d^W, \mathcal {D} \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}, \mathcal {D} | \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}|]$, figures 3 and 4 indicate a reasonably good match between the left-hand side term $\partial /\partial t \left \langle \phi \right \rangle _{{s}}$, shown in red dots, and the sum of the right-hand side terms $\sum {\text {rhs}}$, shown in blue open circles. Such a match is observed for (i) almost all iso-surfaces of $\hat {c} \in (0,1)$ in all three cases, (ii) almost all time instants $0<t<2\tau _F^*$ in the constant-density cases A and B, and (iii) at the fully developed state $t_{\infty }$ in all three cases. Nevertheless, inspection of figures 3 and 4 reveals small mismatches between the left-hand and right-hand side terms.

More specifically, in the turbulent constant-density case B, a small but visible difference between left-hand side and $\sum$rhs is observed during transient evolution, especially for $S_d$ and $S_d^n$ (see figure 3b,d). On the contrary, there is a nearly perfect match for each of the five quantities $\phi$ in the simple shear flow (case A). This difference between cases A and B is associated with the lack of zero-gradient points in case A. As discussed by Yu & Lipatnikov (2019c), the error can be reduced by increasing the number of realizations that are used for sampling transient statistics at later evolution stages, i.e. when an increased amount of zero-gradient points appear on the reaction wave surface disturbed by turbulence. The present data have been obtained using a realization number of $M=170$ to sample transient $c^t$, $M=780$ to sample fully developed $c^s$ in case B and $M=1800$ to sample $c^s$ in case C. Thus, $M$ is significantly larger when statistics are sampled from the fully developed reaction waves. Indeed, in line with the above reasoning, the error is smaller at the fully developed state $t_\infty$; see figure 4 for cases B and C.

The observed reasonable match between the left-hand side and $\sum \text {rhs}$ terms implies that (2.29) for five different $\phi$ can be explored with sufficient confidence by analysing DNS data obtained from turbulent reacting waves both with and without density variation. For instance, a satisfactory match is found in case B, which has a highly complicated reaction zone and a number of zero-gradient points where $\lvert \boldsymbol {\nabla } c \rvert =0$. It should be pointed out that, for (2.29) applied to $\mathcal {D} \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}$, the left-hand side term $\partial /\partial t \left \langle \mathcal {D} \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}} \right \rangle _{{s}}$ is always small (red dots in figure 3g,h), because the averaged mean curvature stays close to zero, i.e. $\left \langle \mathcal {D} \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}} \right \rangle _{{s}} |_{t,\hat {c}} \approx 0$ holds for all time and all iso-surfaces (supplemental figures S7 and S6). This is a consequence of the global geometrical constraint of a statistically planar wave. Such a case configuration justifies an additional investigation using the absolute curvature equation, (2.29) for $\mathcal {D} | \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}|$, because the averaged absolute curvature should be positive in a three-dimensional (3-D) turbulent flow, thus yielding a non-zero left-hand side term $\partial /\partial t \left \langle \mathcal {D}| \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}| \right \rangle _{{s}}$. Figures 3(i) and 3(j) confirm that the absolute curvature equation can also be computed accurately.

4.2. Main equation terms

Concerning overall effects due to the two main factors, turbulence and thermal expansion, a first observation can be made by comparing the range of values that the most relevant equation terms assume in the different DNS cases. For $\phi$ being $S_d$ and its three constituents, the value range is around $[-4, 4]$, $[-40, 40]$ and $[-40, 40]$ in case A, B and C, respectively, cf. the first three rows in figure 4 or 5. Thus, an increase in $Ka$ from 0.9 in case A to 39 in case B results in increasing the range by an order of magnitude. Second, the bottom row in figure 5 shows that the same value range of $[-40, 40]$ is also obtained from the 1-D reference laminar flame L without any turbulence. Therefore, both thermal expansion (case L) and turbulence (case B) yield a comparable increase in the value range with respect to case A. It is worth stressing, however, that dependencies of various terms on $\hat {c}$, plotted for cases B and L in the second and bottom rows, respectively, in figure 5, show substantially different trends. Accordingly, despite comparable value ranges, turbulence and thermal expansion differently affect various terms in the discussed equations. Third, as far as the value range is concerned, thermal expansion and turbulence do not enhance one another in case C, i.e. the same value range of $[-40, 40]$ has been computed. As far as the dependencies of various terms on $\hat {c}$ are concerned, trends shown for cases B and C in figure 5 look similar, but substantially different from trends observed in case L. Accordingly, thermal expansion effects appear to be of less importance when compared with turbulence effects under conditions of the present simulations. Indeed, fourth, figure 4 shows that the effect of thermal expansion does not cause any significant change in the four terms $\left \langle \mathbb {A} \right \rangle _{{s}}$, $\left \langle \mathbb {B} \right \rangle _{{s}}$, $\left \langle \mathbb {C} \right \rangle _{{s}}$ and $\mathcal {U}$, i.e. results computed in the constant-density case B and in the variable-density case C appear to be qualitatively similar.

It may also be noticed that, without heat release, a homogeneous turbulence alone tends to flatten the profile for all (sub-)terms along $\hat {c}$ outside the reaction zone, cf. terms conditioned to $\hat {c}<0.4$ in case B with the corresponding terms in case A in figure 4 or 5. For terms conditioned inside the reaction zone, i.e. $\hat {c}>0.6$, the addition of thermal expansion brings noticeable effects. For instance, the terms $\mathcal {U}(S_d^W)|_{t_\infty }$ and $\left \langle \mathbb {C}(S_d^W) \right \rangle _{{s}} |_{t_\infty }$ are very close in the constant-density case B (figure 4h), but are clearly different in the variable-density case C (figure 4i).

Lines with red plus symbols in figures 3 and 4 show that the magnitude of the ‘direct’ flow effect term, $\mathcal {U}(\phi )|_{\hat {c},t}$, varies in a significantly wider range when compared with terms $\left \langle \mathbb {A} (\phi ) \right \rangle _{{s}} |_{\hat {c},t}$, $\left \langle \mathbb {B} (\phi ) \right \rangle _{{s}} |_{\hat {c},t}$ or $\left \langle \mathbb {C}(\phi ) \right \rangle _{{s}} |_{\hat {c},t}$ for all $\phi$ with the exception of mean curvature in all three DNS cases at different instants. It is worth remembering, however, that the ‘indirect’ flow effect term, $\mathcal {R}(\phi )|_{\hat {c},t_\infty }$, counterbalances $\mathcal {U}(\phi )|_{\hat {c},t_\infty }$ in the fully developed flames, see (2.30), because the left-hand side terms $\partial /\partial t \left \langle \phi \right \rangle _{{s}}$ vanish in this limiting case. As far as developing flames are concerned, the sign of $\mathcal {U}(\phi )|_{\hat {c},t}$ controls the sign of $\partial /\partial t \left \langle \phi \right \rangle _{{s}}$ for various $\phi$ (with the exception of mean curvature). In case A, the magnitudes of the two terms are close at $\hat {c}<0.8$, but $\mathcal {R}(\phi )|_{\hat {c},t_\infty }$ significantly reduces the magnitude of $\partial /\partial t \left \langle \phi \right \rangle _{{s}}$ in case B; see figure 3. A dominant role played by the flow is more evident at the early time immediately after releasing the initial planar wave in the constant-density cases A and B. In fact, any non-zero left-hand side term $\partial /\partial t \left \langle \phi \right \rangle _{{s}}$ can only be kick-started by the initial flow term $\mathcal {U}(\phi )|_{t=0}$ (shown in the supplemental figures), whereas all sub-parts within $\mathcal {R}(\phi )|_{t=0}$ in (2.30) should vanish, because the initial quantities are uniform in space. Outside the reaction zone ($\hat {c}<0.4$) in all three cases, the term $\mathcal {U}(S_d)$ is mainly contributed by its normal diffusion part $\mathcal {U}(S_d^n)$ (second row) while its curvature and reaction parts stay close to zero. Both terms $\mathcal {U}(S_d)|_{\hat {c}}$ and $\mathcal {U}(S_d^n)|_{\hat {c}}$ change sign from positive to negative with an increase in $\hat {c}$. Such sign-flipping behaviour at the fully developed stage is primarily attributed to the sub-term $\left \langle \mathbb {B} _{2u}(S_d^n) \right \rangle _{{s}}$; see red lines in figures 5(b), 5(e) and 5(h).

Regarding the flow curvature contribution term $\mathcal {U}(\mathcal {D} \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}})$ shown in the fourth rows of figures 3 and 4, its value vanishes in constant-density cases A and B which can be explained by the net flow effect on creating positively curved surface elements cancels out with those on negatively curved elements. However, the flow curvature term attains a small negative value in the fully developed variable-density case C (figure 4l) due to non-zero dilatation; see $\left \langle \mathbb {A} _{1u}(\mathcal {D} \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}) \right \rangle _{{s}} |_{t_\infty }$ in figure 5(g). On the contrary, the absolute counterpart $\mathcal {U}(\mathcal {D} | \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}|)|_{\hat {c},t}$ takes a positive value of significant magnitude for all $\hat {c}$ shortly after releasing the initial planar wave (the bottom row in figure 3), reflecting that the flow is very efficient to bend and wrinkle initially planar iso-surfaces, whether in a concave or a convex way. In constant-density cases A and B, which have a constant diffusivity $\mathcal {D}$, the term $\mathcal {U}(\mathcal {D} | \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}|)|_{\hat {c},t}$ shows a flat dependence on $\hat {c}$, suggesting that the effect of bending of surfaces by the flow acts similarly on all different iso-surfaces. The $\hat {c}$-profiles of the fully developed $\mathcal {U}(\mathcal {D} | \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}|)$ are comparable in the constant-density case B and the variable-density case C. This observation indicates, however, that the effectiveness of the flow to bend surfaces is significantly inhibited toward the reaction zone due to density variations in case C. Indeed, $\mathcal {U}(| \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}|)=\mathcal {U}(\mathcal {D} | \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}|)/\mathcal {D}$, with $\mathcal {D}$ being increased by a factor of seven with increasing $\hat {c}$ in case C. It is worth mentioning that, in the turbulent cases B and C, the fully developed term $\mathcal {U}(\mathcal {D}| \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}|)|_{t_\infty }$ remains positive with a significant magnitude, while the same term stays close to zero in case A. This difference appears to indicate that a stationary wave shape adapts itself to a simple shear flow, thus making the curvature almost stationary everywhere.

Now we turn to examining surface-averaged terms in (2.29) for various $\phi$. First, comparison of the relative magnitudes of various terms (i.e. left-hand side, $\left \langle \mathbb {A} \right \rangle _{{s}}$, $\left \langle \mathbb {B} \right \rangle _{{s}}$ and $\left \langle \mathbb {C} \right \rangle _{{s}}$) in the equations for $\mathcal {D} \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}$ and $S_d^W$ shows that, in all three cases A–C, there are only two dominant terms with opposite signs for most of the time. Specifically, in (2.29) for $\mathcal {D} \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}$, the dominant terms are $\left \langle \mathbb {B} \right \rangle _{{s}}$ and $\left \langle \mathbb {C} \right \rangle _{{s}}$, with the left-hand side vanishing and $\left \langle \mathbb {A} \right \rangle _{{s}}$ being small. In (2.29) for $S_d^W$, the dominant terms are $\left \langle \mathbb {A} \right \rangle _{{s}}$ and $\left \langle \mathbb {C} \right \rangle _{{s}}$, with $\left \langle \mathbb {B} \right \rangle _{{s}}$ vanishing and the left-hand side being small; see figures 3(f,h) and 4(g–l).

Regarding the term $\left \langle \mathbb {C}(\phi ) \right \rangle _{{s}}$ for various $\phi$ with the exception of $\phi =S_d^n$, it evolves (from an initial zero in cases A and B) to become negative with substantial magnitude when compared with the other two terms $\left \langle \mathbb {A} ( \phi ) \right \rangle _{{s}}$ and $\left \langle \mathbb {B} ( \phi ) \right \rangle _{{s}}$ in all three cases A–C; see left-pointing triangles in figure 4 and supplemental figures. This observation indicates a significant correlation between $\phi$ and $\mathcal {K}$. Such a correlation does not seem to be surprising at first glance, because the considered $\phi$ are parts of $\mathcal {K}$ according to (2.6). However, $\left \langle \mathbb {C}(S_d^n) \right \rangle _{{s}} |_{t_\infty ,\hat {c}}$ nearly vanishes outside the reaction zone ($\hat {c}<0.4$); see figures 4(e) and 4(f). This result indicates a weak correlation between $\mathcal {K}$ and $S_d^n$, thus implying that the normal diffusion contribution to displacement speed weakly affects the stretch rate outside the reaction zone.

Let us also note that, when compared with the term $\left \langle \mathbb {C}(\mathcal {D} \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}) \right \rangle _{{s}}$, its absolute counterpart $\left \langle \mathbb {C}(\mathcal {D}| \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}|) \right \rangle _{{s}}$ attains a much larger magnitude (primarily negative) during late wave evolution in the turbulent cases B and C. Since $\mathcal {D}\mathcal {K} \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}} =\mathcal {D} a_t \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}} - \mathcal {D} S_d | \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}|^2$, the difference in $\mathcal {D}\mathcal {K} \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}$ and $\mathcal {D} \mathcal {K} | \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}|$ consists solely of the sign of the first term, $\mathcal {D} a_t \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}$ and $\mathcal {D} a_t | \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}|$, respectively, which involves the strain rate $a_t$. Accordingly, the aforementioned difference in the magnitudes of $\left \langle \mathbb {C}(\mathcal {D} \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}) \right \rangle _{{s}}$ and $\left \langle \mathbb {C}(\mathcal {D}| \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}|) \right \rangle _{{s}}$ implies that the surface-averaged strain rate term $\mathcal {D} a_t| \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}|$ is significantly larger than the surface-averaged term $\mathcal {D} S_d | \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}|^2$ in the turbulent cases B and C. On the other hand, $\left \langle \mathbb {C}(\mathcal {D}| \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}|) \right \rangle _{{s}} |_{\hat {c}}$ in case A evolves to become more similar to $\left \langle \mathbb {C}(\mathcal {D} \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}) \right \rangle _{{s}} |_{\hat {c}}$, because the simple shear wave shape evolves to a deep cusp with large (positive) curvature, which dwarfs a minor crest of small (negative) curvature; see the bottom row in supplemental figure S8.

The term $\left \langle \mathbb {A} (\phi ) \right \rangle _{{s}}$, which is associated with evolution of the separation distance between iso-surfaces, flips sign during transient wave evolution for some considered $\phi$; see supplemental figures S1–S3. This observation is associated with transition from thinning of a turbulent reaction wave during an earlier stage of the wave development to rebroadening of the wave during a later stage. This transition is discussed in detail elsewhere (Yu et al. 2019). With exception of a part of $\left \langle \mathbb {A} (S_d^n) \right \rangle _{{s}}$ conditioned to large $\hat {c}$, the fully developed $\left \langle \mathbb {A} (\phi ) \right \rangle _{{s}} |_{\hat {c},t_\infty }$ eventually evolves to become positive in cases A–C for almost all $\hat {c}$ and $\phi$. The fully developed term $\left \langle \mathbb {B} (S_d) \right \rangle _{{s}} |_{t_\infty }$ nearly vanishes outside the reaction zone ($\hat {c}<0.4$) in the constant-density cases A and B. This observation is associated with mutual cancellation of contributions due to a negative normal diffusion term $\left \langle \mathbb {B} (S_d^n) \right \rangle _{{s}} |_{t_\infty }$ and a positive curvature diffusion term $\left \langle \mathbb {B} (\mathcal {D} \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}) \right \rangle _{{s}} |_{t_\infty }$; see figures 4(a), 4(b), 4(d), 4(e), 4(j) and 4(k).

4.3. Fully developed sub-terms in $\left \langle \mathbb {A} \right \rangle _{{s}}$ and $\left \langle \mathbb {B} \right \rangle _{{s}}$

In this section various sub-terms in the two fully developed terms $\left \langle \mathbb {A} \right \rangle _{{s}} |_{t_\infty }$ and $\left \langle \mathbb {B} \right \rangle _{{s}} |_{t_\infty }$ are examined in detail. For brevity, the averaging operator is omitted in expressions for these sub-terms in the present subsection.

In terms of value variation along $\hat {c}$, $\mathbb {B} _{2u}(S_d)$ and $\mathbb {B} _5(S_d)$ (both plotted in red lines in figure 5) are two of the most prominent sub-terms in each of the four cases A, B, C and L. In the first three cases A–C, the term $\mathbb {B} _{2u}(S_d)$ or $\mathbb {B} _5(S_d)$ monotonically decreases or increases along $\hat {c}$, respectively. In case L the two terms vary non-monotonically and reach a high peak inside the reaction zone. Therefore, allowing turbulence to disturb an exothermic reaction planar flame flushes the non-monotonic peaks for both terms. In case L these two terms form a mirror pair (figure 5k,l). In cases A–C variations in these two terms along $\hat {c}$ show opposite trends, thus indicating partial mutual cancellation of the two terms. However, the mirror relation between these two terms does not necessarily hold under all conditions. For instance, the hypothetical configuration considered in appendix B assumes zero flow velocity, with all flow-related sub-terms including $\mathbb {B} _{2u}(S_d)$ vanishing. On the contrary, the term $\mathbb {B} _5(S_d)$ remains significant in that case; see (B12) reduced from (B7).

Comparison of cases C and B shows that thermal expansion affects both aforementioned terms, with the change in $\mathbb {B} _{2u}(S_d)$ being more significant. The difference in $\mathbb {B} _{2u}(S_d)$ between cases C and B is largely contributed by its curvature part, $\mathbb {B} _{2u}(S_d^c)$, representing a correlation between the strain rate and curvature tensors. When only a single effect of either flow turbulence (i.e. case B) or thermal expansion (i.e. case L) is present, the term $\mathbb {B} _{2u}(S_d^c)$ remains small. Accordingly, the large value of this term observed in case C is attributed to the interaction between the two effects.

Let us turn to examining the terms $\mathbb {B} _{1u}(S_d)$ and $\mathbb {B} _{3u}(S_d^c)$. They are both present inside the curvature equation and represent the correlation between a local normal direction of the perturbed wave (${\boldsymbol {n}}$) and a flow-related vector involving a high-order derivative ($\nabla ^2 {\boldsymbol {u}}$ and $\boldsymbol {\nabla } a_n$, respectively). In a random turbulent flow these vectors are unlikely to align with each other. Consequently, the two terms are small in the constant-density cases A and B. On the other hand, in the 1-D reference case L, these two terms form a mirror pair, with a large peak being observed at the reaction zone. Then, in case C, i.e. after allowing turbulence to disturb the planar exothermic flame of case L, these two terms are only moderately modified and still inherit the peak location from case L.

The term $\mathbb {B} _{4}(S_d^c)$ represents tangential diffusion of $S_d$ and has a small magnitude in all cases. On the contrary, the Laplacian diffusion term $\mathbb {B} _{4}(S_d)$ can take non-zero values, with its magnitude being moderately affected by turbulence and thermal expansion. As the turbulence acts to wrinkle iso-surfaces, it is responsible for the positive, flat profile of two comparable sub-terms, i.e. $\mathbb {B} _6(S_d)$ and its part $\mathbb {B} _6(S_d^c)$; see cases B and C in figure 5(f,i). These two terms are barely affected by thermal expansion and vanish in the laminar flame L.

Finally, let us examine sub-terms of $\mathbb {A}$. In the 1-D case L, dilatation due to thermal expansion directly introduces three non-zero terms: positive $\mathbb {A} _{1u}(S_d)$ and $\mathbb {A} _{1u}(S_d^W)$ and the term $\mathbb {A} _{1u}(S_d^n)$ whose sign varies. Indirectly, the thermal expansion induces also three corresponding $\mathbb {A} _2$-terms that mirror the $\mathbb {A} _{1u}$-terms. Note that, due to (3.7), $\mathbb {A} _{1u}(S_d) = S_d^2 \lvert \boldsymbol {\nabla } c \rvert /(c+1/\theta )$ should stay positive everywhere. Under the influence of turbulence, an additional part $\mathbb {A} _{1u}(\mathcal {D} \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}})$ emerges, cf. cases C and L. This part is positive, implying a positive correlation between curvature and dilatation. The full dilatation term $\mathbb {A} _{1u}(S_d)$ is therefore boosted by turbulence, while its two constituents $\mathbb {A} _{1u}(S_d^n)$ and $\mathbb {A} _{1u}(S_d^W)$ barely change.

Comparing cases A and B outside the reaction zone, one can see that the turbulence triggers a negative $\mathbb {A} _2(S_d^c)$ and a positive $\mathbb {A} _2(S_d^n)$ of comparable magnitude, thus resulting in a small value of $\mathbb {A} _2(S_d)$. Comparison of different $\mathbb {A} _2$-terms in cases B and C (figures 5(d) and 5(g), respectively) shows that thermal expansion results in a rather significant change in $\mathbb {A} _2(S_d)$ at large $\hat {c}$, controlled mainly by the reaction term $\mathbb {A} _2(S_d^W)$, whose value is similar to the one in the 1-D case L. In case C the two other terms $\mathbb {A} _2(S_d^c)$ and $\mathbb {A} _2(S_d^n)$ are of the same magnitude as in A and B, but the dependency of $\mathbb {A} _2(S_d^n)$ on $\hat {c}$ is different. The lack of a substantial influence of thermal expansion on the magnitudes of $\mathbb {A} _2(S_d^c)$ and $\mathbb {A} _2(S_d^n)$ is associated with the capability of turbulence to significantly affect these terms under conditions of the present study. However, turbulence does not directly affect $\mathbb {A} _2(S_d^W)$, whereas thermal expansion increases this term, because $S_d^W$ is inversely proportional to the density.

5. Conclusion

Addressed in this work are different terms in the decomposition of the displacement speed, $S_d$, of an iso-scalar surface within a turbulent reacting wave. The decomposition involves the curvature or tangential diffusion $\mathcal {D} \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}$, the normal diffusion $S_d^n$ and the reaction $S_d^W$ terms. Evolution equations for these surface-averaged terms conditioned on the reaction progress variable were derived within a framework developed recently by Yu & Lipatnikov (2019c). All five equations can be written in the same general form, see (2.22), which satisfies a term-wise decomposition relation given by (2.26). More specifically, the derived equations involve three groups of terms responsible for the evolution of the surface-averaged quantity $\left \langle \phi \right \rangle _{{s}}$: (i) term $\left \langle \mathbb {A} (\phi ) \right \rangle _{{s}}$ representing the rate of change in the separation distance between neighbouring iso-surfaces, (ii) term $\left \langle \mathbb {B} (\phi ) \right \rangle _{{s}}$ representing the surface-following rate of change in molecular fluxes, and (iii) term $\left \langle \mathbb {C}(\phi ) \right \rangle _{{s}}$ representing the surface-averaged correlation between fluctuations in the considered quantity and the local stretch rate. In order to asses the effect of the flow, an alternative group of terms $\mathcal {U}(\phi )$ is defined so that it contains all terms that directly involve the flow velocity.

The averaged equations are used to examine a DNS database obtained from turbulent reaction waves characterized by low and moderate Karlovitz numbers. For each surface-averaged evolution equation, the left-hand and right-hand side terms were shown to match reasonably well on all iso-$c$ surfaces and during the entire time interval addressed in the simulations.

The following findings of the present analysis of the DNS data are worth noting.

1. (i) For a turbulent reacting wave at a moderate Karlovitz number, thermal expansion only weakly affects the trends of the main terms in all surface-averaged evolution equations. Further examination of sub-terms reveals that, in addition to introducing a dilatation term $\mathbb {B} _{1u}(\phi )$, thermal expansion also affects the surface-averaged displacement speed and its components through other routes, in particular, through the two opposite-signed sub-terms $\mathbb {B} _{1u}(S_d^c)$ and $\mathbb {B} _{2u}(S_d^c)$, which represent a curvature/dilatation correlation and a strain-rate/curvature-tensor correlation, respectively.

2. (ii) During the early stage of the reaction wave development, turbulent flow controls the evolution of all studied quantities with the exception of the curvature term $\mathcal {D} \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}$. A small magnitude of the total flow term $\mathcal {U}(\mathcal {D} \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}})$ results from mutual cancellation of contributions from flow-bended surface elements with negative or positive curvature.

3. (iii) In (2.29) for $\mathcal {D} \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol {n}}$ and $S_d^W$, there are two dominant terms of opposite signs.

4. (iv) Outside the reaction zone, the correlation between $S_d^n$ and $\mathcal {K}$ nearly vanishes.

5. (v) Among all fully developed surface-averaged sub-terms, $\mathbb {B} _{2u}(S_d)$ and $\mathbb {B} _{5}(S_d)$ are most prominent with large value variation.

Based on the present study, we see several potential directions for future work. While the basic properties of the equations have been established, it remains to be investigated their dependence on various relevant parameters such as Lewis numbers, combustion chemistry, etc. To make a more direct connection with modelling concepts, such as the flame surface density formalism, an extension of the present framework to study the displacement speeds of filtered or averaged quantities can also be of interest. Furthermore, as the present derivation of surface-averaged evolution equations is based on the viewpoint of tracing the evolution of a certain quantity on an iso-surface point, a more systematic examination of such an approach may provide extra insight together with classical Eulerian and Lagrangian approaches.

As a concluding remark regarding potential implications of the present study, we also want to highlight a theoretical analysis of a simplified sub-system of semi-1-D wave evolution over stationary flows; see appendix B. This analysis shows that, under certain conditions, the original governing equation (B1), which contains an arbitrary reaction rate $\mathbb {W}(c)$, can be replaced with a set of (B12) and (B13) for two variables of $S_d$ and $S_d^n$. In this set of equations, the reaction rate function does not explicitly appear and all information pertaining to $\mathbb {W}(c)$ is encoded in the initial conditions.