GIEM Model for Micromixing in Non-Premixed Turbulent

Ind. Eng. Chem. Res. 1998, 37, 2131-2141

2131

The BMC/GIEM Model for Micromixing in Non-Premixed Turbulent Reacting Flows Kuochen Tsai†,‡ and Rodney O. Fox*,§,| Beam Technologies, Inc., Ithaca, New York 14850, and Department of Chemical Engineering, Kansas State University, Manhattan, Kansas 66506-5102

A novel beta mapping closure (BMC)/generalized IEM (GIEM) model is proposed to close the micromixing term in the joint composition probability density function (PDF) transport equation for nonpremixed, multiscalar turbulent reacting flows. Because the beta PDF gives an excellent description of binary inert scalar mixing, it is used to approximate the mixture fraction PDF. The proposed model thus effectively extends presumed PDF methods to handle finite-rate chemistry without introducing ad hoc assumptions about the form of the joint composition PDF. The BMC/GIEM model can be solved efficiently using a particle-based Monte-Carlo code. The model is first compared with the amplitude mapping closure (AMC)/GIEM model and homogeneous reaction-diffusion data for scalar mixing with a two-step reaction. Its application to inhomogeneous flows is demonstrated by simulating a reacting scalar mixing layer with a one-step reaction using a Lagrangian composition PDF (LCPDF) code. The linear-mean-squareestimation (LMSE) model is also evaluated for comparison. The results from the BMC/GIEM model are in excellent agreement with experimental data, whereas the LMSE model gives errors in higher-order statistics. 1. Introduction The numerical simulation of turbulent reacting flows is a major subfield of computational fluid dynamics with many important applications in the chemical industry [Fox, 1996]. Of all the available approaches for treating turbulent reacting flows, probability density function (PDF) methods are particularly attractive because of their ability to treat the chemical reaction terms exactly [O’Brien, 1980; Pope, 1985; Fox, 1996]. Nevertheless, closure of the molecular-mixing (or micromixing) term in the PDF balance equation remains a challenge [Fox, 1996], and none of the currently available micromixing models satisfies all of the desired criteria. In particular, two criteria that must be satisfied to successfully model a single inert scalar are (i) the scalar PDF must relax to a Gaussian form regardless of its initial shape, and (ii) a bounded scalar must remain bounded. The first micromixing model to satisfy both criteria was the amplitude mapping closure (AMC) [Chen et al., 1989]. Other micromixing models, such as the coalescencedispersion (CD) model [Curl, 1963] or the linear-meansquare-estimation (LMSE) model [Dopazo, 1975; O’Brien, 1980] (also referred to as the IEM model [Villermaux, 1991]), do not yield the limiting Gaussian shape. When applied to flows with multiple scalars (i.e., typical reactors in the chemical industry), at least five additional criteria should also be satisfied [Pope, 1995]; these are, (iii) linearity between inert scalars (i.e., elemental balances), (iv) localness of micromixing in * Corresponding author. Telephone: (785) 532-5584. Fax: (785) 532-7372. Email: [email protected]. † Beam Technologies, Inc. ‡ Currently with The Dow Chemical Company, Freeport, TX. § Department of Chemical Engineering. | Written while on sabbatical leave at the E Ä cole Nationale Supe´rieure des Industries Chimiques, Institut National Polytechnique de Lorraine, Nancy, France.

composition space, (v) (possibly time-dependent) boundedness for joint scalar statistics, (vi) appropriate dependence on scalar length scales, and (vii) appropriate dependence on the Reynolds, Schmidt, and Damko¨hler numbers. Criterion (iii) arises from the fact that the transport equations for inert scalars are linear. Thus, any linear combination of transported scalars with equal diffusivity should follow the same transport equation. Criterion (iv) is related to the resolution of the computational grid. Because micromixing occurs at the smaller scales in the flow [Fox, 1996], the scalar values at a particular spatial point should change smoothly and continuously (i.e., not undergo random jumps). Simple “random-pairing” models, such as the CD model, do not distinguish between large and small values of scalar composition; thus, the localness criterion is usually violated. (An example where criterion (iv) is particularly critical is a turbulent diffusion flame far from extinction where the chemical reactions are nearly in local equilibrium. For this example, the simple CD model mixes cold fuel with cold oxidizer, producing an unphysical nonreacting mixture across the flame sheet [Norris and Pope, 1991; Norris and Pope, 1995].) Criterion (v) is associated with possible correlations between scalars. For example, the total mass concentration of any element is always constant during chemical reaction. Stochastic micromixing models, such as the binomial Langevin model [Valinõ and Dopazo, 1991], usually change the scalar composition value of each scalar at random; as a consequence, the physical boundaries in the joint composition space are not preserved. The only micromixing model that satisfies criteria (i)-(v) is again the AMC [Valinõ, 1995]. Unfortunately, its implementation for multiscalar problems is very difficult numerically, and its application to inhomogeneous flows still remains a challenge. Criteria (vi) and (vii) can only be satisfied with higher-order

S0888-5885(97)00589-7 CCC: $15.00 © 1998 American Chemical Society Published on Web 03/06/1998

2132 Ind. Eng. Chem. Res., Vol. 37, No. 6, 1998

