A New Mathematical Framework for Describing Thin-Reaction-Zone Regime of Turbulent Reacting Flows at Low Damköhler Number

Abstract

Recently, Sabelnikov et al. (2019) developed a phenomenological theory of propagation of an infinitely thin reaction sheet, which is adjacent to a mixing layer, in a constant-density turbulent flow in the case of a low Damköhler number. In the cited paper, the theory is also supported by Direct Numerical Simulation data and relevance of such a physical scenario to highly turbulent premixed combustion is argued. The present work aims at complementing the theory with a new mathematical framework that allows for appearance of thick mixing zones adjacent to an infinitely thin reaction sheet. For this purpose, the instantaneous reaction-progress-variable $c(\mathbf{x},t)$ is considered to consist of two qualitatively different zones, that is, (i) mixture of products and reactants, $c(\mathbf{x},t)<1$, where molecular transport plays an important role, and (ii) equilibrium products, $c(\mathbf{x},t)=1$. The two zones are separated by an infinitely thin reaction sheet, where $c(\mathbf{x},t)=1$ and $|\nabla c|$ is fixed in order for the molecular flux into the sheet to yield a constant local consumption velocity equal to the speed of the unperturbed laminar reaction wave. Exact local instantaneous field equations valid in the entire spaceare derived for the conditioned (to the former, mixing, zone) reaction progress variable, its second moment, and instantaneous characteristic functions. Averaging of these equations yields exact, unclosed transport equations for the conditioned reaction-progress-variable moments and Probability Density Function (PDF), as well as a boundary condition for the PDF at the reaction sheet. The closure problem for the derived equations is beyond the scope of the paper.

1. Introduction

Turbulent reacting flows, for example, premixed turbulent combustion, involve highly non-linear, multiscale, unsteady, three-dimensional phenomena, as well as hundreds or thousands of chemical reactions between tens or hundreds species. Accordingly, development of a rigorous theory of turbulent flames does not seem to be feasible in the foreseeable future. In such a situation, applied Computational Fluid Dynamics (CFD) research into turbulent burning in internal combustion engines relies on phenomenological models. However, most such models invoke a number of different simplifications and suffer from the lack of rigor and, sometimes, self-consistency. Nevertheless, there is a seminal exception known as Bray-Moss-Libby (BML) approach, which is widely discussed in turbulent-combustion textbooks [1,2,3] and review papers [4,5,6,7,8]. Based on a single hypothesis that the probability of finding intermediate (between unburned and burned) states of a reacting mixture is much less than unity everywhere in a premixed turbulent flame, because the unburned reactants and burned products are separated by thin zones (flamelets) that the mixture composition variations are localized to; Bray, Moss, and Libby developed a rigorous mathematical framework and obtained several important theoretical results [9,10,11,12,13,14,15,16,17]. Some of these results were earlier published by Prudnikov [18] in Russian, but were poorly known outside USSR. Currently, the BML approach is widely used as the corner stone of various models adopted in CFD research into turbulent flames and as a framework for analyzing experimental data obtained from turbulent flames.

However, the utility of the BML approach for engineering applications is limited, because the aforementioned hypothesis does not hold in highly turbulent flames, which are in the focus of the contemporary research into burning in piston engines or gas turbine combustors. The point is that, in highly turbulent flames, intermediate mixture states are distributed over thick zones by small-scale turbulent eddies and the probability of finding such mixture states is not negligible [7]. While reaction rates vanish in these thick zones, the local mixture composition varies due to molecular diffusion and small-scale turbulent mixing. Accordingly, such zones are known as preheat or mixing zones.

Within the BML framework, the preheat zones are considered to be inherently laminar and thin in weakly turbulent flames, because, under such conditions, the smallest turbulent eddies are too large to penetrate into the preheat zones and disturb them. Such a combustion regime is known as a flamelet regime [1,2,3,4,5,6,7,8]. When turbulence intensity is increased, the scale of the smallest eddies is decreased and they begin penetrating into the preheat zones and broadening them, thus, making the major hypothesis of the BML approach invalid.

Nevertheless, recent experimental and Direct Numerical Simulation (DNS) data reviewed elsewhere [19,20] indicate that, contrary to the mixing zones, there are other zones that remain to be thin even under conditions of intense turbulence. These are the so-called reaction zones, that is, zones that the major chemical reactions and heat release are localized to. In a single-step-chemistry laminar flame, the reaction zone is significantly thinner than the preheat zone [1]. Accordingly, when the length scale of the smallest eddies is decreased, transition from the flamelet combustion regime to a thin-reaction-zone regime [1] occurs. In this regime, the reaction zones are still thin, because the smallest turbulent eddies do not penetrate into them.

The concept of the thin-reaction-zone regime of turbulent burning [1] is based on the Activation Energy Asymptotic (AEA) theory of premixed combustion, which was pioneered by Zel’dovich and Frank-Kamenetskii [21]. Within the framework of the theory, which is also known as ZFK theory and is widely recognized today [1,2,3,22,23], the reaction zone thickness ${\delta }_r$ vanishes asymptotically, that is, ${\delta }_r\to 0$ when the Zel’dovich number $Ze=\Theta (T_b-T_u)/T_b^2$ tends to infinity. Here, $\Theta$ is the activation temperature of the combustion reaction (the theory deals with single-step chemistry), $T_u$ and $T_b$ are the temperatures of unburned reactants and equilibrium adiabatic combustion products, respectively. Moreover, in the asymptotic limit of $Ze\to \infty$, the so-called consumption velocity $u_c$, that is, the local rate of reactant consumption per unit area of the reaction zone, which may be wrinkled and stretched by the local flow, tends to the unperturbed laminar flame speed $S_L$ provided that molecular transport coefficients are equal to each other for the fuel, oxidant, and heat (the Lewis number $Le=1$).

As already noted [19,20], recent experimental and DNS data indicate that such a thin-reaction-zone regime [1] is observed under a wide range of conditions associated with intense turbulence and, for example, with burning in piston engines or gas turbine combustors. Therefore, the regime has been attracting significant amount of attention. In particular, by considering such a regime in the limiting case of the infinitely thin reaction zone (reaction sheet), a phenomenological theory was recently developed [19,24] for predicting the following scaling for turbulent consumption velocity

$$U_T\propto S_{L}R{e}_t^{1/2}\propto u^{\prime }D{a}^{1/2}$$

(1)

in the constant-density and single-step-chemistry case. Here, $U_T$ is the mean rate of consumption of the deficient reactant per unit area of the mean surface of the reaction sheet, $S_L$ is the laminar reaction-wave speed (the laminar flame speed in the case of combustion), $R{e}_t=u^{\prime }L/\nu$ is the turbulent Reynolds number, $Da={\tau }_T/{\tau }_L$ is the Damköhler number, $u^{\prime }$, L, and ${\tau }_T=L/u^{\prime }$ are the rms turbulent velocity, an integral length scale of the turbulence, and its eddy-turn-over time, respectively, ${\tau }_L={\delta }_L/S_L$ and ${\delta }_L=D/S_L$ are the laminar flame time scale and thickness, respectively, $\nu$ and D are the molecular kinematic viscosity of the reacting mixture and molecular diffusivity of the deficient reactant in the mixture, respectively. In the cited paper, Reference [19], the scaling given by Equation (1) was validated by a large amount of DNS data obtained from constant-density single-reaction waves propagating in intense turbulence under conditions of $Da<1$, with the data also indicating thin reaction zones under the simulation conditions. Moreover, at low Damköhler numbers, the same scaling was documented (i) in experiments on reaction-front propagation in aqueous solutions [25] and flames [26] and (ii) in 2D DNSs of constant-density reaction waves [27], 3D DNSs of thermonuclear deflagration [28], and 3D DNSs of lean methane-air and hydrogen-air flames [29].

Accordingly, the regime of thin reaction zones in intense turbulence is of great interest, but the BML framework is not suitable for exploring this regime, because the major BML hypothesis does not hold at low $Da<1$ and large $Ka\gg 1$. The present work aims at developing a mathematical framework for analyzing the thin-reaction-zone regime associated with thick preheat zones, thus, complementing our recent theory [19,24]. Both the BML and present approaches take an advantage of consideration of a local flame zone to be infinitely thin, but these zones are different. The BML approach deals with the asymptotically limiting case of $Da\to \infty$ and $Ka\to 0$ and addresses infinitely thin flamelets, which contain both preheat and reaction zones. The present approach deals with the asymptotically limiting case of $Ze\to \infty$ and addresses infinitely thin reaction zones, but admits thick preheat zones.

In the next section, the considered problem is stated and the invoked simplifications are discussed. Mean and conditioned transport equations are derived in Section 3, followed by derivation of Probability Density Function (PDF) transport equations in Section 4. Conclusions are summarized in Section 5.

2. Statement of the Problem

Let us consider a statistically 1D, planar, single-reaction wave that propagates from right to left along x-axis in homogeneous isotropic turbulence, but does not affect it, because the reaction changes neither the density $\rho$ nor the viscosity $\nu$. Derivation of equations reported in Section 3 and Section 4 can straightforwardly be extended to the case of a variable density, but the present study is restricted to the simplest case of a constant $\rho$ in order to avoid cumbersome equations. We address the case of a single-step chemistry, the Lewis number $Le=1$, and intense turbulence characterized by a high $R{e}_t\gg 1$, a low $Da<1$ and, hence, a high Karlovitz number $Ka=R{e}_t^{1/2}/Da\gg 1$. Under the above assumptions of single-step chemistry and $Le=1$, the state of the mixture in the reaction wave is fully characterized by a single scalar variable [22], for example, a reaction progress variable c, which is equal to zero and unity in fresh reactants and equilibrium products, respectively.

