The transported probability density function (PDF) method is an attractive model for the closure of turbulent mixing and turbulent reactive flows. The transported PDF method suffers from the curse of dimensionality and an efficient numerical implementation of the method has always been a research topic of great importance. The Eulerian Monte Carlo field (EMCF) method (also termed as the stochastic field method) has been proposed as an efficient solution approach for solving the PDF transport equation for two decades. In this work, we revisit the EMCF method and examine its mathematical consistency analytically and numerically. It is found that the EMCF method is not mathematically consistent with the PDF transport equation that the method intends to solve. This creates a serious inconsistency issue and causes uncertainties in the yielded numerical solutions by EMCF. It is imperative to examine the effect of the inconsistency. We evaluate the effect of the inconsistency in a simplified turbulent mixing layer test case. Corrections to remedy the mathematical inconsistency are proposed and examined. The effectiveness of the corrections is demonstrated numerically through convergence testing. The effect of the Reynolds number on the inconsistency is explored. The impact of the inconsistency is also investigated in a thermal wake behind a line source in grid turbulence to assess the importance of the issue in real turbulence problems.

Turbulent reactive flows exist in many energy conversion systems such as furnaces, industrial burners, gasoline engines, gas turbines, and rocket engines. Modeling turbulent mixing and combustion in these systems is quite challenging due to the multi-scale and highly nonlinear nature of the problem. The transported probability density function (PDF) method1–3 is an attractive model for the closure of turbulent mixing and turbulence-chemistry interactions in these problems. The PDF method has been used extensively for the modeling of turbulent reactive flows to study local extinction and re-ignition,4,5 auto-ignition and flame lift-off,6–8 soot formation,9 NOx formation,5,10 etc. The method has also been used in other areas such as cavitating flows,11 atmospheric dispersion,12 groundwater,13 heterogeneous porous media flows,14 and viscoelastic flows.15 

The PDF method suffers from the curse of dimensionality. For real engineering problems involving multiple-scalar transport such as a turbulent combustion system, there can be up to hundreds or even thousands of scalars, resulting in a very high-dimensional turbulent mixing problem for the PDF method. This makes the traditional numerical solution approaches based on the finite-volume or finite-difference methods fail quickly. Different solution approaches are needed for efficiently implementing the PDF method in computational simulations. So far, there are four different types of approaches available for solving the PDF method equations: the Eulerian Monte Carlo particle (EMCP) method16 (also called node-based particle method), the Lagrangian Monte Carlo particle (LMCP) method,1 the Eulerian Monte Carlo field (EMCF) method 17–19 (also commonly termed as the stochastic field method), and the multi-environment PDF with the direct quadrature method of moments (DQMOM) closure.20 The EMCP method was devised by Pope in 1981,16 based on a Monte Carlo method with the particles placed on fixed grid nodes. The EMCP method has been mostly abandoned probably because of its limitation on having only first-order temporal accuracy in numerical integration. In 1985, Pope1 invented the LMCP method which has subsequently become one of the most popular solution approaches for the PDF method. In LMCP, Lagrangian Monte Carlo particles are used as a discrete representation of a PDF function. The mathematical consistency and convergence of the LMCP method have been examined extensively in the past (e.g., Refs. 1 and 21–23). The LMCP method has been used widely in practical simulations.4–7,24,25 The EMCF method17–19 was introduced in the late 1990s as an attractive alternative for solving the PDF transport equation. The method introduces a concept called stochastic fields to discretely represent a PDF function. Although the EMCF method has been used in practical simulations,8,9,11,12,26–28 the examination of the mathematical consistency and convergence of the EMCF method is less comprehensive and the topic is the main focus of this work. Another alternative approach for solving the PDF method is the DQMOM method20 which has also been used in practical simulations.29–32 There have also been nice efforts on assessing the performance difference of the different PDF solution approaches. Akroyd et al.33 compared the EMCF method and the DQMOM method in a constant density tubular reactor. Jaishree and Haworth34 compared the performance of three methods (LMCP, EMCF, and DQMOM) in turbulent non-premixed jet flames with moderate-to-strong turbulence-chemistry interactions. In this work, we will focus only on a theoretical aspect of the EMCF method, i.e., the mathematical consistency of the method.

For the EMCF method, there exist two different formulations. One is based on the Itô stochastic differential equation (SDE) by Valiño17 and the other based on the Stratonovich SDE by Sabel’nikov and Soulard.18 Both forms are statistically equivalent.18 This work uses the Itô form for the discussion of the EMCF method. Less work has been devoted to examining the mathematical consistency and convergence of the EMCF method although the method has existed for two decades. In other words, the mathematical consistency has not been fully examined and verified between the numerical solutions obtained from the EMCF method and the PDF transport equation that the EMCF method intends to solve. The mathematical consistency and convergence of EMCF must be thoroughly examined and verified or quantified if an inconsistency exists before it can be used in model predictions for turbulent mixing and combustion reliably. This work is motivated to examine this mathematical aspect of the EMCF method. Over the years, there have emerged different versions of the EMCF method (in the Itô form) in the literature. Valiño17 developed the original version of the EMCF method. In 2016, Valiño et al.19 pointed out a spurious Wiener term associated with the molecular diffusion in the original EMCF method and introduced a modified EMCF method to remedy the spurious error of the original EMCF. We found, however, that, even with the modification by Valiño et al.,19 the EMCF method is still not fully consistent with the PDF transport equation. The inconsistency issues in the EMCF methods motivated this work, and the objectives of the work are as follows:

  1. To present a mathematical approach to thoroughly examine the mathematical consistency and convergence of the different EMCF methods,

  2. to develop new EMCF methods with a better mathematical consistency,

  3. to conduct numerical convergence tests for the examination and quantification of the mathematical consistency of the EMCF methods, and

  4. to provide an assessment of the impact of the inconsistency of the EMCF methods under experimental turbulence conditions.

The rest of the paper is organized as follows. Section II provides an overview of the transported PDF method and defines criteria for examining the mathematical consistency. Section III reviews the existing EMCF methods and presents a mathematical approach to examine the mathematical consistency of the EMCF methods for solving the transported PDF equation. New EMCF methods are introduced in Sec. III to improve the mathematical consistency. Numerical convergence tests are conducted in an opposed-jet turbulent mixing layer in Sec. IV to confirm the inconsistency issues and to examine their effects. In Sec. V, the performance of the EMCF methods is further examined in an experimental line source dispersion in grid turbulence. Conclusions are drawn in Sec. VI.

We consider a single-scalar field ϕx,t transported by a solenoidal turbulent flow field u = {u1, u2, u3}, ∇ · u = 0. For variable density problems, the discussion is similar, and for simplicity we consider only constant density problems in this work. The conservation equation for ϕ(x, t) is written as

ϕ(x,t)t+uϕ(x,t)=Γϕ(x,t)+S(ϕ),
(1)

where x is the position vector, t is time, Γ is the molecular diffusivity (assumed to be constant for the following discussions), and S(ϕ) is a source term. Treating ϕ(x, t) as a random variable, we can define a single-point single-time PDF of ϕ(x, t) as fϕ(ψ; x, t), where ψ is the sample space variable of ϕ(x, t). [The PDF discussed in this work is placed in the context of Reynolds-averaged Navier-Stokes (RANS).35 In large-eddy simulations (LES),35 the PDF is sometimes termed as the filtered-density function (FDF),36 and the issues discussed are applicable to LES as well.] The PDF fϕ(ψ; x, t) can be expressed as the mathematical expectation of a Dirac delta function1 δ(ψϕ(x, t)),

fϕ(ψ;x,t)=δ(ψϕ(x,t)).
(2)

By using the standard technique of Pope,1 we can derive the transport equation for fϕ(ψ; x, t) as

fϕ(ψ;x,t)t+u|ϕ=ψfϕ(ψ;x,t)  =Γfϕ(ψ;x,t)2ψ2    ×Γϕϕϕ=ψfϕ(ψ;x,t)    S(ψ)fϕ(ψ;x,t)ψ,
(3)

where ⟨·|ϕ = ψ⟩ denotes conditional mathematical expectations. The following equations and models can be used to close the PDF transport equation (3):

u|ϕ=ψ=u+u|ϕ=ψ,
(4)
u|ϕ=ψfϕ(ψ;x,t)=Γtfϕ(ψ;x,t),
(5)
ϕϕϕ=ψ=ϕϕ+ϕϕϕ=ψ,
(6)
ψΓϕϕϕ=ψfϕ(ψ;x,t)=ωψϕfϕ(ψ;x,t),
(7)

