An extended quadrature method of moments using the $\beta $ kernel density function ($\beta $-EQMOM) is used to approximate solutions to the evolution equation for univariate and bivariate composition probability distribution functions (PDFs) of a passive scalar for binary and ternary mixing. The key element of interest is the molecular mixing term, which is described using the Fokker–Planck (FP) molecular mixing model. The direct numerical simulations (DNSs) of Eswaran and Pope [“Direct numerical simulations of the turbulent mixing of a passive scalar,” Phys. Fluids **31**, 506 (1988)] and the amplitude mapping closure (AMC) of Pope [“Mapping closures for turbulent mixing and reaction,” Theor. Comput. Fluid Dyn. **2**, 255 (1991)] are taken as reference solutions to establish the accuracy of the FP model in the case of binary mixing. The DNSs of Juneja and Pope [“A DNS study of turbulent mixing of two passive scalars,” Phys. Fluids **8**, 2161 (1996)] are used to validate the results obtained for ternary mixing. Simulations are performed with both the conditional scalar dissipation rate (CSDR) proposed by Fox [*Computational Methods for Turbulent Reacting Flows* (Cambridge University Press, 2003)] and the CSDR from AMC, with the scalar dissipation rate provided as input and obtained from the DNS. Using scalar moments up to fourth order, the ability of the FP model to capture the evolution of the shape of the PDF, important in turbulent mixing problems, is demonstrated. Compared to the widely used assumed $\beta $-PDF model [S. S. Girimaji, “Assumed *β*-pdf model for turbulent mixing: Validation and extension to multiple scalar mixing,” Combust. Sci. Technol. **78**, 177 (1991)], the $\beta $-EQMOM solution to the FP model more accurately describes the initial mixing process with a relatively small increase in computational cost.

## I. INTRODUCTION

The spatial and temporal evolution of a passive scalar in a turbulent mixing process can be described using an evolution equation for the probability density function (PDF) of the system composition.^{1,2} Solutions of the equation for the composition PDF provide a one-point statistical description of the behavior of the scalar composition field, during its evolution from the unmixed to the mixed state.^{3–9} Due to the importance of turbulent mixing in industrial and environmental applications, several methods have been developed to obtain such solutions.

Before summarizing the most common approaches used to obtain the composition PDF from its evolution equation, it is worth reminding that the most accurate approach for modeling turbulent mixing is direct numerical simulation (DNS),^{9–13} which consists in directly discretizing the physical space and solving the transport equation for the passive scalar along with the Navier-Stokes equation.^{1} However, this method is computationally prohibitive when applied to practical problems, due to the requirements in terms of spatial resolution, which needs to explicitly resolve the smallest of the Kolmogorov and Batchelor scales. Other strategies include large-eddy simulation, such as those discussed in Refs. 14–16. In the context of PDF methods, which aim at finding approximate solutions for the evolution equation of the composition PDF, Lagrangian PDF methods^{1,2,17,18} are computationally more affordable than DNS, but still computationally intensive for many practical applications at the industrial and environmental scales. Stochastic field methods were proposed,^{19–27} in which the evolution of the composition PDF is approximated using a set of Eulerian stochastic fields, defined over the entire computational domain. The spatio-temporal evolution of these fields is regulated by stochastic partial differential equations (PDEs). A number of PDEs equal to the product between the number of composition variables and the number of stochastic fields need to be solved. For practical applications, Reynolds-averaged Navier-Stokes (RANS) simulations are particularly attractive to describe turbulent mixing processes due to their low computational cost. In RANS methods, the evolution equation of the composition PDF is closed by introducing molecular mixing models, which are defined in terms of one-point turbulence statistics. As it will be illustrated later, these models typically require the scalar dissipation rate to be modeled, which may not be a trivial task in general non-homogeneous problems.

In many applications, a presumed $\beta $-PDF model is employed to describe the mixture fraction,^{6} which requires the solution of two RANS equations: one for the mixture-fraction mean and one for its variance. To go beyond the $\beta $-PDF model, the Fokker–Planck (FP) model^{2,28,29} can be used to close the molecular mixing term in the composition PDF transport equation. The resulting nonlinear FP equation must then be solved, which is a nontrivial task.

The Fokker-Planck equation was developed by Fokker^{30} and Planck,^{31} to describe the time evolution of the velocity PDF of a particle under the effect of random motions and drag forces. Kolmogoroff^{32} also obtained the equation independently, while Smoluchowski^{33} obtained the equivalent of the Fokker-Planck equation for the particle position in Brownian motion.

