# Comparison of Presumed PDF Models of Turbulent Flames

Authors: Chen Huang; Andrei N. Lipatnikov
Journal: Journal of Combustion
 (2012)

---

## Abstract

Over the past years, the use of a presumed probability density function (PDF) for combustion progress variable or/and mixture fraction has been becoming more and more popular approach to average reaction rates in premixed and partially premixed turbulent flames. Commonly invoked for this purpose is a beta-function PDF or a combination of Dirac delta functions, with the parameters of the two PDFs being determined based on the values of their first and second moments computed by integrating proper balance equations. Because the choice of any of the above PDFs appears to be totally arbitrary as far as underlying physics of turbulent combustion is concerned, the use of such PDFs implies weak sensitivity of the key averaged quantities to the PDF shape. The present work is aimed at testing this implicit assumption by comparing mean heat release rates, burning velocities, and so forth, averaged by invoking the aforementioned PDFs, with all other things being equal. Results calculated in the premixed case show substantial sensitivity of the mean heat release rate to the shape of presumed combustion-progress-variable PDF, thus, putting the approach into question. To the contrary, the use of a presumed mixture-fraction PDF appears to be a sufficiently reasonable simplification for modeling the influence of fluctuations in the mixture fraction on the mean burning velocity provided that the mixture composition varies within flammability limits.


---

## Body

## 1.  Introduction