in which ⟨·⟩ denotes Reynolds averaging, the prime “′” denotes fluctuations, Γt is the turbulent eddy diffusivity, ω is the scalar mixing frequency, a gradient-diffusion transport model35 is used in Eq. (5), and the interaction by exchange with the mean (IEM) mixing model37 is used to close the transport of PDF caused by the conditional scalar dissipation in Eq. (7). After substituting Eqs. (4)–(7) into Eq. (3), we obtained the following modeled PDF equation in a closed form:

fϕ(ψ;x,t)t+ufϕ(ψ;x,t)   =Γt+Γfϕ(ψ;x,t)     2ψ2Γϕϕfϕ(ψ;x,t)     +ωψϕfϕ(ψ;x,t)ψS(ψ)fϕ(ψ;x,t)ψ.
(8)

With proper initial and boundary conditions, the modeled PDF Eq. (8) can be solved mathematically to obtain the solution of the PDF fϕ(ψ; x, t). An analytical solution for fϕ(ψ; x, t) can hardly be found in practical problems and a numerical solution approach is always needed. The EMCF methods17–19 as the solution approach for solving Eq. (8) numerically are the main focus of this paper. Particularly, we will examine in detail the mathematical and numerical consistency between the modeled PDF Eq. (8) and the solutions obtained from the EMCF methods. The details of the EMCF methods will be reviewed in Sec. III.

To examine the mathematical consistency of the EMCF methods, we require a set of criteria to assess the issue. For the PDF method, it is difficult to directly examine the mathematical consistency of the solution of the PDF itself. Instead a weak approach21,23 is often used by examining the mathematical consistency of the moments of scalars obtained from the EMCF methods. The mth central moments (m = 1, 2, …, ) of the scalar ϕm is defined as

ϕm=ψmfϕ(ψ;x,t)dψ.
(9)

Following  Appendix A, we can derive the following moment equation from the modeled PDF Eq. (8):

ϕmt+uϕm   =(Γt+Γ)ϕmm(m1)Γϕϕϕm2     mωϕmϕϕm1+mϕm1S(ϕ).
(10)

The yielded moment equation by the EMCF methods needs to reproduce the moment equation (10) to be fully consistent with the modeled PDF Eq. (8). The mathematical consistency of the EMCF methods through the moment equation (10) will be examined in detail in Sec. III. Practically, it is infeasible to examine an infinite number of moments in Eq. (10), and in this work we mainly focus on the first two moments [m = 1 and 2 in Eq. (10)], which are written as

ϕt+uϕ=(Γt+Γ)ϕ+S(ϕ),
(11)
ϕ2t+uϕ2=(Γt+Γ)ϕ22Γϕϕ2ωϕ2ϕϕ+2ϕS(ϕ).
(12)

The transport equation of the scalar variance define as ϕ2=ϕ2ϕ2 can be obtained from Eqs. (11) and (12) as follows:

ϕ2t+uϕ2=(Γt+Γ)ϕ2+2Γtϕϕ2ωϕ2+2ϕS,
(13)

where ϕS=ϕSϕS. The mean equation (11) and the variance equation (13) will be used as the mathematical consistency criteria for the discussion of the EMCF methods in Sec. III.

In the EMCF methods, the concept of Monte Carlo fields (or stochastic fields) is introduced to provide a discrete representation of the PDF fϕ(ψ; x, t) to be solved. The evolution of the fields is governed by a set of stochastic partial differential equations (SPDEs). The solutions of these SPDEs approximate the solution of the modeled PDF transport equation (8). Existing EMCF methods are reviewed below with their mathematical consistency thoroughly examined. New versions of EMCF are developed in Sec. III C to improve the mathematical consistency.

The EMCF methods define a Monte Carlo field ϕ*(x, t) with its sample space variable ψ*. The PDF of ϕ* is defined as, similar to Eq. (2),

fϕ*(ψ*;x,t)=δ(ψ*ϕ*(x,t)).
(14)

The PDF fϕ*(ψ*;x,t) can be approximated by Nf discrete Monte Carlo fields ϕnf(x,t), (nf = 1, …, Nf),

fϕ*(ψ*;x,t)1Nfnf=1Nfδψ*ϕnf(x,t).
(15)

When Nf, the approximated PDF in Eq. (15) becomes the true PDF in Eq. (14). The main question addressed in this paper is to find the mathematical equivalence between fϕ*(ψ*;x,t) obtained from the EMCF methods in Eq. (15) when Nf and fϕ(ψ; x, t) governed by Eq. (8).

In the original EMCF method by Valiño,17 the SPDEs for the Nf Monte Carlo fields ϕnf(x,t) are written as

dϕnf(x,t)=uϕnf(x,t)dt+(Γt+Γ)ϕnf(x,t)dt+(2Γt+2Γ)1/2ϕnf(x,t)dWnf(t)ωϕnf(x,t)ϕ*(x,t)Nfdt+Sϕnf(x,t)dt,
(16)

where Wnf(t) is a vector-valued Wiener process with dWnf=0 and dWnfdWnf=Idt and I is a unit matrix. We denote this version of EMCF as the EMCF-O method.

In order to examine the mathematical consistency between EMCF-O in Eq. (16) and the modeled PDF transport equation (8), we need to derive the equivalent moment transport equation from Eq. (16) which can then be compared directly with the moment equation (10) derived from the PDF Eq. (8) for consistency checking. Deriving the moment equation from Eq. (16) is not trivial and has not been reported in any previous studies. We present a new approach for this derivation in  Appendix B and the derived moment equation is

ϕ*mt+uϕ*m=(Γt+Γ)ϕ*mmωϕ*mϕ*ϕ*m1+mϕ*m1Sϕ*.
(17)

With the moments equation (17) obtained from the EMCF-O method, we can readily examine the mathematical consistency of the method. By directly comparing Eq. (17) with the moment equation (10) from the modeled PDF transport equation in Sec. II, we can immediately see the inconsistency of the EMCF-O method since the two moment equations are not in an identical form. In other words, the solutions from the EMCF-O method do not converge to the true solution admitted by the modeled PDF transport equation (8) even when dt → 0 (dt can be viewed as a finite time step size used in a numerical simulation) and Nf. The inconsistency of the EMCF-O method can be further examined by focusing on the first two moments (the mean and variance) by setting m = 1 or 2 in Eq. (17). The yielded mean and variance transport equations are as follows:

ϕ*t+uϕ*=(Γt+Γ)ϕ*+Sϕ*,
(18)
ϕ*2t+uϕ*2   =(Γt+Γ)ϕ*2+2Γtϕ*ϕ*     +2Γϕ*ϕ*T62ωϕ*2+2ϕ*S.
(19)

By comparing the transport equations for the mean ⟨ϕ⟩ in (11) and for ϕ* in (18), we can see that they are in an identical form and hence yield identical solutions if the same initial and boundary conditions are imposed. This indicates that the EMCF-O method is indeed able to yield a consistent solution for the scalar mean. For the scalar variance, by comparing Eq. (19) with Eq. (13), we see an extra term T6 in (19) in the EMCF-O method, which indicates that the EMCF-O method cannot yield a consistent solution to the scalar variance. The term T6 in Eq. (19) is non-negative and is caused by molecular diffusion to produce a spurious variance production. Since it is caused by molecular diffusion, the significance of the term T6 can be questioned generally in the RANS modeling of high-Re turbulence problems, and the importance can be argued only in low-Re number turbulence (such as wall turbulence) and near laminar-flow regions.19 However, in LES of turbulence, even for high-Re problems, the significance of the term T6 can still be dramatic because the magnitude of T6 can be comparable to the variance production term [the second term on the right-hand side of Eq. (19)] since Γ and Γt can be in a similar magnitude.38 In general, the mathematical inconsistency of the EMCF-O method appears in all the moments with m ≥ 2.

Recently, Valiño et al.19 recognized the inconsistency of the EMCF-O method and proposed a modified version of the EMCF method to remedy the issue. Although Valiño et al.19 has reported the same inconsistency issue of EMCF-O, their approach is very different from ours. The approach we developed in Sec. III A is more general and can be used to get to the root of the mathematical inconsistency of all EMCF methods including the modified EMCF method. The modified EMCF method is written as19 

dϕnf(x,t)=uϕnf(x,t)dt+(Γt+Γ)ϕnf(x,t)dt+(2Γt)1/2ϕnf(x,t)dWnf(t)ωϕnf(x,t)ϕ*(x,t)Nfdt+Sϕnf(x,t)dt.
(20)

We denote this version of the EMCF method as the EMCF-M method. In comparison with the EMCF-O method in Eq. (16), Eq. (20) has a slightly different form for the Wiener term with the molecular diffusivity removed. This was intended to eliminate the spurious variance production in the EMCF-O method. The full consistency of the EMCF-M method with the modeled PDF transport equation (8), however, has not been examined and will be done below.

