Application of the Fokker-Planck molecular mixing model to turbulent scalar mixing using moment methods

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.¹ 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 $β$-PDF model is employed to describe the mixture fraction,⁶ which requires the solution of two RANS equations: one for the mixture-fraction mean and one for its variance. To go beyond the $β$-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³⁰ and Planck,³¹ to describe the time evolution of the velocity PDF of a particle under the effect of random motions and drag forces. Kolmogoroff³² also obtained the equation independently, while Smoluchowski³³ 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.³¹ As a consequence, numerical methods have been developed to find approximate solutions in cases of general interest. Spencer and Bergman³⁴ and Kumar and Narayanan³⁵ 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.³⁶ 

In this work, the extended quadrature method of moments (EQMOM) with $β$ kernel density functions³⁷ 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³⁸ for problems involving aerosols, and later applied to chemical engineering problems by Marchisio et al.³⁹ In the EQMOM approach,³⁷ 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 $β$ 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⁴⁰ 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 $β$-EQMOM to the composition PDF equation, closed with the FP model, will be referred to as the $β$-EQMOM-FP model.

The primary objectives of the work are twofold. First, the version of the FP molecular mixing model proposed by Fox² is validated against DNS for binary⁹ and ternary⁴⁵ 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),⁴⁶ including at very short times, where the $β$-PDF model yields poor predictions. Second, the $β$-EQMOM approach is shown to yield accurate approximations for the composition PDF equation using as few as four moments (i.e., $⟨ϕ⟩$⁠, $⟨ϕ2⟩$⁠, $⟨ϕ3⟩$⁠, $⟨ϕ4⟩$⁠) for binary mixing. This result opens the possibility of applying the $β$-EQMOM-FP model in practical computational fluid dynamic (CFD) simulations of complex flows as an alternative to the widely used assumed $β$-PDF model. In addition, we demonstrate that the $β$-EQMOM-FP model can be used for ternary mixing problems for which assumed PDF models are not available.² 

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 $β$-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 $β$-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⁹ using the conditional scalar dissipation rate (CSDR) proposed by Fox,² and Sec. IV D presents the corresponding results obtained using the CSDR from AMC. Section IV E concludes with the application of the $β$-EQMOM-FP approach to the ternary mixing case of Juneja and Pope.⁴⁵ 

The evolution equation of the composition PDF $fϕψ;x,t$ for nonreacting, passive scalars, assuming constant density, reads^1,2

where repeated indices imply summation. The ith composition variable is denoted by $ϕi$⁠, and its phase-space component by $ψi$⁠. Fluctuations about the mean are denoted with a prime: $ϕi′$⁠. 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 $Γi$⁠. For statistically homogeneous systems, $∇2ϕ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 $𝜀ij=(ΓiΓj)1∕2∇ϕi′⋅∇ϕj′$ is the scalar covariance dissipation rate (SDR) for component i, j, $ϕi$ is the scalar mean for component i, and $ϕj′ϕk′−1$ is the i, j component of the inverse covariance matrix of the scalars.² In applications, it is almost always assumed that the time scales $1∕τii=∑j𝜀ijϕj′ϕi′−1$ are the same for all components i. Moreover, it is also assumed that the scalars are independent so that $∑j𝜀ijϕj′ϕk′−1=0$ when $i≠k$⁠.^1,2 While the IEM model has been widely used, in its original formulation,⁴⁸ 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 $Γi$ are equal, the conditional diffusion can be related to the joint CSDR $𝜖ij|ψ=(ΓiΓj)1∕2∇ϕi′⋅∇ϕj′|ψ$ by² 

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

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

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ϕ=fϕψ;t$ and $Rψ;t$ is the source term due to micromixing. The models for $Rψ;t$ using the IEM and FP models are, respectively,

where the SDR $𝜀ϕ$ is the expected value of the CSDR, which we indicate with $𝜖ϕ|ψ$⁠. The form of the CSDR must be provided by the user to close the FP model.² In Sec. III A, we consider two possible closures based on the $β$ PDF² and AMC.⁴⁶ 

In the case of a bivariate composition PDF for the mixture-fraction vector $ξ$ with phase-space vector $ζ$⁠,² the statistically homogeneous version of Eq. (1) with the mixing term closed using the FP model yields

where $𝜀ij$ is the expected value of the joint CSDR $𝜖ij|ζ1,ζ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² proposed the following closure:

where