micromixing models that include a description of the scalar length scales as well as the scalars [Fox, 1995 and 1997]. Most such models (e.g., the Fokker-Planck (FP) closure [Fox, 1994]) are more complicated than the LSME or CD models and thus difficult to apply to inhomogeneous flows. In summary, given the difficulty for any single micromixing model to satisfy all seven criteria, most practical applications of PDF methods to date have opted for the LSME model because of its well-known computational efficiency [Fox, 1996]. Nevertheless, a potentially attractive alternative micromixing model (the generalized IEM model) has been recently developed and successfully applied to describe nonpremixed, multiscalar homogeneous reacting flows [Tsai and Fox, 1995]. In that study, the micromixing term of the mixture fraction [a required input to the generalized IEM (GIEM) model] was found from an explicit analytical solution of the AMC [Gao, 1991]. The combined AMC/GIEM model preserves all the aforementioned criteria except (vi) and (vii). Unfortunately, because the inverse mapping function that maps the surrogate scalar composition field onto the Gaussian reference field is difficult to implement numerically, the AMC part of the model cannot be easily implemented for real (inhomogeneous) flows. Because application of the GIEM model is independent of the AMC, we propose in this work an alternative mapping function based on the beta PDF for the mixture fraction [Fox, 1996] that is applicable to inhomogeneous flows. The novel beta mapping closure (BMC)/GIEM model not only yields an asymptotic joint Gaussian PDF for inert scalars, but, because each scalar has the same relaxation rate as the mixture fraction, it also fulfills criterion (iv). Like the AMC, the BMC is known to yield a good approximation for nonpremixed binary mixing; indeed, it has been shown to yield a better description of direct numerical simulation (DNS) data for the initial stages of mixing [Girimaji, 1991]. More importantly, because both the beta PDF and its cumulative distribution function (CDF) can be accurately approximated using efficient numerical algorithms [Press et al., 1990], the difficult inverse mapping problem associated with the AMC solution is completely avoided. The BMC/ GIEM model thus yields a computationally tractable description of micromixing for binary, inhomogeneous multiscalar turbulent reacting flows. However, because the BMC/GIEM model is only intended to close the micromixing term in the PDF balance equation at the one-point level, the difficulties associated with modeling the scalar length scales are not addressed here [i.e., criterion (vi)]. Such information can be obtained from other well-established closures (e.g., equilibrium turbulence models), the spectral relaxation model [Fox, 1995 and 1997], or, as is done in this work, directly from DNS/experimental data. To demonstrate the applicability of the BMC/GIEM model to inhomogeneous flows, here we apply it to simulate a reacting scalar mixing layer experiment [Bilger et al., 1991] in conjunction with a Lagrangian composition PDF code (LCPDF); [Tsai and Fox, 1994, Fox, 1996]. (Note that the experiment has been repeated [Li et al., 1992] to correct errors in the original data. The same data set has been employed in an earlier modeling study [Wang and Tarbell, 1993].) The LCPDF code is a modified version of the Lagrangian velocity-composition PDF code PDF2DS [Pope, 1992]

and, for single-phase incompressible flows, has been shown to produce nearly identical results [Tsai and Fox, 1994]. The reacting scalar mixing-layer experiment contains several important features suitable for model validation, including data for the means and variances of the mixture fraction and reacting scalars, and for the skewness and kurtosis of the mixture fraction. Moreover, the scalar mixing layer setup allows us to avoid uncertainties in modeling low Reynolds number flows caused by wall boundary layers. Both the BMC/GIEM and LMSE models have been integrated into the LCPDF code. By comparing the simulation results with DNS/ experimental data, the advantages of the combined BMC/GIEM model over the LMSE model can be verified. The remainder of this paper is organized as follows: the theoretical background of the BMC/GIEM model is given in Section 2; the LCPDF code and the Monte-Carlo algorithm used to implement the new multiscalar micromixing model is reviewed in Section 3; the simulation results and comparisons between the BMC/GIEM and the LMSE models, and DNS/experimental data are discussed in Section 4; and conclusions are drawn in Section 5. 2. The BMC/GIEM Model Our earlier study [Tsai and Fox, 1995] suggested that the application of the GIEM model to inhomogeneous flows could be achieved by using an inverse function closure based on the AMC. However, computational considerations favor a more practical algorithm based on the BMC. Both approaches are reviewed below. A. The Amplitude Mapping Closure. The AMC/ GIEM model starts from the following equality for the transport of the mixture fraction PDF, f(ξ; t), in constant density homogeneous turbulence:

∂2 ∂f ∂ ) [〈D∇2ζ|ξ〉f] ) - 2[〈D(∇ζ)2|ξ〉f] ∂t ∂ξ ∂ξ

(1)

where ζ is the mixture fraction, ξ is the corresponding composition variable in phase space, D∇2ζ is the molecular diffusion term, D(∇ζ)2 is the scalar dissipation term, and the symbol 〈‚|ξ〉 denotes the expected value conditioned on the event {ζ(x,t) ) ξ}. By extending the LMSE model, the micromixing term 〈D∇2ζ|ξ〉 can be modeled for binary mixing as follows:

〈D∇2ζ|ξ〉 )

1 ∂ [〈D∇ζ)2|ξ〉f] ) C(ξ;t) [β(t) - ξ] (2) f(ξ;t) ∂ξ

where C(ξ;t) is a function of ξ (instead of a constant as in the LMSE model) and β(t) is the zero-flux point in phase space [Tsai and Fox, 1995]. This formulation has two advantages over other micromixing models: first, the closure function has the same form when extended to multiple scalars, thus preserving criterion (ii); and second, with the proper choice of C(ξ;t), the asymptotic form of the mixture fraction PDF will be Gaussian. Thus, for nonpremixed binary mixing, C(ξ;t) can be expressed as follows:

C(ξ;t) )

∂ 1 [〈D∇ζ)2|ξ〉f] [β(t) - ξ] f(ξ;t) ∂ξ

with the zero-flux point β(t) defined by eq 4:

(3)

Ind. Eng. Chem. Res., Vol. 37, No. 6, 1998 2133

∂ [〈D∇ζ)2|ξ〉f]|ξ)β ) 0 ∂ξ

(4)

Using the AMC, 〈D∇ζ)2|ξ〉 and f(ξ;t) have the following forms [O’Brien and Jiang, 1991]:

〈D∇ζ)2|ξ〉 ) F(t) exp{-2[erf-1(2ξ - 1)]2}

(5)

and

f(ξ;t) ) fG(ξ0)

(

-1

(6)

where fG(ξ0) is the Gaussian density function of the reference field ξ0, with zero mean and unit variance:

X(ξ0;t) ) ξ ) erf

(

ξ0e-τ(t) - C0

)

x2 - 2e-2τ(t)

(7)

where τ(t) is the normalized time, C0 is determined by the area under the left-hand peak of the initial PDF f(ξ;0):

a0 )

1 ∫ x2π

C0 -x2/2 e -∞

dx

(8)

and F(t) is a function of time that can be obtained by the following relation:

∫01exp{-2[erf-1(2ξ - 1)]2} dξ ) 〈D(∇ζ)2〉 ) 〈φ〉

F(t)

(9) As noted in the Introduction, the scalar dissipation 〈φ〉 can be determined from other one-point models, such as the k- or Reynolds-stress model, or from an explicit transport equation [Fox, 1995 and 1997]. In practice, it is not necessary to evaluate C(ξ;t) through eq 3. Because f(ξ;t) is known at all times (eq 6), integrating both sides of the first equality in eq 1 yields:

∂ ∂t