Following the same mathematical approach in Sec. III A, we can readily derive the transport equation for the moments ϕ*m from Eq. (20) as

ϕ*mt+uϕ*m   =(Γt+Γ)ϕ*mm(m1)Γϕ*ϕ*ϕ*m2     mωϕ*mϕ*ϕ*m1+mϕ*m1Sϕ*.
(21)

Comparing this moment equation with Eq. (10) from the PDF method, we see immediately that the two equations are also not in an identical form. In particular, the second terms on the right-hand side of both equations are different. This indicates evidently that, even with the modification made by Valiño et al.,19 the EMCF-M method is still not mathematically consistent with the PDF equation. To further examine this inconsistency, we next focus on the first two moments. By setting m = 1 in Eq. (21), we obtain the mean scalar equation which is the same as Eq. (18) [also the same as Eq. (11) from the modeled PDF equation]. This indicates that the EMCF-M method is also consistent for yielding the correct transport of the scalar mean, just like the EMCF-O method. This observation is not surprising since the Wiener term in Eq. (20) does not affect the scalar mean. The scalar variance equation by the EMCF-M method can be derived from Eq. (21) as

ϕ*2t+uϕ*2   =(Γt+Γ)ϕ*2+2Γtϕ*ϕ*     ×2Γϕ*ϕ*T72ωϕ*2+2ϕ*S,
(22)

where ϕ*ϕ*=ϕ*ϕ*ϕ*ϕ*. The modification in the EMCF-M method is indeed effective to eliminate the spurious variance production term [T6 in Eq. (19)]. However, the modification also adversely introduces the term T7 in Eq. (22) which does not exist in the variance transport equation (13) from the modeled PDF equation, making the EMCF-M method also mathematically inconsistent for the scalar variance transport. The term T7 is non-positive, and hence it introduces a spurious dissipation effect to the variance transport, as opposed to the spurious production in the EMCF-O method. As a result, it is expected that the EMCF-M method tends to under-predict the variance because of the spurious dissipation, while the EMCF-O method tends to over-predict the variance. Theoretically, the spurious dissipation T7 represents the numerical dissipation resolved by the smooth Monte Carlo fields. In practical simulations, the computational grid used is unable to resolve all the turbulent fluctuations and hence the smooth Monte Carlo fields can only resolve little scalar gradient. For low Re turbulence, the physical scalar dissipation [the fourth term on the right-hand side of Eq. (22)] weakens which makes the effect of the spurious variance dissipation T7 in Eq. (22) more evident. In the limit of Re, the effect of the spurious variance dissipation is expected to vanish. The inconsistency in the EMCF-M method also appears in all moments with m ≥ 2, just like the EMCF-O method.

A new mathematical approach is presented to rigorously examine the mathematical consistency of the different versions of the EMCF methods. Both versions are found to be inconsistent with the modeled PDF transport equation. The effect of the inconsistency observed in EMCF-O and EMCF-M has not been evaluated in actual simulations. A comparative assessment of the two methods is also missing to provide a guideline for model selection in future studies. These will be done in Secs. IV and V.

The developed mathematical approach in Sec. III A can also be used to explore the possibility of developing fully consistent EMCF methods. Unfortunately, based on our analysis, it appears to be very difficult if not impossible to develop fully consistent EMCF methods for solving the modeled PDF transport equation (8). Thus, the inconsistency discussed above inevitably becomes part of the sources of errors in the EMCF methods, and it is important to have the inconsistency errors thoroughly examined and quantified. As long as the inconsistency error is dominated by other errors such as the physical model error involved in the IEM mixing, the EMCF method is still a valuable solution approach for efficiently solving the modeled PDF transport equation. On the other hand, if we can somewhat minimize the inconsistency of the EMCF method, that will also be valuable. This minimization of the inconsistency is attempted here.

We consider a corrected EMCF method written as follows:

dϕnf(x,t)=uϕnf(x,t)dt+(Γt+Γ)ϕnf(x,t)dt+(2Γt+2Γc)1/2ϕnf(x,t)dWnf(t)(ω+ωc)ϕnf(x,t)ϕ*(x,t)Nfdt+Sϕnf(x,t)dt,
(23)

in which two new parameters, Γc and ωc, are added to minimize the inconsistency in the EMCF methods. We choose to add a correction to the Wiener term [the third term on the right-hand side of Eq. (23)] and the mixing term (the fourth term on the right-hand side) to ensure that the field mean is not affected since EMCF-O and EMCF-M can yield the field mean correctly. A criterion for the inconsistency minimization will be established below, which will be used to specify the parameters Γc and ωc. Following the same mathematical approach for deriving the moment equation in Sec. III A, we can derive the moment equation from the corrected EMCF method in Eq. (23) as follows:

ϕ*mt+uϕ*m   =(Γt+Γ)ϕ*mm(m1)(ΓΓc)     ×ϕ*ϕ*ϕ*m2m(ω+ωc)     ×ϕ*mϕ*ϕ*m1+mϕ*m1Sϕ*.
(24)

For this moment equation to be identical to Eq. (10) obtained directly from the modeled PDF Eq. (8), we find the following constraint equation for Γc and ωc:

m(m1)(ΓΓc)ϕ*ϕ*ϕ*m2   mωcϕ*mϕ*ϕ*m1      =m(m1)Γϕ*ϕ*ϕ*m2,
(25)

for any value of m. This equation is over-constrained for determining the two model parameters Γc and ωc. It is noted that the constraint equation holds when m = 1 (for the scalar mean) for any choice of Γc and ωc. If we just correct the inconsistency for the variance (m = 2), there is an infinite number of specifications for Γc and ωc. Two example (and probably the simplest) specifications are shown as follows:

Γc=Γωc=Γϕ*ϕ*ϕ*2ϕ*2,
(26)

or

Γc=Γ1ϕ*ϕ*ϕ*ϕ*ωc=0.
(27)

We limit our discussions on the corrected EMCF method for variance to the above two specifications. The consistency of the variance transport is chosen for the criterion of minimizing the inconsistency of the EMCF method. This cannot correct the inconsistency for high-order moments, but achieving consistency for low-order moments is valuable since they are of most practical interest. We can additionally use the consistency of the third-order moment (m = 3) to provide a unique specification of the two correction parameters. By solving Eq. (25) with m = 2 and 3, we can readily obtain

Γc=Γϕ*ϕ*ϕ*ϕ*ϕ*3ϕ*ϕ*2Γϕ*ϕ*2ϕ*ϕ*2ϕ*2ϕ*2ϕ*ϕ*ϕ*3ϕ*ϕ*2ϕ*2ϕ*2ϕ*ϕ*2ωc=Γϕ*ϕ*2ϕ*ϕ*2ϕ*ϕ*Γϕ*ϕ*ϕ*ϕ*ϕ*ϕ*2ϕ*ϕ*ϕ*3ϕ*ϕ*2ϕ*2ϕ*2ϕ*ϕ*2,
(28)

which appears very complicated. For easy discussion of the corrected EMCF methods, we name them the EMCF-C1 method in Eq. (23) with Eq. (26), the EMCF-C2 method with Eq. (27), and the EMCF-C3 method with Eq. (28). We choose to focus on the first two moments following general practice for turbulence modeling, and hence will consider only EMCF-C1 and EMCF-C2 in the later discussions. This correction strategy is mathematical in nature and it does not yield mathematical consistency for all the moments (or the PDF). Nevertheless, it is still useful to have versions of EMCF that are superior over the existing ones with mathematical consistency for the transport of scalar variance, especially in the context where the scalar mean and variance are of most interest as in many engineering problems. In the current work, the model correction is biased toward lower-order moments. Errors are still present in the high-order moments and the probability distribution which are responsible for important physics such as turbulent intermittency and low-probability events like localized flame extinction in combustion. Theoretically, it is possible to introduce more correction parameters into Eq. (23) to match more higher-order moments, but it will add significant complexity to the EMCF method and hence is not attempted here. The EMCF-C1 and EMCF-C2 models discussed above are for constant density flow problems. Extension of these corrected EMCF models to variable density and reactive flow problems is straightforward by incorporating density-weighted averaging. The corrected models in Eqs. (26) and (27) can cause numerical difficulties when the denominators become zero, i.e., the variance is zero or the fields are uniform. In the numerical implementation of the models, we simply add a tiny positive number on the order of 10−15 to the denominators to avoid the numerical difficulties.