When simulating combustion of a turbulent gas mixture, characterized by large magnitude of temperature fluctuations, the key challenge consists of averaging (or filtering if large eddy simulation is concerned) quantities that depend nonlinearly on the temperature, for example, reaction rates, densityρ, and so forth. To resolve the problem, fresh reactants and equilibrium combustion products are often assumed to be separated by thin, inherently laminar, self-propagating layers (sometimes called flamelets) that are wrinkled and stretched by turbulent eddies [1–3]. Within the framework of such a paradigm, the probability of finding intermediate (between unburned and burned) states of the mixture is assumed to be much less than unity and turbulent flame propagation is addressed by solving a balance equation for the Favre-averaged combustion progress variable c~≡ρc¯/ρ¯, which is equal to zero (unity) in unburned (burned) mixture and is often associated with properly normalized temperature. Accordingly, the probability of finding combustion products is equal to the Reynolds-averaged combustion progress variable c¯ and the Reynolds-averaged density ρ¯ is evaluated as follows [4]:
(1)ρ¯=ρu(1-c¯)+ρbc¯,(2)ρbc¯=ρ¯c~,
where subscripts u and b designate unburned and burned mixture, respectively. To close the c~-equation, the mean mass rate W¯ of product creation is modeled by considering the influence of turbulence on flamelet surface area and internal structure. Alternatively, a G equation [5] is also invoked to simulate premixed turbulent combustion within the framework of the same paradigm of thin flamelets. In the latter case, turbulent flame speed St should be modeled.While significant progress was already made by invoking the above approaches to simulate premixed turbulent flames, certain important problems have not yet been resolved. For instance, first, the influence of turbulent fluctuations in equivalence ratioF on W¯ and St should be thoroughly studied, because the composition of unburned reactants is inhomogeneous in many internal combustion engines developed currently to respond environmental challenges. Second, to simulate emissions from the engines, not only a single global rate W but also rates of pollutant formation should be averaged invoking a chemical mechanism of combustion. Third, in sufficiently intense turbulence, the paradigm of thin flamelets does not work properly [1–3] and, in particular, (2) may be wrong.Among methods developed to resolve the above problems, the so-called presumed probability density function (PDF) approach has been becoming more and more popular tool over past years. The approach consists of specifying a shape of a Favre PDFρ(ξ)P(ξ)/ρ¯ for a random variable ξ associated with either c or mixture fraction f, which is commonly used instead of F in numerical simulations of partially premixed or nonpremixed flames [5], because variations in f are limited by unity. If the shape of P~(ξ) is specified, the PDF may be determined by using the first ξ~, second ρ(ξ)ξ′′2¯/ρ¯, and eventually higher moments ρ(ξ)ξ′′n¯/ρ¯, with these moments being computed by solving proper balance equations closed by invoking P~(ξ). Such an approach was already used for evaluating W¯=∫01‍W(c)P(c)dc [6–10], or ∫01‍W¯(c~,f)P(f)df and ρb¯=∫01‍ρb(f)P(f)df [11–13], or S¯t=∫01‍St(f)P(f)df and ρb¯ [14–19], or W¯=∬01W(c,f)P(c,f)dcdf and ρ¯=∬01ρ(c,f)P(c,f)dcdf [20–25]. Note that the use of the G-equation to simulate partially premixed burning [14–19] implies that (1) holds provided that ρu and ρb are substituted with ρu(f)¯ and ρb(f)¯, respectively, with the dependence of the unburned density on the mixture fraction being commonly disregarded. In the following, the density will be normalized using the unburned gas density, that is, ρu=1.It is worth also stressing that the well-known Bray-Moss-Libby (BML) PDF, that is, the sum of two Dirac delta functions associated withc=0 and 1 [4], does not allow us to average W(c), because W(c=0)=W(c=1)=0. Accordingly, the following discussion is solely restricted to PDFs presumed in the range of 0<ξ<1, while the BML PDF, which does not address the intermediate range of 0<c<1, is beyond the scope of the present paper.In the vast majority of the cited papers, two basically different shapes of PDFs were presumed. First, a beta-function PDF given by the following expression(3)P(ξ,x,t)=Γ(a+b)Γ(a)Γ(b)ξa(x,t)-1(1-ξ)b(x,t)-1
was introduced into the combustion literature by Janicka and Kollmann [26] who simulated a turbulent nonpremixed flame and considered P(f,x,t), that is, ξ=f in (3). Later, a beta-function PDF P(c,x,t), that is, (3) with ξ=c, was applied to premixed and partially premixed flames in a couple of papers [6–13, 15, 17, 20, 21]. Here, t is time, vector x is associated with spatial coordinates, and a ratio of gamma functions Γ(z)≡∫0∞‍ηz-1e-ηdη is used in order to satisfy the normalizing constraint of ∫01‍P(ξ)dξ=1.Equation (3) involves two unknown parameters a(x,t) and b(x,t), which are commonly evaluated using the following expressions:
(4)a(x,t)=ξ¯(g-1-1),b(x,t)=(1-ξ¯)(g-1-1),
where
(5)g≡ξ′2¯ξ¯(1-ξ¯)
is a segregation factor. To do so, the first ξ¯(x,t) and second ξ′2¯(x,t) moments are determined by solving proper balance equations. If Favre-averaged balance equations are integrated to compute ξ~(x,t) and ρξ′′2¯(x,t), then (i) the right hand side (RHS) of (3) is used to approximate the Favre PDF
(6)P~(ξ,x,t)≡ρ(ξ)P(ξ,x,t)ρ¯(x,t)
and (ii) ξ¯(x,t) and ξ′2¯(x,t) are substituted with ξ~(x,t) and ρξ′′2¯(x,t)/ρ¯(x,t), respectively, in (4) and (5).The shape of the beta-function PDF is very flexible and looks like either a Gaussian function or the sum of two Dirac delta functions atg≪1 and 1-g≪1, respectively, see Figure 1. For instance, if a<1 (or b<1), then the PDF given by (3) tends to infinity and resembles Dirac delta function at ξ→0 (or ξ→1). However, contrary to Dirac delta function, substitution of (3) into the integral ∫0ϵ‍P(ξ)dξ with ϵ→0 yields vanishing probability of finding ξ=0, for example, vanishing probability of finding unburned mixture in the case of premixed combustion (ξ=c). Indeed, if c¯≪1 and g differs markedly from zero and unity, then, a≪1, Γ(a+b)≈Γ(b), and the above integral is approximately equal to ϵa/[aΓ(a)]. Moreover, if a≪1, then, Γ(a)≈a-1, because Γ(z+1)=zΓ(z) and Γ(1)=1. Therefore, in the considered example, the integral ∫0ϵ‍P(c)dc≈ϵa vanishes as ϵ→0, that is, the probability of finding unburned mixture (in the case of premixed combustion) or oxidizer (in the case of nonpremixed burning) vanishes. Similarly, one can easily show that the probability of finding burned mixture (in the case of premixed combustion) or fuel (in the case of nonpremixed burning) also vanishes if the beta-function PDF is invoked. This basic limitation of the approach should be borne in mind.Figure 1
Beta-function PDFs calculated using (3) for (1)ξ¯=0.5 and g=0.05, (2)ξ¯=0.5 and g=0.15, (3)ξ¯=0.5 and g=0.40, (4)ξ¯=0.1 and g=0.15, (5)ξ¯=0.8 and g=0.25.Second, another presumed PDF, that is, a combination of Dirac delta functionsδ(ξ)(7)P(c,f,x,t)=α(x,t)δ[c-c1(x,t)]δ[f-f1(x,t)]+[1-α(x,t)]δ[c-c2(x,t)]δ[f-f2(x,t)]
was introduced by Libby and Williams [27] to simulate partially premixed flames. The second term on the RHS is multiplied with factor [1-α(x,t)] in order to satisfy the normalizing constraint of ∬01P(c,f)dcdf=1.The use of (7) substantially simplifies averaging reaction rates, for example,
(8)W¯(x,t)=α(x,t)W[ξ1(x,t)]+[1-α(x,t)]W[ξ2(x,t)]
if we consider a double-delta-function PDF
(9)P(ξ,x,t)=α(x,t)δ[ξ1(x,t)-ξ]+[1-α(x,t)]δ[ξ2(x,t)-ξ]
for a single random variable ξ. However, evaluation of the parameters α(x,t), ξ1(x,t), ξ2(x,t) is more difficult, because the first two moments ξ¯(x,t) and ξ′2¯(x,t) are not sufficient to determine the three unknown parameters and the system of three nonlinear algebraic equations
(10)ξ¯(x,t)=α(x,t)ξ1(x,t)+[1-α(x,t)]ξ2(x,t),ξ′2¯(x,t)=ξ2¯(x,t)-ξ¯2(x,t)(11)=α(x,t)ξ12(x,t)+[1-α(x,t)]ξ22(x,t)-ξ¯2(x,t),
and, for example,
(12)ξ′3¯(x,t)=ξ3¯(x,t)-3ξ¯(x,t)ξ′2¯(x,t)-ξ¯3(x,t)=α(x,t)ξ13(x,t)+[1-α(x,t)]ξ23(x,t)-3ξ¯(x,t)ξ′2¯(x,t)-ξ¯3(x,t)
should be solved for this purpose. To avoid such numerical complications, Ribert et al. [22] have proposed a simplified model that yields
(13)ξ1(x,t)=ξ¯(x,t)[1-g1/2(x,t)],ξ2(x,t)=ξ¯(x,t)+[1-ξ¯(x,t)]g1/2(x,t),α(x,t)=1-ξ¯(x,t)
in the case of a single random variable ξ. Note that the Reynolds-averaged moments ξn¯(x,t) used in the above equations may be substituted with the Favre-averaged moments ρξn¯(x,t)/ρ¯(x,t) provided that P(ξ,x,t) is simultaneously substituted with the Favre PDF defined by (6).In the following, we will call the presumed PDF given by (7) or (9) and (10)–(12) the 2δ-PDF, while the presumed PDF given by (7) or (9) and (13) will be called the R-2δ-PDF, by referring to the paper by Ribert et al. [22]. Moreover, the dependencies of the above moments, of the parameters of presumed PDFs, and so forth on time t and spatial coordinates x will not be specified for the sake of brevity.As far as premixed or partially premixed turbulent combustion is concerned, neither the beta-function shape of the presumed PDF forξ=c (or ξ=f) nor (7) has been substantiated by basic physical arguments and, to the best of the present authors knowledge, the main reasons for invoking either (3) or (7) consist of (i) numerical efficiency and (ii) simplicity of implementation. Such reasons are certainly of importance, but only if the mean reaction rate and other key mean quantities are mainly controlled by the first moments of a presumed PDF, but are weakly sensitive to its shape. If, however, two different presumed PDFs characterized by equal ξ¯ and by equal ξ′2¯ yield substantially different W¯, or S¯t, or ρ¯, then, the use of such PDFs does not seem to be a solid predictive approach until the shape of one of them is justified by basic arguments or by extensive experimental or DNS study.The present work is aimed at (i) investigating the sensitivity ofW¯,S¯t, and ρ¯ to the shape of presumed PDFs P(c) and P(f) and (ii) at identifying conditions under that the sensitivity may be considered to be sufficiently weak and the use of the studied presumed PDFs appears to be justified.


