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.

## I. INTRODUCTION

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) method^{1–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} NO_{x} 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) method^{16} (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, Pope^{1} 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 method^{17–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 method^{20} 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 Haworth^{34} 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ño^{17} 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ño^{17} 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:

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

to develop new EMCF methods with a better mathematical consistency,

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

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.

## II. TRANSPORTED PDF METHOD

We consider a single-scalar field $\varphi x,t$ transported by a solenoidal turbulent flow field **u** = {*u*_{1}, *u*_{2}, *u*_{3}}, ∇ · **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

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 function^{1} *δ*(*ψ* − *ϕ*(**x**, *t*)),

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

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

in which ⟨·⟩ denotes Reynolds averaging, the prime “′” denotes fluctuations, Γ_{t} is the turbulent eddy diffusivity, *ω* is the scalar mixing frequency, a gradient-diffusion transport model^{35} is used in Eq. (5), and the interaction by exchange with the mean (IEM) mixing model^{37} 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:

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 methods^{17–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 approach^{21,23} is often used by examining the mathematical consistency of the moments of scalars obtained from the EMCF methods. The *m*th central moments (*m* = 1, 2, …, *∞*) of the scalar $\varphi m$ is defined as

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

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

The transport equation of the scalar variance define as $\varphi \u20322=\varphi 2\u2212\u27e8\varphi \u27e92$ can be obtained from Eqs. (11) and (12) as follows:

## III. CONSISTENCY OF EULERIAN MONTE CARLO FIELDS METHOD

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.

### A. Eulerian Monte Carlo fields method: The original form

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),

The PDF $f\varphi *(\psi *;x,t)$ can be approximated by *N*_{f} discrete Monte Carlo fields $\varphi nf(x,t)$, (*n*_{f} = 1, …, *N*_{f}),

When *N*_{f} → *∞*, 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\varphi *(\psi *;x,t)$ obtained from the EMCF methods in Eq. (15) when *N*_{f} → *∞* and *f*_{ϕ}(*ψ*; **x**, *t*) governed by Eq. (8).

In the original EMCF method by Valiño,^{17} the SPDEs for the *N*_{f} Monte Carlo fields $\varphi nf(x,t)$ are written as

where $Wnf(t)$ is a vector-valued Wiener process with $\u27e8dWnf\u27e9=0$ and $\u27e8dWnfdWnf\u27e9=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

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 *N*_{f} → *∞*. 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:

By comparing the transport equations for the mean ⟨*ϕ*⟩ in (11) and for $\varphi *$ 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.

### B. Eulerian Monte Carlo fields method: The modified form

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 as^{19}

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 $\varphi *m$ from Eq. (20) as

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

where $\u2207\varphi *\u2032\u22c5\u2207\varphi *\u2032=\u2207\varphi *\u22c5\u2207\varphi *\u2212\u2207\varphi *\u22c5\u2207\varphi *$. 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.

### C. Eulerian Monte Carlo fields method: Consistency correction

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:

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:

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}:

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:

or

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

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.

## IV. CONVERGENCE STUDY OF EULERIAN MONTE CARLO FIELD METHOD

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. Opposed-jet turbulent mixing layer test case

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, $\u016b$(*x*) = −*a*_{s}*x*, where *a*_{s} 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.

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 $\u2212L,\u2009L$ with *L* = 0.04 m, the strain rate *a*_{s} = 1.0 s^{−1}, the kinematic viscosity *ν* = 1.65 × 10^{−5} m^{2}/s, the turbulent kinetic energy (homogeneous) *k* = 5 × 10^{−4} m^{2}/s^{2}, 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 scale^{40} varies over a range between [20, 400] in the discussions below. With these parameters, the turbulent kinetic energy dissipation rate is computed^{40} as $\epsilon $ = 2*k*^{2}/(3*νRe*), the molecular diffusivity Γ = *ν*/*σ*, and the turbulent diffusivity Γ_{t} = *ν*_{t}/*σ*_{t}. The turbulent viscosity *ν*_{t} is modeled following the *k*–$\epsilon $ turbulence model^{35} as *ν*_{t} = *C*_{μ}*k*^{2}/$\epsilon $ with *C*_{μ} = 0.09. The mixing frequency *ω* is commonly modeled^{1} as *ω* = *c*_{ϕ}$\epsilon $/(2*k*) with *c*_{ϕ} = 2.0.

The governing equations for the scalar mean ⟨*ϕ*⟩ in (11), the variance $\varphi \u20322$ 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,

where *E* is the total error in the prediction of a statistical result such as ⟨*ϕ*⟩ from EMCF, *E*_{m} is the model error (*E*_{m} = 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 *N*_{f} (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 *N*_{f} and it scales^{22,42} as $Nf\u22121$. In the following numerical convergence testing, we specify a set of baseline numerical parameters, Δ*x*_{b} = 8 × 10^{−3} m, Δ*t*_{b} = 0.08 s, and *N*_{f,b} = 4. Uniform grids are used and the baseline grid size Δ*x*_{b} corresponds to ten uniform grid cells in the computational domain [−*L*, *L*]. The baseline time step size Δ*t*_{b} is specified by satisfying the Courant-Friedrichs-Lewy (CFL) condition, $\Delta tb=CFL\xd7min\Delta xb/(Las),\Delta xb2/(\Gamma +\Gamma t),1/(2\omega )$ with CFL = 0.4. Once we have a baseline case, we then vary the numerical parameters according to $\Delta x2=\Delta xb2/M$, Δ*t* = Δ*t*_{b}/*M*, and *N*_{f} = *N*_{f,b}*M* with *M* being a free parameter such that the error terms in Eq. (29) change in a proportional manner,

When *M* → *∞*, the numerical errors tend to zero and *E* → *E*_{m}.

The mixing layer test case is statistically stationary for which a time-averaging strategy^{42} 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

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.

### B. Convergence results and discussion

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 *N*_{f} = 400). Figure 2 shows the centerline profiles of the predicted mean ⟨*ϕ*⟩ and variance $\varphi \u20322$ and their relative errors *E*_{⟨ϕ⟩}/⟨*ϕ*⟩_{FV,x=0} and $E\varphi \u20322\u27e8\varphi \u20322\u27e9FV,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.

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\u27e8\varphi \u20322\u27e9x=0\u27e8\varphi \u20322\u27e9FV,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 *E*_{m} ≠ 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.

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.

### C. Model performance under different Reynolds numbers

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, $\Delta t=CFL\xd7min\Delta xb/(Las),\Delta xb2/(\Gamma +\Gamma t),1/\omega $ with CFL = 0.2, and *N*_{f} = 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\u27e8\varphi \u20322\u27e9,max\u27e8\varphi \u20322\u27e9FV,x=0$ ,where the subscript “max” denotes the maximum error in terms of the error magnitude, against the different values of *Re* ∈ [20, 400].

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.

## V. MODEL ASSESSMENT IN A THERMAL WAKE FROM A LINE SOURCE IN GRID TURBULENCE

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 Warhaft^{43} to further examine the EMCF methods in terms of the model performance in real turbulent mixing problems in this section.

### A. Thermal wake behind a line source in grid turbulence

A classic line source dispersion experiment in grid turbulence^{43} 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 *M*_{G} = 2.54 × 10^{−2} m. The turbulent kinetic energy at *x*/*M*_{G} = 1 is measured^{43–45} to be *k*_{G} = 8.13 m^{2}/s^{2}. A heated wire is placed at *x*_{0}/*M*_{G} = 52 to form a thermal line source. The wire is small in diameter *d*_{0} and hence its disturbance to the mean flow is assumed to be small. Two different cases are considered with *d*_{0} = 1.27 × 10^{−4} m and *d*_{0} = 2.5 × 10^{−5} m. The turbulence at the line source location *x* = *x*_{0} 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} m^{2}/s.

### B. Modeling approaches

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

In the EMCF modeling of scalar transport and mixing, we obtain the turbulent kinetic energy from the experiment^{43} as follows^{44,45}:

where *k*_{G} = 8.13 m^{2}/s^{2} and *m* = 1.4. The turbulent kinetic energy dissipation rate $\epsilon $(*t*) is then obtained from

The turbulent diffusivity Γ_{t} is modeled^{47} as

where *d*(*t*) is the characteristic width of the mean field,^{46}

in which $\u27e8v\u20322\u27e9t=0$ is the variance of the velocity in the y direction at *t* = 0 (or *x* = *x*_{0}), $r=m3C0/2\u22121/2+1,s=m3C0/2+1/2\u22121$, and the constant^{44,46}*C*_{0} = 3.0. The mixing frequency *ω* is modeled following Sawford:^{44}

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* = *x*_{0}) when the initial mean and variance of the scalar are specified as

The Monte Carlo fields are initialized identically at *t* = 0, *ϕ*_{i}(0) = ⟨*ϕ*⟩_{t=0} (*i* = 1, …, *N*_{f}). The initial computational domain in *y* spans between [−*y*_{L}, *y*_{L}] = [−15*d*_{0}, 15*d*_{0}]. A uniform grid is used and the initial grid size Δ*x*/*d*_{0} = 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 $\Delta t=CFL\xd7min\Delta x2/(\Gamma +\Gamma t),1/\omega $ with CFL = 0.25. The number of fields used in the simulations is *N*_{f} = 400.

### C. Results and discussion

Figure 6 shows the time history of the characteristic width of the thermal wake *d*(*t*)/*M*_{G} and the scalar root mean square (RMS) $\varphi RMS,y=0=\u27e8\varphi \u20322\u27e9y=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 *d*_{0}. 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.

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.

## VI. CONCLUSION

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.

## ACKNOWLEDGMENTS

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.

### APPENDIX A: DERIVATION OF MOMENT EQUATION FROM THE PDF TRANSPORT EQUATION

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

By using integration by parts and the property^{1} *f*_{ϕ}(*ψ* = ±*∞*; **x**, *t*) = 0, we obtain the three terms $T1$, $T2$, and $T3$ in Eq. (A1) as

### APPENDIX B: DERIVATION OF MOMENTS EQUATION FROM THE EMCF FIELD EQUATION

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

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 *m*th moment, we take the *m*th power of both sides of Eq. (B1) to yield

in which all the higher-order terms are grouped into $R\u223cOdt3/2$. The half-order *dt*^{1/2} arises from the Wiener term^{1} since $dWnf=\eta \u22c5dt1/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 $\u22c5Nf$ over all the Monte Carlo fields as

where *Q* is an arbitrary function. In the limit of *N*_{f} → *∞*, the field-averaging operator converges to the Reynolds averaging

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

where the statistical independence between $\varphi *m\u22121(x,t)\u2207\varphi *$$(x,t)$ and *d***W**^{*}(*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

in which “:” denotes the double dot product, the statistical independence between $\varphi *m\u22122(x,t)\u2207\varphi *(x,t)\u2207\varphi *(x,t)$ and *d***W**^{*}(*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 $limdt\u21920R/dt=0$ since *R* contains all high-order terms and

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.

## REFERENCES

_{2}/N

_{2}jet flame in a vitiated coflow

_{4}/H

_{2}flames under the MILD combustion condition