In summary, we presented a mathematical approach to examine the mathematical consistency of the different versions of the EMCF methods. The two existing versions, the original method (EMCF-O) and the modified method (EMCF-M), are found to be not fully consistent with the modeled PDF transport equation. This finding is important to enable a full understanding of the existing EMCF methods. It also warrants a need for the quantification of the effect of the inconsistency and further improvement to the method. A few corrected EMCF methods (EMCF-C1, EMCF-C2, and EMCF-3) are proposed to partially eliminate the mathematical inconsistency in the existing EMCF methods. The impact of the mathematical inconsistency in actual simulations and the effectiveness of the consistency correction will be examined in Sec. IV.

In this section, we conduct a numerical convergence study of the different EMCF methods discussed in Sec. III to understand and to quantify the effect of the mathematical inconsistency of the EMCF methods. An opposed-jet turbulent mixing layer test case is employed as the numerical convergence test case. The purpose of the numerical convergence tests is to confirm the mathematical inconsistency of the EMCF methods observed in Sec. III as well as to examine the effect of the inconsistency. The effect of different Re numbers on the model inconsistency will also be evaluated.

A highly simplified opposed-jet turbulent mixing layer test case, illustrated in Fig. 1, is employed for the convergence study of the EMCF methods. Homogeneous turbulence is assumed in the whole mixing field, and the mean flow is assumed to follow a potential flow, ū(x) = −asx, where as is the mean strain rate. The statistics of scalars (or the PDF of scalars) are assumed to be independent of the radial direction, following the assumption commonly used in laminar opposed-jet flame simulations.39 The validity of these assumptions is not important since the purpose of the opposed-jet turbulent mixing layer test case is to have a simplified problem for examining the different numerical approaches. With the above assumptions, the mixing layer problem becomes statistically one dimensional along the axis so that the modeled PDF transport equation (8) and the moment equations (11) and (13) are highly simplified for the testing. A constant-density flow with a conserved scalar [the source term S(ϕ) = 0] is assumed.

FIG. 1.

Sketch of an opposed-jet turbulent mixing layer.

FIG. 1.

Sketch of an opposed-jet turbulent mixing layer.

Close modal

In the mixing layer simulations, a single scalar ϕ is considered and its statistics are a function of x only. The one-dimensional computational domain is set to be L,L with L = 0.04 m, the strain rate as = 1.0 s−1, the kinematic viscosity ν = 1.65 × 10−5 m2/s, the turbulent kinetic energy (homogeneous) k = 5 × 10−4 m2/s2, the molecular Schmidt number σ = 0.7, and the turbulent Schmidt number σt = 0.9. The turbulent Reynolds number based on the fluctuating velocity scale and the integral length scale40 varies over a range between [20, 400] in the discussions below. With these parameters, the turbulent kinetic energy dissipation rate is computed40 as ε = 2k2/(3νRe), the molecular diffusivity Γ = ν/σ, and the turbulent diffusivity Γt = νt/σt. The turbulent viscosity νt is modeled following the kε turbulence model35 as νt = Cμk2/ε with Cμ = 0.09. The mixing frequency ω is commonly modeled1 as ω = cϕε/(2k) with cϕ = 2.0.

The governing equations for the scalar mean ⟨ϕ⟩ in (11), the variance ϕ2 in (13), and the Monte Carlo field equations in (16), (20), and (23) are in a similar format and can be solved by using standard finite volume (FV) methods.41 A second-order upwind scheme is used for discretizing the convection terms, and central difference schemes are used for all the other gradient terms. A first-order Euler scheme is used to advance the equations in time until reaching the statistically stationary state. The mean and variance equations (11) and (13) are solved to provide model-consistent results for examining the mathematical consistency and convergence of the different EMCF methods. Only the first two scalar moments are used for the convergence testing since it is difficult for the EMCF methods to achieve consistency for higher-order moments as discussed in Sec. III. With the above numerical discretization, the error involved in the statistical results of the EMCF methods can be generally approximated as, based on the Taylor series analysis,

E=Em+cΔxΔx2+cΔtΔt+cNfNf1,
(29)

where E is the total error in the prediction of a statistical result such as ⟨ϕ⟩ from EMCF, Em is the model error (Em = 0 if the model is mathematically consistent), and cΔx, cΔt, and cNf are constants in the numerical errors that are independent of the grid size Δx, the time step size Δt, and Nf (the number of the Monte Carlo fields). The last term in Eq. (29) represents the numerical bias involved in EMCF due to the use of a finite value of Nf and it scales22,42 as Nf1. In the following numerical convergence testing, we specify a set of baseline numerical parameters, Δxb = 8 × 10−3 m, Δtb = 0.08 s, and Nf,b = 4. Uniform grids are used and the baseline grid size Δxb corresponds to ten uniform grid cells in the computational domain [−L, L]. The baseline time step size Δtb is specified by satisfying the Courant-Friedrichs-Lewy (CFL) condition, Δtb=CFL×minΔxb/(Las),Δxb2/(Γ+Γt),1/(2ω) with CFL = 0.4. Once we have a baseline case, we then vary the numerical parameters according to Δx2=Δxb2/M, Δt = Δtb/M, and Nf = Nf,bM with M being a free parameter such that the error terms in Eq. (29) change in a proportional manner,

E=Em+cΔxΔxb2+cΔtΔtb+cNfNf,b1M.
(30)

When M, the numerical errors tend to zero and EEm.

The mixing layer test case is statistically stationary for which a time-averaging strategy42 can be used to reduce the statistical errors involved in the EMCF methods. The time-averaged statistics are collected over about 1 000 000 time steps in the simulations to reduce the statistical errors. The total error in the statistical results in Eq. (30) is estimated from

EQ=QQFV,
(31)

where Q denote a statistical result such as ⟨ϕ⟩ obtained from the EMCF methods and the subscript “FV” denotes that the results are obtained by directly computing the mean or variance equations in (11) and (13) by using the FV method. To reduce the numerical errors in the FV results, we use a very fine grid for the FV simulations with Δx = 4 × 10−5 m (2000 uniform grid cells) and Δt = 4 × 10−6 s. Second-order accuracy in both space and time is used in the FV simulations.

A small turbulent Reynolds number Re = 20 is used in the convergence tests in the following Sec. IV B to focus on the mathematical consistency of the different EMCF methods. It is noted that the identified model inconsistency in Sec. III is closely related to molecular diffusion and hence examining the inconsistency under low Re is a wise choice to isolate its effect from other physics such as turbulent transport and other errors such as grid discretization errors. The sensitivity of the model inconsistency to Re will be examined later in Sec. IV C.

We first examine the mixing layer simulation results for a specified set of numerical parameters with M = 100 (Δx = 8 × 10−4 m, Δt = 8 × 10−4 s, and Nf = 400). Figure 2 shows the centerline profiles of the predicted mean ⟨ϕ⟩ and variance ϕ2 and their relative errors Eϕ/⟨ϕFV,x=0 and Eϕ2ϕ2FV,x=0 according to Eq. (31) against x/L. The subscript “FV, x = 0” denotes the centerline FV results at x = 0. For the predictions of ⟨ϕ⟩, the four different EMCF methods (EMCF-O, EMCF-M, EMCF-C1, and EMCF-C2) yield very close results that are in excellent agreement with the FV results. The relative errors Eϕ/⟨ϕFV,x=0 estimated from Eq. (31) are also very small, on the order of O(10−3). The slight difference in the errors for the predictions of the mean ⟨ϕ⟩ with the different EMCF methods is caused by statistical errors arising from the Wiener terms in the field equations. The results in Fig. 2 are consistent with the fact that all the EMCF methods in Sec. III yield consistently the same transport equation for ⟨ϕ⟩ as that from the modeled PDF transport equation (8). For the predictions of ⟨ϕ2⟩ in Fig. 2, the results from the different models are quite different from each other. The EMCF-O method significantly over-predicts ⟨ϕ2⟩ by about 48% relative to the FV results at x = 0. This over-prediction confirms the mathematical inconsistency in the EMCF-O method discussed in Sec. III A, in particular, the identified spurious variance production in the EMCF-O method [the term T6 in Eq. (19)] that can cause over-prediction of the scalar variance. The EMCF-M method under-predicts ⟨ϕ2⟩ by about 9.3% relative to the FV result at x = 0 in Fig. 2. This is consistent with the identified spurious variance dissipation in Sec. III B for the EMCF-M method. In comparison, the two corrected models EMCF-C1 and EMCF-C2 in Sec. III C yield good agreement with the FV results (relative error about 1.2%), which supports that both EMCF-C1 and EMCF-C2 methods are mathematically consistent for the predictions of variance as discussed in Sec. III C.

FIG. 2.