Let us study a reaction whose rate W depends on c in the extremely nonlinear manner, that is, $W\left[c\right(\mathbf{x},t\left)\right]$ vanishes outside a very thin reaction zone, whose thickness ${\delta }_r$ is much less than both the Kolmogorov length scale ${\eta }_K=LR{e}_t^{-3/4}$ and the laminar wave thickness ${\delta }_L$. In the asymptotically limiting case of $Ze\to \infty$ [21,22], the reaction zone degenerates to a reaction surface $c({\mathbf{x}}_r,t)=1$, which separates the mixture of reactants and products, where $0\le c(\mathbf{x},t)<1$, from the equilibrium products, where $c(\mathbf{x},t)=1$ and $\mathbf{x}\ne {\mathbf{x}}_r$. In the following, (i) the mixture of the reactants and the products is designated with subscript u and is called reactants for brevity, that is, $0\le c<1$ in the reactants and (ii) the equilibrium products are designated with subscript b.

In the reactant region, $\nabla {\mathbf{u}}_u=0$ and a common diffusion equation without source term holds, that is,

$$\begin{array}{c}\hfill \frac{\partial c_u}{\partial t}+{\mathbf{u}}_u·\nabla c_u=D\Delta c_u\end{array}$$

(2)

subject to the to the following Dirichlet and Neumann boundary conditions [22,30,31]

$$\begin{array}{c}\hfill c_u({\mathbf{x}}_{\mathbf{r}},t)=1,|\mathbf{n}·\nabla c_u{|}_r={|\nabla c_u|}_r={\left|\frac{\partial c_u}{\partial n}\right|}_r=\frac{S_L}{D}=\frac{1}{{\delta }_L}\end{array}$$

(3)

on the reactant side of the reaction sheet. Initial conditions read

$$\begin{array}{c}\hfill c_u(\mathbf{x},t=0)=c_0\left(\mathbf{x}\right),c_0(x\to \infty )\to 0,\end{array}$$

(4)

where, without loss of generality, $c_0\left(\mathbf{x}\right)$ may be set using the following ZFK solution [22]

$$\begin{array}{c}\hfill c_0(x\ge 0)=exp\left(-\frac{x}{{\delta }_L}\right),c_0(x<0)=1\end{array}$$

(5)

for the preheat zone of an unperturbed laminar flame whose reaction zone is located at $x=0$. Henceforth, subscript r designates quantities or differential operators taken at the reactant side of the reaction sheet, ${\mathbf{n}}_u=-(\nabla c_u/|\nabla c_u{\left|\right)}_r$ is the unit vector normal to the reaction sheet and pointing to the reactants, and n is spatial distance counted from the reaction sheet along the ${\mathbf{n}}_u$-direction. Equation (3) warrants that the reactant flux $D|\partial c_u{/\partial n|}_r$ towards the reaction sheet is equal to the rate $S_L$ of the reactant consumption per unit sheet area.

The product region, where $c_b(\mathbf{x},t)=1$ is separated from the reactant region by a moving boundary (reaction sheet) whose coordinates are not known beforehand. The initial boundary value problem (IBVP) set above, that is, Equations (2)–(5), is a free boundary problem. A region where the problem should be solved is unknown in advance and has to be found as a part of the solution. Propagation of the free boundary (reaction sheet) can be found due to an additional boundary condition stated on the reaction sheet, that is, the Neumann condition given by the second set of equalities in Equation (3). Such problems appear in different scientific areas, for example, the Stefan problem is a classical free boundary problem [32].

The reaction sheet $c({\mathbf{x}}_r,t)=1$ propagates with respect to the local flow at a displacement speed [1]

$$\begin{array}{c}\hfill S_d^r=D(\Delta c_u/|\nabla c_u{\left|\right)}_r=D{\delta }_L{(\Delta c_u)}_r.\end{array}$$

(6)

Accordingly, motion of the reaction sheet can be tracked using $S_d^r$ and the normal component ${\mathbf{u}}_u({\mathbf{x}}_r,t)·{\mathbf{n}}_u({\mathbf{x}}_r,t)$ of the flow field ${\mathbf{u}}_u(\mathbf{x},t)$, see Equation (10) discussed later.

In an inhomogeneous flow, $S_d^r$ can significantly differ from $S_L$ [31]. For example, if term $D{(\Delta c_u)}_r$ is rewritten in the spherical coordinate system, the displacement speed $S_d^r$ involves an extra term whose magnitude $2D/R_r$ [22] is inversely proportional to the curvature radius $R_r$ of the reaction zone. Accordingly, if $R_r={\delta }_L$ and ${(\partial c_u/\partial r)}_r<0$ (the curvature center in products), the extra term overwhelms $S_L$ and makes $S_d^r$ negative. Strong variations in $S_d^r$ can also be caused by local velocity gradients even in a planar case [31]. In the following theoretical analysis, the reaction sheet is not tracked and Equation (6) is not used.

The problem stated above differs fundamentally from the classical problem of front propagation in a turbulent medium, for example, see a paper by Mayo and Kerstein [33] and references quoted therein. The point is that molecular mixing, that is, the term on the Right Hand Side (RHS) of Equation (2), is not directly addressed in the latter case. Accordingly, the latter problem is associated with $L\gg {\delta }_L$ and $Da\gg 1$ (the BML framework), whereas the present paper addresses the case of a low $Da$. It is worth stressing that molecular mixing smooths out small-scale wrinkles of reaction-zone surface, generated by turbulent eddies, and, therefore, significantly reduces turbulent consumption velocity $U_T$. Indeed, if the Schmidt number $Sc=\mathrm{O}\left(1\right)$ and small-scale turbulent eddies are assumed to be able to wrinkle the reaction surface so that the local curvature radius $R_r$ of a wrinkle is on the order of the Kolmogorov length scale ${\eta }_K$, then the aforementioned mixing contribution $2D/{\eta }_K$ to the displacement speed $S_d^r$ is locally comparable with the Kolmogorov velocity $v_K=u^{\prime }R{e}_t^{-1/4}$ and is much larger than $S_L$ at $Ka\gg 1$. Moreover, the mixing-controlled displacement speed $2D/R_r$ locally affects the wrinkle even after dissipation of short-living Kolmogorov eddies that created the wrinkle. Consequently, the wrinkle is rapidly smoothed out. A recent DNS study [34,35] does show that small-scale (when compared to ${\delta }_L$) wrinkles of a reaction-zone surface are efficiently smoothed out by molecular mixing, with this effect significantly reducing $U_T$ when compared to a linear dependence of $U_T\propto u^{\prime }$ simulated in the case of front propagation in the statistically same turbulence [36].

3. Local, Instantaneous, Mean, and Conditioned Transport Equations

Within the framework of the concept discussed above, the instantaneous reaction-progress- variable field consists of two qualitatively different zones; (i) “reactants”, that is, mixture of products and reactants, where $c_u(\mathbf{x},t)<1$ and molecular transport plays an important role, and (ii) equilibrium products, where $c_b(\mathbf{x},t)=1$. These two zones are separated by a moving infinitely thin reaction surface, where $c_u(\mathbf{x},t)$ and $|\nabla c_u{|}_r$ is fixed, see Equation (3). Derivation of conditionally averaged equations for such a field is not a trivial task. To elucidate relevant difficulties, let us briefly present the zone-conditioning methodology developed in the context of turbulent/non-turbulent interfaces [37,38,39,40] and two-phase flows [41]. It is based on the Navier-Stokes and scalar transport equations, which involve molecular viscous stress tensor and molecular diffusion flux, respectively. Therefore, the velocity and scalar fields, as well as their first and second spatial derivatives, are considered to be continuous in the entire space. Under such common conditions, the methodology consists of the following three steps [37,38,39,40]: (i) introduction of the indicator (characteristic) functions for each region (e.g., turbulent and non/turbulent flows), (ii) multiplication of the governing equations, which are valid in the entire space, with the indicator functions, and (iii) subsequent rearrangement and averaging of various terms in the obtained equations. Finally, the methodology yields transport equations for conditioned quantities, each of them is averaged over a single region of the flow, that is, conditioned to a single state of the fluid.