Analytical solutions for FP equations are possible only in a limited set of cases.^{31} As a consequence, numerical methods have been developed to find approximate solutions in cases of general interest. Spencer and Bergman^{34} and Kumar and Narayanan^{35} used the finite-difference method to directly discretize the spatial and composition spaces. As with Lagrangian PDF methods, due to the transient nature of mixing problems and to the dimensionality of the composition PDF, this approach still carries a significant computational cost, which makes it unattractive for many applications. As a consequence, alternative approaches have been developed to compute moments of the composition PDF, rather than to capture its entire evolution. These methods are based on the idea that the quantities of interest in most applications can be related to the calculated moments. However, when applied to reactive scalars with nonlinear chemical source terms, integrals with respect to the unknown composition PDF must be evaluated.^{2,14} For this purpose, an efficient and accurate method for reconstructing the composition PDF from its moments is required.^{36}

In this work, the extended quadrature method of moments (EQMOM) with $\beta $ kernel density functions^{37} is adopted to approximate solutions to the statistically homogeneous composition PDF equation, which is closed using the FP model for molecular mixing. EQMOM is an extension of the quadrature method of moments (QMOM), initially proposed by McGraw^{38} for problems involving aerosols, and later applied to chemical engineering problems by Marchisio *et al.*^{39} In the EQMOM approach,^{37} the composition PDF is approximated by a discrete sum of positive kernel density functions (KDF), instead of the Dirac delta functions used in QMOM. This allows a continuous reconstruction of the PDF to be obtained from a finite set of integer moments. The choice of the functional form of the kernel density function is based on the support of the PDF, whose moments are to be reconstructed, and on the explicit knowledge of the recurrence relation of the polynomials orthogonal to such a kernel density function. In particular, if the PDF is defined on the compact support [0, 1], the $\beta $ density function is an appropriate choice, while for a PDF defined on the positive real like, the gamma or log-normal density functions would be adequate.

It is worth noting that, differently from the direct QMOM^{40} used in multi-environment PDF models,^{2,41,42} which directly computes the quadrature weights and abscissae associated with the moments of the composition PDF, EQMOM solves for the moments of the PDF, which are conserved quantities in the calculation. For spatially inhomogeneous cases, this feature allows a realizable finite-volume scheme to be formulated, as illustrated in the work of Vikas *et al.*^{43,44} While we consider homogeneous problems in this case, a brief discussion of the extension to non-homogeneous problems is provided in the Conclusions. Hereinafter, the turbulent mixing model resulting from the application of $\beta $-EQMOM to the composition PDF equation, closed with the FP model, will be referred to as the $\beta $-EQMOM-FP model.

The primary objectives of the work are twofold. First, the version of the FP molecular mixing model proposed by Fox^{2} is validated against DNS for binary^{9} and ternary^{45} mixing. To the best of our knowledge, this step has not been previously reported for the FP model. In doing so, we demonstrate that the FP model for binary mixing is essentially identical to the amplitude mapping closure (AMC),^{46} including at very short times, where the $\beta $-PDF model yields poor predictions. Second, the $\beta $-EQMOM approach is shown to yield accurate approximations for the composition PDF equation using as few as four moments (i.e., $\u27e8\varphi \u27e9$, $\u27e8\varphi 2\u27e9$, $\u27e8\varphi 3\u27e9$, $\u27e8\varphi 4\u27e9$) for binary mixing. This result opens the possibility of applying the $\beta $-EQMOM-FP model in practical computational fluid dynamic (CFD) simulations of complex flows as an alternative to the widely used assumed $\beta $-PDF model. In addition, we demonstrate that the $\beta $-EQMOM-FP model can be used for ternary mixing problems for which assumed PDF models are not available.^{2}

The remainder of the paper is structured as follows. In Sec. II the evolution equation for the composition PDF is reported, and the FP molecular mixing model is summarized. The $\beta $-EQMOM approach is introduced in Sec. IV, while the moment transport equations and the corresponding source terms originating from the equation for the composition PDF are presented in Sec. III. The AMC, used as a reference to validate the FP model, is summarized in Sec. V. Next, the $\beta $-EQMOM-FP approach is validated in Sec. VI. In particular, Sec. VI C shows the validation results for the binary mixing cases of Eswaran and Pope^{9} using the conditional scalar dissipation rate (CSDR) proposed by Fox,^{2} and Sec. IV D presents the corresponding results obtained using the CSDR from AMC. Section IV E concludes with the application of the $\beta $-EQMOM-FP approach to the ternary mixing case of Juneja and Pope.^{45}

## II. EVOLUTION EQUATION FOR THE COMPOSITION PDF

The evolution equation of the composition PDF $f\varphi \psi ;x,t$ for nonreacting, passive scalars, assuming constant density, reads^{1,2}

where repeated indices imply summation. The *i*th composition variable is denoted by $\varphi i$, and its phase-space component by $\psi i$. Fluctuations about the mean are denoted with a prime: $\varphi i\u2032$. In this work, we consider only statistically homogeneous cases for which the convection in physical space is null. The molecular diffusivity of species *i* is denoted by $\Gamma i$. For statistically homogeneous systems, $\u22072\varphi i$ is null, leaving only the first term on the right-hand side of (1). In order to close this term, a molecular mixing model is needed. Two models will be considered in this work: the interaction-by-exchange-with-the-mean (IEM) model^{47,48} and the FP model.^{2,28,29,49}