The predicted centerline profiles of the scalar mean ⟨ϕ⟩ and variance ⟨ϕ2⟩ and their normalized errors Eϕ/⟨ϕFV,x=0 and Eϕ2/ϕ2FV,x=0 against x/L in the opposed-jet turbulent mixing layer by using the different models (FV, EMCF-O, EMCF-M, EMCF-C1, and EMCF-C2).

FIG. 2.

The predicted centerline profiles of the scalar mean ⟨ϕ⟩ and variance ⟨ϕ2⟩ and their normalized errors Eϕ/⟨ϕFV,x=0 and Eϕ2/ϕ2FV,x=0 against x/L in the opposed-jet turbulent mixing layer by using the different models (FV, EMCF-O, EMCF-M, EMCF-C1, and EMCF-C2).

Close modal

A numerical convergence test of the different EMCF methods is conducted next with respect to the different values of M between 1 and 100. The relative errors (in magnitude) Eϕ2x=0ϕ2FV,x=0 obtained by using the different EMCF methods at the center of the mixing layer x = 0 are shown against M in Fig. 3. Only the convergence of the variance is considered since all the models are fully consistent for the predictions of ⟨ϕ⟩. From the figure, we can evidently see the non-convergent behavior of EMCF-O and EMCF-M when M increases, i.e., non-zero model error Em ≠ 0 when M due to the spurious variance production/dissipation involved in the models as discussed in Sec. III. On the contrary, EMCF-C1 and EMCF-C2 yield consistent convergence behavior over the whole range of M that has been examined, which again supports the consistency of the corrected EMCF for the variance predictions.

FIG. 3.

The relative prediction errors (in magnitude) Eϕ2x=0ϕ2FV,x=0 at the center of the mixing layer x = 0 against the values of M in the opposed-jet turbulent mixing layer by using the different EMCF methods. The symbols are the results by EMCF-O (red squares), EMCF-M (green circles), EMCF-C1 (blue diamonds), and EMCF-C2 (magenta triangles). The lines are curve fits according to Eq. (30).

FIG. 3.

The relative prediction errors (in magnitude) Eϕ2x=0ϕ2FV,x=0 at the center of the mixing layer x = 0 against the values of M in the opposed-jet turbulent mixing layer by using the different EMCF methods. The symbols are the results by EMCF-O (red squares), EMCF-M (green circles), EMCF-C1 (blue diamonds), and EMCF-C2 (magenta triangles). The lines are curve fits according to Eq. (30).

Close modal

In summary, through the numerical convergence tests of the turbulent mixing layer problem, we are able to verify the different convergence properties of the different EMCF methods discussed in Sec. III. These tests are done with a relatively low turbulent Re = 20 that does not represent real problems. To understand the effect of the model inconsistency under real turbulence conditions, we continue to examine the turbulent mixing layer test case under different Re numbers in Sec. IV C.

To examine the effect of Re on the performance of the EMCF methods, we specify the numerical parameters in the simulations of the opposed-jet mixing layer test case as Δx = L/750, Δt=CFL×minΔxb/(Las),Δxb2/(Γ+Γt),1/ω with CFL = 0.2, and Nf = 500 in the mixing layer test case and vary the turbulent Re number up to 400. The domain size [−L, L] is adaptively chosen to sufficiently cover the mixing layer under different Re. All the other parameters used in the simulations are the same as those used in Sec. IV B. Figure 4 shows the relative prediction errors for the scalar variance, Eϕ2,maxϕ2FV,x=0 ,where the subscript “max” denotes the maximum error in terms of the error magnitude, against the different values of Re ∈ [20, 400].

FIG. 4.

The relative prediction errors Eϕ2,maxϕ2FV,x=0 against the values of Re in the opposed-jet turbulent mixing layer by using the different EMCF methods. The symbols are the results by EMCF-O (red squares), EMCF-M (green circles), EMCF-C1 (blue diamonds), and EMCF-C2 (magenta triangles).

FIG. 4.

The relative prediction errors Eϕ2,maxϕ2FV,x=0 against the values of Re in the opposed-jet turbulent mixing layer by using the different EMCF methods. The symbols are the results by EMCF-O (red squares), EMCF-M (green circles), EMCF-C1 (blue diamonds), and EMCF-C2 (magenta triangles).

Close modal

From Fig. 4, we can evidently see that the relative error of the over-predicted variance by EMCF-O reduces sharply from about 48% at Re = 20 to about 2% at Re = 400. This confirms that the mathematical inconsistency in EMCF-O is purely a spurious molecular diffusion effect which tends to diminish for high Re number problems. This makes the model inconsistency in EMCF-O less concerning when the model is used in RANS modeling of practical applications where high Re is typically observed. This kind of model inconsistency due to molecular diffusion also widely exists in the turbulence modeling of multi-scalar transport like combustion, e.g., assuming equal molecular diffusivities.1,40 Although it has not been tested in this work, the importance of the mathematical inconsistency in EMCF-O is generally anticipated in LES since in the LES context the turbulent Re number associated with the unresolved turbulence scales is generally small.

The EMCF-M method in Fig. 4 yields negative errors (under-prediction of scalar variance) for the prediction of variance with the different values of Re which is consistent with the spurious variance dissipation involved in the model. The magnitudes of the errors increase slightly and reach a maximum in the magnitude of 13% at Re = 80, and then decrease slowly to 10% at Re = 400, as opposed to 2% in the other models at the same Re. This decrease in error with Re also supports the fact that the inconsistency in EMCF-M is a spurious molecular diffusion effect which tends to zero at high Re. However, it appears that the model error in EMCF-M decreases at a much slower rate than that in EMCF-O when Re increases. Although the intention of the EMCF-M method was to correct the spurious variance production in the original EMCF-O,19 it does so by unintentionally introducing a spurious variance dissipation that appears to be more serious at high Re.

The two new EMCF methods with corrections in Sec. III C, EMCF-C1 and EMCF-C2, yield errors mostly independent of Re, as shown in Fig. 4. The magnitudes of the errors are about 2%, which are expected to be purely numerical errors due to the use of finite values of numerical parameters in the simulations. The results confirm the effectiveness of the model correction proposed in Sec. III C for the spurious production/dissipation involved in EMCF-O and EMCF-M.

In summary, both EMCF-O and EMCF-M are inconsistent with the modeled PDF equation they intend to solve in terms of the variance predictions. The spurious variance production/dissipation involved in both EMCF-O and EMCF-M is expected to be purely a molecular diffusion effect and vanish when Re. Hence, despite the inconsistency, both models are still applicable to RANS modeling of many engineering problems where high Re is observed. In the LES context, however, the significance of the mathematical inconsistency can be more severe based on the same Re number argument (defined based on unresolved velocity scales). The EMCF-M method is found to be not effective to correct the inconsistency in EMCF-O as intended and the model error decays at a much slower rate than that by EMCF-O when Re increases in the opposed-jet turbulent mixing layer. The two corrected models, EMCF-C1 and EMCF-C2, are shown to be consistent (for the first two moments) numerically.

The opposed-jet turbulent mixing layer test case used in Sec. IV is a highly simplified mixing problem with many assumptions involved such as the assumed homogeneous turbulence without decaying. These assumptions are logically necessary in order to have an isolated system for a focused study of the mathematical inconsistency of the EMCF methods. To extend the examination of the EMCF methods to realistic problems, we use an experimental thermal wake behind a line source in grid turbulence by Warhaft43 to further examine the EMCF methods in terms of the model performance in real turbulent mixing problems in this section.

A classic line source dispersion experiment in grid turbulence43 is chosen as the test case, as illustrated in Fig. 5. The mean flow speed of a wind tunnel is U = 7 m/s. The turbulence is generated by a grid at x = 0 with grid size MG = 2.54 × 10−2 m. The turbulent kinetic energy at x/MG = 1 is measured43–45 to be kG = 8.13 m2/s2. A heated wire is placed at x0/MG = 52 to form a thermal line source. The wire is small in diameter d0 and hence its disturbance to the mean flow is assumed to be small. Two different cases are considered with d0 = 1.27 × 10−4 m and d0 = 2.5 × 10−5 m. The turbulence at the line source location x = x0 has a Reynolds number Re = 129 and decays following the mean flow. The temperature excess (denoted as ϕ) produced by the line source is considered to be a conserved scalar in the wake region behind the line source. The molecular (thermal) diffusivity is assumed to be uniform in the flow field and is approximated to be Γ = 2.1 × 10−5 m2/s.

FIG. 5.

Sketch of the experiment of the thermal wake behind a line source in grid turbulence by Warhaft.43 

FIG. 5.

Sketch of the experiment of the thermal wake behind a line source in grid turbulence by Warhaft.43 