∫0ξf(x;t) dx ) -〈D∇2ζ|ξ〉 f(ξ;t)

(10)

Then, applying eq 2 to the right-hand side of 10, C(ξ;t) can be derived as follows:

∂ 1 C(ξ;t) ) [ξ - β(t)] f(ξ;t) ∂t

∫0 f(x;t) dx ξ

(11)

Using the equality in eq 12:

f(ξ;t) dξ ) fG(ξ0) dξ0

(12)

Equation 11 can be rewritten as follows:

C(ξ;t) )

∂ 1 [ξ - β(t)] f(ξ;t) ∂t

∫-∞ξ (ξ)fG(x) dx 0

(13)

B. The Beta Mapping Closure. Although theoretically sound, eq 13 requires the evaluation of ξ0(ξ) ) erf-1(ξ) and because ξ0 has a range of (-∞, ∞), the numerical algorithm cannot be easily implemented. Fortunately, because the beta PDF gives an excellent approximation for the mixture fraction PDF, instead of using the AMC solution, f(ξ;t) can be approximated by a beta PDF. With this assumption, C(ξ;t) can be found

1 ∂ I (a(t),b(t)) [ξ - β(t)] fβ(ξ;t) ∂t ξ

C(ξ;t) )

(14)

where fβ(ξ;t) is the beta PDF and Iξ(a(t),b(t)) is the beta CDF. The incomplete beta function is defined by

Iξ(a,b) )

)

∂X(ξ0;t) ∂ξ0

from

∫0ξsa-1(1 - s)b-1 ds

1 B(a,b)

(a, b > 0)

(15)

In eq 15, the variable t has been dropped to simplify the notation, and B(a,b) is given by eq 16:

B(a,b) )

∫01sa-1(1 - s)b-1 ds

(16)

where a and b are related to the mean (〈ζ〉) and variance (〈ζ′2〉) through eqs 17 and 18:

a ) 〈ζ〉

(

〈ζ〉(1 - 〈ζ〉)

(

b ) (1 - 〈ζ〉)

〈ζ′2〉

〈ζ〉(1 - 〈ζ〉) 〈ζ′2〉

)

-1

(17)

)

-1

(18)

Note that a and b are implicit functions of t that can be related to the scalar dissipation through eq 19: 2 ∂a ∂a ∂〈ζ′ 〉 ∂a ∂〈ζ〉 ∂a (19) ) + ) -2〈φ〉 ∂t ∂〈ζ′2〉 ∂t ∂t ∂〈ζ〉 ∂〈ζ′2〉

by using the fact that micromixing leaves the mean mixture fraction unchanged (i.e., ∂〈ζ〉/∂t ) 0). C. The GIEM Model. For multiple scalars, the PDF balance equation for molecular diffusion in homogeneous turbulence [O’Brien, 1980] can be written as follows:

∂f ∂ [〈D∇2φk|ψ〉f] )∂t ∂ψk

(20)

where φk are the reacting scalars, ψ is the vector of ψk, and f(ψ;t) is the joint PDF of ψk. (The repeated index k implies summation.) The micromixing term on the right-hand side can be closed by the LMSE model as follows:

〈D∇2φk|ψ〉 )

Cφ (〈φ 〉 - ψk) τφ k

(21)

where Cφ is the mechanical-to-scalar timescale ratio and τφ is the micromixing timescale [Fox, 1996]. As noted in the Introduction, although this model is easy to apply and fulfills the linearity criterion, it fails to relax the inert scalar PDF to an asymptotic Gaussian form [criterion (i)]. The problem arises from the fact that τφ is, in general, not a constant, which suggests the following alternative representation of Cφ/τφ:

C(ξ;t) )

〈 |〉 Cφ ξ τφ

(22)

where C(ξ;t) (i.e., eq 13 or 14) can now be regarded as the local micromixing rate conditioned on the mixture fraction. Assuming (as in eq 21) that all reacting scalars have the same micromixing rate, the GIEM model for multiscalar Fickian diffusion can be formulated as follows:


〈D∇2φk|ψ〉 ) 〈D∇2φk|ξ, ψk〉 ) C(ξ;t)[βk(t) - ψk]

(23)

Using the fact that ∂〈φk〉/∂t ) -〈D∇2φk〉 ) 0 during micromixing, βk(t) must satisfy eq 24:

βk(t) )

〈C(ζ;t)φk〉 〈C(ζ;t)〉

3. Monte-Carlo Implementation Assuming constant density, the transport equation for the joint composition PDF f(ψ;x,t) in an inhomogeneous turbulent reacting flow can be written as follows [Pope, 1985]:

∂ ∂ ∂f + [S (ψ)f] ) + 〈uj〉 ∂t ∂xj ∂ψk k ∂ ∂ ∂ (〈u′|ψ〉f) [〈D∇2φk|ψ〉f] - [〈D∇2ζ|ψ〉f] (25) ∂xj j ∂ψk ∂ξ where uj is the turbulent velocity field and Sk(ψ) is the chemical source term for φk. (The repeated indices j and k imply summation.) Note that the terms on the lefthand side are closed, which includes transport by the mean velocity and the chemical source term [Fox, 1996]. The three terms on the right-hand side represent the probability fluxes due to velocity fluctuations, and micromixing of the reacting scalars and the mixture fraction, respectively. These terms are closed by the gradient-diffusion assumption and BMC/GIEM model as follows:

∂f ∂xj

〈D∇2φk|ψ〉 ) C(ξ;x,t) [βk(x,t) - ψk] 〈D∇2ζ|ψ〉 )

∂ 1 I (a(x,t),b(x,t)) fβ(ξ;x,t) ∂t ξ

1/2 dx* ) (〈u〉 + ∇ΓT/F)x* dt + (2ΓT/F)x* dW(t) (29)

(24)

where the expected value is with respect to the joint composition PDF. The combination of eqs 14 and 23 (i.e., the BMC/ GIEM model) offers several advantages: • The linear relationship between the reacting scalars is guaranteed. • The asymptotic PDF of the mixture fraction is Gaussian. • Because C(ξ;t) is deterministic, the localness requirement for multiscalar micromixing models is also guaranteed, which ensures the correct behavior of the joint statistics between the mixture fraction and reacting scalars. • Because both fβ(ξ;t) and Iξ(a,b) can be efficiently evaluated numerically [Press et al., 1990], the conditional micromixing rate C(ξ;t) can be computed for inhomogeneous flows at a reasonable cost.