The IEM model is the simplest mixing model as it assumes a linear relaxation of the concentration towards its mean.^{1,2} The generalized IEM expression for the conditional diffusion is

where $\mathit{\epsilon}ij=(\Gamma i\Gamma j)1\u22152\u2207\varphi i\u2032\u22c5\u2207\varphi j\u2032$ is the scalar covariance dissipation rate (SDR) for component *i*, *j*, $\varphi i$ is the scalar mean for component *i*, and $\varphi j\u2032\varphi k\u2032\u22121$ is the *i*, *j* component of the inverse covariance matrix of the scalars.^{2} In applications, it is almost always assumed that the time scales $1\u2215\tau ii=\u2211j\mathit{\epsilon}ij\varphi j\u2032\varphi i\u2032\u22121$ are the same for all components *i*. Moreover, it is also assumed that the scalars are independent so that $\u2211j\mathit{\epsilon}ij\varphi j\u2032\varphi k\u2032\u22121=0$ when $i\u2260k$.^{1,2} While the IEM model has been widely used, in its original formulation,^{48} it does not account for differential diffusion (see Ref. 50 for an extension including this effect), and it assumes that all the scalars mix with a single time scale. As a consequence, in the absence of gradients in the scalar mean, the shape of the scalar PDF does not evolve with time.

In order to capture the change in shape, Fox^{28,29} introduced the FP molecular mixing model, as a possible solution to the limitations of the IEM model. We observe that when the $\Gamma i$ are equal, the conditional diffusion can be related to the joint CSDR $\mathit{\epsilon}ij|\psi =(\Gamma i\Gamma j)1\u22152\u2207\varphi i\u2032\u22c5\u2207\varphi j\u2032|\psi $ by^{2}

However, Eq. (3) leads to negative diffusion in phase space, as discussed by Pope.^{1} To avoid this problem, in the FP model the conditional diffusion is modeled as^{2}

where *c*_{FP} is a positive constant that controls the diffusive relaxation rate for the shape of the PDF. The generalized IEM model in (2) is then used to close the first term on the right-hand side of Eq. (4),^{2}

This formulation yields the generalized IEM model in the case of *c*_{FP} = 0. Here, *c*_{FP} = 10 is chosen to ensure an accurate reproduction of the results obtained in comparison with AMC (see Sec. IV B for details).

Eliminating the terms representing transport in physical space, the evolution equation for a univariate composition PDF can be written as

where $f\varphi =f\varphi \psi ;t$ and $R\psi ;t$ is the source term due to micromixing. The models for $R\psi ;t$ using the IEM and FP models are, respectively,

where the SDR $\mathit{\epsilon}\varphi $ is the expected value of the CSDR, which we indicate with $\mathit{\epsilon}\varphi |\psi $. The form of the CSDR must be provided by the user to close the FP model.^{2} In Sec. III A, we consider two possible closures based on the $\beta $ PDF^{2} and AMC.^{46}

In the case of a bivariate composition PDF for the mixture-fraction vector $\xi $ with phase-space vector $\zeta $,^{2} the statistically homogeneous version of Eq. (1) with the mixing term closed using the FP model yields

where $\mathit{\epsilon}ij$ is the expected value of the joint CSDR $\mathit{\epsilon}ij|\zeta 1,\zeta 2$. The components of the inverse covariance matrix are defined by

As in the univariate case, a closure is required for the joint CSDR. Fox^{2} proposed the following closure:

where

As described in the work of Fox,^{2} this closure ensures that the joint PDF found by solving (9) is nonzero only inside the phase-space triangle defined by $0\u2264\zeta i$ and $0\u2264\zeta 1+\zeta 2\u22641$ (i.e., the phase-space diffusive flux is null in the direction normal to the boundary). The three SDR components in (12) must be supplied by a separate model for the scalar covariance dissipation rate. By definition, the $2\xd72$ matrix with components $\mathit{\epsilon}ij$ has to be symmetric and non-negative. In Sec. IV E we relate $\xi $ to the two scalars ($\varphi 1,\varphi 2$) used in the ternary mixing case of Juneja and Pope,^{45} for which the covariance and the cross SDR (i.e., $\mathit{\epsilon}12$) are null.

## III. TRANSPORT EQUATIONS FOR THE MOMENTS OF THE COMPOSITION PDF

Transport equations for the moments are now obtained for the univariate and bivariate cases considered in the present work.

### A. Univariate case

The moment of order *m* of the univariate PDF $f\varphi \psi ;t$ is defined as

Multiplying both sides of Eq. (6) by $\psi m$ and integrating with respect to $\psi $, we obtain

Observing that differentiation and integration commute on the left-hand side of Eq. (14), and indicating with *R*_{m}(*t*) the moment of order *m* of $R(\psi ,t)$,