Close modal

The line source problem is commonly formulated into a transient problem by following a moving frame x(t),44–46 

x(t)=x0+Ut.
(32)

In the EMCF modeling of scalar transport and mixing, we obtain the turbulent kinetic energy from the experiment43 as follows44,45:

k(t)=kGx(t)MGm,
(33)

where kG = 8.13 m2/s2 and m = 1.4. The turbulent kinetic energy dissipation rate ε(t) is then obtained from

ε(t)=ddtk(t).
(34)

The turbulent diffusivity Γt is modeled47 as

Γt(t)=12ddtd2(t)Γ,
(35)

where d(t) is the characteristic width of the mean field,46 

d2(t)=d02+2Γt+2v2t=0x0U2×(1+Ut/x0)rsr(rs)+(1+Ut/x0)srs1s(rs),
(36)

in which v2t=0 is the variance of the velocity in the y direction at t = 0 (or x = x0), r=m3C0/21/2+1,s=m3C0/2+1/21, and the constant44,46C0 = 3.0. The mixing frequency ω is modeled following Sawford:44 

ω1=0.7t1+tanhln(Ut/x0)+4.52.
(37)

With the above models and modeling parameters, the EMCF model equations in Sec. III are in closed forms and can be solved numerically. The simulations start at t = 0 (or x = x0) when the initial mean and variance of the scalar are specified as

ϕt=0=expy22d02,
(38)
ϕ2t=0=0.
(39)

The Monte Carlo fields are initialized identically at t = 0, ϕi(0) = ⟨ϕt=0 (i = 1, …, Nf). The initial computational domain in y spans between [−yL, yL] = [−15d0, 15d0]. A uniform grid is used and the initial grid size Δx/d0 = 0.018 75 (corresponding to 1600 grid cells). The computational domain adaptively expands in time to have at least 320 grid cells inside the line source wake. The time step size Δt is determined from the CFL condition Δt=CFL×minΔx2/(Γ+Γt),1/ω with CFL = 0.25. The number of fields used in the simulations is Nf = 400.

Figure 6 shows the time history of the characteristic width of the thermal wake d(t)/MG and the scalar root mean square (RMS) ϕRMS,y=0=ϕ2y=01/2 at y = 0 in the thermal wake behind a line source. The characteristic width of the thermal wake d(t) in Eq. (36) represents the predictions of the scalar mean ⟨ϕ⟩ and shows excellent agreement with the experiments for both experimental cases with different d0. No obvious difference is observed between the FV and the different EMCF results, which is consistent with the fact that all the EMCF methods are consistent for the scalar mean. For ϕRMS, EMCF-O consistently over-predicts ϕRMS along the centerline at y = 0, which is consistent with the identified spurious variance production in Sec. III A and also the results shown in Sec. IV B. The EMCF-M method consistently under-predicts ϕRMS along the centerline expectedly due to its spurious variance dissipation. The two new EMCF methods, EMCF-C1 and EMCF-C2, show overall best agreement with the experiments. These results demonstrate the apparent effect of the mathematical inconsistency of the different EMCF methods on the model predictions of real turbulent mixing.

FIG. 6.

Time history of the characteristic width of the thermal wake d(t)/MG (left subplot) and the scalar RMS ϕRMS,y=0 at y = 0 (right subplot) in the thermal wake behind a line source. The symbols (red solid squares and diamonds) denote two different experimental conditions by Warhaft,43 the circles are predictions by the FV method and the lines are the results of the different EMCF methods in Sec. III, EMCF-O (red solid lines), EMCF-M (green dashed lines), EMCF-C1 (blue dash-dotted lines), and EMCF-C2 (magenta dotted lines).

FIG. 6.

Time history of the characteristic width of the thermal wake d(t)/MG (left subplot) and the scalar RMS ϕRMS,y=0 at y = 0 (right subplot) in the thermal wake behind a line source. The symbols (red solid squares and diamonds) denote two different experimental conditions by Warhaft,43 the circles are predictions by the FV method and the lines are the results of the different EMCF methods in Sec. III, EMCF-O (red solid lines), EMCF-M (green dashed lines), EMCF-C1 (blue dash-dotted lines), and EMCF-C2 (magenta dotted lines).

Close modal

Figure 7 further examines the profiles of the normalized scalar RMS, ϕRMS/ϕRMS,y=0 against y/d(t) in the thermal wake behind a line source. Due to the normalization by using the centerline RMS ϕRMS,y=0, the y axis in Fig. 7 is always one at y = 0. All the model results are in reasonable agreement with the experimental data. The FV results are in excellent agreement with the experimental data as shown in the figure. The three models, EMCF-O, EMCF-C1, and EMCF-C2, show good agreement with the results as well (for the normalized scalar RMS). The EMCF-M method shows the most discrepancy from the experimental data due to the spurious variance dissipation which decays slowly when Re increases as shown in Fig. 4. The modification introduced in Ref. 19, unfortunately, makes the EMCF method even worse compared to the original version for this line source dispersion problem.

FIG. 7.

The profiles of the normalized scalar RMS, ϕRMS/ϕRMS,y=0 against y/d(t) in the thermal wake behind a line source. The symbols (red diamonds and squares) denote the experimental data,43 the circles are the FV results, and the lines are the predictions by the different EMCF methods in Sec. III, EMCF-O (red solid lines), EMCF-M (green dashed lines), EMCF-C1 (blue dash-dotted lines), and EMCF-C2 (magenta dotted lines).

FIG. 7.

The profiles of the normalized scalar RMS, ϕRMS/ϕRMS,y=0 against y/d(t) in the thermal wake behind a line source. The symbols (red diamonds and squares) denote the experimental data,43 the circles are the FV results, and the lines are the predictions by the different EMCF methods in Sec. III, EMCF-O (red solid lines), EMCF-M (green dashed lines), EMCF-C1 (blue dash-dotted lines), and EMCF-C2 (magenta dotted lines).

Close modal

In summary, a thermal wake behind a line source experiment by Warhaft43 is briefly examined in this section by using the different EMCF methods. The results support all the conclusions drawn from the numerical convergence study of the idealized turbulent mixing layer test case in Sec. IV.

This paper presents a mathematical approach to thoroughly examine the mathematical consistency of the EMCF methods as an efficient approach for solving the transported PDF equation. Two existing versions of the EMCF methods (EMCF-O and EMCF-M) are examined. The EMCF-O method is found to be consistent only for yielding the scalar mean and to contain a spurious variance production in the model. This spurious variance production is mainly caused by molecular diffusion and it vanishes in the limit of infinite Re. The EMCF-M method is also found to be consistent only for yielding the scalar mean and contains a spurious various dissipation. This spurious variance dissipation is also found to be caused by molecular diffusion and is expected to reduce when Re increases. However, based on this study, we find that the spurious variance dissipation reduces at a much slower rate when Re increases than the spurious variance production in the EMCF-O method. Corrected EMCF methods are introduced to eliminate the spurious variance production/dissipation for the variance transport. The mathematical inconsistency of the methods is confirmed numerically by conducting numerical convergence tests in an opposed-jet turbulent mixing layer. An experimental thermal wake behind a line source in grid turbulence is also examined briefly to further confirm the findings. The mathematical inconsistency found in the existing EMCF methods is expected to be marginally important in RANS modeling where low Re-number regions (e.g., near wall flow) are encountered. In LES, however, the importance of the mathematical inconsistency can be more evident because the Re number associated with the unresolved velocity scale can be small when a reasonable LES resolution is used.

This paper is based upon work supported by the National Science Foundation under Grant No. CBET-1336075. This research was supported in part through computational resources provided by Information Technology at Purdue University, West Lafayette, Indiana.

By multiplying the modeled PDF Eq. (8) by ψm and then integrating over the entire sample space ψ ∈ (−, ), we obtain the integrated moments equation as

ϕmt+uϕm=(Γt+Γ)ϕm   ψm2ψ2Γϕϕfϕ(ψ;x,t)dψT1   +ψmωψϕfϕ(ψ;x,t)ψdψT2   ψmS(ψ)fϕ(ψ;x,t)ψdψT3.
(A1)

By using integration by parts and the property1 fϕ(ψ = ±; x, t) = 0, we obtain the three terms T1, T2, and T3 in Eq. (A1) as

T1=ψm2ψ2Γϕϕfϕ(ψ;x,t)dψ=m(m1)Γϕϕϕm2,
(A2)
T2=ψmωψϕfϕ(ψ;x,t)ψdψ=mωϕmϕϕm1,
(A3)
T3=ψmS(ψ)fϕ(ψ;x,t)ψdψ=mϕm1S(ϕ).
(A4)