〈u′j|ψ〉f ) -ΓT(x,t)

equation [Gardiner, 1990]. The details can be found elsewhere [Pope, 1985] and are not repeated here. The resultant SDE for the position of a Lagrangian particle x*(t) reads as follows:

where the subscript x* denotes the spatial location where the coefficients are evaluated. In the MonteCarlo simulation, dW(t) is approximated by η(∆t)1/2, where η is a standardized joint normal random vector. Likewise, the scalar values in a particle obey eq 30:

dψk ) C(ξ;x*,t) [βk(x*,t) - ψk] + Sk(ψ) dt

Note that the micromixing rate C(ξ;x*,t) and βk(x*,t) are functions of the particle location x*(t), thereby coupling eq 30 to eq 29. Similarly, the SDEs are coupled to the turbulent velocity field through the dependence on 〈u〉 and ΓT. Using time splitting [Pope, 1985], the micromixing and reaction terms in eq 30 are computed separately. Thus, the change in the scalars due to micromixing can be updated using eq 14 (with the x*-dependence dropped to simplify the notation):

ψk,t+∆t ) ψk,t +

(28)

where ΓT(x,t) is the turbulent diffusivity (discussed later), and the micromixing rate C(ξ;x,t) is found by applying eq 14 at every point x in the flow. Because the dimension of f(ψ;x,t) increases linearly with the number of reacting scalars, the direct computation of eq 25 is intractable. Instead, a set of Lagrangian stochastic differential equations (SDEs) is formulated [Fox, 1996]. The relationship between the SDEs and eq 25 is derived using the Fokker-Plank

∆ξ(t)

[βk(t) - ψk,t] (31)

[β(t) - ξt]f(ξ;t)

with

∆ξ(t) ) Iξ(a(t+∆t),b(t+∆t)) - Iξ(a(t),b(t))

(32)

Furthermore, Iξ(a,b) can be approximated by a continuous fraction [Press et al., 1990] with an error of 1.0 × 10-7 everywhere on the interval ξ ∈ (0,1). Values of a and b can be found from eq 19 with the mixture fraction variances 〈ζ′2〉t and 〈ζ′2〉t+∆t found before and after the micromixing step with the known scalar dissipation rate 〈φ〉. However, initially, ξ only has values of 0 and 1, which produces a Heaviside function for Iξ(a,b), and the continuous fraction approximation cannot be used. To overcome this difficulty, the concept of the AMC can be employed to construct a uniform reference field ξ0 (instead of a Gaussian reference field used in the AMC) through I-1 ξ (a,b) on the interval [0,1]. Thus, the updated mixture fraction can be found from eq 33:

(26) (27)

(30)

(a(t+∆t),b(t+∆t)) ξt+∆t ) Iξ-1 0

(33)

using the following probability equality:

ξ0 )

∫0ξ ds ) ∫0ξfβ(s;t) ds ) Iξ(a(t),b(t)) 0

(34)

In other words, the difficulty of evaluating the Heaviside function at t ) 0 can be avoided by using eqs 33 and 34 at t ) ∆t. Because Iξ(a,b) can be efficiently approximated, the root search for ξ in eq 34, which is needed for evaluating ∆ξ (eq 32), can be implemented using standard techniques (e.g., Brent’s method [Press et al., 1990] is adopted in this study). Note that 〈ζ〉 and 〈ζ′2〉 change between micromixing and reaction steps due to turbulent advection (i.e., eq 29).


Figure 1. Normalized conditional dissipation predicted by the AMC and BMC at different values of 〈ζ′2〉 with 〈ζ〉 ) 0.9.

Figure 2. Scalar PDFs predicted by the BMC at different values of 〈ζ′2〉 with 〈ζ〉 ) 0.9.

The numerical procedure just described for the micromixing step is repeated at every point (i.e., computational cell) in the flow and can be summarized as follows: 1. Compute the mixture fraction variance 〈ζ′2〉t and 〈ζ′2〉t+∆t with the known scalar dissipation 〈φ〉 (i.e., found from DNS/experimental data for the mixture fraction). 2. Compute a(t), b(t), a(t + ∆t) and b(t + ∆t) using eq 19. 3. Compute β(t) by finding the root of ∆ξ(t) ) 0. 4. Update ξ with eqs 33 and 34. 5. Compute βk(t) using eq 24. 6. Update ψk with eq 31. It should be noted that the aforementioned implementation is in complete agreement with the BMC outlined in Section 2. Both models apply the principle that the CDFs in the reference and the physical fields are equal. The difficulty associated with a “singular” PDF arises only because of the finite samples used in the numerical simulation. The chemical source term Sk(ψ) is a function of the local composition vector ψ and can be analytically solved for the one-step bilinear reaction studied in the reacting scalar mixing layer experiment. More details concerning the numerical implementation of the Monte-Carlo simulations can be found elsewhere [Pope, 1985; Tsai and Fox, 1994; Fox, 1996], including recent advances in computing the chemical source term for complex reactions [Pope, 1997].

solution, which has the same functional form at all time, as shown by eq 9. Similar results are found using the BMC. Although the curves are slightly shifted to the right for the BMC, both models predict the same concave form. Although a constant shape for this function has been argued to be incorrect [O’Brien and Sahay, 1992] (because the end values of the scalar concentration will disappear gradually, leading to a delta function at the mean), this argument should indeed be taken carefully because once the probability of a scalar value is zero, its dissipation rate can no longer be defined. Thus, the shape of the conditional dissipation function should not be an important issue as long as the relaxation process of the scalar PDF can be correctly approximated. Shown in Figure 2 are the corresponding PDFs of the mixture fraction, which follow the beta PDF given by eq 6 with known mean and variance. Note that the PDFs are still far from Gaussian for the range of variances shown. With the knowledge of the mixture fraction PDF, the evolution of the conditional diffusion term, 〈D∇2ζ|ξ〉, can be derived using eq 14 and is shown in Figure 3. In contrast to the LMSE model, which gives a linear function with a constant slope, the BMC generates a nonlinear function with zero values at 0 and 1 and with the zero flux point β(t) moving from 0.5 (t ) 0) towards 〈ζ〉 ) 0.9 as t increases. (The zero-flux points are the abscissas of the crossing points between the 〈D∇2ζ|ξ〉 curves and the horizontal line 〈D∇2ζ|ξ〉 ) 0.) This behavior mimics binary mixing results found by the DNS [Leonard and Hill, 1991] and thus greatly enhances the accuracy of the model predictions. Homogeneous Reacting Flow. For completeness, the BMC/GIEM model is first compared with the AMC/ GIEM model for the following two-step reaction in a homogeneous reactor:

4. Results and Discussion Before applying the BMC/GIEM model to multiple reacting scalars, we first compare the prediction by BMC of the conditional scalar dissipation to the AMC solution (eq 5). Using the equality between the second and the third terms in eq 1, the BMC yields the following expression:

〈D(∇ζ)2|ξ〉 )

1 ∂ I (a,b) ds ∫0ξf (s;t) ∂t s

(35)

k1

A + B 98 R k2

β

B + R 98 S

Shown in Figure 1 is the predicted scalar conditional dissipation at different values of 〈ζ′2〉 (i.e., in place of t). The initial PDF is a double-delta function with 〈ζ〉 ) 0.9 and 〈ζ′2〉 ) 0.09. The solid line shows the AMC

with k1 ) 100 and k2 ) 10. (The reaction rate units have no effect on the model comparison because the micromixing rate is computed directly from the mixture fraction variance in the homogeneous reactor [Tsai and

(36)


Figure 3. Normalized conditional Laplacians predicted by the BMC at different values of 〈ζ′2〉 with 〈ζ〉 ) 0.9.

Fox, 1995]. Nevertheless, because the scalar dissipation rate falls quickly as a function of time for this flow [Fox, 1994], both reaction rates remain well above the mixing rate as time increases.) The initial concentrations of A and B are both set at 1.0 with a volume ratio of 1.0. This reaction is known to be sensitive to micromixing; thus, the predicted mean concentrations are also sensitive to the choice of the micromixing model. A one-dimensional, reaction-diffusion equation is solved with randomly distributed slabs of A and B to mimic the behavior of the two-step reaction in a homogeneous reactor without turbulent stretching [Muzzio and Ottino, 1989; Fox, 1994; Tsai and Fox, 1995]. Similar concepts can also be found in the linear-eddy model [Kerstein, 1988]. The initial distribution of the slab thickness is assumed to be Gaussian. This system has two advantages over DNS: first, a higher resolution can be achieved without limiting the reaction rates (i.e., Damko¨hler numbers); and second, the initial distributions of A and B can be made completely segregated. Moreover, as we have demonstrated elsewhere [Tsai et al., 1996], turbulent stretching only affects the scalar variance decay rate and has no direct influence on the joint scalar PDF. Figures 4-6 show the evolution of mean concentrations for all species as found from the reaction-diffusion equation, and the AMC/GIEM, BMC/GIEM and LMSE models, respectively. It can be seen that the results of the AMC/GIEM and BMC/GIEM models are virtually identical and that both give excellent agreement with the results from the reaction-diffusion equation. However, this result is not surprising because both models predict similar local micromixing rates and mixture fraction PDFs. In contrast, the predictions of the LMSE model deviate significantly from the results from the reaction-diffusion equation. (Deviations are even larger when higher volume ratios are employed.) It should be noted that such homogeneous systems present a highly rigorous test for micromixing models. Indeed, as demonstrated by the reacting scalar mixing layer, turbulent advection a plays major role in determining the shape of the joint scalar PDF with inhomogeneous flows, thereby reducing the direct effect of micromixing. Reacting Scalar Mixing Layer. The performance of BMC/GIEM model for inhomogeneous turbulent

Figure 4. Evolution of scalar means 〈φA〉 (+), 〈φB〉 (*), 〈φR〉 (]), and 〈φS〉 (4) predicted by the one-dimensional, reaction-diffusion (RD) equation and the AMC/GIEM model.

Figure 5. Same as in Figure 4, but with data predicted by the one-dimensional, reaction-diffusion (RD) equation and the BMC/ GIEM model.

reacting flows has been validated using data for a reacting scalar mixing layer [Bilger et al., 1991; Li et al., 1992] with the following one-step, irreversible reaction: k1

NO + O3 98 NO2 + O2 + 200 kJ/mol

(37)

where the rate constant k1 is equal to 0.37 ppm-1 s-1 at 20 °C [Chameides et al., 1977]. The second experiment [Li et al., 1992] repeated the original experiment with some modifications to the experiment setup and a corrected computer code, which reduced the velocity results by a factor of x2. Both experiments studied the reacting scalar mixing layer in grid turbulence using nonpremixed, dilute nitric oxide (NO) and ozone (O3) as reactants. The mixture fraction is defined as [Bilger et al., 1991] follows:

ζ)

φA - φB + φB2 φA1 + φB2

(38)

where φA is the NO concentration, φB is the O3 concen-


variances near x ) 0. However, because this developing stage is relatively short and the scalar mixing layer is asymptotically self similar, this assumption only slightly affects the simulation results. The values of u′ and were fitted at x/M ) 12 and 21 for Re ) 8273 (11 700 [Bilger et al., 1991] before the computer code correction [Li et al., 1992]). Here Re is defined as 〈u〉M/ν with ν the kinematic viscosity and M ) 320 mm. The turbulent diffusion coefficient, ΓT, which is invoked in the gradient-diffusion assumption, is related to the eddy viscosity by [Jones and Launder, 1972]:

ΓT )

Figure 6. Same as in Figure 4, but with data predicted by the one-dimensional, reaction-diffusion (RD) equation and the LMSE model.

Figure 7. Schematic diagram for the reactive scalar mixing layer.

tration, φA1 is the NO concentration in stream 1, and φB2 is the O3 concentration in stream 2, as shown in Figure 7. The velocity field can be well approximated by assuming a constant mean velocity 〈u〉 ()0.389 m/s), and the turbulence intensity u′ decays as [Bilger et al., 1991; Bilger, 1993] follows:

( ) u′ 〈u〉

2

()

) c1

x M

-m

(39)

where m ) 1.3. Thus, the timescale τu for the turbulent kinetic energy k decays as follows:

〈u〉 1 ) ) c2m τu k x

(40)

where M is the grid spacing. The “standard” model for the scalar dissipation is

〈u〉 2 〈ζ′ 〉 x

〈φ〉 ) aTm

(41)