the transport equations for the moments of $f\varphi \psi ;t$ are obtained,

The expression for the source term obtained from the IEM model is

where $\varphi \u20322=\varphi 2\u2212\varphi 2$. Note that this expression is closed for the integer moments up to order $m\u22650$.

The form of the source term *R*_{m,FP} for the FP model depends on the choice for the functional form of the CSDR. As this choice essentially determines the predicted shape of the scalar PDF, it is the most important part of the FP model. Based on the $\beta $ PDF with support [0, 1], Fox^{2} proposed

which leads to a source term for the moment transport equations given by

It is important to note that (19) is closed for the moments up to order $m\u22650$. (The cases *m* = 0, 1 are trivial.) In other words, the predicted values of the moments will not depend on how the scalar PDF is reconstructed from the moments (e.g., using $\beta $-EQMOM). The case of a scalar PDF with support [−1, 1] is treated by applying a linear transformation to map the support to [0, 1], in which case the numerator of (18) is $1\u2212\psi 2$.^{2} In Sec. V we introduce an alternative choice to (18) based on the AMC.^{46}

### B. Bivariate case

The moments of the joint mixture-fraction PDF $f\xi \zeta ;t$, in which $\zeta =(\zeta 1,\zeta 2)\u2208\Omega ={0\u2264\zeta i,\u20090\u2264\zeta 1+\zeta 2\u22641}$, are defined as

Multiplying both sides of Eq. (9) by $\zeta 1m\zeta 2n$, and integrating with respect to $\zeta 1$ and $\zeta 2$, yields the transport equations for the moments of $f\xi \zeta 1,\zeta 2;t$,

where $Rmnt$ is the moment of the source term

For the FP model with the joint CSDR closure in (11), this yields

where $\xi i\u2032\xi j\u2032\u22121$ are the inverse components defined in (10). If we denote the maximum moment order by *K* (i.e., $0\u2264n+m\u2264K$), it can be observed that (23) is closed for the set of bivariate moments up to order *K*. Thus, as in the univariate case in (19), the FP model solution for the bivariate moments does not depend on how the joint mixture-fraction PDF is reconstructed.

If the time scales $1\u2215\tau i=\u2211j\mathit{\epsilon}ij\xi j\u2032\xi i\u2032\u22121$ are assumed to be equal to $\tau \varphi $, then $\mathit{\epsilon}12=\xi 1\u2032\xi 2\u2032\u2215\tau \varphi $ in the standard model. The time scale $\tau \varphi $ must then be related to the turbulence time scales to close the FP model. In this work, we either take $\tau \varphi $ directly from the DNS data or simply use it as a fixed parameter to rescale the time *t* in the moment evolution equations.

## IV. EXTENDED QUADRATURE METHOD OF MOMENTS

The $\beta $-EQMOM^{36,37,51–53} is used in this study to compute an approximate solution to the univariate PDF transport equation. In this method, a univariate PDF is approximated by a weighted sum of non-negative KDFs $\delta \sigma \psi ,\psi i$,

where $w\alpha $ and $\varphi \alpha $ are, respectively, the weights and abscissae associated with the Gaussian quadrature computed from the first 2N transported moments. The right-hand side of Eq. (24) contains 2N + 1 unknowns: N weights $w1,w2,\u2026,wN$, N abscissae $\varphi 1,\varphi 2,\u2026,\varphi N$, and the parameter $\sigma $, assumed to be common for all KDFs.^{37} The value of $\sigma $ is determined by enforcing that the last moment of a set with an odd number of moments (a total of 2N + 1 moments) is preserved by the approximated PDF.^{53}

The choice of the KDF depends on the support of the PDF to be reconstructed and on the ease of computation of the orthogonal polynomials associated with the KDF. Previous works on the validation of the assumed PDF approach in turbulent mixing^{54,55} compared the $\beta $-PDF with the mixture-fraction PDF reconstructed from DNS in both homogeneous flow and turbulent shear flow. In general, it was found that the assumed $\beta $-PDF gives a reasonable approximation of the mixture-fraction PDF, especially when the mixture-fraction variance is much less than one. We will then use the $\beta $ KDF^{37} to formulate the $\beta $-EQMOM closures used in this work. In theory, any other KDF with support [0, 1] could be used for this purpose.^{37}

The $\beta $ KDF is defined as

where $\lambda \alpha =\varphi \alpha \u2215\sigma $, $\mu \alpha =1\u2212\varphi \alpha \u2215\sigma $, and

represents the beta function.

A linear transformation allows this distribution to be defined on the arbitrary compact support $a,b$.^{56} The integer moment of order *k* of the $\beta $ KDF is

with *M*_{0} = 1. Thus, the integer moments of the scalar PDF to be approximated as in (24) can be written as

where

Equation (28) defines a lower triangular system of equations to find the 2N + 1 unknowns, i.e., N weights, N abscissae, and $\sigma $,