## 2.  Method of Research

The following analysis is based on comparison of the results of averaging dependencies ofw(c,f),ρ(c),sL(f),sL1/2(f), and ρb(f) obtained by invoking presumed PDFs given by either (3) or (7), with all other things being equal. Here, the combustion progress variable c=(T-Tu)/[Tb(f)-Tu] is equal to the normalized temperature, the burned temperature Tb depends on the mixture fraction defined as follows f=(F-Fmin)/(Fmax-Fmin), where Fmin and Fmax are presumed limits of fluctuations in the equivalence ratio, w(c,f)≡W(c,f)/Ω(f) is the heat release rate W(c,f) normalized using its maximum value Ω=max{W(c)} calculated by varying the combustion progress variable, sL(f)=SL(f)/max{SL(f)} is the laminar flame speed SL(f) normalized using its maximum value max{SL(f)} calculated by varying the mixture fraction, and the densities ρ and ρb(f) are normalized using the unburned gas density. The normalized laminar flame speed sL and its square root are addressed, because variations in turbulent flame speed St with mixture composition are mainly controlled by the dependence of SL on the equivalence ratio and numerous experimental data reviewed elsewhere [28] indicate that St∝SLq, with q being close to 0.5.Dependencies to be averaged are shown in Figure2. They have been computed by running the PREMIX code [32] of the Chemkin-II package [33] and using a semidetailed (111 species and 616 reactions) chemical mechanism for gasoline-like fuel mixture (iso-octane, toluene, and n-heptane in volumetric proportions of 55%, 35%, and 10%, resp.) developed by Golovitchev et al. [34, 35]. Dependencies of SL on the equivalence ratio, pressure, and temperature, yielded by this mechanism, are compared with available experimental data [29–31] in Figure 3, while results of validating the mechanism against published data on ignition delay times are reported elsewhere [35].Data for averaging invoking presumed PDFs. (a) Normalized heat release ratesw(c) versus combustion progress variable c, calculated for various gasoline-air mixtures, with the equivalence ratio being specified in legends. (b) Normalized laminar flame speed sL (solid curve) and maximum heat release rate Ω (dashed curve) versus equivalence ratio F.
(a)
(b)Validation of the chemical mechanism. Symbols show experimental data, and curves have been computed. (a) Laminar flame speed versus equivalence ratio.(1)Tu=298K, data by Metghalchi and Keck [29]. (2)Tu=350K, data by Metghalchi and Keck [29]. (3)Tu=353K, data by Zhao et al. [30]. (4)Tu=500K, data by Zhao et al. [30]. (b) Laminar flame speeds versus pressure at three different equivalence ratios specified in legends. Tu=373K. Data by Jerzembeck et al. [31].
(a)
(b)We will firstly consider premixed flames. In this case, the equivalence ratio is stationary and spatially uniform and thef-dependent part of a presumed PDF is reduced to a single Dirac delta function at a given f0, for example, f1=f2=f0 in (7). Accordingly, we will compare the normalized mean heat release rates
(14)w¯(c¯,c′2¯,f0)=∫01‍w(c,f0)P(c,c¯,c′2¯)dc,
the normalized densities
(15)ρ¯(c¯,c′2¯,f0)=∫01‍ρ(c,f0)P(c,c¯,c′2¯)dc,
and the normalized turbulent burning velocities
(16)ut=δt∫-∞∞‍w¯dζ=δt∫01‍w¯(dc¯dζ)-1dc¯
obtained for various given first c¯ and second c′2¯ moments by substituting either (i) (3)-(4), or (ii) the 2δ-PDF, or (iii) the R-2δ-PDF into the integrals on the RHSs, with ξ being associated with c in all the three cases. Here,
(17)ρ=11+τc,τ=ρb-1-1 is the heat release parameter, ζ=(z-zf)/δt is the normalized spatial distance z counted along the normal to the mean flame brush, which has a presumed thickness δt, which was not varied in the present simulations, zf is associated with c¯=0.5, and the well-known complementary-error-function profile [28]
(18)c¯=1-12erfc⁡(πζ)
is invoked to calculate the last integral on the RHS of (16).When addressing fluctuations in the mixture fraction, we will compare the mean quantities(19)q¯(f¯,f′2¯)=∫01‍q(f)P(f,f¯,f′2¯)df
obtained for various given first f¯ and second f′2¯ moments by substituting either (i) (3)-(4), or (ii) the 2δ-PDF, or (iii) the R-2δ-PDF into the integral on the RHS, with ξ being associated with f in all the three cases. Here, q is equal to either ρb(f), or to sL(f), or to sL1/2(f), or to the normalized maximum (for various c) heat release rate ω(f)=Ω(f)/max{Ω(f)}, where max{Ω(f)} is the maximum value of Ω(f) determined by varying f, while Ω(f) is the maximum value of W(f,c) calculated by varying only c.Note that averagingw(c)-curves implies that turbulence does not change the dependence of the normalized heat release rate on the normalized temperature c in a premixed flame. Moreover, averaging sL(f), sL1/2(f), and ω(f) implies that a partially premixed turbulent flame is associated with an ensemble of premixed flames characterized by various equivalence ratios.Because the use of the 2δ-PDF requires knowledge of the third moment ξ′3¯, it was evaluated invoking the beta-function PDF P(ξ,ξ¯,ξ′2¯) and was subsequently set to determine the 2δ-PDF. Thus, when comparing quantities averaged invoking either the beta-function PDF or the 2δ-PDF, both functions were characterized not only by the same first and the same second moments, but also by the same third moments. To the contrary, when the R-2δ-PDF was invoked, the third moment yielded by it could differ significantly from ξ′3¯ calculated using the beta-function PDF with the same ξ¯ and the same ξ′2¯. However, obtained numerical results indicate that the beta-function and R-2δ PDFs yield almost equal ξ′3¯ for the same ξ¯ and the same g if g≥0.5 (cf. curves 1 and 2 in Figure 4(c) and note that curves 1 and 3 are indistinguishable). When g is decreased, the difference in ξ′3¯ given by the two PDFs is increased (see Figures 4(a) and 4(b)).Dependencies ofξ′3¯ on ξ¯, calculated for the beta-function PDF (curves 1), the R-2δ-PDF (curves 2), and the 2δ-PDF (curves 3). (a) g=0.05, (b) g=0.25, (c) g=0.5.
(a)
g
=
0.05
(b)
g
=
0.25
(c)
g
=
0.5