where aT is the mechanical-to-scalar timescale ratio (usually assumed to be ∼1.0 [Pope, 1985]). Because both u′ and have been measured for the reacting scalar mixing-layer experiment [Bilger et al., 1991], the values of c1 and c2 have been estimated from the experimental data. Their values at x ) 0 were taken as the cellaverage within the first computational cell to avoid singularity. This procedure may result in a different spreading rate of the scalar mixing layer and the scalar

0.09k2 σφ

(42)

where σφ is the turbulent Schmidt number and has a value of 0.35 for the reacting scalar mixing layer experiment [Bilger et al., 1991]. The remaining undetermined coefficient is the mechanical-to-scalar timescale ratio aT (i.e., Cφ). Because modeling the exact details of the grid turbulence is not possible, this value must be adjusted for the model to predict reasonable values of the mixture fraction variances. The effect of varying Cφ is discussed later. Shown in Figure 7 is the schematic diagram for the experiment and the PDF simulation. The simulation domain has a dimension of 7 × 4 m2 and a grid size of 101 × 41. The grid size was doubled and no significant differences were found. The width of the simulation domain was determined by giving enough distance for ζ to reach 0 and 1 smoothly and ensures that the spread of the scalar mixing layer is not affected by use of a finite grid. Thus, the upper and lower boundary conditions were set to have zero gradients in the Y direction for the scalars, which gives reflection boundary conditions for the notional particles [Tsai and Fox, 1994]. Simulations with a width of 3 m have also been attempted, and no changes in the scalar mixing layer thickness were found. The number of notional particles in each computational cell has been varied from about 100 to 400, and it was found that the mean and higherorder statistics are insensitive to the change in particle numbers. Thus, a particle number of 100 (approximately) in each cell was adopted for all the PDF simulations. Note that the higher-order statistics, including the variance, skewness, and kurtosis, were computed by accumulating particle values over 50 time steps after the solution had reached a steady state. Thus, each cell had a sample size of 5000, which was found to be sufficient. The profile in Figure 8 is of the mean mixture fraction across the scalar mixing layer predicted by the BMC/ GIEM and LMSE models, and the experimental data [Li et al., 1992] at x/M ) 21, with the initial concentrations of NO and O2 at 4.08 and 3.85 ppm, respectively. The scalar mixing layer thickness δ, which is defined as the distance from 〈ζ〉 ) 0.1 to 〈ζ〉 ) 0.9, has a value of 1.4 m measured from the experimental data and a value of 1.0 m found from the simulation. This difference resulted from the choice of the turbulent Schmidt number, σφ. Because the value of σφ adopted in the simulation is from the experimental data measured at x/M ) 16 and 21, its value at earlier stages of the scalar mixing layer is unknown; thus, the expansion rate is underestimated. Fortunately, because the scaled concentration reaches a self-similar state at early stages,


Figure 8. Profile of the mean mixture fraction 〈ζ〉 across the scalar mixing layer with σφ ) 0.08 and Cφ ) 0.8, where ] denotes the averaged experimental data [Li et al., 1992] and the vertical bars indicate their range.

this value does not affect the profiles of the higher-order moments across the scalar mixing layer, which has been verified by comparing the profiles at different downstream locations from the simulations (not shown). In order for the PDF simulation to produce the same δ as measured by the experiments, σφ must be set at 0.08 (i.e., much lower than the “standard” value). It can be seen that both micromixing models produce similar results for the mean values and approximate the experimental data closely. Predicting the profile of the mixture fraction root mean square (rms), defined as

ζrms ) 〈(ζ - 〈ζ〉)2〉1/2

(43)

is not as straightforward as predicting the mean value. As found in earlier studies [Fox, 1995; Tsai and Fox, 1994], the mechanical-to-scalar timescale ratio, Cφ, which determines the scalar variance decay rate relative to the energy dissipation rate, is highly flow dependent. Because the PDF simulation cannot simulate the early stages of the scalar mixing layer correctly for the reasons mentioned earlier, the value of Cφ is subject to modification. After some trial and error, it was found that both the BMC/GIEM and LMSE models predict the correct profile for 〈ζ′2〉 with Cφ) 0.8, as shown in Figure 9. It should be noted that because both micromixing models produce the same evolution equation for 〈ζ′2〉, this result also demonstrates the accuracy of the BMC/ GIEM model. Although the PDF of the mixture fraction can be described very closely by a beta PDF, the scalar dissipation rate cannot be reproduced exactly due to the finite sample size in each computational cell. Without any correction, the BMC for C(ξ;t) (eq 14) produces an average error of 3% for the scalar dissipation rate. If a high-accuracy solution is desired, the following errorcorrecting algorithm can be applied to the micromixing step to minimize the statistical error:

Figure 9. Profile of the mixture fraction rms ζrms across the scalar mixing layer with σφ and Cφ the same as in Figure 8.

1. Calculate the updated and correct mixture frac′2 〉 using eq 31 and ∂〈ζ′2〉/ tion variances 〈ζt+∆t ′2 〉 and 〈ζexact ∂t ) -2〈φ〉, respectively. ′2 〉 2. Multiply 〈ζexact ′2 〉 by the ratio [1 - (〈ζt+∆t 2 2 〈ζexact ′ 〉)/〈ζexact ′ 〉] when evaluating a(t + ∆t) and b(t + ∆t) in eq 31 for the next time step. 3. Repeat steps 1 and 2 for every time step. Because the solution will reach a steady state, the aforementioned procedure is able to reduce the statistical error to near 0.5% with a sample size of 100. However, because the error is already small without any correction, this procedure does not have a significant effect on the overall results. The profiles of 〈φA〉 and 〈φB〉 are shown in Figure 10. Both the LMSE and BMC/GIEM models give excellent agreement with the experimental data and both profiles are insensitive to the value of Cφ. This result indicates that for simple (i.e., one-step reactions) isothermal turbulent reacting flows, the LMSE model is able to predict satisfactory results for the mean values of reacting species. The profiles of the normalized rms values of φA and φB, defined by eqs 44 and 45, are shown in Figure 11:

γ′A )

〈(φA - 〈φA〉)2〉1/2 φA1

(44)

γ′B )

〈(φB - 〈φB〉)2〉1/2 φB2

(45)

Again, both micromixing models give reasonable approximations to the experimental data, which is not surprising because the profiles of ζrms closely match the experimental data and because turbulent advection has a strong effect on the scalar PDFs in the reacting scalar mixing layer. However, the difference between these two micromixing models can be seen from the profiles of higher-order statistics, which are more sensitive to the tails of the