where the non-negative coefficients $\chi k$ depend only on $\sigma $, and $Mk*=\u2211\alpha =1Nw\alpha \varphi \alpha k$. The improved version^{53} of the iterative procedure of Yuan *et al.*^{37} is used to compute the values of the quadrature weights, abscissae, and of the parameter $\sigma $ from the 2N + 1 transported scalar moments (1, $\varphi $, $\varphi 2$, $\u2026$, $\varphi 2N$). Note that fixing N = 1 reduces $\beta $-EQMOM to the assumed $\beta $-PDF approximation.^{6}

As illustrated in Sec. IV E, the bivariate cases considered in this work are reduced to the application of the univariate $\beta $-EQMOM procedure for each of the scalars, argument of the joint composition PDF, by means of a linear mapping. This is possible because only the marginal composition PDFs, and not the joint composition PDF, are reconstructed in the bivariate case concerning ternary mixing.

## V. MAPPING CLOSURE

Mapping closures were initially developed by Kraichnan^{57,58} and Chen *et al.*^{59} and were applied to model turbulent scalar mixing. Pope^{46} used the same approach to study the evolution of an inert scalar in isotropic turbulence decaying from a PDF initially set to be constituted by a double delta function. He compared the results obtained from the solution of the mapping closure with those obtained from the DNS of Eswaran and Pope^{9} and demonstrated the agreement between the two models. An analytical expression for the PDF is^{46}

where

As pointed out by Pope,^{46} the time *t* enters into (31) only through the scalar variance. In other words, for a fixed scalar variance the shape of the PDF does not depend on *t*. Similar behavior was observed by Eswaran and Pope.^{9} In the literature, the AMC in (31) is considered to be the best model for binary mixing in the case of homogeneous turbulence, and hence it has often been used for validating other molecular mixing models.^{17,18}

Using the AMC, the CSDR can be expressed as^{2}

where

The expression for *R*_{m} in Eq. (16) obtained from the AMC model is

where *N*_{N}(*t*) is found by substituting the quadrature representation of the PDF into Eq. (34),^{37}

with $Ni\u2265N+1$. The secondary weights $wij$ and abscissae $\psi ij$ are the Gauss-Jacobi quadrature corresponding to the $\beta $ KDF, which can be computed using the methodology described in the work of Yuan *et al.*^{37} The value of N_{i} is chosen large enough that the numerical value of $NN$ is insensitive to N_{i}. Compared to (19), the choice of the AMC CSDR leads to an unclosed moment system due to (34). We are thus interested to know which choice (18) or (33) performs best in the FP model.

## VI. RESULTS AND DISCUSSION

In order to investigate its capability to capture the known shape of the PDF, and to determine the required number of KDFs, $\beta $-EQMOM is applied first to reconstruct the composition PDF using the moments from AMC. After establishing this capability, a set of simulations is then performed using $\beta $-EQMOM-FP to determine the effect of the model constant *c*_{FP}, in order to match the AMC PDF. Finally, the $\beta $-EQMOM-FP model is validated against the results of the DNS of binary mixing reported in the work of Eswaran and Pope,^{9} and of ternary mixing reported in the work of Juneja and Pope.^{45}

### A. Verification of $\beta $-EQMOM

The first step is to verify that the $\beta $-EQMOM is able to accurately reconstruct the PDF obtained with AMC in Sec. V and to determine the number of KDFs required to achieve this result. The reader should note that this step is of pure numerical interest, and does not involve the FP closure model, but simply the $\beta $-EQMOM procedure. Moments of the AMC PDF are computed and provided as input to $\beta $-EQMOM, and the number N of KDFs is varied to determine its effect on the shape of the reconstructed PDF. Figure 1 shows the results obtained using from one to four KDFs. It is possible to observe that with N = 1, the agreement at the initial time is poor, but $\beta $-EQMOM closely matches the AMC PDF at later times. As noted earlier, the latter is expected because both the AMC and the $\beta $ PDF approach a Gaussian PDF for large times. Results with $N\u22652$ show that the reconstructed PDFs agree closely with the AMC PDFs at all times. These results demonstrate that $\beta $-EQMOM is capable of reconstructing the shape of the PDFs of interest in the binary mixing problems given a known set of moments as input. Moreover, as the difference between N = 2 and 3 is negligible, it is possible to conclude that for binary mixing only two KDFs are required to accurately reproduce the AMC PDF. Thus, CFD models for turbulent mixing can make use of the FP model by solving for as few as four mixture-fraction moments: $\u27e8\xi \u27e9$, $\xi 2$, $\xi 3$, $\xi 4$. Compared to the assumed $\beta $-PDF model, the two additional moments are needed to capture the shape of the PDF when the variance is near its maximum value (i.e., the PDF is close to two Dirac delta functions) as seen in Fig. 1.

### B. Choice of the value of the diffusive relaxation rate constant