As far as the problem studied in the present paper is concerned, application of such a well-established methodology [37,38,39,40] to the reaction-progress-variable field $c(\mathbf{x},t)$ defined in the entire space is not straightforward, because the first and second spatial derivatives of the $c(\mathbf{x},t)$-field are discontinuous at the reaction sheet. Indeed, the considered IBVP, see Equations (2)–(5), is solely stated in the reactant region, rather than the entire space. Consequently, the field $c_u(\mathbf{x},t)$ found by solving this IBVP has continuous second spatial derivatives only within the reactant region. After extension to the entire space, the reaction progress variable $c(\mathbf{x},t)$, which is defined in the entire space and is equal to $c_u(\mathbf{x},t)$ in the reactant region, conserves the continuity, but its gradient $\nabla c(\mathbf{x},t)$ is discontinuous at the reaction surface. Indeed, at the reactant side of the surface ${|\nabla c|}_r={|\nabla c_u|}_r$ is given by Equation (3), whereas $|\nabla c|=|\nabla c_b|=0$ at the product side. Therefore, the second derivative ${\partial }^{2}c/\partial n^2$ of $c(\mathbf{x},t)$, taken in the direction normal to the reaction surface, is not a regular function, but is a distribution (generalized function) [42], because the derivative is proportional to Dirac delta function with the support at the reaction surface. Consequently, first, the Laplacian $\Delta c$ involves not only a regular term, but also a term proportional to Dirac delta function. Second, to counterbalance the singular term in the Laplacian, another singular term proportional to the same Dirac delta function should appear in the transport equation for the reaction progress variable extended to the entire space, as will be shown later. Accordingly, the straightforward application of the classical methodology [37,38,39,40], which was developed for fields whose first and second spatial derivatives are continuous in the entire space, to the problem studied here results in appearance of a product of the discontinuous Heaviside function and the singular Dirac delta function, but such a product is not defined in the theory of generalized functions [42]. The present authors are not aware on extension of the discussed methodology to the considered problem, while such an extension may be admissible.

A different approach to solving the problem was developed by Kataoka [43], who studied two-phase flows using indicator functions, but did not multiply the governing equations with the indicator functions. Therefore, in Kataoka’s approach, a product of a discontinuous Heaviside function and a singular Dirac delta function does not appear [43]. Accordingly, it is Kataoka’s approach that is applied here to derive transport equations that describe a thick reaction wave with infinitely thin reaction zone in intense ($Da\ll 1$) turbulence.

Within the framework of that approach, various scalar or vector field quantities $q(\mathbf{x},t)$, for example, velocity, pressure, temperature, or species concentrations, which are continuous or bounded at the reaction surface, are defined in the entire space as follows

$$\begin{array}{c}\hfill q(\mathbf{x},t)={\chi }_u(\mathbf{x},t)q_u(\mathbf{x},t)+{\chi }_b(\mathbf{x},t)q_b(\mathbf{x},t),\end{array}$$

(7)

where subscripts u and b refer to the reactant and product regions, respectively. If the field quantity $q(\mathbf{x},t)$ is discontinuous at the reaction surface, the values of $q_u(\mathbf{x},t)$ and $q_b(\mathbf{x},t)$ are taken at the reactant and product sides of the reaction surface, respectively.The indicator functions ${\chi }_b(\mathbf{x},t)$ and ${\chi }_u(\mathbf{x},t)$ are defined as follows

$$\begin{array}{c}\hfill {\chi }_b(\mathbf{x},t)=H\left[f(\mathbf{x},t)\right],{\chi }_u(\mathbf{x},t)=1-{\chi }_b(\mathbf{x},t)\end{array}$$

(8)

using Heaviside function $H\left[f\right(\mathbf{x},t\left)\right]$ and invoking an arbitrary non-dimensional, continuous, and smooth function $f(\mathbf{x},t)$ that (i) is positive in the products, (ii) negative in the reactants ($0\le c_u<1$), and (iii) satisfies the following condition

$$\begin{array}{c}\hfill f({\mathbf{x}}_r,t)=0\end{array}$$

(9)

on the reaction surface ${\mathbf{x}}_r(\mathbf{x},t)$. The propagation of the surface $f({\mathbf{x}}_r,t)=0$ can be tracked by solving the following kinematic equation

$$\begin{array}{c}\hfill \frac{\partial f}{\partial t}+\nabla ·\left({\mathbf{u}}_{u}f\right)=S_d^r|\nabla f|,\end{array}$$

(10)

where the displacement speed is defined by Equation (6).

The function $f(\mathbf{x},t)$ somehow resembles G-function used to track flame surface [1,44] within the framework of level set methods, but there are fundamental differences. For instance, the level set methods commonly invoke an analytical expression for a displacement speed required to track the surface of $G(\mathbf{x},t)=0$, but, within the framework of the studied problem, $S_d^r$ should be evaluated solving Equations (2) and (6) and is significantly affected by molecular mixing ahead of the reaction surface, as discussed in the previous section.

Taking spatial derivatives of Equation (8), we have

$$\begin{array}{c}\hfill \nabla {\chi }_u=-\nabla {\chi }_b=-\delta \left[f(\mathbf{x},t)\right]\nabla f={\mathbf{n}}_u|\nabla f|\delta \left[f(\mathbf{x},t)\right]={\mathbf{n}}_u\Sigma (\mathbf{x},t),\end{array}$$

(11)

where $\delta \left[f\right(\mathbf{x},t\left)\right]$ is Dirac delta function, the unit normal vector ${\mathbf{n}}_u=-{(\nabla f/|\nabla f\left|\right)}_r=-(\nabla c_u/|\nabla c_u{\left|\right)}_r$ points to the reacting mixture, and the instantaneous reaction-surface density $\Sigma (\mathbf{x},t)$ is equal to

$$\begin{array}{c}\hfill \Sigma (\mathbf{x},t)=|\nabla {\chi }_u|=|\nabla f|\delta \left[f(\mathbf{x},t)\right]\end{array}$$

(12)

using Equation (9). Henceforth, dependence of the instantaneous quantity q on $\mathbf{x}$ and t is often skipped for brevity.

It should be stressed that (i) all field quantities are continuous at the reaction surface and (ii) the decomposition defined by Equation (7) holds both in constant density and variable density flows Decomposition of the spatial derivatives (of arbitrary order) of the field quantity $q(\mathbf{x},t)$ can be obtained by differentiating Equation (7). However, in a general case, differentiation of Equation (7) may yield additional terms when compared to straightforward application of Equation (7) to the derivative referred to. For example, let us consider the decompositions for the gradient $\nabla q$ and the Laplacian $\Delta q$ of the scalar field $q(\mathbf{x},t)$. Taking gradient of Equation (7) yields

$$\begin{array}{c}\hfill \nabla q=\nabla ({\chi }_u{q}_u+{\chi }_b{q}_b)=q_u\nabla {\chi }_u+{\chi }_u\nabla q_u+q_b\nabla {\chi }_b+{\chi }_b\nabla q_b\\ \hfill =(q_b-q_u)\nabla f\delta \left[f(\mathbf{x},t)\right]+{\chi }_u\nabla q_u+{\chi }_b\nabla q_b={\chi }_u\nabla q_u+{\chi }_b\nabla q_b.\end{array}$$

(13)

The term with Dirac delta function vanishes, because the scalar $q(\mathbf{x},t)$ is continuous at the reaction sheet. Therefore, both differentiation of Equation (7) and its straightforward application to $\nabla q$ yield the same result even if $\nabla q$ may be discontinuous at the reaction surface.

On the contrary, as far as the Laplacian $\Delta q$ is concerned, taking divergence of Equation (13) results in

$$\begin{array}{c}\hfill \Delta q=\nabla ·\nabla q=\nabla ·({\chi }_u\nabla q_u+{\chi }_b\nabla q_b)=\nabla {\chi }_u·\nabla q_u+{\chi }_u\Delta q_u+\nabla {\chi }_b·\nabla q_b+{\chi }_b\Delta q_b\\ \hfill =(\nabla q_b-\nabla q_u)·\nabla f\delta \left[f(\mathbf{x},t)\right]+{\chi }_u\Delta q_u+{\chi }_b\Delta q_b.\end{array}$$

(14)

Equation (14) reduces to the result of the straightforward application of Equation (7) to the Laplacian $\Delta q$ only if the gradient $\nabla q(\mathbf{x},t)$ is continuous at the reaction surface, that is, $\nabla q_u({\mathbf{x}}_r,t)=\nabla q_b({\mathbf{x}}_r,t)$.

Let us apply Equations (7), (13) and (14) to the reaction-progress-variable field $c(\mathbf{x},t)$ extended to the entire space. Using $c_b=1$, we arrive at

$$\begin{array}{c}\hfill c={\chi }_u{c}_u+{\chi }_b{c}_b={\chi }_u{c}_u+{\chi }_b,\end{array}$$

(15)

$$\begin{array}{c}\hfill \nabla c={\chi }_u\nabla c_u+{\chi }_b\nabla c_b={\chi }_u\nabla c_u,\end{array}$$

(16)

that is, Equation (7) holds for the discontinuous (but bounded) gradient $\nabla c$, and

$$\begin{array}{c}\hfill \Delta c=-\frac{1}{{\delta }_L}|\nabla f|\delta \left(f\right)+{\chi }_u\Delta c_u=-\frac{\Sigma }{{\delta }_L}+{\chi }_u\Delta c_u.\end{array}$$

(17)

The first term on the RHS of Equation (17) results from the first term on the RHS of the last equality in Equation (14), because (i) $\nabla c_b=0$, (ii) $-\nabla c_u·\nabla f\delta \left[f(\mathbf{x},t)\right]=\nabla c_u·{\mathbf{n}}_u|\nabla f|\delta \left[f(\mathbf{x},t)\right]=\nabla c_u·{\mathbf{n}}_u\Sigma$ due to Equation (12), and (iii) at the reaction sheet, $\nabla c_u·{\mathbf{n}}_u=-{|\nabla c_u|}_r={\delta }_L^{-1}$ due to Equation (3). Straightforward differentiation of Equation (16) yields the same result. Indeed,