Substituting Eqs. (A2)–(A4) to Eq. (A1), we obtain the moment equation (10) derived from the modeled PDF Eq. (8).

To derive the moment equation from the EMCF field equation, we take an Euler step with an infinitesimal time step size dt for each Monte Carlo field from time t in Eq. (16) so that the updated fields at time t + dt become

ϕnf(x,t+dt)=ϕnf(x,t)uϕnf(x,t)dt+(Γt+Γ)ϕnf(x,t)dt+(2Γt+2Γ)1/2ϕnf(x,t)dWnf(t)ωϕnf(x,t)ϕ*(x,t)Nfdt+Sϕnf(x,t)dt.
(B1)

This is equivalent to integrating equation (16) from time t to t + dt. Since dt is infinitely small, the right-hand side of Eq. (16) can be frozen at time t in Eq. (B1). To find the mth moment, we take the mth power of both sides of Eq. (B1) to yield

ϕnfm(x,t+dt)=ϕnf(x,t)uϕnf(x,t)dt+(Γt+Γ)ϕnf(x,t)dt+(2Γt+2Γ)1/2ϕnf(x,t)dWnf(t)ωϕnf(x,t)ϕ*(x,t)Nfdt+Sϕnf(x,t)dtm=ϕnfm(x,t)mϕnfm1(x,t)uϕnf(x,t)dt+mϕnfm1(x,t)(Γt+Γ)ϕnf(x,t)dt+mϕnfm1(x,t)(2Γt+2Γ)1/2ϕnf(x,t)dWnf(t)+m(m1)2ϕnfm2(x,t)(2Γt+2Γ)ϕnf(x,t)dWnf(t)2mϕnfm1(x,t)ωϕnf(x,t)ϕ*(x,t)Nfdt+mϕnfm1(x,t)Sϕnf(x,t)dt+R=ϕnfm(x,t)uϕnfm(x,t)dt+(Γt+Γ)ϕnfm(x,t)dtm(m1)ϕnfm2(x,t)(Γt+Γ)ϕnf(x,t)ϕnf(x,t)dt+mϕnfm1(x,t)(2Γt+2Γ)1/2ϕnf(x,t)dWnf(t)+m(m1)ϕnfm2(x,t)(Γt+Γ)ϕnf(x,t)dWnf(t)2mϕnfm1(x,t)ωϕnf(x,t)ϕ*(x,t)Nfdt+mϕnfm1(x,t)Sϕnf(x,t)dt+R,
(B2)

in which all the higher-order terms are grouped into ROdt3/2. The half-order dt1/2 arises from the Wiener term1 since dWnf=ηdt1/2, where η is a vector with each element being a Gaussian random number with zero mean and unit variance.

Now we introduce a field-averaging operator Nf over all the Monte Carlo fields as

Qϕ*Nf=1Nfnf=1NfQϕnf,
(B3)

where Q is an arbitrary function. In the limit of Nf, the field-averaging operator converges to the Reynolds averaging

limNfQϕ*Nf=Qϕ*.
(B4)

Applying the field-averaging operator in Eqs. (B3) to (B2), we have

ϕ*m(x,t+dt)Nf=ϕ*m(x,t)Nfuϕ*m(x,t)Nfdt+(Γt+Γ)ϕ*m(x,t)Nfdtm(m1)(Γt+Γ)ϕ*m2(x,t)ϕ*(x,t)ϕ*(x,t)Nfdt+m(2Γt+2Γ)1/2ϕ*m1(x,t)ϕ*(x,t)dW*(t)Nf+m(m1)(Γt+Γ)ϕ*m2(x,t)ϕ*(x,t)dW*(t)2Nfmωϕ*m(x,t)Nfϕ*(x,t)Nfϕ*m1(x,t)Nfdt+mϕ*m1(x,t)Sϕ*(x,t)Nfdt+RNf.
(B5)

By taking the limit of Nf in Eq. (B5) and using Eq. (B4), we obtain

ϕ*m(x,t+dt)=ϕ*m(x,t)uϕ*m(x,t)dt+(Γt+Γ)ϕ*m(x,t)dtm(m1)(Γt+Γ)ϕ*m2(x,t)ϕ*(x,t)ϕ*(x,t)dt+m(2Γt+2Γ)1/2ϕ*m1(x,t)ϕ*(x,t)dW*(t)T4+m(m1)(Γt+Γ)ϕ*m2(x,t)ϕ*(x,t)dW*(t)2T5mωϕ*m(x,t)ϕ*(x,t)ϕ*m1(x,t)dt+mϕ*m1(x,t)Sϕ*(x,t)dt+R.
(B6)

The term T4 in Eq. (B6) is simplified to

T4=m(2Γt+2Γ)1/2ϕ*m1(x,t)ϕ*(x,t)dW*(t)=m(2Γt+2Γ)1/2ϕ*m1(x,t)ϕ*(x,t)dW*(t)=0,
(B7)

where the statistical independence between ϕ*m1(x,t)ϕ*(x,t) and dW*(t) has been used in the second line of the equation, and dW*(t)=0 by definition has been used in the third line. The term T5 in Eq. (B6) is reduced to

T5=m(m1)(Γt+Γ)ϕ*m2(x,t)ϕ*(x,t)dW*(t)2=m(m1)(Γt+Γ)ϕ*m2(x,t)ϕ*(x,t)ϕ*(x,t):dW*(t)dW*(t)=m(m1)(Γt+Γ)ϕ*m2(x,t)ϕ*(x,t)ϕ*(x,t)dt,
(B8)

in which “:” denotes the double dot product, the statistical independence between ϕ*m2(x,t)ϕ*(x,t)ϕ*(x,t) and dW*(t) has been used in the second line of the equation, and dW*(t)dW*(t)=Idt has been used in the third line of the equation. Substituting Eqs. (B7) and (B8) into Eq. (B6), we obtain the moment equation (17) in which we have used limdt0R/dt=0 since R contains all high-order terms and

ϕ*m(x,t)t=limdt0ϕ*m(x,t+dt)ϕ*m(x,t)dt.
(B9)

Hereby we have presented a mathematical approach for a rigorous examination of the mathematical consistency of the EMCF methods through the moment equations. The derivation of the moment equation from the SPDEs is new and has not been reported before.