Figure 10. Profiles of 〈φA〉 (4) and 〈φB〉 (3) across the reactive scalar mixing layer predicted by the LMSE and BMC/GIEM models with σφ and Cφ the same as in Figure 8. Symbols show the experimental data [Bilger et al., 1991] without the range bars.

Figure 11. Profiles of the normalized rms values γ′A (4) and γ′B (4) across the reactive scalar mixing layer predicted by the LMSE and BMC/GIEM models with σφ and Cφ the same as in Figure 8.

scalar PDF. The skewness and kurtosis of the mixture fraction, defined by

S)

〈(ζ - 〈ζ〉)3〉 ζ3rms

(46)

and

K)

〈(ζ - 〈ζ〉)4〉 ζ4rms

(47)

Figure 12. Profile of mixture fraction skewness S across the scalar mixing layer predicted by the LMSE and BMC/GIEM models with σφ and Cφ the same as in Figure 8.

Figure 13. Profile of mixture fraction kurtosis K across the scalar mixing layer predicted by the LMSE and BMC/GIEM models with σφ and Cφ the same as in Figure 8.

were extensively studied with a wide range of experiments including different initial concentrations and different flow velocities [Bilger et al., 1991; Li et al., 1992]. Except the uncharacteristic kinks found in one experiment [Bilger et al., 1991], which were later corrected [Li et al., 1992], the experiments produced consistent profiles for both quantities. Shown in Figure 12 are the profiles of the mixture fraction skewness. Because the LMSE model preserves the shape of the scalar PDF, it predicts higher values of skewness at y/δ < 0.5 and lower values at y/δ > 0.5 compared with the experimental data. Similar behavior is also found in


Figure 14. PDFs of ζ, φA, φB, and φR predicted by the LMSE and BMC/GIEM models at y/δ ) 0.4 with σφ and Cφ the same as in Figure 8.

Cφ have been fixed. Thus, the predictions of skewness and kurtosis can only be affected by the choice of micromixing model. This situation should be of importance for other more sensitive systems, such as the multistep reaction discussed earlier, or nonpremixed reactions with high-volume ratios [Fox, 1996]. The pair-wise scatter plots of the joint statistics between the mixture fraction and reacting scalars predicted by the LMSE and BMC/GIEM models at y/δ ) 0 are shown in Figure 15. It can be seen that the tails predicted by the BMC/GIEM model are longer than those predicted by the LMSE model. This phenomenon also reflects the wider marginal scalar PDFs predicted by the BMC/GIEM model, as shown in Figure 14. The scatter plots also show an important difference between the GIEM and the conditional moment closure (CMC) [Bilger, 1993]. Our earlier study [Tsai and Fox, 1995] demonstrated that the AMC/GIEM model produces results similar to the CMC model in homogeneous flows because both models predict a one-dimensional support for the pair-wise joint PDFs. However, for inhomogeneous flows, the BMC/GIEM model yields considerable scatter caused by turbulent advection as shown in Figure 15. In contrast, because the CMC models the reacting scalar compositions as functions of the mixture fraction, such scattering cannot be predicted, even when turbulent advection plays a major role. 5. Conclusions

Figure 15. Pair-wise scatter plots of the mixture fraction and reactive scalars predicted by the LMSE (top) and BMC/GIEM (bottom) models, respectively.

the profiles of the mixture fraction kurtosis predicted by the LMSE model, as shown in Figure 13. In contrast, the BMC/GIEM model consistently predicts satisfactory values for these two quantities. As mentioned earlier, this difference stems from the fact that the LMSE model does not relax the PDF shape. This phenomenon is best illustrated by Figure 14 which shows the marginal PDFs of ζ, φA, φB, and φR at y/δ ) 0.4. It can be clearly seen that the BMC/GIEM model predicts wider scalar PDFs than the LMSE model. Note that the PDF simulation has no adjustable parameters once the values of σφ and

A novel micromixing model for nonpremixed, multiscalar chemical reactions in turbulent flows has been proposed and validated. Because the BMC is able to relax the mixture fraction PDF to a Gaussian form for nonpremixed initial conditions using the easily computable incomplete beta function, it has been combined with the GIEM model to extend the latter to inhomogeneous flows. Unlike earlier micromixing models that are either unable to relax the scalar PDF to a Gaussian form (e.g., the LMSE and CD models) or cannot be easily applied to inhomogeneous flows (e.g., the AMC), the BMC/GIEM model provides an attractive alternative micromixing model for nonpremixed, finite rate, turbulent reacting flows. Other criteria for modeling the micromixing term in the PDF formulation for multiscalar reactions include localness, boundedness, and linearity. Although the binomial Langevin model [Valinõ and Dopazo, 1991] can be applied to multiscalar reactions in inhomogeneous flows, it violates all of these criteria because each scalar composition follows its own stochastic path. The Euclidean Minimum Spanning Trees (EMST) model [Pope, 1995, Masri et al., 1996] provides a solution by assuming that the particle compositions evolve by interactions with neighboring particles, thereby guaranteeing the localness and boundedness criteria. Unfortunately, the asymptotic PDF of the EMST model can be nonGaussian, and the linearity criterion is violated. The BMC/GIEM model is the first micromixing model that preserves all of the advantages of the AMC, but can still be applied to inhomogeneous turbulent flows with nonpremixed (binary), multiscalar, finite rate chemical reactions. Because many chemical reactors involve nonpremixed initial conditions [Fox, 1996], the BMC/GIEM model should thus provide a valuable alternative to presumed PDF methods that require chemical equilibrium assumptions. However, more experimental data with mixing-sensitive reactions are