As described in the work of Fox,² this closure ensures that the joint PDF found by solving (9) is nonzero only inside the phase-space triangle defined by $0≤ζi$ and $0≤ζ1+ζ2≤1$ (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×2$ matrix with components $𝜀ij$ has to be symmetric and non-negative. In Sec. IV E we relate $ξ$ to the two scalars (⁠$ϕ1,ϕ2$⁠) used in the ternary mixing case of Juneja and Pope,⁴⁵ for which the covariance and the cross SDR (i.e., $𝜀12$⁠) are null.

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

The moment of order m of the univariate PDF $fϕψ;t$ is defined as

Multiplying both sides of Eq. (6) by $ψm$ and integrating with respect to $ψ$⁠, 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(ψ,t)$⁠,

the transport equations for the moments of $fϕψ;t$ are obtained,

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

where $ϕ′2=ϕ2−ϕ2$⁠. Note that this expression is closed for the integer moments up to order $m≥0$⁠.

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 $β$ PDF with support [0, 1], Fox² 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≥0$⁠. (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 $β$-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−ψ2$⁠.² In Sec. V we introduce an alternative choice to (18) based on the AMC.⁴⁶ 

The moments of the joint mixture-fraction PDF $fξζ;t$⁠, in which $ζ=(ζ1,ζ2)∈Ω={0≤ζi,0≤ζ1+ζ2≤1}$⁠, are defined as

Multiplying both sides of Eq. (9) by $ζ1mζ2n$⁠, and integrating with respect to $ζ1$ and $ζ2$⁠, yields the transport equations for the moments of $fξζ1,ζ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 $ξi′ξj′−1$ are the inverse components defined in (10). If we denote the maximum moment order by K (i.e., $0≤n+m≤K$⁠), 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∕τi=∑j𝜀ijξj′ξi′−1$ are assumed to be equal to $τϕ$⁠, then $𝜀12=ξ1′ξ2′∕τϕ$ in the standard model. The time scale $τϕ$ must then be related to the turbulence time scales to close the FP model. In this work, we either take $τϕ$ directly from the DNS data or simply use it as a fixed parameter to rescale the time t in the moment evolution equations.

The $β$-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 $δσψ,ψi$⁠,

where $wα$ and $ϕα$ 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,…,wN$⁠, N abscissae $ϕ1,ϕ2,…,ϕN$⁠, and the parameter $σ$⁠, assumed to be common for all KDFs.³⁷ The value of $σ$ 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.⁵³ 

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 $β$-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 $β$-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 $β$ KDF³⁷ to formulate the $β$-EQMOM closures used in this work. In theory, any other KDF with support [0, 1] could be used for this purpose.³⁷ 

The $β$ KDF is defined as

where $λα=ϕα∕σ$⁠, $μα=1−ϕα∕σ$⁠, and

represents the beta function.

A linear transformation allows this distribution to be defined on the arbitrary compact support $a,b$⁠.⁵⁶ The integer moment of order k of the $β$ KDF is

with M₀ = 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 $σ$⁠,

where the non-negative coefficients $χk$ depend only on $σ$⁠, and $Mk*=∑α=1Nwαϕαk$⁠. The improved version⁵³ of the iterative procedure of Yuan et al.³⁷ is used to compute the values of the quadrature weights, abscissae, and of the parameter $σ$ from the 2N + 1 transported scalar moments (1, $ϕ$⁠, $ϕ2$⁠, $…$⁠, $ϕ2N$⁠). Note that fixing N = 1 reduces $β$-EQMOM to the assumed $β$-PDF approximation.⁶ 

As illustrated in Sec. IV E, the bivariate cases considered in this work are reduced to the application of the univariate $β$-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.

Mapping closures were initially developed by Kraichnan^57,58 and Chen et al.⁵⁹ and were applied to model turbulent scalar mixing. Pope⁴⁶ 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⁹ and demonstrated the agreement between the two models. An analytical expression for the PDF is⁴⁶ 

where

As pointed out by Pope,⁴⁶ 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.⁹ 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² 

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

with $Ni≥N+1$⁠. The secondary weights $wij$ and abscissae $ψij$ are the Gauss-Jacobi quadrature corresponding to the $β$ KDF, which can be computed using the methodology described in the work of Yuan et al.³⁷ 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.

In order to investigate its capability to capture the known shape of the PDF, and to determine the required number of KDFs, $β$-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 $β$-EQMOM-FP to determine the effect of the model constant c_FP, in order to match the AMC PDF. Finally, the $β$-EQMOM-FP model is validated against the results of the DNS of binary mixing reported in the work of Eswaran and Pope,⁹ and of ternary mixing reported in the work of Juneja and Pope.⁴⁵ 

The first step is to verify that the $β$-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 $β$-EQMOM procedure. Moments of the AMC PDF are computed and provided as input to $β$-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 $β$-EQMOM closely matches the AMC PDF at later times. As noted earlier, the latter is expected because both the AMC and the $β$ PDF approach a Gaussian PDF for large times. Results with $N≥2$ show that the reconstructed PDFs agree closely with the AMC PDFs at all times. These results demonstrate that $β$-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: $⟨ξ⟩$⁠, $ξ2$⁠, $ξ3$⁠, $ξ4$⁠. Compared to the assumed $β$-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.

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.⁶⁰ Here, as in MMC, it is necessary to consider two micromixing time scales: $τϕ$ and $τϕ∕cFP$⁠. 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 $τϕ$⁠.

It is worth noting that the moments $ϕ$ and $ϕ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=∞$ 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⁶⁰ 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 $τϕ∝k∕𝜀$ where k is the turbulent kinetic energy and $𝜀$ its dissipation, if the Schmidt number is greater than one, then the diffusion time scale would yield $cFP∝k∕𝜀ν$⁠, where $ν$ is the kinematic viscosity of the fluid.

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 $β$-EQMOM-FP model against the DNS results of Eswaran and Pope.⁹ 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 [$𝜀ϕ(t)$] obtained from the DNS Eswaran and Pope,⁹ for the three cases, is provided as input to the $β$-EQMOM-FP model. The evolution of $𝜀ϕ$ as a function of the dimensionless time $τ$ 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 $ϕm$ for $m∈{1,2,3,4}$⁠. The parameter $cFP=10$ was chosen to reproduce the AMC behavior as described in Sec. IV B.

Figure 5 shows the scalar PDF reported by Eswaran and Pope⁹ compared to the results obtained using $β$-EQMOM-FP at different values of the scalar variance $Φ$⁠, for case 1. Two KDFs (N = 2) were used to obtain these results. It is apparent that the results provided by $β$-EQMOM-FP after the first few time steps differ from those of the DNS.⁹ This difference, however, can be explained by noting that the initial condition defined in the DNS simulations⁹ is not an exact pair of Dirac delta functions, due to the difficulties of generating this configuration in a discretized space. However, with $β$-EQMOM-FP, as with AMC, it is possible to have an exact initial definition of the scalar PDF. Further advancing the $β$-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 $β$-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 $β$-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.

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≥2N+1$ secondary nodes³⁷ 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.

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 $β$-EQMOM-FP model is capable of describing such cases, we consider the ternary mixing problem studied by Juneja and Pope,⁴⁵ who reported the evolution of the joint PDF $fϕ(ψ1,ψ2;t)$ of the two scalars. We will use their results for the marginal PDFs to validate the results obtained with $β$-EQMOM-FP, starting from the bivariate mixture-fraction model introduced in Sec. III B.

Juneja and Pope⁴⁵ 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 $ϕ1$ and $ϕ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 $ψ1$ and $ψ2$⁠. Additionally, the two scalars are uncorrelated so that $ϕ1ϕ2=0$⁠,⁴⁵ and hence the cross SDR is also zero when written in terms of $ϕ$⁠.

In order to solve the problem using the mixture-fraction PDF, the composition phase space $ψ$ is mapped linearly onto the mixture-fraction phase space $ζ$⁠, so that the matrix for the joint CSDR introduced in the work of Fox,² defined in terms of $ζ$⁠, 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 $ζ$ and $ψ$ is found

Note that using these transformation it is straightforward to write the moments $ϕ1mϕ2n$ in terms of the moments $ξ1mξ2n$⁠. Also note that $ξ1=ξ2=1∕3$⁠, and the mixture-fraction variances are $ξ1′2=13ϕ12+19ϕ22$ and $ξ2′2=49ϕ22$⁠, while the mixture-fraction covariance $ξ1′ξ2′=−29ϕ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 $ξ$ to those defined in terms of $ϕ$ by using the second part of (38). The transformed joint scalar dissipation matrix can be written as

where $𝜀ϕ11=ϕ12τϕ$ and $𝜀ϕ22=ϕ22τϕ$ are the SDR components for $ϕ$⁠, and the scalar mixing time $τϕ$ is found from the DNS. Note that (39) implies that $∑j𝜀ijξj′ξk′−1=1τϕδi,k$ in (9), which is just the single scalar time scale model used in the IEM model.² 

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

By expanding the right-hand side of (40), the univariate moments needed for $β$-EQMOM-FP reconstruction of the two marginal PDFs for $ϕ1$ and $ϕ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,