## 3. Results: Combustion Progress Variable PDF

Figures5 and 6 clearly show that the mean heat release rate may be very sensitive to the shape of presumed PDF even if two different PDFs are characterized by the same first, the same second, and the same third moments. Such sensitivity is weakly pronounced only at low g (see Figure 5(a)), but is substantial even at g=0.25 (see Figure 5(b)) and becomes extremely strong as g→1, see Figure 5(c).Dependencies of the normalized mean heat release rate on the Reynolds-averaged combustion progress variable, calculated for the beta-function PDF (curves 1), the R-2δ-PDF (curves 2), and the 2δ-PDF (curves 3). F=1.0. (a) g=0.05, (b) g=0.25, (c) g=0.9.
(a)
g
=
0.05
(b)
g
=
0.25
(c)
g
=
0.9Figure 6
Dependencies of the normalized mean heat release rate on the Reynolds-averaged combustion progress variable, calculated for various equivalence ratios specified in legends andg=0.5 using the beta-function PDF (thin curves) and the 2δ-PDF (bold curves).This trend is caused by the highly nonlinear dependence of the heat release rate on the combustion progress variable, whereas for a weakly nonlinear function likep0+p1ξ+p2ξ2+p3ξ3, its mean value controlled by the first three moments should be the same for the two PDFs. Because the nonlinearity of w(c) is more pronounced in lean and rich mixtures as compared with near-stoichiometric ones (cf. the half-widths of the computed curves plotted in Figure 2(a)), the influence of the shape of a presumed PDF on w¯(c¯) is enhanced by an increase in |1-F|.Figure5(b) indicates that the use of the δ-PDFs may yield two peaks in the obtained dependencies of w¯(c¯). This observation is explained in Figure 7. If the 2δ-PDF is invoked and c¯ is sufficiently large, then (10)–(12) may not have a solution consistent with c2≤1 (thin dashed curve does not reach c¯=1 in Figure 7). At slightly lower c¯, the parameter c1 rapidly grows with c¯ (see thin dashed curve) and reaches a value close to cm=0.685 that is associated with the maximum w(c) for F=1, see solid curve in Figure 2(a). Accordingly, the first term on the RHS of (8) is significantly increased by c¯ as c1→cm and the mean rate w¯ calculated using this equation is also increased by c¯ yielding the second peak of the w¯(c¯)-curve. If the R-2δ-PDF is invoked (see solid curves in Figure 7), then a similar second peak is also observed (see curve 2 in Figure 5(b)), but the effect is less pronounced, because c1 is substantially less than cm and w(c1)≪max{w(c)}. Accordingly, an increase in w(c1) by c¯ on the RHS of (8) is less pronounced and counteracted by a decrease in α as c¯ tends to unity, see (13). Moreover, when the R-2δ-PDF was invoked, we obtained 0≤c1<c2≤1 in all studied cases, contrary to the 2δ-PDF.Figure 7
Dependencies of the parametersc1 (thin curves) and c2 (bold curves) on the Reynolds-averaged combustion progress variable, calculated for F=1.0 and g=0.5 using the R-2δ-PDF (curves 1) and the 2δ-PDF (curves 2).Due to the sensitivity of the mean heat release rate to the shape of a presumed PDF, the normalized burning velocityut is also sensitive to the PDF shape, with the effect being most pronounced in rich mixtures, see Figure 8. It is worth noting that the observed decrease in the normalized ut with g does not mean a decrease in the turbulent burning velocity with the segregation factor, because the mean turbulent flame brush thickness, used to evaluate ut, may also depend on g. For the goals of the present work (i.e., comparison of various presumed PDFs), eventual dependence of δt on g is of minor importance and was disregarded, but it should be taken into account in order to investigate the influence of the segregation factor on the burning rate integral.Figure 8
Dependencies of the normalized burning rate integralut on the segregation factor c′2¯/[c¯(1-c¯)], calculated for various equivalence ratios specified in legends using the beta-function PDF (thin curves) and the 2δ-PDF (bold curves).Figure8 also shows that the dependencies of ut(g) calculated invoking the 2δ-PDF (bold curves) are nonmonotonic and may have two local maxima. This behavior is explained in Figure 9(a), which indicates that an increase in the segregation factor results in (i) widening the computed w¯(c¯)-curves and (ii) decreasing their maxima. The former trend is associated with a weaker dependence of c2 on c¯ at a larger g, see bold curves in Figure 9(b) and note that the second term dominates on the RHS of (8) under typical conditions. At moderate g and small (large) c¯, the parameter c2 is lower (higher) than the value of cm=0.868 associated the peak w(c) for the considered mixture (F=1.6). Therefore, there is a c¯ such that c2(c¯)=cm, with this c¯ decreasing when the segregation factor is increased, see bold curves in Figure 9(b). Accordingly, the parameter α evaluated for this particular c¯ increases, see Figure 9(c); that is, the second term on the RHS of (8) decreases, thus, reducing max{w¯(c¯)}, which is roughly equal to [1-α(c¯=cm)]w(cm). However, at a larger  g, this trend is overwhelmed by the widening of the computed w¯(c¯)-curves and ut∝∫-∞∞‍w¯dζ is increased by g. To the contrary, if the segregation factor is sufficiently large, for example, g=0.8, see dotted-dashed curves in Figure 9, then c2>cm even for low c¯ and a further increase in g reduces not only (1-α), but also w(c2) on the RHS of (8). Accordingly, a decrease in max{w¯(c¯)} with increasing g becomes much more pronounced, overwhelms the widening of the computed w¯(c¯)-curves, and results in decreasing ut.Dependencies of (a) the normalized mean heat release rate, (b) parametersc1 (thin curves) and c2 (bold curves), and (c) parameter α on the Reynolds-averaged combustion progress variable, calculated invoking the 2δ-PDF for F=1.6 and different segregation factors specified in legends.
(a)
(b)
(c)The above results, particularly Figures5, 6, and 8, imply that the use of a presumed PDF P(c) in order to average heat release rate and evaluate turbulent burning velocity is a flawed approach unless the shape of the PDF is supported by solid physical arguments or by a wide set of experimental or DNS data. The approach could be useful only in the case of low g<0.2 (cf. thin and bold curves in Figure 8), but a typical premixed turbulent flame is characterized by a larger g, as reviewed elsewhere, see [36, Section 3.1].As far as the mean density is concerned, the studied presumed PDFs yield similar dependencies ofρ¯(c¯) at various g (cf. thin curves 1, 2, and 3 in Figure 10), because the dependence of ρ on c, given by (17) is weakly nonlinear. However, at moderate g, the obtained ρ¯(c¯)-curves 1–3 differ substantially from curves 4 calculated using the BML Equations (1) and (2). When the segregation factor tends to unity or zero, curves 1–4 become close to each other (not shown in Figure 10). The above difference is controlled by the difference between bold c~(c¯)-curves 1–3, obtained invoking the studied PDFs, and bold curves 4, computed using (2). It is worth remembering that if (17) is averaged, then (1) holds for any g and any P(c). Indeed, on the one hand, the Favre-averaged quantity ρ-1 is equal to 1+τc~ due to (17), but, on the other hand, the Favre-averaged ρ-1 is equal to ρ¯-1 by definition of a Favre-averaged quantity q~≡ρq¯/ρ¯.Dependencies of the normalized mean density (thin curves) and Favre-averaged combustion progress variablec~ (bold curves) on the Reynolds-averaged combustion progress variable, calculated invoking the beta-function PDF (curves 1), the R-2δ-PDF (curves 2), the 2δ-PDF (curves 3), and the BML equations (curves 4) for F=1.0 and (a) g=0.25, (b) g=0.5, and (c) g=0.75.
(a)
g
=
0.25
(b)
g
=
0.5
(c)
g
=
0.75As far as each of the two presumed PDFs investigated here is concerned separately, the sum of Dirac delta functions atc1 and c2 does not seem to be a proper approach to average w(c) for the following reasons.First, Figure11 shows dependencies of the parameters of the 2δ-PDF (bold curves) and R-2δ-PDF (thin curves, which are very close to the counterpart bold curves) on c¯, calculated in the flamelet regime of premixed turbulent combustion (g=0.9). In this case, c1≪1 (dashed curves) and 1-c2≪1 (dotted-dashed curves). Therefore, both w(c1) and w(c2) are small (see solid line in Figure 2(a)) and the use of the 2δ-PDFs strongly underestimates the mean heat release rate (see Figure 5(c)) and turbulent burning velocity, compare thin and bold solid curves in Figure 8.Figure 11
Dependencies of the parameters of the R-2δ (thin curves) and 2δ (bold curves) PDFs on the Reynolds-averaged combustion progress variable, calculated for g=0.9 and F=1.0.  (1)α, (2)c1, (3)c2.Second, using (9) with ξ=c, one can easily obtain
(20)c¯=αc1+(1-α)c2,c¯(1-c¯)=[αc1+(1-α)c2][1-αc1-(1-α)c2],c′2¯=α(c1-c¯)2+(1-α)(c2-c¯)2=α(1-α)2(c1-c2)2+α2(1-α)(c1-c2)2=α(1-α)(c1-c2)2.
In the limit case of g→1 and, hence, c′2¯→c¯(1-c¯), (20) results in
(21)α⟶c2(1-c2)c2(1-c2)-c1(1-c1).
Because c1(1-c1)≥0, the RHS of (21) may be within the limits of [0,1] if either (i) c1=0 and c2<1, or (ii) c1>0 and c2=1, or (iii) c1=0 and c2=1. For brevity, we consider c1≤c2 here. In cases (i) and (ii), α=1 and α=0, respectively, and the first equality in  (20) does not hold. In case (iii), w(c1)=w(c2)=0 and (8) yields w¯=0. Thus, the PDF given by (9) cannot be used to average w(c) in the flamelet regime of premixed turbulent combustion.The beta-function PDFP(c) does not seem to be well tailored for simulating the flamelet regime either. Indeed, (4) yields a≪1 and b≪1 in the case of 1-g≪1. Accordingly, even a small error in computing c¯ and/or c′2¯ by solving proper balance equations can result in substantial errors in calculating the small parameters a and b and, therefore, substantial errors in evaluating w¯ by invoking the beta-function PDF P(c). Too small parameters a and b cannot be determined accurately by subtracting one finite quantity from another finite quantity if the two quantities are calculated even with small errors, which are inevitable, not only for numerical reasons, but also and mainly due to limitations of physical models invoked to close the relevant balance equations. A possible way of resolving the problem could consist of replacing a balance equation for c′2¯ with a balance equation for 1-g [37].Finally, it is worth noting the following interesting feature of the use of a presumed beta-function PDFP(c). In order to determine the shape of such a PDF, one has to know its first and second moments. The two moments are commonly evaluated by solving proper balance equations, with the Favre-averaged equations being addressed. Accordingly, the PDF shape should be found based on known values of c~ and ρc′′2¯/ρ¯. This could easily be done if the RHS of (3) is assumed to approximate the Favre PDF P~=ρP/ρ¯. In such a case, ρcn¯/ρ¯ and w~ evaluated using P~(c) are equal to cn¯ and w¯, respectively, calculated invoking P(c) provided that both P~(c) and P(c) are approximated by the same beta function.However, the above substitution ofP(c) with P~(c) substantially affects the shape of computed w¯(c¯)-curves. To show this effect, we compared the following two ways of averaging the heat release rate. First, we specified c¯ and gR≡c′2¯/[c¯(1-c¯)] and calculated w¯(c¯)-curve by varying c¯, keeping gR constant, and using (3), (4), and (14). Second, we specified c~ and gF≡ρc′′2¯/[ρ¯c~(1-c~)] and computed w¯(c~)-curve by varying c~, keeping gF constant and equal to the above gR, and using the following equations:
(22)w¯=∫01w(c)P(c)dc=∫01w(c)ρ¯P~(c)ρdc=Γ(a+b)Γ(a)Γ(b)∫01w(c)ρ¯ρca-1(1-c)b-1dc,a=c~(gF-1-1),b=(1-c~)(gF-1-1).
Subsequently, we transformed the latter w¯(c~)-curve to w¯(c¯)-curve using
(23)c¯(c~,gF)=Γ(a+b)Γ(a)Γ(b)∫01‍cρ¯ρca-1(1-c)b-1dc.
Results are shown in Figure 12 in thin and bold curves, respectively. While the w¯(c¯)-curves computed invoking P(c,c¯,c′2¯) look symmetrical with respect to c¯=0.5 and reach peak values around c¯=0.5, in line with numerous experimental studies [38–43] of flame surface density (FSD); bold curves calculated using the Favre PDF are strongly nonsymmetrical with the peaks being significantly shifted to the trailing edge of the flame brush. The shift (i) is mainly associated with the fact that the Reynolds-averaged c¯ is significantly larger than the Favre-averaged c~ in the largest part of a mean turbulent flame brush and (ii) holds for other magnitudes of the segregation factor, see Figure 13.Figure 12
Dependencies of the normalized mean heat release rate on the Reynolds-averaged combustion progress variable, calculated forg=0.9 and various equivalence ratios, specified in legends, by invoking the beta-function PDF. Thin and bold curves are associated with P(c,c¯,c′2¯) and P~(c,c~,ρc′′2¯/ρ¯), respectively.Figure 13
Dependencies of the Reynolds-averaged (curves 1) and Favre-averaged (curves 2) combustion progress variables, associated with the maximum mean heat release rate, on the segregation factorg. F=1.0. Thin and bold curves are associated with P(c,c¯,c′2¯) and P~(c,c~,ρc′′2¯/ρ¯), respectively.These numerical data and the aforementioned experimental results that indicate almost symmetrical dependencies of the FSD onc¯ favor the use of P(c,c¯,c′2¯) at least at larger g. However, such a method is more complicated, because it requires evaluation of c¯ and c′2¯ based on computed c~ and ρc′′2¯/ρ¯.