$$\begin{array}{c}\hfill \Delta c=\nabla {\chi }_u\nabla c_u+{\chi }_u\Delta c_u=\nabla c_u·{\mathbf{n}}_u\Sigma +{\chi }_u\Delta c_u\end{array}$$

(18)

using Equation (11). As already shown, see items (ii) and (iii) below Equation (17), $\nabla c_u·{\mathbf{n}}_u\Sigma ={\delta }_L^{-1}\Sigma$. Thus, Equation (17) is recovered.

On the contrary, the straightforward application of Equation (7) to the Laplacian $\Delta c$ yields a wrong equation

$$\begin{array}{c}\hfill \Delta c={\chi }_u\Delta c_u+{\chi }_b\Delta c_b={\chi }_u\Delta c_u,\end{array}$$

(19)

which does not involve the singular term $\nabla {\chi }_u\nabla c_u=-{\delta }_L^{-1}|\nabla f|\delta \left(f\right)=-{\delta }_L^{-1}\Sigma$. The point is that the contribution from the infinitesimally thin reaction surface associated with the discontinuous gradient of the reaction progress variable and, hence, infinitely large Laplacian $\Delta c$ is overlooked in Equation (19). On the contrary, Equation (19) holds for a $c(\mathbf{x},t)$-field whose spatial derivatives are continuous at a reaction surface.

Ensemble-averaged Equation (7) reads

$$\begin{array}{c}\hfill \overline{q}(\mathbf{x},t)={\gamma }_u(\mathbf{x},t){\langle q\rangle }_u(\mathbf{x},t)+{\gamma }_b(\mathbf{x},t){\langle q\rangle }_b(\mathbf{x},t),\end{array}$$

(20)

where ${\gamma }_u(\mathbf{x},t)\equiv \overline{{\chi }_u(\mathbf{x},t)}$ and ${\gamma }_b(\mathbf{x},t)\equiv \overline{{\chi }_b(\mathbf{x},t)}=1-{\gamma }_u(\mathbf{x},t)$ are probabilities of finding the reactants and products, respectively, and the conditioned quantities ${\langle q\rangle }_u$ and ${\langle q\rangle }_b$ are defined as follows

$$\begin{array}{c}\hfill {\gamma }_b(\mathbf{x},t){\langle q\rangle }_b(\mathbf{x},t)\equiv \overline{{\chi }_b(\mathbf{x},t)q_b(\mathbf{x},t)},{\gamma }_u(\mathbf{x},t){\langle q\rangle }_u(\mathbf{x},t)\equiv \overline{{\chi }_u(\mathbf{x},t)q_u(\mathbf{x},t)}.\end{array}$$

(21)

Transport equation for the indicator function ${\chi }_u(\mathbf{x},t)$ reads

$$\begin{array}{c}\hfill \frac{\partial {\chi }_u}{\partial t}+\nabla ·\left({\mathbf{u}}_u{\chi }_u\right)=-\delta \left(f\right)\left[\frac{\partial f}{\partial t}+\nabla ·\left({\mathbf{u}}_{u}f\right)\right]\\ \hfill =-S_d^r|\nabla f|\delta \left(f\right)=-S_d^r\Sigma \end{array}$$

(22)

using Equations (10) and (12). Transport equation for the indicator function ${\chi }_b(\mathbf{x},t)$ reads

$$\begin{array}{c}\hfill \frac{\partial {\chi }_b}{\partial t}+\nabla ·\left({\mathbf{u}}_b{\chi }_b\right)=\delta \left(f\right)\left[\frac{\partial f}{\partial t}+\nabla ·\left({\mathbf{u}}_{b}f\right)\right]=S_d^r\Sigma ,\end{array}$$

(23)

with ${\mathbf{u}}_b={\mathbf{u}}_u$ at the reaction sheet.

Ensemble-averaged Equations (22) and (23) read

$$\begin{array}{c}\hfill \frac{\partial {\gamma }_u}{\partial t}+\nabla ·\left({\gamma }_u{\langle \mathbf{u}\rangle }_u\right)=-\overline{S_d^r\Sigma }\equiv -{\langle S_d\rangle }_r\overline{\Sigma },\frac{\partial {\gamma }_b}{\partial t}+\nabla ·\left({\gamma }_b{\langle \mathbf{u}\rangle }_b\right)={\langle S_d\rangle }_r\overline{\Sigma },\end{array}$$

(24)

where $\overline{\Sigma }$ is the mean reaction-surface density and ${\langle S_d\rangle }_r$ designates the displacement speed averaged over the reaction sheet. At first glance, Equation (24) looks similar to transport equations for probabilities of finding unburned reactants ($c_u=0$) and equilibrium products, respectively, derived within the BML framework earlier [45,46,47,48]. However, there is an important difference between the studied and flamelet regimes, because $0\le c<1$ in the reactants and, therefore, ${\langle c\rangle }_u\ne 1$ and ${\langle c^2\rangle }_u\ne 1$ in the present case. Accordingly, field equations for ${\langle c\rangle }_u$ and ${\langle c^2\rangle }_u$ are required.

They can be obtained as follows. First,

$$\begin{array}{c}\hfill \frac{\partial {\chi }_u{c}_u}{\partial t}+\nabla ·\left({\chi }_u{\mathbf{u}}_u{c}_u\right)={\chi }_u\left[\frac{\partial c_u}{\partial t}+\nabla ·\left({\mathbf{u}}_u{c}_u\right)\right]+c_u\left[\frac{\partial {\chi }_u}{\partial t}+\nabla ·\left({\mathbf{u}}_u{\chi }_u\right)\right]\\ \hfill ={\chi }_{u}D\Delta c_u-S_d^r\Sigma \end{array}$$

(25)

using Equation (22) and substituting $c_u$ with $c_u=1$ in the last term on the RHS, because $\Sigma$ involves Dirac delta function, see Equation (12). Since

$$\begin{array}{c}\hfill {\chi }_{u}D\Delta c_u=D\nabla ·\left({\chi }_u\nabla c_u\right)-D\nabla c_u·\nabla {\chi }_u\\ \hfill =D\nabla ·\left({\chi }_u\nabla c_u\right)-D{(\nabla c_u)}_r·\mathbf{n}\Sigma =D\nabla ·\left({\chi }_u\nabla c_u\right)+S_L\Sigma \end{array}$$

(26)

by virtue of Equations (3) and (11), Equation (25) also reads

$$\begin{array}{c}\hfill \frac{\partial {\chi }_u{c}_u}{\partial t}+\nabla ·\left({\chi }_u{\mathbf{u}}_u{c}_u\right)=D\nabla ·\left({\chi }_u\nabla c_u\right)+(S_L-S_d^r)\Sigma .\end{array}$$

(27)

Application of a similar method to the field equation for $c_u^2$ yields

$$\begin{array}{c}\hfill \frac{\partial {\chi }_u{c}_u^2}{\partial t}+\nabla ·\left({\chi }_u{\mathbf{u}}_u{c}_u^2\right)=2{\chi }_{u}D{c}_u\Delta c_u-S_d^r\Sigma \end{array}$$

(28)

and

$$\begin{array}{c}\hfill \frac{\partial {\chi }_u{c}_u^2}{\partial t}+\nabla ·\left({\chi }_u{\mathbf{u}}_u{c}_u^2\right)=D\nabla ·\left({\chi }_u\nabla c_u^2\right)-S_d^r\Sigma +2S_L\Sigma -2{\chi }_u{N}_u,\end{array}$$

(29)

where $N_u=D{(\nabla c_u)}^2$ is scalar dissipation in the mixture of reactants and products.

Second, ensemble-averaged Equations (25) and (27) read

$$\begin{array}{c}\hfill \frac{\partial {\gamma }_u{\langle c\rangle }_u}{\partial t}+\nabla ·\left({\gamma }_u{\langle \mathbf{u}c\rangle }_u\right)={\gamma }_{u}D{\langle \Delta c\rangle }_u-{\langle S_d\rangle }_r\overline{\Sigma }\end{array}$$

(30)

and

$$\begin{array}{c}\hfill \frac{\partial {\gamma }_u{\langle c\rangle }_u}{\partial t}+\nabla ·\left({\gamma }_u{\langle \mathbf{u}c\rangle }_u\right)=D\nabla ·\left({\gamma }_u{\langle \nabla c\rangle }_u\right)+(S_L-{\langle S_d\rangle }_r)\overline{\Sigma },\end{array}$$

(31)

respectively. Similarly, ensemble-averaged Equations (28) and (29) read

$$\begin{array}{c}\hfill \frac{\partial {\gamma }_u{\langle c^2\rangle }_u}{\partial t}+\nabla ·\left({\gamma }_u{\langle \mathbf{u}{c}^2\rangle }_u\right)=2{\gamma }_{u}D{\langle c\Delta c\rangle }_u-{\langle S_d\rangle }_r\overline{\Sigma },\end{array}$$

(32)

and

$$\begin{array}{c}\hfill \frac{\partial {\gamma }_u{\langle c^2\rangle }_u}{\partial t}+\nabla ·\left({\gamma }_u{\langle \mathbf{u}{c}^2\rangle }_u\right)=D\nabla ·\left({\gamma }_u{\langle \nabla c^2\rangle }_u\right)-{\langle S_d\rangle }_r\overline{\Sigma }+2S_L\overline{\Sigma }-2{\gamma }_u{\langle N\rangle }_u,\end{array}$$