needed to validate the BMC/GIEM model for other applications (e.g., combustion). For example, recent natural-gas diffusion flame experiments [Nooren et al., 1997] are good candidates for future validation because PDF simulations have shown that the LMSE model performs poorly relative to the AMC and CD models. Acknowledgment The authors thank Prof. S. B. Pope for allowing us to access his FORTRAN code PDF2DS. This work was supported by the National Science Foundation under grant CTS-9158124 (PYI Award), and The Dow Chemical Company. The second author also acknowledges the Institut National Polytechnique de Lorraine, France, for sabbatical leave support during the final phase of the project. Finally, we thank Prof. P. Givi for bringing the second reacting mixing layer experiment to our attention. Literature Cited Bilger, R. W. Conditional Moment Closure for Turbulent Reacting Flow. Phys. Fluids A 1993, 5, 436. Bilger, R. W.; Saetran, L. R.; Krishnamoorthy, L. V. Reaction in a Scalar Mixing Layer. J. Fluid Mech. 1991, 233, 211. Chameides, W. L.; Stedman, D. H. Tropospheric Ozone Coupling Transport and Photochemistry. J. Geophys. Res. 1977, 82, 1787. Chen, H.; Chen, S.; Kraichnan, R. H. Probability Distribution of a Stochastically Advected Scalar Field. Phys. Rev. Lett. 1989, 63, 2657. Curl, R. L. Dispersed Phase Mixing I. Theory and Effects in Simple Reactors. AIChE J. 1963, 9, 175. Dopazo, C. Probability Density Function Approach for a Turbulent Axisymmetric Heated Jet Centerline Evolution. Phys. Fluids 1975, 18, 397. Fox, R. O. Improved Fokker-Planck Model for the Joint Scalar, Scalar Gradient PDF. Phys. Fluids 1994, 6, 334. Fox, R. O. The Spectral Relaxation Model of the Scalar Dissipation Rate in Homogeneous Turbulence. Phys. Fluids 1995, 7, 1082. Fox, R. O. Computational Methods for Turbulent Reacting Flows in the Chemical Process Industry. Rev. Inst. Franç ais du Pe´ trole 1996, 51, 215. Fox, R. O. The Lagrangian Spectral Relaxation Model of the Scalar Dissipation in Homogeneous Turbulence. Phys. Fluids 1997, 9, 2364. Gao, F. An Analytical Solution for the Scalar Probability Density Function in Homogeneous Turbulence. Phys. Fluids A 1991, 3, 511. Gardiner, C. W. Handbook of Stochastic Methods, 2nd ed.; Springer: New York, 1990. Girimaji, S. S. Assumed β-PDF Model for Turbulent Mixing: Validation and Extension to Multiple Scalar Mixing. Combust. Sci. Technol. 1991, 78, 177. Jones, W. P.; Launder, B. E. The Prediction of Laminarization with a Two-Equation Model of Turbulence. Intl. J. Heat Mass Transfer 1972, 15, 301. Kerstein, A. R. A Linear Eddy Model of Turbulent Scalar Transport and Mixing. Combust. Sci. Technol. 1988, 60, 391.

Leonard, A. D.; Hill, J. C. Scalar Dissipation and Mixing in Turbulent Reacting Flows. Phys. Fluids A 1991, 3, 1286. Li, J. D.; Brown, R. J.; Bilger, R. W. Experimental Study of a Scalar Mixing Layer Using Reactive and Passive Scalars, 11th Australian Fluid Mechanics Conference, Tasmania, 1992. Masri, A. R.; Subramaniam, S.; Pope, S. B. A Mixing Model to Improve the PDF Simulation of Turbulent Diffusion Flames. In Proc. 24th Symp. (Intl.) on Combustion (Naples); The Combustion Institute: Pittsburgh, 1996. Muzzio, F. J.; Ottino, J. M. Evolution of a Lamellar System with Diffusion and Reaction: A Scaling Approach. Phys. Rev. Lett. 1989, 63, 47. Nooren, P. A.; Wouters, H. A.; Peeters, T. W. J.; Roekaerts, D.; Maas, U.; Schmidt, D. Monte-Carlo PDF Modelling of a Turbulent Natural-Gas Diffusion Flame. Combust. Theory Modelling 1997, 1, 79. Norris, A. T.; Pope, S. B. Turbulent Mixing Model Based on Ordered Pairing. Combust. Flame 1991, 83, 27. Norris, A. T.; Pope, S. B. Modeling of Extinction in Turbulent Diffusion Flames by the Velocity-Dissipation-Composition PDF Method. Combust. Flame 1995, 100, 211. O’Brien, E. E. In Turbulent Reactive Flows; Libby, P. A., Williams, F. A., Eds.; Springer: Berlin, 1980. O’Brien, E. E.; Jiang, T.-L. The Conditional Dissipation Rate of an Initially Binary Scalar in Homogeneous Turbulence. Phys. Fluids A 1991, 3, 3121. O’Brien. E. E.; Sahay, A. Asymptotic Behavior of the Amplitude Mapping Closure. Phys. Fluids A 1992, 4, 1773. Pope, S. B. PDF Methods for Turbulent Reactive Flows. Prog. Energy Combust. Sci. 1985, 11, 119. Pope, S. B. PDF2DS. Unpublished FORTRAN code, 1992. Pope, S. B. The PDF Approach to Reactive Turbulent Mixing, AIChE Annual Meeting, Miami, November 1995; Paper 122A. Pope, S. B. Computationally Efficient Implementation of Combustion Chemistry Using In Situ Adaptive Tabulation. Combust. Theory Modelling 1997, 1, 41. Press, W. H.; Flannery, B. P.; Teukolsky, S. A.; Vetterling, W. T. Numerical Recipes in C; Cambridge: New York, 1990. Tsai, K.; Fox, R. O. PDF Simulation of a Turbulent Series-Parallel Reaction in an Axisymmetric Reactor. Chem. Eng. Sci. 1994, 49, 5141. Tsai, K.; Fox, R. O. The Generalized IEM Model for Reactive Scalar Mixing. Phys. Fluids 1995, 7, 2820. Tsai, K.; Chakrabarti, M.; Fox, R. O.; Hill, J. C. Assessment of PDF Micromixing Models Using DNS Data for a Two-Step Reaction, APS Div. Fluid Dynamics Annual Meeting, Syracuse, November 1996. Valinõ, L. Multiscalar Mapping Closure for Mixing in Homogeneous Turbulence. Phys. Fluids 1995, 7, 144. Valinõ, L.; Dopazo, C. A Binomial Langevin Model for Turbulent Mixing. Phys. Fluids A 1991, 3, 3034. Villermaux, J. Mixing Effects on Complex Chemical Reactions in a Stirred Reactor. Rev. Chem. Eng. 1991, 7, 51. Wang, K.; Tarbell, J. M. Closure Models for Turbulent Reacting Flows with a Nonhomogeneous Concentration Field. Chem. Eng. Sci. 1993, 48, 3907.

Received for review August 25, 1997 Revised manuscript received November 6, 1997 Accepted November 6, 1997 IE970589J

GIEM Model for Micromixing in Non-Premixed Turbulent

Recommend Documents