The constant *c*_{FP}, which appears in the FP model, controls the diffusive relaxation rate of the PDF. In this section, we obtain solutions with different values of *c*_{FP} and compare the results with those obtained with the AMC. Figure 2 shows that as *c*_{FP} increases, the distributions are comparable to those obtained with AMC. In particular, a value of *c*_{FP} = 10 or higher yields results in good agreement with AMC. For smaller values of *c*_{FP}, the PDF is closer to that predicted with the IEM model. As *c*_{FP} fixes the time scale for diffusive relaxation, it is reminiscent of the “generation time” used in the multiple mapping closure (MMC) as described by Sundaram *et al.*^{60} Here, as in MMC, it is necessary to consider two micromixing time scales: $\tau \varphi $ and $\tau \varphi \u2215cFP$. Thus, if the FP model is used for a reactive scalar, a second Damköhler number will arise to describe the time scale ratio between diffusive relaxation and chemical reactions. It can be anticipated that this new Damköhler number will better characterize ignition/extinction phenomena than the one based on integral-scale mixing $\tau \varphi $.

It is worth noting that the moments $\varphi $ and $\varphi 2$ do not depend on *c*_{FP}. Thus, changing its value only affects moments of orders higher than two. In general, in order to recover the AMC PDF, *c*_{FP} must be set large enough that it no longer affects these higher-order moments. In fact, if the only requirement were to recover the AMC PDF, then we can set $cFP=\u221e$ and use algebraic expressions for the third- and fourth-order moments. These observations are confirmed in Fig. 3, which shows for large times the scaled fourth-order central moment relaxes to 3, which is the value for the Gaussian distribution. The curves for different *c*_{FP} overlap for values higher than 5. For pure mixing problems like the ones considered in this work, the value of *c*_{FP} can be set arbitrarily large. However, when additional physics are included in the PDF transport equation, a more precise definition will be needed for fixing its value relative to other time scales. The theoretical estimates developed for the MMC^{60} may be helpful in this respect. Nonetheless, as the AMC PDF holds for purely diffusive problems,^{28,29} one cannot link *a priori* the value of *c*_{FP} to the turbulence time scales without considering the length scales of molecular diffusion. For fully developed turbulent flows for which $\tau \varphi \u221dk\u2215\mathit{\epsilon}$ where *k* is the turbulent kinetic energy and $\mathit{\epsilon}$ its dissipation, if the Schmidt number is greater than one, then the diffusion time scale would yield $cFP\u221dk\u2215\mathit{\epsilon}\nu $, where $\nu $ is the kinematic viscosity of the fluid.

### C. Validation of $\beta $-EQMOM-FP model against DNS of binary mixing

Three cases of mixing in homogeneous isotropic turbulence, with different time evolutions of the scalar dissipation rate (Fig. 4), are now considered to validate the $\beta $-EQMOM-FP model against the DNS results of Eswaran and Pope.^{9} In these cases, the initial scalar fields constitute of two Dirac delta functions, located at points −1 and 1 in composition space, to represent a binary mixing problem, characterized by the presence of the scalar at two different states. The time-dependent SDR [$\mathit{\epsilon}\varphi (t)$] obtained from the DNS Eswaran and Pope,^{9} for the three cases, is provided as input to the $\beta $-EQMOM-FP model. The evolution of $\mathit{\epsilon}\varphi $ as a function of the dimensionless time $\tau $ is shown in Fig. 4 for the three cases considered here. The SDR is combined with the moment source term in (19) to solve for the closed system of moments $\varphi m$ for $m\u2208{1,2,3,4}$. The parameter $cFP\u2009=\u200910$ was chosen to reproduce the AMC behavior as described in Sec. IV B.

Figure 5 shows the scalar PDF reported by Eswaran and Pope^{9} compared to the results obtained using $\beta $-EQMOM-FP at different values of the scalar variance $\Phi $, for case 1. Two KDFs (N = 2) were used to obtain these results. It is apparent that the results provided by $\beta $-EQMOM-FP after the first few time steps differ from those of the DNS.^{9} This difference, however, can be explained by noting that the initial condition defined in the DNS simulations^{9} is not an exact pair of Dirac delta functions, due to the difficulties of generating this configuration in a discretized space. However, with $\beta $-EQMOM-FP, as with AMC, it is possible to have an exact initial definition of the scalar PDF. Further advancing the $\beta $-EQMOM-FP solution in time shows that the model predicts the evolution towards the well-mixed condition in reasonable agreement with the DNS results and converges towards a Gaussian PDF [Fig. 5(e)].

Figures 6 and 7 show the results obtained for Cases 2 and 3, respectively. In both cases, the $\beta $-EQMOM-FP model satisfactorily predicts the evolution of the PDF, and its relaxation to a Gaussian distribution, as it is observed in the DNS. The satisfactory performance of the $\beta $-EQMOM-FP model for binary mixing is a direct reflection of its agreement with AMC as seen in Sec. IV B. We show next that the quality of the predictions depends on the choice of the CSDR.