(33)

respectively. Here, ${\langle \nabla c\rangle }_u$, ${\langle \nabla c^2\rangle }_u$, ${\langle \Delta c\rangle }_u$, ${\langle c\Delta c\rangle }_u$, ${\langle \mathbf{u}c\rangle }_u$, ${\langle \mathbf{u}{c}^2\rangle }_u$, and ${\langle N\rangle }_u$ are defined substituting $q_u=\nabla c_u$, $q_u=\nabla c_u^2$, $q_u=\Delta c_u$, $q_u=c_u\Delta c_u$, $q_u=\mathbf{u}{c}_u$, $q_u=\mathbf{u}{c}_u^2$, and $q_u=N_u$, respectively, into Equation (21).

Summation of Equation (23) with Equation (25) or (27) yields

$$\begin{array}{c}\hfill \frac{\partial c}{\partial t}+\nabla ·\left(\mathbf{u}c\right)={\chi }_{u}D\Delta c_u=D\nabla ·\left({\chi }_u\nabla c_u\right)+S_L\Sigma =D\nabla ·\left({\chi }_u\nabla c_u\right)+W\end{array}$$

(34)

or

$$\begin{array}{c}\hfill \frac{\partial c}{\partial t}+\nabla ·\left(\mathbf{u}c\right)={\chi }_{u}D\Delta c_u=D\Delta c+S_L\Sigma =D\Delta c+W,\end{array}$$

(35)

respectively, using Equations (16) and (12). The singular term $S_L\Sigma$ is the rate W of product creation, which involves Dirac delta function within the present mathematical framework. If differentiation is treated in the sense of generalized functions, Equation (35) reduces to Equation (2) in the reactant region and to the boundary conditions given by Equation (3) at the reaction sheet $f(\mathbf{x},t)=0$.

Similarly, summation of Equation (23) with Equation (28) or (29)

$$\begin{array}{c}\hfill \frac{\partial c^2}{\partial t}+\nabla ·\left(\mathbf{u}{c}^2\right)=2{\chi }_{u}D{c}_u\Delta c_u\end{array}$$

(36)

or

$$\begin{array}{c}\hfill \frac{\partial c^2}{\partial t}+\nabla ·\left(\mathbf{u}{c}^2\right)=D\nabla ·\left({\chi }_u\nabla c_u^2\right)+2S_L\Sigma -2{\chi }_u{N}_u,\end{array}$$

(37)

respectively. The latter equation reads

$$\begin{array}{c}\hfill \frac{\partial c^2}{\partial t}+\nabla ·\left(\mathbf{u}{c}^2\right)=D\Delta c^2+2S_L\Sigma -2N,\end{array}$$

(38)

because $\nabla c_b=N_b=0$, $\nabla c^2={\chi }_u\nabla c_u^2$, and $N={\chi }_u{N}_u$.

Ensemble-averaged Equations (34)–(38) read

$$\begin{array}{c}\hfill \frac{\partial \overline{c}}{\partial t}+\nabla ·\overline{\mathbf{u}c}={\gamma }_{u}D{\langle \Delta c\rangle }_u=D\nabla ·\left({\gamma }_u{\langle \nabla c\rangle }_u\right)+S_L\overline{\Sigma }.=D\Delta \overline{c}+S_L\overline{\Sigma },\end{array}$$

(39)

$$\begin{array}{c}\hfill \frac{\partial \overline{c^2}}{\partial t}+\nabla ·\overline{\mathbf{u}{c}^2}=2{\gamma }_{u}D{\langle c\Delta c\rangle }_u=D\nabla ·\left({\gamma }_u{\langle \nabla c^2\rangle }_u\right)+2S_L\overline{\Sigma }-2{\gamma }_u{\langle N\rangle }_u.\\ \hfill =D\Delta \overline{c^2}+2S_L\overline{\Sigma }-2\overline{N}.\end{array}$$

(40)

The term $S_L\overline{\Sigma }$ on the RHS of Equation (39) is the mean rate $\overline{W}$ of product creation within the present mathematical framework.

The instantaneous bulk consumption rate Q is equal to either $\rho \partial {\chi }_b/\partial t$ or $\rho \partial c/\partial t$ integrated over a volume V that is sufficiently large in order for the convection and molecular fluxes of c to vanish at the volume boundaries. For instance, in a statistically planar 1D case and in the coordinate framework attached to the mean flow (i.e., $\overline{\mathbf{u}}=0$), V is a thick layer that envelopes the entire reaction wave and is parallel to the mean position of the reaction surface. Integration of Equations (23) and (34) yields

$$\begin{array}{c}\hfill \frac{Q}{\rho }={\int }_V\frac{\partial {\chi }_b}{\partial t}\mathrm{d}V={\int }_V{S}_d^r\Sigma \mathrm{d}V,\frac{Q}{\rho }={\int }_V\frac{\partial c}{\partial t}\mathrm{d}V={\int }_V{S}_L\Sigma \mathrm{d}V,\end{array}$$

(41)

respectively, where $A_r={\int }_V\Sigma \mathrm{d}V$ is the reaction-surface area. Ensemble-averaged Equation (41) reads

$$\begin{array}{c}\hfill \frac{\overline{Q}}{\rho }={\int }_V\frac{\partial {\gamma }_b}{\partial t}\mathrm{d}V={\int }_V{\langle S_d\rangle }_r\overline{\Sigma }\mathrm{d}V,\frac{\overline{Q}}{\rho }={\int }_V\frac{\partial \overline{c}}{\partial t}\mathrm{d}V={\int }_V{S}_L\overline{\Sigma }\mathrm{d}V.\end{array}$$

(42)

In particular, in the statistically planar 1D case and in the coordinate framework attached to the mean flow, Equations (41) and (42) read

$$\begin{array}{c}\hfill U_T={\int }_{-\infty }^{\infty }{S}_d^r\Sigma \mathrm{d}x=S_L{\int }_{-\infty }^{\infty }\Sigma \mathrm{d}x\end{array}$$

(43)

and

$$\begin{array}{c}\hfill \overline{U_T}={\int }_{-\infty }^{\infty }{\langle S_d\rangle }_r\overline{\Sigma }\mathrm{d}x=S_L{\int }_{-\infty }^{\infty }\overline{\Sigma }\mathrm{d}x,\end{array}$$

(44)

respectively.

Thus, in spite of the fact that local differences in the displacement speed $S_d^r$ and consumption velocity $S_L$ may be very large and even the product of $S_d^r{S}_L$ may be negative, as discussed in Section 2, the two speeds multiplied with $\Sigma$ and integrated over the reaction-wave volume are equal to one another. In other words, positive and negative fluctuations in $S_d^r-S_L$ completely balance one another after the volume integration. Therefore, in spite of the well-known fact that the local displacement speed cannot characterize the local consumption velocity, the volume-integrated $S_d^r\Sigma$ or ${\langle S_d\rangle }_r\overline{\Sigma }$ properly characterizes the volume-integrated consumption velocity. It is worth remembering, however, that the local ${|\nabla c|}_r$ is given by Equation (3) and does not fluctuate in the studied case of an infinitely thin reaction zone. On the contrary, in the case of a reaction zone of a finite thickness, $|\nabla c|$ conditioned to the reaction zone may fluctuate and these fluctuations may correlate with fluctuations in the local displacement speed. Therefore, the volume-integrated $S_d^r\Sigma$ or ${\langle S_d\rangle }_r\overline{\Sigma }$ may differ from the volume-integrated $S_L\Sigma$ or $S_L\overline{\Sigma }$ in such a case.

Finally, it is worth noting that

$$\begin{array}{c}\hfill {\int }_V\frac{\partial {\chi }_u{c}_u}{\partial t}\mathrm{d}V={\int }_V(S_L-S_d^r)\Sigma \mathrm{d}V=0\end{array}$$

(45)

by virtue of Equations (27) and (41). Consequently,

$$\begin{array}{c}\hfill {\int }_V\frac{\partial {\gamma }_u{\langle c\rangle }_u}{\partial t}\mathrm{d}V={\int }_V(S_L-{\langle S_d\rangle }_r)\overline{\Sigma }\mathrm{d}V=0,\end{array}$$

(46)

see also Equations (31) and (42). The Left Hand Side (LHS) integral vanishes, because consumption of reactants not only reduces the probability ${\gamma }_u$, but also appropriately increases the value of ${\langle c\rangle }_u$.

4. Transport Equations for PDF of Reaction Progress Variable

To begin, let us recall that there are two methods for deriving transport equations for PDF. Such equations can be derived using either characteristic functions [23,49,50] or fine-grained PDF’s [51,52]. Both methods require adaptation to be applied to the problem of propagation of a reaction wave with the infinitely thin reaction zone, that is, the singular behavior of the Laplacian $\Delta c$ of the reaction progress variable $c(\mathbf{x},t)$ extended to the entire space, see Equation (17).

Here, the PDF transport equations for the reaction progress variable are derived adapting the former method, which was pioneered by Kuznetsov [50] and was discussed in detail elsewhere [23]. In particular, the method of instantaneous characteristic functions is adapted following the approach developed by Kataoka [43] and used already in Section 3.

By setting $q=exp\left(i\lambda c\right)$ in Equation (7), the following instantaneous characteristic functions