## 4.  Results: Mixture Fraction PDF

Figure14 shows that two very different presumed mixture-fraction PDFs (beta function and the sum of two Dirac delta functions) yield sufficiently close dependencies of the mean normalized density ρ¯b of burned gas, the mean normalized maximum heat release rate ω¯, and sL1/2¯ on the mean equivalence ratio F¯ even if gf≡f′2¯/[f¯(1-f¯)] is as large as 0.4. It is worth remembering that, contrary to fluctuations in the combustion progress variable, which may be characterized by g→1, fluctuations in mixture fraction are characterized by significantly lower gf in a typical stratified turbulent flame.Dependencies of the normalized mean densityρ¯b of burned gas (a), the normalized maximum heat release rate ω¯ (b), and the square root sL1/2¯ of the normalized laminar flame speed (c) on the mean equivalence ratio F¯, calculated either invoking the beta-function (curves 1), the R-2δ (curves 2), and 2δ (curves 3) PDFs, or neglecting fluctuations in the mixture composition (curves 4). Fmin=0.5 and Fmax=1.9,gf=0.4.
(a)
(b)
(c)The difference between results obtained by invoking the considered presumed PDFs is decreased whengf is decreased, for example, see Figure 15 and note that a similar trend is observed for other simulated quantities q={ρb,ω,sL,sL1/2}. Results computed using the 2δ-PDF are not shown here and in the following, because they are close to results calculated using the R-2δ-PDF.Figure 15
Dependencies of the normalized mean laminar flame speeds¯L on the mean equivalence ratio F¯, calculated invoking the beta-function (thin curves) and R-2δ (bold curves) PDFs for different gf specified in legends. Fmin=0.5 and Fmax=1.9.It is worth also stressing that the values ofρb(F¯), ω(F¯), and sL1/2(F¯), obtained by neglecting fluctuations in mixture composition, differ substantially from counterpart quantities averaged by invoking the presumed PDFs (cf. curves 4 and 1-3, respectively, in Figure 14 or unity with the maxima of curves in Figure 15). Therefore, the above results (sufficiently weak dependence of q(f)¯ on the shape of presumed PDF P(f), but significant difference in q(f)¯ and q(f¯)) support invoking a presumed PDF P(f) for simulating stratified turbulent combustion if fluctuations in mixture fraction are weak or moderate, for example, g≤0.4.However, such a conclusion is correct only if the equivalence ratio always fluctuates within flammability limitsFl<F<Fr. If the amplitude of fluctuations in F is increased so that inflammable mixture composition may be observed with a finite probability pq=∫FminFl‍P(F)dF+∫FrFmax‍P(F)dF, then computed q(f)¯ is more sensitive to the shape of presumed PDF. For instance, Figure 16 indicates that a decrease in Fmin from 0.5 to 0.4 and an increase in Fmax from 1.9 to 2.2 make sL1/2¯ evaluated in statistically lean or rich mixtures sensitive to the shape of presumed PDF (cf. thin and bold solid curves). In the present simulations, the lean and rich flammability limits are arbitrarily set to be equal to 0.5 and 2.0, respectively. Such a simplification appears to be reasonable in order to gain insight into the influence of PDF shape on averaged quantities. However, it is worth remembering that local flammability limits in a partially premixed turbulent flame are affected by local gradients of mixture fraction, stretch rate, transient effects, and so forth and a model capable for predicting these local phenomena has not yet been elaborated.Figure 16
Dependencies of the square rootsL1/2¯ of the normalized mean laminar flame speed on the mean equivalence ratio F¯, calculated invoking the beta-function (thin curves) and the R-2δ (bold curves) PDFs. g=0.2. (1)Fmin=0.5 and Fmax=1.9, (2)Fmin=0.4 and Fmax=2.2, (3)Fmin=0.3 and Fmax=2.5, (4)Fmin=0.1 and Fmax=4.0.If the sum of two Dirac delta functions is invoked, thensL1/2¯ drops at certain F¯1<1 and F¯2>1 (see bold solid line), because both f1 and f2 are increased by f¯ (this trend is shown for ξ=c in Figure 11) and f1 (or f2) reaches the lean (or rich) flammability limit when F¯=F¯1 (or F¯=F¯2). Accordingly, if F¯<F¯1 (or F¯>F¯2), then f1 (or f2) is outside the lean (or rich) flammability limit, hence, sL(f1)=0 [or sL(f2)=0], and, consequently, sL1/2¯ is small. Further extension of the range of fluctuations in the mixture composition makes the discussed effects more pronounced (cf. thin and bold dashed curves in Figure 16), and the dependencies of sL1/2¯ on F¯, calculated by invoking the beta-function and R-2δ PDFs, differ substantially from one another if, for example, Fmin=0.1 and Fmax=4.0 (cf. thin and bold dotted curves).Figure17 shows that the probability pf=1-pq of finding flammable mixture is very sensitive to the shape of presumed PDF P(f). If the beta-function PDF is invoked, pf depends smoothly on F¯ (thin curves), whereas the probability calculated using the sum of two Dirac delta functions drops sharply (bold curves) as f1 or/and f2 moves beyond the flammable range fl<f<fr.Figure 17
Dependencies of the probability of finding flammable mixture composition on the mean equivalence ratioF¯, calculated invoking the beta-function (thin curves) and the R-2δ (bold curves) PDFs. g=0.2. (1)Fmin=0.4 and Fmax=2.2, (2)Fmin=0.3 and Fmax=2,5, (3)Fmin=0.1 and Fmax=4.0.