### D. $\beta $-EQMOM-FP model of binary mixing with the AMC CSDR

The solution for the three cases considered in Sec. VI C using the AMC CSDR found by solving (35) is considered here. Two KDFs (N = 2) and eight secondary nodes (N_{i} = 5) were used to study the evolution of the PDF with the AMC CSDR. Figure 8 shows the results obtained using the AMC CSDR for Cases 1 and 2. Using the AMC CSDR, the PDF relaxes to a Gaussian distribution as expected. However, the agreement with the DNS and AMC is not as good as with the CSDR given in (18). It is worth mentioning that increasing the number of KDFs does not improve the predicted shape of the PDF. Quadrature methods are accurate for $Ni\u22652N+1$ secondary nodes^{37} and increasing N_{i} does not affect the results. Different values for $cFP$ were tested, and smaller values make the result closer to the IEM model while values of 10 or higher do not improve the prediction of the PDF. Using AMC CSDR, the moment transport equations must be closed using the quadrature representation of the PDF. This increases the cost of computation in the sense that an inversion is required at each time step to compute the source terms.

### E. Validation of $\beta $-EQMOM-FP model against DNS results for ternary mixing

In realistic problems, such as in chemical reacting flows, it is not infrequent to have reactants introduced in more than two inlets. In order to show that the $\beta $-EQMOM-FP model is capable of describing such cases, we consider the ternary mixing problem studied by Juneja and Pope,^{45} who reported the evolution of the joint PDF $f\varphi (\psi 1,\psi 2;t)$ of the two scalars. We will use their results for the marginal PDFs to validate the results obtained with $\beta $-EQMOM-FP, starting from the bivariate mixture-fraction model introduced in Sec. III B.

Juneja and Pope^{45} considered two scalars with an initial joint PDF constituted by a triple delta function, corresponding to blobs of fluids at three distinct states. Figure 9(b) shows the bivariate phase space. The three initial states correspond to the vertices of an equilateral triangle, and have equal probability. Consequently, the means $\varphi 1$ and $\varphi 2$ are zero, which ensures the three-way symmetry of the initial state in composition phase space. The setup is designed so that the bivariate PDF relaxes to zero for both $\psi 1$ and $\psi 2$. Additionally, the two scalars are uncorrelated so that $\varphi 1\varphi 2=0$,^{45} and hence the cross SDR is also zero when written in terms of $\varphi $.

In order to solve the problem using the mixture-fraction PDF, the composition phase space $\psi $ is mapped linearly onto the mixture-fraction phase space $\zeta $, so that the matrix for the joint CSDR introduced in the work of Fox,^{2} defined in terms of $\zeta $, can be used. A graphical representation of this mapping is illustrated in Fig. 9. The linear mapping is defined by the following linear algebraic equation:

After simplifying, the relationship between the components of $\zeta $ and $\psi $ is found

Note that using these transformation it is straightforward to write the moments $\varphi 1m\varphi 2n$ in terms of the moments $\xi 1m\xi 2n$. Also note that $\xi 1=\xi 2=1\u22153$, and the mixture-fraction variances are $\xi 1\u20322=13\varphi 12+19\varphi 22$ and $\xi 2\u20322=49\varphi 22$, while the mixture-fraction covariance $\xi 1\u2032\xi 2\u2032=\u221229\varphi 22$ is negative for this example.

In order to employ the FP model for the mixture-fraction moments in (23), we must relate the SDR components defined in terms of $\xi $ to those defined in terms of $\varphi $ by using the second part of (38). The transformed joint scalar dissipation matrix can be written as

where $\mathit{\epsilon}\varphi 11\u2009=\u2009\varphi 12\tau \varphi $ and $\mathit{\epsilon}\varphi 22\u2009=\u2009\varphi 22\tau \varphi $ are the SDR components for $\varphi $, and the scalar mixing time $\tau \varphi $ is found from the DNS. Note that (39) implies that $\u2211j\mathit{\epsilon}ij\xi j\u2032\xi k\u2032\u22121=1\tau \varphi \delta i,k$ in (9), which is just the single scalar time scale model used in the IEM model.^{2}

Univariate moments defined over the space $\psi $ are obtained from the bivariate moments over $\zeta $ by means of Eq. (38), leading to

By expanding the right-hand side of (40), the univariate moments needed for $\beta $-EQMOM-FP reconstruction of the two marginal PDFs for $\varphi 1$ and $\varphi 2$ can be related to the bivariate moments of mixture fraction. In this study N = 3 KDFs are used to reconstruct each of the two marginal PDFs. The number of primary quadrature nodes was chosen to be able to exactly represent the initial condition, which is constituted by three Dirac delta distributions. As a consequence of this choice of N, 27 joint mixture-fraction moments up to sixth order are computed,