$$\begin{array}{c}\hfill \phi \left(\lambda \right)={\chi }_u{\phi }_u\left(\lambda \right)+{\chi }_b{\phi }_b\left(\lambda \right),{\phi }_u\left(\lambda \right)=exp\left(i\lambda c_u\right),{\phi }_b\left(\lambda \right)=exp\left(i\lambda \right)\end{array}$$

(47)

are introduced. Then,

$$\begin{array}{c}\hfill \frac{\partial {\chi }_u{\phi }_u}{\partial t}+\nabla ·\left({\mathbf{u}}_u{\chi }_u{\phi }_u\right)={\chi }_u\left[\frac{\partial {\phi }_u}{\partial t}+\nabla ·\left({\mathbf{u}}_u{\phi }_u\right)\right]+{\phi }_u\left[\frac{\partial {\chi }_u}{\partial t}+\nabla ·\left({\mathbf{u}}_u{\chi }_u\right)\right]\\ \hfill =i\lambda {\chi }_u{\phi }_u\left[\frac{\partial c_u}{\partial t}+\nabla ·\left({\mathbf{u}}_u{c}_u\right)\right]-{\phi }_u{S}_d^r\Sigma =i\lambda {\chi }_u{\phi }_{u}D\Delta c_u-exp\left(i\lambda \right)S_d^r\Sigma \end{array}$$

(48)

and

$$\begin{array}{c}\hfill \frac{\partial {\chi }_b{\phi }_b}{\partial t}+\nabla ·\left({\mathbf{u}}_u{\chi }_b{\phi }_b\right)={\phi }_b\left[\frac{\partial {\chi }_b}{\partial t}+\nabla ·\left({\mathbf{u}}_u{\chi }_b\right)\right]=exp\left(i\lambda \right)S_d^r\Sigma ,\end{array}$$

(49)

see Equations (22) and (23). Note that ${\phi }_u=exp\left(i\lambda c_u\right)$ is substituted with $exp\left(i\lambda \right)$ in the last term on the RHS of Equation (48), because this term involves $\Sigma =|\nabla f|\delta \left(f\right)$ and $c=1$ if $f=0$. Summation of Equations (48) and (49) yields

$$\begin{array}{c}\hfill \frac{\partial \phi }{\partial t}+\nabla ·\left(\mathbf{u}\phi \right)=i\lambda {\chi }_u{\phi }_{u}D\Delta c_u.\end{array}$$

(50)

If $\lambda =0$, Equations (48) and (49) reduce to Equations (22) and (23), respectively. Similarly, Equations (25) and (34) can be obtained by taking the first derivatives of Equations (48) and (50) with respect to $i\lambda$ and, subsequently, setting $\lambda =0$. Furthermore, Equations (28) and (36) can be obtained by taking the second derivatives of Equations (48) and (50) with respect to $i\lambda$ and, subsequently, setting $\lambda =0$.

Using Equation (26) and

$$\begin{array}{c}\hfill i\lambda {\phi }_{u}D\nabla ·\left({\chi }_u\nabla c_u\right)=i\lambda D\nabla ·\left({\chi }_u{\phi }_u\nabla c_u\right)-i\lambda D{\chi }_u\nabla c_u·\nabla {\phi }_u\\ \hfill =i\lambda D\nabla ·\left({\chi }_u{\phi }_u\nabla c_u\right)-{\left(i\lambda \right)}^{2}D{\chi }_u{\phi }_u\nabla c_u·\nabla c_u,\end{array}$$

(51)

the RHS of Equation (48) reduces to

$$\begin{array}{c}\hfill i\lambda {\phi }_{u}D\nabla ·\left({\chi }_u\nabla c_u\right)=i\lambda D\nabla ·\left({\chi }_u{\phi }_u\nabla c_u\right)+{\lambda }^2{\chi }_u{\phi }_u{N}_u.\end{array}$$

(52)

Then, Equations (48) and (50) read

$$\begin{array}{c}\hfill \frac{\partial {\chi }_u{\phi }_u}{\partial t}+\nabla ·\left({\mathbf{u}}_u{\chi }_u{\phi }_u\right)\\ \hfill =i\lambda D\nabla ·\left({\chi }_u{\phi }_u\nabla c_u\right)+{\lambda }^2{\chi }_u{\phi }_u{N}_u+i\lambda exp\left(i\lambda \right)S_L\Sigma -exp\left(i\lambda \right)S_d^r\Sigma \end{array}$$

(53)

and

$$\begin{array}{c}\hfill \frac{\partial \phi }{\partial t}+\nabla ·\left(\mathbf{u}\phi \right)=i\lambda D\nabla ·\left({\chi }_u{\phi }_u\nabla c_u\right)+{\lambda }^2{\chi }_u{\phi }_u{N}_u+i\lambda exp\left(i\lambda \right)S_L\Sigma .\end{array}$$

(54)

Again, if $\lambda =0$, Equation (53) reduces to Equation (22). Similarly, Equations (27) and (34) can be obtained by taking the first derivatives of Equations (53) and (54) with respect to $i\lambda$ and, subsequently, setting $\lambda =0$. Furthermore, Equations (29) and (36) can be obtained by taking the second derivatives of Equations (53) and (54) with respect to $i\lambda$ and, subsequently, setting $\lambda =0$.

The first term on the RHS of Equation (53) can be rewritten as follows

$$\begin{array}{c}\hfill i\lambda D\nabla ·\left({\chi }_u{\phi }_u\nabla c_u\right)=D\nabla ·\left({\chi }_{u}i\lambda {\phi }_u\nabla c_u\right)=D\nabla ·\left({\chi }_u\nabla {\phi }_u\right).\end{array}$$

(55)

Substitution of Equations (55) into the RHS’s of Equations (53) and (54) yields

$$\begin{array}{c}\hfill \frac{\partial {\chi }_u{\phi }_u}{\partial t}+\nabla ·\left({\mathbf{u}}_u{\chi }_u{\phi }_u\right)\\ \hfill =D\nabla ·\left({\chi }_u\nabla {\phi }_u\right)+{\lambda }^2{\chi }_u{\phi }_u{N}_u+i\lambda exp\left(i\lambda \right)S_L\Sigma -exp\left(i\lambda \right)S_d^r\Sigma \end{array}$$

(56)

and

$$\begin{array}{c}\hfill \frac{\partial \phi }{\partial t}+\nabla ·\left(\mathbf{u}\phi \right)=D{\nabla }^2\phi +{\lambda }^2{\chi }_u{\phi }_u{N}_u+i\lambda exp\left(i\lambda \right)S_L\Sigma .\end{array}$$

(57)

Ensemble-averaged Equations (50), (54) and (57) read

$$\begin{array}{c}\hfill \frac{\partial \overline{\phi }}{\partial t}+\nabla ·\overline{\mathbf{u}\phi }=i\lambda D\overline{{\chi }_u{\phi }_u\Delta c_u},\end{array}$$

(58)

$$\begin{array}{c}\hfill \frac{\partial \overline{\phi }}{\partial t}+\nabla ·\overline{\mathbf{u}\phi }=i\lambda D\nabla ·\overline{{\chi }_u{\phi }_u\nabla c_u}+{\lambda }^2\overline{{\chi }_u{\phi }_u{N}_u}+i\lambda exp\left(i\lambda \right)S_L\overline{\Sigma },\end{array}$$

(59)

and

$$\begin{array}{c}\hfill \frac{\partial \overline{\phi }}{\partial t}+\nabla ·\overline{\mathbf{u}\phi }=D{\nabla }^2\overline{\phi }+{\lambda }^2\overline{{\chi }_u{\phi }_u{N}_u}+i\lambda exp\left(i\lambda \right)S_L\overline{\Sigma },\end{array}$$

(60)

respectively. Here,

$$\begin{array}{c}\hfill \overline{\mathbf{u}\phi }=\overline{{\chi }_u{\phi }_u{\mathbf{u}}_u}+\overline{{\chi }_b{\phi }_b{\mathbf{u}}_b}={\gamma }_u{\langle \phi \mathbf{u}\rangle }_u+{\gamma }_b{\langle exp\left(i\lambda \right)\mathbf{u}\rangle }_b.\end{array}$$

(61)

In the considered case, the reaction-progress-variable PDF may be decomposed as follows

$$\begin{array}{c}\hfill P(c,\mathbf{x},t)={\gamma }_u(\mathbf{x},t)P_u(c,\mathbf{x},t)+{\gamma }_b(\mathbf{x},t)\delta (c-1)\end{array}$$

(62)

and the following transport equation

$$\begin{array}{c}\hfill \frac{\partial P}{\partial t}+\nabla ·\left({\langle \mathbf{u}\rangle }_{c}P\right)=-\frac{\partial }{\partial c}\left[{\gamma }_{u}D{\langle \Delta c\rangle }_{u,c}{P}_u\right]\end{array}$$

(63)

can be obtained by applying the inverse Fourier transformation to Equation (58). Here,

$$\begin{array}{c}\hfill {\langle \mathbf{u}\rangle }_{c}P={\langle \mathbf{u}\rangle }_{u,c}{\gamma }_u{P}_u+{\langle \mathbf{u}\rangle }_b{\gamma }_b\delta (c-1),\end{array}$$

(64)

by virtue of Equation (62) and subscript $u,c$ in ${\langle q\rangle }_{u,c}$ designates value of a quantity q conditioned to a particular $c_u$ in the reactants.