## 5.  Conclusions

The use of a presumed combustion-progress-variable PDFP(c) for averaging heat release rate in a premixed turbulent flame characterized by sufficiently large c′2¯/[c¯(1-c¯)]>0.4 does not seem to be justified, because the mean rate is sensitive to the shape of P(c) even if first c¯, second c′2¯, and third c′3¯ moments are kept constant when changing the PDF shape.Becauseq(f)¯ differs substantially from q(f¯) and q(f)¯ depends sufficiently weak on the shape of a presumed mixture-fraction PDF P(f,f¯,f′2¯), with f¯ and f′2¯ being kept constant and f′2¯/[f¯(1-f¯)]<0.4, the use of P(f) for simulating the influence of fluctuations in mixture composition on mean quantities appears to be a reasonable approach to modeling stratified turbulent combustion, but only if the mixture composition always fluctuates within flammability limits. If the probability of finding inflammable mixture is substantial, the mean burning rate is sensitive to the PDF shape.Although results computed by invoking the beta-function PDF given by (3) look often more reasonable than results obtained by using the Dirac-delta-function PDF given by (7), this trend puts the latter PDF into question, but does not validate the former PDF.Ifg>0.5, then the R-2δ-PDF model by Ribert et al. [22] yields the third moment approximately equal to the third moment calculated invoking the beta-function PDF provided that the two PDFs are characterized by the same first and the same second moments.



---

*Source: 564621-2012-07-31.xml*