The initial conditions for this moment set are easily found from the known bivariate PDF. Recall from Sec. III B that the 27 moment transport equations are closed. Thus, their values do not depend on how the PDFs are reconstructed from the moments. For the results presented here, the diffusive relaxation parameter is again $cFP=10$.

Figure 10 shows the results for the marginal PDFs of $\varphi 1$ and Fig. 11 for $\varphi 2$ obtained with $\beta $-EQMOM-FP, compared to the DNS.^{45} The marginal PDFs are evolving and finally converge to Gaussian distributions. The marginal PDF of $\varphi 2$ obtained using $\beta $-EQMOM-FP with unequal initial conditions is in very good agreement with that obtained by DNS. As the scalar $\varphi 2$ represents binary mixing with unequal volumes, this result demonstrates the ability of the FP model to treat such cases. For $\varphi 1$, $\beta $-EQMOM-FP captures the peak in the middle for the marginal PDF; however, it is not as sharp as the one obtained by the DNS. This is a result of the $\beta $-EQMOM reconstruction methodology not accounting for the fact that the marginal PDF is a projection of the bivariate PDF for the correlated variables $\xi 1$ and $\xi 2$. Thus, in order to improve the reconstruction, it would be necessary to reconstruct the bivariate mixture-fraction PDF using a method such as CQMOM.^{61,62} However, before this can be done in a general way, further work will be needed to extend CQMOM to bounded convex domains in phase space such as those in Fig. 9. Given that binary mixing required moments up to fourth order, it is likely that a bivariate CQMOM using fourth-order moments will be successful for most applications, thereby reducing the number of moments in (41) that would have to be computed for ternary mixing.

## VII. CONCLUSIONS

For the first time, the extended quadrature method of moments with $\beta $ KDFs has been applied to find approximate solutions for the composition PDF using the FP molecular mixing model.^{2,28} This approach allows any form of the conditional scalar dissipation rate $\u27e8\mathit{\epsilon}\varphi |\psi \u27e9$ to be used in the calculations. Differently from models relying on a $\beta $-PDF to approximate the mixture-fraction PDF, it is able to represent any form of the scalar PDF, due to its capability of reproducing the moments of arbitrary distributions. For example, here we have shown that the mixture-fraction PDF can be better captured using moments up to fourth order.

Two mixing cases were used to validate the $\beta $-EQMOM-FP against results obtained with the DNS.^{9,45} For binary mixing, the initial scalar PDF is specified as two delta functions at different locations in composition space. Results show that the $\beta $-EQMOM-FP predicts the evolution of scalar PDF in agreement with AMC^{46} and the DNS results in the cases examined. For ternary mixing, the capability of the $\beta $-EQMOM-FP to predict the time evolution of the moments of the bivariate scalar PDF was demonstrated considering the DNS of Juneja and Pope.^{45} Here, the $\beta $-EQMOM-FP approach predicted the evolution of the marginal PDFs towards Gaussian. However, for the true ternary mixing problem with three delta functions for the univariate PDF, the middle peak was not as well captured. Nonetheless, the applicability of the $\beta $-EQMOM-FP approach to solve turbulent mixing problems has been demonstrated for the first time.

Finally, the closure for the conditional scalar dissipation rate in the FP model was shown to have a strong effect on the predictions, and it was found that the simplest closure produces the best results. This conclusion was initially verified by reconstructing the PDF obtained with known moments from the AMC. In this manner, it was demonstrated how the diffusive relaxation time scale in the FP model affects the solution.

Future work should focus on applying quadrature-based moment methods (QBMM) with the PDF transport equation to investigate turbulent reacting scalars. In principle, once the PDF equation has been closed (e.g., using the FP molecular mixing model), its solution can be approximated using QBMM.

The extension of the proposed approach to non-homogeneous problems can be achieved by keeping the spatial transport terms in the moment transport equations obtained from (1). These will consist of a diffusion and an advection term, in physical space, for each moment. The solution of the moment transport equations obtained in this way has to be performed with a numerical scheme that ensures the realizability of the transported moment set. As discussed in Ref. 44, the discretization of the moment diffusion term does not present difficulties, if a traditional second-order finite-volume scheme is used. More difficulties are faced when dealing with the discretization of the moment advection term, which is known to be problematic,^{63} because it may compromise the realizability of the set of transported moments, when using numerical schemes other than the first order upwind. In the framework of quadrature-based moment methods, the closure of the advection term can be achieved using a kinetic-based flux.^{64–66} In these schemes, the PDF reconstructed by means of the quadrature approximation is used to close the convective term of the moment transport equations, as illustrated in Refs. 65 and 66. These schemes were extended to higher order of discretization in Refs. 43 and 67, where the technical details of the discretization of the moment transport equations are provided.

## ACKNOWLEDGMENTS

The authors gratefully acknowledge the support of the US National Science Foundation under the SI^{2}–SSE Award No. NSF–ACI 1440443.