Since Equation (63) involves Dirac delta function, it should be considered to be an equation for generalized functions defined when $-\infty <c<\infty$. Such a generalized-function analysis of inert turbulent mixing of substances u and b allowed Kuznetsov [49] to obtain (i) a transport equation for continuous (smooth) PDF $P_u\left(c\right)$, (ii) transport equations for the probabilities ${\gamma }_u$ and ${\gamma }_b$, and (iii) a boundary condition for the PDF $P_u\left(c\right)$ at $c=1$. In the following, this method is applied and briefly described. The reader interested in details is referred to a book by Kuznetsov and Sabelnikov [23].

Let us (i) multiply Equation (63) with an arbitrary good, for example, analytic, function $F\left(c\right)$, (ii) integrate the obtained equation from $-\infty$ to ∞ in the c-space in order to avoid Dirac delta function $\delta (c-1)$, and, then, (iii) substitute Equations (62) and (64) into the integrated equation. Following rules of operations with generalized functions [42], we arrive at

$$\begin{array}{c}\hfill \frac{\partial }{\partial t}\left[{\gamma }_u{\int }_0^1{P}_{u}F\mathrm{d}c+{\gamma }_{b}F\left(1\right)\right]+\nabla ·\left[{\int }_0^1{\langle \mathbf{u}\rangle }_{u,c}{\gamma }_u{P}_{u}F\mathrm{d}c+{\langle \mathbf{u}\rangle }_b{\gamma }_{b}F\left(1\right)\right]\\ \hfill =-{\gamma }_{u}D{\int }_{-\infty }^{\infty }F\frac{\partial }{\partial c}\left({\langle \Delta c\rangle }_{u,c}{P}_u\right)\mathrm{d}c={\gamma }_{u}D{\int }_{-\infty }^{\infty }{\langle \Delta c\rangle }_{u,c}{P}_u\frac{\mathrm{d}F}{\mathrm{d}c}\mathrm{d}c\\ \hfill ={\gamma }_{u}D{\int }_0^1{\langle \Delta c\rangle }_{u,c}{P}_u\frac{\mathrm{d}F}{\mathrm{d}c}\mathrm{d}c\\ \hfill ={\gamma }_{u}D{P}_{u,c=1}F\left(1\right){\langle \Delta c\rangle }_{u,c=1}-{\gamma }_{u}D{\int }_0^{1}F\frac{\partial }{\partial c}\left({\langle \Delta c\rangle }_{u,c}{P}_u\right)\mathrm{d}c.\end{array}$$

(65)

Since the function F is arbitrary, the sum of all terms that involve $F\left(1\right)$ should vanish, as well as the sum of all terms that involve $F\left(c\right)$ within integrals. Accordingly, the following two equations

$$\begin{array}{c}\hfill \frac{\partial {\gamma }_b}{\partial t}+\nabla ·\left({\langle \mathbf{u}\rangle }_b{\gamma }_b\right)={\gamma }_{u}D{P}_{u,c=1}{\langle \Delta c\rangle }_{u,c=1},\end{array}$$

(66)

$$\begin{array}{c}\hfill \frac{\partial {\gamma }_u{P}_u}{\partial t}+\nabla ·\left({\langle \mathbf{u}\rangle }_{u,c}{\gamma }_u{P}_u\right)=-\frac{\partial }{\partial c}\left({\gamma }_{u}D{\langle \Delta c\rangle }_{u,c}{P}_u\right)\end{array}$$

(67)

should hold. Finally, integration of Equation (67) from $c=0$ to $c=1$ yields the following transport equation for the probability ${\gamma }_u$

$$\begin{array}{c}\hfill \frac{\partial {\gamma }_u}{\partial t}+\nabla ·\left({\langle \mathbf{u}\rangle }_u{\gamma }_u\right)=-{\gamma }_{u}D{\langle \Delta c\rangle }_{u,c=1}{P}_{u,c=1},\end{array}$$

(68)

which is written in a form alternative to Equation (24). Equation (68) can also be obtained using Equation (7) for $q=\mathbf{u}$, a constraint of ${\gamma }_u+{\gamma }_b=1$, Equation (66), and the ensemble-averaged continuity equation of $\nabla ·\overline{\mathbf{u}}=0$.

Application of the inverse Fourier transformation to Equations (59) and (60) yields the following generalized-function equation

$$\begin{array}{c}\hfill \frac{\partial P}{\partial t}+\nabla ·\left({\langle \mathbf{u}\rangle }_{c}P\right)\\ \hfill =-\frac{\partial }{\partial c}\nabla ·\left({\gamma }_{u}D{\langle \nabla c\rangle }_{u,c}{P}_u\right)-\frac{{\partial }^2}{\partial c^2}\left({\gamma }_u{\langle N\rangle }_{u,c}{P}_u\right)-S_L\overline{\Sigma }\frac{\mathrm{d}}{\mathrm{d}c}\delta (c-1)\end{array}$$

(69)

and

$$\begin{array}{c}\hfill \frac{\partial P}{\partial t}+\nabla ·\left({\langle \mathbf{u}\rangle }_{c}P\right)\\ \hfill =-D\Delta \left({\gamma }_u{P}_u\right)-\frac{{\partial }^2}{\partial c^2}\left({\gamma }_u{\langle N\rangle }_{u,c}{P}_u\right)-S_L\overline{\Sigma }\frac{\mathrm{d}}{\mathrm{d}c}\delta (c-1),\end{array}$$

(70)

respectively. Adapting the method already used to transform Equations (63)–(69), we arrive at

$$\begin{array}{c}\hfill \frac{\partial }{\partial t}\left[{\gamma }_u{\int }_0^1{P}_{u}F\mathrm{d}c+{\gamma }_{b}F\left(1\right)\right]+\nabla ·\left[{\int }_0^1{\langle \mathbf{u}\rangle }_{u,c}{\gamma }_u{P}_{u}F\mathrm{d}c+{\langle \mathbf{u}\rangle }_b{\gamma }_{b}F\left(1\right)\right]\\ \hfill =\nabla ·\left[{\gamma }_u{P}_{u,c=1}D{\langle \nabla c\rangle }_{u,c=1}F\left(1\right)\right]-\nabla ·{\int }_0^{1}F\frac{\partial }{\partial c}\left({\gamma }_u{P}_{u}D{\langle \nabla c\rangle }_{u,c}\right)\mathrm{d}c\\ \hfill -{\gamma }_u{\langle N\rangle }_{u,c=1}{P}_{u,c=1}{\left.\frac{\mathrm{d}F}{\mathrm{d}c}\right|}_{c=1}+{\left[\frac{\partial }{\partial c}\left({\gamma }_u{\langle N\rangle }_{u,c}{P}_u\right)\right]}_{c=1}F\left(1\right)\\ \hfill -{\int }_0^{1}F\frac{{\partial }^2}{\partial c^2}\left({\gamma }_u{\langle N\rangle }_{u,c}{P}_u\right)\mathrm{d}c+S_L\overline{\Sigma }{\left.\frac{\mathrm{d}F}{\mathrm{d}c}\right|}_{c=1},\end{array}$$

(71)

because

$$\begin{array}{c}\hfill -\nabla ·{\int }_{-\infty }^{\infty }F\frac{\partial }{\partial c}\left({\gamma }_u{P}_{u}D{\langle \nabla c\rangle }_{u,c}\right)\mathrm{d}c=\nabla ·{\int }_0^1{\gamma }_u{P}_{u}D{\langle \nabla c\rangle }_{u,c}\frac{\mathrm{d}F}{\mathrm{d}c}\mathrm{d}c\\ \hfill =\nabla ·\left[{\gamma }_u{P}_{u,c=1}D{\langle \nabla c\rangle }_{u,c=1}F\left(1\right)\right]-\nabla ·{\int }_0^{1}F\frac{\partial }{\partial c}\left({\gamma }_u{P}_{u}D{\langle \nabla c\rangle }_{u,c}\right)\mathrm{d}c,\end{array}$$

(72)

$$\begin{array}{c}\hfill -{\int }_{-\infty }^{\infty }F\frac{{\partial }^2}{\partial c^2}\left({\gamma }_u{\langle N\rangle }_{u,c}{P}_u\right)\mathrm{d}c=-{\int }_0^1{\gamma }_u{\langle N\rangle }_{u,c}{P}_u\frac{{\mathrm{d}}^{2}F}{\mathrm{d}{c}^2}\mathrm{d}c\\ \hfill =-{\gamma }_u{\langle N\rangle }_{u,c=1}{P}_{u,c=1}{\left.\frac{\mathrm{d}F}{\mathrm{d}c}\right|}_{c=1}+{\int }_0^1\frac{\mathrm{d}F}{\mathrm{d}c}\frac{\partial }{\partial c}\left({\gamma }_u{\langle N\rangle }_{u,c}{P}_u\right)\mathrm{d}c\\ \hfill =-{\gamma }_u{\langle N\rangle }_{u,c=1}{P}_{u,c=1}{\left.\frac{\mathrm{d}F}{\mathrm{d}c}\right|}_{c=1}+{\left[\frac{\partial }{\partial c}\left({\gamma }_u{\langle N\rangle }_{u,c}{P}_u\right)\right]}_{c=1}F\left(1\right)\\ \hfill -{\int }_0^{1}F\frac{{\partial }^2}{\partial c^2}\left({\gamma }_u{\langle N\rangle }_{u,c}{P}_u\right)\mathrm{d}c,\end{array}$$