1.
S. B.
Pope
, “
PDF methods for turbulent reactive flows
,”
Prog. Energy Combust. Sci.
11
,
119
192
(
1985
).
2.
C.
Dopazo
, “
Recent developments in PDF methods
,” in
Turbulent Reacting Flows
, edited by
F. A. W. P. A.
Libby
(
Academic Press
,
1994
), pp.
376
457
.
3.
D. C.
Haworth
, “
Progress in probability density function methods for turbulent reacting flows
,”
Prog. Energy Combust. Sci.
36
,
168
259
(
2010
).
4.
J.
Xu
and
S. B.
Pope
, “
PDF calculations of turbulent nonpremixed flames with local extinction
,”
Combust. Flame
123
,
281
307
(
2000
).
5.
R. P.
Lindstedt
,
S. A.
Louloudi
, and
E. M.
Vaos
, “
Joint scalar probability density function modeling of pollutant formation in piloted turbulent jet diffusion flames with comprehensive chemistry
,”
Proc. Combust. Inst.
28
,
149
156
(
2000
).
6.
R.
Cabra
,
T.
Myhrvold
,
J. Y.
Chen
,
R. W.
Dibble
,
A. N.
Karpetis
, and
R. S.
Barlow
, “
Simultaneous laser Raman-Rayleigh-LIF measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow
,”
Proc. Combust. Inst.
29
,
1881
1888
(
2002
).
7.
R. R.
Cao
,
S. B.
Pope
, and
A. R.
Masri
, “
Turbulent lifted flames in a vitiated coflow investigated using joint PDF calculations
,”
Combust. Flame
142
,
438
453
(
2005
).
8.
W. P.
Jones
and
S.
Navarro-Martinez
, “
Study of hydrogen auto-ignition in a turbulent air co-flow using a large eddy simulation approach
,”
Comput. Fluids
37
,
802
808
(
2008
).
9.
F.
Sewerin
and
S.
Rigopoulos
, “
An LES-PBE-PDF approach for predicting the soot particle size distribution in turbulent flames
,”
Combust. Flame
189
,
62
76
(
2018
).
10.
G.
Bulat
,
W. P.
Jones
, and
A. J.
Marquis
, “
NO and CO formation in an industrial gas-turbine combustion chamber using LES with the Eulerian sub-grid PDF method
,”
Combust. Flame
161
,
1804
1825
(
2014
).
11.
J.
Dumond
,
F.
Magagnato
, and
A.
Class
, “
Stochastic-field cavitation model
,”
Phys. Fluids
25
,
073302
(
2013
).
12.
M.
Cassiani
,
J. F.
Vinuesa
,
S.
Galmarini
, and
B.
Denby
, “
Stochastic fields method for sub-grid scale emission heterogeneity in mesoscale atmospheric dispersion models
,”
Atmos. Chem. Phys.
10
,
267
277
(
2010
).
13.
N.
Suciu
,
L.
Schüler
,
S.
Attinger
, and
P.
Knabner
, “
Towards a filtered density function approach for reactive transport in groundwater
,”
Adv. Water Resour.
90
,
83
98
(
2016
).
14.
D. W.
Meyer
,
P.
Jenny
, and
H. A.
Tchelepi
, “
A joint velocity-concentration PDF method for tracer flow in heterogeneous porous media
,”
Water Resour. Res.
46
,
W12522
, https://doi.org/10.1029/2010wr009450 (
2010
).
15.
H. C.
Öttinger
,
B. H. A. A.
van den Brule
, and
M. A.
Hulsen
, “
Brownian configuration fields and variance reduced CONNFFESSIT
,”
J. Non-Newtonian Fluid Mech.
70
,
255
261
(
1997
).
16.
S. B.
Pope
, “
A Monte Carlo method for the PDF equations of turbulent reactive flow
,”
Combust. Sci. Technol.
25
,
159
174
(
1981
).
17.
L.
Valiño
, “
A field Monte Carlo formulation for calculating the probability density function of a single scalar in a turbulent flow
,”
Flow, Turbul. Combust.
60
,
157
172
(
1998
).
18.
V.
Sabel’nikov
and
O.
Soulard
, “
Rapidly decorrelating velocity-field model as a tool for solving one-point Fokker-Planck equations for probability density functions of turbulent reactive scalars
,”
Phys. Rev. E
72
,
016301
(
2005
).
19.
L.
Valiño
,
R.
Mustata
, and
K. B.
Letaief
, “
Consistent behavior of Eulerian Monte Carlo fields at low Reynolds numbers
,”
Flow, Turbul. Combust.
96
,
503
512
(
2016
).
20.
R. O.
Fox
,
Computational Models for Turbulent Reacting Flows
(
Cambridge University Press
,
2003
).
21.
S. B.
Pope
, “
Particle method for turbulent flows: Integration of stochastic model equations
,”
J. Comput. Phys.
117
,
332
349
(
1995
).
22.
J.
Xu
and
S. B.
Pope
, “
Assessment of numerical accuracy of PDF/Monte-Carlo methods for turbulent reactive flows
,”
J. Comput. Phys.
152
,
192
230
(
1999
).
23.
H.
Wang
,
P. P.
Popov
, and
S. B.
Pope
, “
Weak second-order splitting schemes for Lagrangian Monte Carlo particle methods for the composition PDF/FDF transport equations
,”
J. Comput. Phys.
229
,
1852
1878
(
2010
).
24.
H.
Wang
and
P.
Zhang
, “
A unified view of pilot stabilized turbulent jet flames for model assessment across different combustion regimes
,”
Proc. Combust. Inst.
36
,
1693
1703
(
2017
).
25.
P.
Zhang
,
A. R.
Masri
, and
H.
Wang
, “
Studies of the flow and turbulence fields in a turbulent pulsed jet flame using LES/PDF
,”
Combust. Theory Modell.
21
,
897
924
(
2017
).
26.
A.
Garmory
and
E.
Mastorakos
, “
Aerosol nucleation and growth in a turbulent jet using the stochastic fields method
,”
Chem. Eng. Sci.
63
,
4078
4089
(
2008
).
27.
R.
Mustata
,
L.
Valiño
,
C.
Jiménez
,
W. P.
Jones
, and
S.
Bondi
, “
A probability density function Eulerian Monte Carlo field method for large eddy simulations: Application to a turbulent piloted methane/air diffusion flame (Sandia D)
,”
Combust. Flame
145
,
88
104
(
2006
).
28.
A.
Avdić
,
G.
Kuenne
,
F.
di Mare
, and
J.
Janicka
, “
LES combustion modeling using the Eulerian stochastic field method coupled with tabulated chemistry
,”
Combust. Flame
175
,
201
219
(
2017
).
29.
Q.
Tang
,
W.
Zhao
,
M.
Bockelie
, and
R. O.
Fox
, “
Multi-environment probability density function method for modelling turbulent combustion using realistic chemical kinetics
,”
Combust. Theory Modell.
11
,
889
907
(
2007
).
30.
J.
Lee
,
S.
Jeon
, and
Y.
Kim
, “
Multi-environment probability density function approach for turbulent CH4/H2 flames under the MILD combustion condition
,”
Combust. Flame
162
,
1464
1476
(
2015
).
31.
J.
Lee
and
Y.
Kim
, “
DQMOM based PDF transport modeling for turbulent lifted nitrogen-diluted hydrogen jet flame with autoignition
,”
Int. J. Hydrogen Energy
37
,
18498
18508
(
2012
).
32.
H.
Koo
,
P.
Donde
, and
V.
Raman
, “
A quadrature-based LES/transported probability density function approach for modeling supersonic combustion
,”
Proc. Combust. Inst.
33
,
2203
2210
(
2011
).
33.
J.
Akroyd
,
A. J.
Smith
,
L. R.
McGlashan
, and
M.
Kraft
, “
Comparison of the stochastic fields method and DQMoM-IEM as turbulent reaction closures
,”
Chem. Eng. Sci.
65
,
5429
5441
(
2010
).
34.
J.
Jaishree
and
D. C.
Haworth
, “
Comparisons of Lagrangian and Eulerian PDF methods in simulations of non-premixed turbulent jet flames with moderate-to-strong turbulence-chemistry interactions
,”
Combust. Theory Modell.
16
,
435
463
(
2012
).
35.
S. B.
Pope
,
Turbulent Flows
(
Cambridge University Press
,
Cambridge
,
2000
), p.
806
.
36.
S. B.
Pope
, “
Computations of turbulent combustion: Progress and challenges
,”
Symp. (Int.) Combust.
23
,
591
612
(
1990
).
37.
C.
Dopazo
and
E. E.
O’Brien
, “
An approach to the autoignition of a turbulent mixture
,”
Acta Astronaut.
1
,
1239
1266
(
1974
).
38.
K. A.
Kemenov
,
H.
Wang
, and
S. B.
Pope
, “
Turbulence resolution scale dependence in large-eddy simulations of a jet flame
,”
Flow, Turbul. Combust.
88
,
529
561
(
2012
).
39.
G.
Dixon-Lewis
,
T.
David
,
P. H.
Gaskell
,
S.
Fukutani
,
H.
Jinno
,
J. A.
Miller
,
R. J.
Kee
,
M. D.
Smooke
,
N.
Peters
,
E.
Effelsberg
,
J.
Warnatz
, and
F.
Behrendt
, “
Calculation of the structure and extinction limit of a methane-air counterflow diffusion flame in the forward stagnation region of a porous cylinder
,”
Symp. (Int.) Combust.
20
,
1893
1904
(
1984
).
40.
H.
Wang
, “
Consistent flamelet modeling of differential molecular diffusion for turbulent non-premixed flames
,”
Phys. Fluids
28
,
035102
(
2016
).
41.
J. D.
Anderson
,
Computational Fluid Dynamics: The Basics with Applications
(
McGraw-Hill
,
1985
).
42.
H.
Wang
and
S. B.
Pope
, “
Time-averaging strategies in the finite-volume/particle hybrid algorithm for the joint PDF equation of turbulent reactive flows
,”
Combust. Theory Modell.
12
,
529
544
(
2008
).
43.
Z.
Warhaft
, “
The interference of thermal fields from line sources in grid turbulence
,”
J. Fluid Mech.
144
,
363
387
(
1984
).
44.
B. L.
Sawford
, “
Micro-mixing modelling of scalar fluctuations for plumes in homogeneous turbulence
,”
Flow, Turbul. Combust.
72
,
133
160
(
2004
).
45.
S.
Viswanathan
and
S. B.
Pope
, “
Turbulent dispersion from line sources in grid turbulence
,”
Phys. Fluids
20
,
101514
(
2008
).
46.
M. S.
Anand
and
S. B.
Pope
, “
Diffusion behind a line source in grid turbulence
,” in
Turbulent Shear Flows
(
Springer
,
1985
), Vol. 4, pp.
46
61
.
47.
J. W.
Deardorff
, “
Closure of second-and third-moment rate equations for diffusion in homogeneous turbulence
,”
Phys. Fluids
21
,
525
530
(
1978
).