(73)

and

$$\begin{array}{c}\hfill -S_L\overline{\Sigma }{\int }_{-\infty }^{\infty }F\frac{\mathrm{d}}{\mathrm{d}c}\delta (c-1)\mathrm{d}c=S_L\overline{\Sigma }{\int }_{-\infty }^{\infty }\delta (c-1)\frac{\mathrm{d}F}{\mathrm{d}c}\mathrm{d}c=S_L\overline{\Sigma }{\left.\frac{\mathrm{d}F}{\mathrm{d}c}\right|}_{c=1}.\end{array}$$

(74)

Again, the sum of all terms that involve ${(\mathrm{d}F/\mathrm{d}c)}_{c=1}$ should vanish, the sum of all terms that involve $F\left(1\right)$ should also vanish, as well as the sum of all terms that involve $F\left(c\right)$ within integrals. Accordingly, the following three equations

$$\begin{array}{c}\hfill {\gamma }_u{\langle N\rangle }_{u,c=1}{P}_{u,c=1}=S_L\overline{\Sigma },\end{array}$$

(75)

$$\begin{array}{c}\hfill \frac{\partial {\gamma }_b}{\partial t}+\nabla ·\left({\langle \mathbf{u}\rangle }_b{\gamma }_b\right)=D\nabla ·\left({\gamma }_u{P}_{u,c=1}{\langle \nabla c\rangle }_{u,c=1}\right)+{\left[\frac{\partial }{\partial c}\left({\gamma }_u{\langle N\rangle }_{u,c}{P}_u\right)\right]}_{c=1},\end{array}$$

(76)

and

$$\begin{array}{c}\hfill \frac{\partial {\gamma }_u{P}_u}{\partial t}+\nabla ·\left({\langle \mathbf{u}\rangle }_{u,c}{\gamma }_u{P}_u\right)=-D\frac{\partial }{\partial c}\nabla ·\left({\gamma }_u{P}_u{\langle \nabla c\rangle }_{u,c}\right)-{\gamma }_u\frac{{\partial }^2}{\partial c^2}\left({\langle N\rangle }_{u,c}{P}_u\right)\end{array}$$

(77)

should hold. Moreover, the following equation

$$\begin{array}{c}\hfill \frac{\partial {\gamma }_u}{\partial t}+\nabla ·\left({\langle \mathbf{u}\rangle }_u{\gamma }_u\right)=-D\nabla ·\left({\gamma }_u{P}_{u,c=1}{\langle \nabla c\rangle }_{u,c=1}\right)-{\left[\frac{\partial }{\partial c}\left({\gamma }_u{\langle N\rangle }_{u,c}{P}_u\right)\right]}_{c=1}\end{array}$$

(78)

results from Equation (76) using Equation (7) for $q=\mathbf{u}$, a constraint of ${\gamma }_u+{\gamma }_b=1$, and the ensemble-averaged continuity equation of $\nabla ·\overline{\mathbf{u}}=0$.

Equations (3) and (75) result in

$$\begin{array}{c}\hfill {\gamma }_u{P}_{u,c=1}={\delta }_L\overline{\Sigma }.\end{array}$$

(79)

Equation (76) can be rewritten in another form

$$\begin{array}{c}\hfill \frac{\partial {\gamma }_b}{\partial t}+\nabla ·\left({\langle \mathbf{u}\rangle }_b{\gamma }_b\right)=-S_L\nabla ·\left({\gamma }_u{P}_{u,c=1}{\langle \mathbf{n}\rangle }_{u,c=1}\right)+{\left[\frac{\partial }{\partial c}\left({\gamma }_u{\langle N\rangle }_{u,c}{P}_u\right)\right]}_{c=1},\end{array}$$

(80)

because

$$\begin{array}{c}\hfill {\langle \nabla c\rangle }_{u,c=1}=-\langle |\nabla c_u{|\mathbf{n}\rangle }_{u,c=1}=-\frac{1}{{\delta }_L}{\langle \mathbf{n}\rangle }_{u,c=1}.\end{array}$$

(81)

Comparison of (24) and (80) results in

$$\begin{array}{c}\hfill -S_L\nabla ·\left({\gamma }_u{P}_{u,c=1}{\langle \mathbf{n}\rangle }_{u,c=1}\right)+{\left[\frac{\partial }{\partial c}\left({\langle N\rangle }_{u,c}{P}_u\right)\right]}_{c=1}={\langle S_d\rangle }_r\overline{\Sigma }.\end{array}$$

(82)

It is of interest to note that integration of (79) yields

$$\begin{array}{c}\hfill U_T={\int }_{-\infty }^{\infty }{S}_L\overline{\Sigma }\mathrm{d}x=\frac{1}{{\tau }_L}{\int }_{-\infty }^{\infty }{\gamma }_u{P}_{u,c=1}\mathrm{d}x=\frac{{\delta }_T}{{\tau }_L}{\int }_{-\infty }^{\infty }{\gamma }_u{P}_{u,c=1}\mathrm{d}\zeta ,\end{array}$$

(83)

where $\zeta =x/{\delta }_T$ is normalized distance. Consequently, using (1) and the following scaling

$$\begin{array}{c}\hfill {\delta }_T\propto {\delta }_L{\left(ScR{e}_t\right)}^{1/2}\propto LD{a}^{-1/2}\end{array}$$

(84)

obtained for the considered combustion regime earlier [19], we arrive at

$$\begin{array}{c}\hfill \frac{Da}{{\tau }_T}\propto \frac{1}{{\tau }_L}{\int }_{-\infty }^{\infty }{\gamma }_u{P}_{u,c=1}\mathrm{d}\zeta ,\end{array}$$

(85)

that is, by virtue of the definition of $Da={\tau }_T/{\tau }_L$, the integral on the RHS of Equation (85) does not depend on $Da$, $Ka$, or $R{e}_t$ in the considered regime of turbulent wave propagation.

Finally, application of the method adapted above to Equation (70) instead of Equation (69) results in

$$\begin{array}{c}\hfill \frac{\partial {\gamma }_b}{\partial t}+\nabla ·\left({\langle \mathbf{u}\rangle }_b{\gamma }_b\right)=D\Delta {\gamma }_b+{\gamma }_u{\left[\frac{\partial }{\partial c}\left({\langle N\rangle }_{u,c}{P}_u\right)\right]}_{c=1},\end{array}$$

(86)

$$\begin{array}{c}\hfill \frac{\partial {\gamma }_u{P}_u}{\partial t}+\nabla ·\left({\langle \mathbf{u}\rangle }_{u,c}{\gamma }_u{P}_u\right)=D\Delta \left({\gamma }_u{P}_u\right)-{\gamma }_u\frac{{\partial }^2}{\partial c^2}\left({\langle N\rangle }_{u,c}{P}_u\right),\end{array}$$

(87)

and

$$\begin{array}{c}\hfill \frac{\partial {\gamma }_u}{\partial t}+\nabla ·\left({\langle \mathbf{u}\rangle }_u{\gamma }_u\right)=D\Delta {\gamma }_u-{\gamma }_u{\left[\frac{\partial }{\partial c}\left({\langle N\rangle }_{u,c}{P}_u\right)\right]}_{c=1}.\end{array}$$

(88)

5. Conclusions

A new mathematical framework for describing thin-reaction-zone regime of propagation of a single-reaction wave in a constant-density turbulent flow at low Damköhler numbers is developed. In this regime, the molecular mixing plays a vital role, because it smooths out small-scale wrinkles of reaction-zone surface, generated by turbulent eddies, and, therefore, significantly reduces the reaction-zone surface area and turbulent consumption velocity. This regime differs fundamentally from the classical problem of front propagation in a turbulent medium at large Damköhler numbers, because the molecular mixing is of the secondary importance and reaction-progress-variable PDF $P\left(c\right)$ involves two Dirac delta functions in the latter case. On the contrary, in the studied case of a low $Da$, the PDF $P\left(c\right)$ contains a single Dirac delta function and is formally similar to the scalar PDF in the case of passive scalar mixing in the fully developed region of a turbulent jet.

The major results of the present work are as follows.

First, exact, unclosed transport equations for the unconditioned, see Equations (35) and (36), and conditioned (to $c<1$), see Equations (26)–(29), first and second-order moments of the reaction-progress-variable moments are derived.

Second, exact, unclosed transport equations for the unconditioned, see Equations (59), (65) and (66), and conditioned, see Equations (62)–(64), (72)–(74) and (82)–(84)), reaction-progress-variable PDFs are derived.

Third, for the conditioned PDF, the boundary condition at the reaction sheet is derived in two different forms, see Equations (71) and (75).

Forth, a constraint on the PDF conditioned to the reactant side of the reaction sheet and integrated over space is established, see Equation (81).

It should be stressed that the closure problem for the derived equations is beyond the scope of the present study.

Finally, when applying the obtained results to more complicated problems such as premixed turbulent combustion, it is worth remembering important effects, for example, thermal expansion, preferential diffusion, complex chemistry, or/and heat losses, that vanish within the framework of the simplified problem addressed in the present work. As far as the thermal expansion effects are concerned, the developed mathematical framework can easily be extended to allow for density variations.