Comparison of Presumed PDF Models of Turbulent Flames

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.

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 𝑐̃≡𝜌𝑐/𝜌, 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 𝑐 and the Reynolds-averaged density 𝜌 is evaluated as follows [4]: 𝜌=𝜌_𝑢𝜌1−𝑐+𝜌_𝑏𝑐,(1)_𝑏𝑐=𝜌𝑐̃,(2) where subscripts 𝑢 and 𝑏 designate unburned and burned mixture, respectively. To close the 𝑐̃-equation, the mean mass rate 𝑊 of product creation is modeled by considering the influence of turbulence on flamelet surface area and internal structure. Alternatively, a 𝐺 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 𝑆_𝑡 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 ratio 𝐹 on 𝑊 and 𝑆_𝑡 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 𝑊 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 𝜌(𝜉)𝑃(𝜉)/𝜌 for a random variable 𝜉 associated with either 𝑐 or mixture fraction 𝑓, which is commonly used instead of 𝐹 in numerical simulations of partially premixed or nonpremixed flames [5], because variations in 𝑓 are limited by unity. If the shape of 𝑃(𝜉) is specified, the PDF may be determined by using the first 𝜉̃, second 𝜌(𝜉)𝜉^2/𝜌, and eventually higher moments 𝜌(𝜉)𝜉^𝑛/𝜌, with these moments being computed by solving proper balance equations closed by invoking 𝑃(𝜉). Such an approach was already used for evaluating 𝑊=∫¹⁰𝑊(𝑐)𝑃(𝑐)𝑑𝑐 [6–10], or ∫¹⁰𝑊(𝑐̃,𝑓)𝑃(𝑓)𝑑𝑓 and 𝜌_𝑏=∫¹⁰𝜌_𝑏(𝑓)𝑃(𝑓)𝑑𝑓 [11–13], or 𝑆_𝑡=∫¹⁰𝑆_𝑡(𝑓)𝑃(𝑓)𝑑𝑓 and 𝜌_𝑏 [14–19], or 𝑊=∬¹⁰𝑊(𝑐,𝑓)𝑃(𝑐,𝑓)𝑑𝑐𝑑𝑓 and 𝜌=∬⁰¹𝜌(𝑐,𝑓)𝑃(𝑐,𝑓)𝑑𝑐𝑑𝑓 [20–25]. Note that the use of the 𝐺-equation to simulate partially premixed burning [14–19] implies that (1) holds provided that 𝜌_𝑢 and 𝜌_𝑏 are substituted with 𝜌_𝑢(𝑓) and 𝜌_𝑏(𝑓), 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, 𝜌_𝑢=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 with 𝑐=0 and 1 [4], does not allow us to average 𝑊(𝑐), because 𝑊(𝑐=0)=𝑊(𝑐=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<𝑐<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 𝑃(𝜉,𝐱,𝑡)=Γ(𝑎+𝑏)Γ(𝑎)Γ(𝑏)𝜉^{𝑎(𝐱,𝑡)−1}(1−𝜉)^{𝑏(𝐱,𝑡)−1}(3) was introduced into the combustion literature by Janicka and Kollmann [26] who simulated a turbulent nonpremixed flame and considered 𝑃(𝑓,𝐱,𝑡), that is, 𝜉=𝑓 in (3). Later, a beta-function PDF 𝑃(𝑐,𝐱,𝑡), that is, (3) with 𝜉=𝑐, was applied to premixed and partially premixed flames in a couple of papers [6–13, 15, 17, 20, 21]. Here, 𝑡 is time, vector 𝐱 is associated with spatial coordinates, and a ratio of gamma functions Γ(𝑧)≡∫^∞0𝜂^𝑧−1𝑒^−𝜂𝑑𝜂 is used in order to satisfy the normalizing constraint of ∫¹⁰𝑃(𝜉)𝑑𝜉=1.

Equation (3) involves two unknown parameters 𝑎(𝐱,𝑡) and 𝑏(𝐱,𝑡), which are commonly evaluated using the following expressions: 𝑎(𝐱,𝑡)=𝜉𝑔⁻¹−1,𝑏(𝐱,𝑡)=1−𝜉𝑔⁻¹−1,(4) where 𝑔≡𝜉^2𝜉1−𝜉(5) is a segregation factor. To do so, the first 𝜉(𝐱,𝑡) and second 𝜉^2(𝐱,𝑡) moments are determined by solving proper balance equations. If Favre-averaged balance equations are integrated to compute 𝜉̃(𝐱,𝑡) and 𝜌𝜉^2(𝐱,𝑡), then (i) the right hand side (RHS) of (3) is used to approximate the Favre PDF 𝑃(𝜉,𝐱,𝑡)≡𝜌(𝜉)𝑃(𝜉,𝐱,𝑡)𝜌(𝐱,𝑡)(6) and (ii) 𝜉(𝐱,𝑡) and 𝜉^2(𝐱,𝑡) are substituted with 𝜉̃(𝐱,𝑡) and 𝜌𝜉^2(𝐱,𝑡)/𝜌(𝐱,𝑡), 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 at 𝑔≪1 and 1−𝑔≪1, respectively, see Figure 1. For instance, if 𝑎<1 (or 𝑏<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𝑃(𝜉)𝑑𝜉 with 𝜖→0 yields vanishing probability of finding 𝜉=0, for example, vanishing probability of finding unburned mixture in the case of premixed combustion (𝜉=𝑐). Indeed, if 𝑐≪1 and 𝑔 differs markedly from zero and unity, then, 𝑎≪1, Γ(𝑎+𝑏)≈Γ(𝑏), and the above integral is approximately equal to 𝜖^𝑎/[𝑎Γ(𝑎)]. Moreover, if 𝑎≪1, then, Γ(𝑎)≈𝑎⁻¹, because Γ(𝑧+1)=𝑧Γ(𝑧) and Γ(1)=1. Therefore, in the considered example, the integral ∫^𝜖0𝑃(𝑐)𝑑𝑐≈𝜖^𝑎 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 𝑔=0.05, (2)𝜉=0.5 and 𝑔=0.15, (3)𝜉=0.5 and 𝑔=0.40, (4)𝜉=0.1 and 𝑔=0.15, (5)𝜉=0.8 and 𝑔=0.25.

Second, another presumed PDF, that is, a combination of Dirac delta functions 𝛿(𝜉)𝑃(𝑐,𝑓,𝐱,𝑡)=+𝛼[(1𝐱,−𝑡)𝛼𝛿(𝐱𝑐,−𝑡)]𝑐𝛿₁(𝐱𝑐,−𝑡)𝑐₂𝛿𝑓−𝑓₁(𝐱,𝑡)(𝐱,𝑡)𝛿𝑓−𝑓₂(𝐱,𝑡)(7) was introduced by Libby and Williams [27] to simulate partially premixed flames. The second term on the RHS is multiplied with factor [1−𝛼(𝐱,𝑡)] in order to satisfy the normalizing constraint of ∬¹⁰𝑃(𝑐,𝑓)𝑑𝑐𝑑𝑓=1.

The use of (7) substantially simplifies averaging reaction rates, for example, 𝑊(𝐱,𝑡)=𝛼(𝐱,𝑡)𝑊𝜉₁(𝐱,𝑡)+[1−𝛼(𝐱,𝑡)]𝑊𝜉₂(𝐱,𝑡)(8) if we consider a double-delta-function PDF 𝑃(𝜉,𝐱,𝑡)=+𝛼[(1𝐱,−𝑡)𝛼𝛿(𝐱𝜉,₁𝑡()𝐱],𝛿𝑡)𝜉−₂𝜉(𝐱,𝑡)−𝜉(9) for a single random variable 𝜉. However, evaluation of the parameters 𝛼(𝐱,𝑡), 𝜉₁(𝐱,𝑡), 𝜉₂(𝐱,𝑡) is more difficult, because the first two moments 𝜉(𝐱,𝑡) and 𝜉^2(𝐱,𝑡) are not sufficient to determine the three unknown parameters and the system of three nonlinear algebraic equations 𝜉(𝐱,𝑡)=𝛼(𝐱,𝑡)𝜉₁(𝐱,𝑡)+[1−𝛼(𝐱,𝑡)]𝜉₂(𝐱,𝑡),𝜉^2(𝐱,𝑡)==𝛼𝜉²((𝐱𝐱,,𝑡𝑡))𝜉−₂₁𝜉²(𝐱,𝑡)(10)(𝐱,𝑡)+[1−𝛼(𝐱,𝑡)]𝜉²²(𝐱,𝑡)−𝜉²(𝐱,𝑡),(11) and, for example, 𝜉^3(𝐱,𝑡)==𝛼𝜉³((𝐱𝐱,,𝑡𝑡))𝜉−₃₁3𝜉(𝐱,𝑡)𝜉′²(𝐱,𝑡)−𝜉³(𝐱,𝑡)−3(𝐱,𝑡)+[1−𝛼(𝐱,𝑡)]𝜉³²(𝐱,𝑡)𝜉(𝐱,𝑡)𝜉′²(𝐱,𝑡)−𝜉³(𝐱,𝑡)(12) should be solved for this purpose. To avoid such numerical complications, Ribert et al. [22] have proposed a simplified model that yields 𝜉𝜉₁(𝐱,𝑡)=𝜉(𝐱,𝑡)1−𝑔^1/2(𝐱,𝑡)₂,(𝐱,𝑡)=𝜉𝛼(𝐱(𝐱,𝑡,)𝑡)+=11−−𝜉(𝐱,𝑡)𝑔^1/2(𝐱,𝑡),𝜉(𝐱,𝑡)(13) in the case of a single random variable 𝜉. Note that the Reynolds-averaged moments 𝜉^𝑛(𝐱,𝑡) used in the above equations may be substituted with the Favre-averaged moments 𝜌𝜉^𝑛(𝐱,𝑡)/𝜌(𝐱,𝑡) provided that 𝑃(𝜉,𝐱,𝑡) 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 𝑡 and spatial coordinates 𝐱 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 𝜉=𝑐 (or 𝜉=𝑓) 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 𝑊, or 𝑆_𝑡, 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 of 𝑊,𝑆_𝑡, and 𝜌 to the shape of presumed PDFs 𝑃(𝑐) and 𝑃(𝑓) 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 of 𝑤(𝑐,𝑓),𝜌(𝑐),𝑠_𝐿(𝑓),𝑠^𝐿1/2(𝑓), and 𝜌_𝑏(𝑓) obtained by invoking presumed PDFs given by either (3) or (7), with all other things being equal. Here, the combustion progress variable 𝑐=(𝑇−𝑇_𝑢)/[𝑇_𝑏(𝑓)−𝑇_𝑢] is equal to the normalized temperature, the burned temperature 𝑇_𝑏 depends on the mixture fraction defined as follows 𝑓=(𝐹−𝐹_min)/(𝐹_max−𝐹_min), where 𝐹_min and 𝐹_max are presumed limits of fluctuations in the equivalence ratio, 𝑤(𝑐,𝑓)≡𝑊(𝑐,𝑓)/Ω(𝑓) is the heat release rate 𝑊(𝑐,𝑓) normalized using its maximum value Ω=max{𝑊(𝑐)} calculated by varying the combustion progress variable, 𝑠_𝐿(𝑓)=𝑆_𝐿(𝑓)/max{𝑆_𝐿(𝑓)} is the laminar flame speed 𝑆_𝐿(𝑓) normalized using its maximum value max{𝑆_𝐿(𝑓)} calculated by varying the mixture fraction, and the densities 𝜌 and 𝜌_𝑏(𝑓) are normalized using the unburned gas density. The normalized laminar flame speed 𝑠_𝐿 and its square root are addressed, because variations in turbulent flame speed 𝑆_𝑡 with mixture composition are mainly controlled by the dependence of 𝑆_𝐿 on the equivalence ratio and numerous experimental data reviewed elsewhere [28] indicate that 𝑆_𝑡∝𝑆^𝑞𝐿, with 𝑞 being close to 0.5.

Dependencies to be averaged are shown in Figure 2. 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 (𝑖𝑠𝑜-octane, toluene, and 𝑛-heptane in volumetric proportions of 55%, 35%, and 10%, resp.) developed by Golovitchev et al. [34, 35]. Dependencies of 𝑆_𝐿 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].

(a)

(b)

Figure 2 

Data for averaging invoking presumed PDFs. (a) Normalized heat release rates 𝑤(𝑐) versus combustion progress variable 𝑐, calculated for various gasoline-air mixtures, with the equivalence ratio being specified in legends. (b) Normalized laminar flame speed 𝑠𝐿 (solid curve) and maximum heat release rate Ω (dashed curve) versus equivalence ratio 𝐹.

(a)

(b)

Figure 3 

Validation of the chemical mechanism. Symbols show experimental data, and curves have been computed. (a) Laminar flame speed versus equivalence ratio. (1)𝑇𝑢=298K, data by Metghalchi and Keck [29]. (2)𝑇𝑢=350K, data by Metghalchi and Keck [29]. (3)𝑇𝑢=353K, data by Zhao et al. [30]. (4)𝑇𝑢=500K, data by Zhao et al. [30]. (b) Laminar flame speeds versus pressure at three different equivalence ratios specified in legends. 𝑇𝑢=373K. Data by Jerzembeck et al. [31].

We will firstly consider premixed flames. In this case, the equivalence ratio is stationary and spatially uniform and the 𝑓-dependent part of a presumed PDF is reduced to a single Dirac delta function at a given 𝑓₀, for example, 𝑓₁=𝑓₂=𝑓₀ in (7). Accordingly, we will compare the normalized mean heat release rates 𝑤𝑐,𝑐^2,𝑓₀=¹⁰𝑤𝑐,𝑓₀𝑃𝑐,𝑐,𝑐^2𝑑𝑐,(14) the normalized densities 𝜌𝑐,𝑐^2,𝑓₀=¹⁰𝜌𝑐,𝑓₀𝑃𝑐,𝑐,𝑐′²𝑑𝑐,(15) and the normalized turbulent burning velocities 𝑢_𝑡=𝛿_𝑡^∞−∞𝑤𝑑𝜁=𝛿_𝑡¹⁰𝑤𝑑𝑐𝑑𝜁⁻¹𝑑𝑐(16) obtained for various given first 𝑐 and second 𝑐^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 𝑐 in all the three cases. Here, 𝜌=11+𝜏𝑐,(17)𝜏=𝜌^𝑏−1−1 is the heat release parameter, 𝜁=(𝑧−𝑧_𝑓)/𝛿_𝑡 is the normalized spatial distance 𝑧 counted along the normal to the mean flame brush, which has a presumed thickness 𝛿_𝑡, which was not varied in the present simulations, 𝑧_𝑓 is associated with 𝑐=0.5, and the well-known complementary-error-function profile [28] 𝑐=1−12erfc√𝜋𝜁(18) 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 𝑞𝑓,𝑓′²=¹⁰𝑞(𝑓)𝑃𝑓,𝑓,𝑓^2𝑑𝑓(19) obtained for various given first 𝑓 and second 𝑓′² 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 𝑓 in all the three cases. Here, 𝑞 is equal to either 𝜌_𝑏(𝑓), or to 𝑠_𝐿(𝑓), or to 𝑠^𝐿1/2(𝑓), or to the normalized maximum (for various 𝑐) heat release rate 𝜔(𝑓)=Ω(𝑓)/max{Ω(𝑓)}, where max{Ω(𝑓)} is the maximum value of Ω(𝑓) determined by varying 𝑓, while Ω(𝑓) is the maximum value of 𝑊(𝑓,𝑐) calculated by varying only 𝑐.

Note that averaging 𝑤(𝑐)-curves implies that turbulence does not change the dependence of the normalized heat release rate on the normalized temperature 𝑐 in a premixed flame. Moreover, averaging 𝑠_𝐿(𝑓), 𝑠^𝐿1/2(𝑓), and 𝜔(𝑓) 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 𝑃(𝜉,𝜉,𝜉^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 𝜉′³ for the same 𝜉 and the same 𝑔 if 𝑔≥0.5 (cf. curves 1 and 2 in Figure 4(c) and note that curves 1 and 3 are indistinguishable). When 𝑔 is decreased, the difference in 𝜉^3 given by the two PDFs is increased (see Figures 4(a) and 4(b)).

(a) 𝑔=0.05

(b) 𝑔=0.25

(c) 𝑔=0.5

Figure 4 

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) 𝑔=0.05, (b) 𝑔=0.25, (c) 𝑔=0.5.

3. Results: Combustion Progress Variable PDF

Figures 5 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 𝑔 (see Figure 5(a)), but is substantial even at 𝑔=0.25 (see Figure 5(b)) and becomes extremely strong as 𝑔→1, see Figure 5(c).

(a) 𝑔=0.05

(b) 𝑔=0.25

(c) 𝑔=0.9

Figure 5 

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). 𝐹=1.0. (a) 𝑔=0.05, (b) 𝑔=0.25, (c) 𝑔=0.9.

Figure 6 

Dependencies of the normalized mean heat release rate on the Reynolds-averaged combustion progress variable, calculated for various equivalence ratios specified in legends and 𝑔=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 like 𝑝₀+𝑝₁𝜉+𝑝₂𝜉²+𝑝₃𝜉³, its mean value controlled by the first three moments should be the same for the two PDFs. Because the nonlinearity of 𝑤(𝑐) 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 𝑤(𝑐) is enhanced by an increase in |1−𝐹|.

Figure 5(b) indicates that the use of the 𝛿-PDFs may yield two peaks in the obtained dependencies of 𝑤(𝑐). This observation is explained in Figure 7. If the 2𝛿-PDF is invoked and 𝑐 is sufficiently large, then (10)–(12) may not have a solution consistent with 𝑐₂≤1 (thin dashed curve does not reach 𝑐=1 in Figure 7). At slightly lower 𝑐, the parameter 𝑐₁ rapidly grows with 𝑐 (see thin dashed curve) and reaches a value close to 𝑐_𝑚=0.685 that is associated with the maximum 𝑤(𝑐) for 𝐹=1, see solid curve in Figure 2(a). Accordingly, the first term on the RHS of (8) is significantly increased by 𝑐 as 𝑐₁→𝑐_𝑚 and the mean rate 𝑤 calculated using this equation is also increased by 𝑐 yielding the second peak of the 𝑤(𝑐)-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 𝑐₁ is substantially less than 𝑐_𝑚 and 𝑤(𝑐₁)≪max{𝑤(𝑐)}. Accordingly, an increase in 𝑤(𝑐₁) by 𝑐 on the RHS of (8) is less pronounced and counteracted by a decrease in 𝛼 as 𝑐 tends to unity, see (13). Moreover, when the R-2𝛿-PDF was invoked, we obtained 0≤𝑐₁<𝑐₂≤1 in all studied cases, contrary to the 2𝛿-PDF.

Figure 7 

Dependencies of the parameters 𝑐1 (thin curves) and 𝑐2 (bold curves) on the Reynolds-averaged combustion progress variable, calculated for 𝐹=1.0 and 𝑔=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 velocity 𝑢_𝑡 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 𝑢_𝑡 with 𝑔 does not mean a decrease in the turbulent burning velocity with the segregation factor, because the mean turbulent flame brush thickness, used to evaluate 𝑢_𝑡, may also depend on 𝑔. For the goals of the present work (i.e., comparison of various presumed PDFs), eventual dependence of 𝛿_𝑡 on 𝑔 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 integral 𝑢𝑡 on the segregation factor 𝑐2/[𝑐(1−𝑐)], calculated for various equivalence ratios specified in legends using the beta-function PDF (thin curves) and the 2𝛿-PDF (bold curves).

Figure 8 also shows that the dependencies of 𝑢_𝑡(𝑔) 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 𝑤(𝑐)-curves and (ii) decreasing their maxima. The former trend is associated with a weaker dependence of 𝑐₂ on 𝑐 at a larger 𝑔, see bold curves in Figure 9(b) and note that the second term dominates on the RHS of (8) under typical conditions. At moderate 𝑔 and small (large) 𝑐, the parameter 𝑐₂ is lower (higher) than the value of 𝑐_𝑚=0.868 associated the peak 𝑤(𝑐) for the considered mixture (𝐹=1.6). Therefore, there is a 𝑐 such that 𝑐₂(𝑐)=𝑐_𝑚, with this 𝑐 decreasing when the segregation factor is increased, see bold curves in Figure 9(b). Accordingly, the parameter 𝛼 evaluated for this particular 𝑐 increases, see Figure 9(c); that is, the second term on the RHS of (8) decreases, thus, reducing max{𝑤(𝑐)}, which is roughly equal to [1−𝛼(𝑐=𝑐_𝑚)]𝑤(𝑐_𝑚). However, at a larger 𝑔, this trend is overwhelmed by the widening of the computed 𝑤(𝑐)-curves and 𝑢_𝑡∝∫^∞−∞𝑤𝑑𝜁 is increased by 𝑔. To the contrary, if the segregation factor is sufficiently large, for example, 𝑔=0.8, see dotted-dashed curves in Figure 9, then 𝑐₂>𝑐_𝑚 even for low 𝑐 and a further increase in 𝑔 reduces not only (1−𝛼), but also 𝑤(𝑐₂) on the RHS of (8). Accordingly, a decrease in max{𝑤(𝑐)} with increasing 𝑔 becomes much more pronounced, overwhelms the widening of the computed 𝑤(𝑐)-curves, and results in decreasing 𝑢_𝑡.

(a)

(b)

(c)

Figure 9 

Dependencies of (a) the normalized mean heat release rate, (b) parameters 𝑐1 (thin curves) and 𝑐2 (bold curves), and (c) parameter 𝛼 on the Reynolds-averaged combustion progress variable, calculated invoking the 2𝛿-PDF for 𝐹=1.6 and different segregation factors specified in legends.

The above results, particularly Figures 5, 6, and 8, imply that the use of a presumed PDF 𝑃(𝑐) 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 𝑔<0.2 (cf. thin and bold curves in Figure 8), but a typical premixed turbulent flame is characterized by a larger 𝑔, as reviewed elsewhere, see [36, Section 3.1].

As far as the mean density is concerned, the studied presumed PDFs yield similar dependencies of 𝜌(𝑐) at various 𝑔 (cf. thin curves 1, 2, and 3 in Figure 10), because the dependence of 𝜌 on 𝑐, given by (17) is weakly nonlinear. However, at moderate 𝑔, the obtained 𝜌(𝑐)-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 𝑐̃(𝑐)-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 𝑔 and any 𝑃(𝑐). Indeed, on the one hand, the Favre-averaged quantity 𝜌⁻¹ is equal to 1+𝜏𝑐̃ due to (17), but, on the other hand, the Favre-averaged 𝜌⁻¹ is equal to 𝜌⁻¹ by definition of a Favre-averaged quantity 𝑞̃≡𝜌𝑞/𝜌.

(a) 𝑔=0.25

(b) 𝑔=0.5

(c) 𝑔=0.75

Figure 10 

Dependencies of the normalized mean density (thin curves) and Favre-averaged combustion progress variable 𝑐̃ (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 𝐹=1.0 and (a) 𝑔=0.25, (b) 𝑔=0.5, and (c) 𝑔=0.75.

As far as each of the two presumed PDFs investigated here is concerned separately, the sum of Dirac delta functions at 𝑐₁ and 𝑐₂ does not seem to be a proper approach to average 𝑤(𝑐) for the following reasons.

First, Figure 11 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 𝑐, calculated in the flamelet regime of premixed turbulent combustion (𝑔=0.9). In this case, 𝑐₁≪1 (dashed curves) and 1−𝑐₂≪1 (dotted-dashed curves). Therefore, both 𝑤(𝑐₁) and 𝑤(𝑐₂) 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 𝑔=0.9 and 𝐹=1.0. (1)𝛼, (2)𝑐1, (3)𝑐2.

Second, using (9) with 𝜉=𝑐, one can easily obtain 𝑐=𝛼𝑐₁+(1−𝛼)𝑐₂,𝑐1−𝑐=𝛼𝑐₁+(1−𝛼)𝑐₂1−𝛼𝑐₁−(1−𝛼)𝑐₂,𝑐^2=𝛼𝑐₁−𝑐²+(1−𝛼)𝑐₂−𝑐²=𝛼(1−𝛼)²𝑐₁−𝑐₂²+𝛼²(1−𝛼)𝑐₁−𝑐₂²=𝛼(1−𝛼)𝑐₁−𝑐₂².(20) In the limit case of 𝑔→1 and, hence, 𝑐′²→𝑐(1−𝑐), (20) results in 𝛼⟶𝑐²1−𝑐₂𝑐₂1−𝑐₂−𝑐₁1−𝑐₁.(21) Because 𝑐₁(1−𝑐₁)≥0, the RHS of (21) may be within the limits of [0,1] if either (i) 𝑐₁=0 and 𝑐₂<1, or (ii) 𝑐₁>0 and 𝑐₂=1, or (iii) 𝑐₁=0 and 𝑐₂=1. For brevity, we consider 𝑐₁≤𝑐₂ here. In cases (i) and (ii), 𝛼=1 and 𝛼=0, respectively, and the first equality in (20) does not hold. In case (iii), 𝑤(𝑐₁)=𝑤(𝑐₂)=0 and (8) yields 𝑤=0. Thus, the PDF given by (9) cannot be used to average 𝑤(𝑐) in the flamelet regime of premixed turbulent combustion.

The beta-function PDF 𝑃(𝑐) does not seem to be well tailored for simulating the flamelet regime either. Indeed, (4) yields 𝑎≪1 and 𝑏≪1 in the case of 1−𝑔≪1. Accordingly, even a small error in computing 𝑐 and/or 𝑐^2 by solving proper balance equations can result in substantial errors in calculating the small parameters 𝑎 and 𝑏 and, therefore, substantial errors in evaluating 𝑤 by invoking the beta-function PDF 𝑃(𝑐). Too small parameters 𝑎 and 𝑏 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 𝑐^2 with a balance equation for 1−𝑔 [37].

Finally, it is worth noting the following interesting feature of the use of a presumed beta-function PDF 𝑃(𝑐). 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 𝑐̃ and 𝜌𝑐^2/𝜌. This could easily be done if the RHS of (3) is assumed to approximate the Favre PDF 𝑃=𝜌𝑃/𝜌. In such a case, 𝜌𝑐^𝑛/𝜌 and 𝑤 evaluated using 𝑃(𝑐) are equal to 𝑐^𝑛 and 𝑤, respectively, calculated invoking 𝑃(𝑐) provided that both 𝑃(𝑐) and 𝑃(𝑐) are approximated by the same beta function.

However, the above substitution of 𝑃(𝑐) with 𝑃(𝑐) substantially affects the shape of computed 𝑤(𝑐)-curves. To show this effect, we compared the following two ways of averaging the heat release rate. First, we specified 𝑐 and 𝑔_𝑅≡𝑐^2/[𝑐(1−𝑐)] and calculated 𝑤(𝑐)-curve by varying 𝑐, keeping 𝑔_𝑅 constant, and using (3), (4), and (14). Second, we specified 𝑐̃ and 𝑔_𝐹≡𝜌𝑐^2/[𝜌𝑐̃(1−𝑐̃)] and computed 𝑤(𝑐̃)-curve by varying 𝑐̃, keeping 𝑔_𝐹 constant and equal to the above 𝑔_𝑅, and using the following equations: 𝑤==∫∫¹⁰𝑤(𝑐)𝑃(𝑐)𝑑𝑐₁₀𝑤(𝑐)𝜌𝑃(𝑐)=Γ(𝑎+𝑏)𝜌𝑑𝑐Γ(𝑎)Γ(𝑏)∫¹⁰𝑤(𝑐)𝜌𝜌𝑐^𝑎−1𝑎=𝑐̃𝑔_𝐹−1(1−𝑐)^𝑏−1𝑑𝑐,−1,𝑏=(1−𝑐̃)𝑔^𝐹−1−1.(22) Subsequently, we transformed the latter 𝑤(𝑐̃)-curve to 𝑤(𝑐)-curve using 𝑐𝑐̃,𝑔_𝐹=Γ(𝑎+𝑏)Γ(𝑎)Γ(𝑏)¹⁰𝑐𝜌𝜌𝑐^𝑎−1(1−𝑐)^𝑏−1𝑑𝑐.(23) Results are shown in Figure 12 in thin and bold curves, respectively. While the 𝑤(𝑐)-curves computed invoking 𝑃(𝑐,𝑐,𝑐^2) look symmetrical with respect to 𝑐=0.5 and reach peak values around 𝑐=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 𝑐 is significantly larger than the Favre-averaged 𝑐̃ 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 for 𝑔=0.9 and various equivalence ratios, specified in legends, by invoking the beta-function PDF. Thin and bold curves are associated with 𝑃(𝑐,𝑐,𝑐2) and 𝑃(𝑐,𝑐̃,𝜌𝑐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 factor 𝑔. 𝐹=1.0. Thin and bold curves are associated with 𝑃(𝑐,𝑐,𝑐2) and 𝑃(𝑐,𝑐̃,𝜌𝑐2/𝜌), respectively.

These numerical data and the aforementioned experimental results that indicate almost symmetrical dependencies of the FSD on 𝑐 favor the use of 𝑃(𝑐,𝑐,𝑐^2) at least at larger 𝑔. However, such a method is more complicated, because it requires evaluation of 𝑐 and 𝑐^2 based on computed 𝑐̃ and 𝜌𝑐^2/𝜌.

4. Results: Mixture Fraction PDF

Figure 14 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 𝜌_𝑏 of burned gas, the mean normalized maximum heat release rate 𝜔, and 𝑠^𝐿1/2 on the mean equivalence ratio 𝐹 even if 𝑔_𝑓≡𝑓^2/[𝑓(1−𝑓)] is as large as 0.4. It is worth remembering that, contrary to fluctuations in the combustion progress variable, which may be characterized by 𝑔→1, fluctuations in mixture fraction are characterized by significantly lower 𝑔_𝑓 in a typical stratified turbulent flame.

(a)

(b)

(c)

Figure 14 

Dependencies of the normalized mean density 𝜌𝑏 of burned gas (a), the normalized maximum heat release rate 𝜔 (b), and the square root 𝑠𝐿1/2 of the normalized laminar flame speed (c) on the mean equivalence ratio 𝐹, 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). 𝐹min=0.5 and 𝐹max=1.9,𝑔𝑓=0.4.

The difference between results obtained by invoking the considered presumed PDFs is decreased when 𝑔_𝑓 is decreased, for example, see Figure 15 and note that a similar trend is observed for other simulated quantities 𝑞={𝜌_𝑏,𝜔,𝑠_𝐿,𝑠^𝐿1/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 speed 𝑠𝐿 on the mean equivalence ratio 𝐹, calculated invoking the beta-function (thin curves) and R-2𝛿 (bold curves) PDFs for different 𝑔𝑓 specified in legends. 𝐹min=0.5 and 𝐹max=1.9.

It is worth also stressing that the values of 𝜌_𝑏(𝐹), 𝜔(𝐹), and 𝑠^𝐿1/2(𝐹), 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 𝑞(𝑓) on the shape of presumed PDF 𝑃(𝑓), but significant difference in 𝑞(𝑓) and 𝑞(𝑓)) support invoking a presumed PDF 𝑃(𝑓) for simulating stratified turbulent combustion if fluctuations in mixture fraction are weak or moderate, for example, 𝑔≤0.4.

However, such a conclusion is correct only if the equivalence ratio always fluctuates within flammability limits 𝐹_𝑙<𝐹<𝐹_𝑟. If the amplitude of fluctuations in 𝐹 is increased so that inflammable mixture composition may be observed with a finite probability 𝑝_𝑞=∫^𝐹𝑙_𝐹min𝑃(𝐹)𝑑𝐹+∫^𝐹max_𝐹𝑟𝑃(𝐹)𝑑𝐹, then computed 𝑞(𝑓) is more sensitive to the shape of presumed PDF. For instance, Figure 16 indicates that a decrease in 𝐹_min from 0.5 to 0.4 and an increase in 𝐹_max from 1.9 to 2.2 make 𝑠^𝐿1/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 root 𝑠𝐿1/2 of the normalized mean laminar flame speed on the mean equivalence ratio 𝐹, calculated invoking the beta-function (thin curves) and the R-2𝛿 (bold curves) PDFs. 𝑔=0.2. (1)𝐹min=0.5 and 𝐹max=1.9, (2)𝐹min=0.4 and 𝐹max=2.2, (3)𝐹min=0.3 and 𝐹max=2.5, (4)𝐹min=0.1 and 𝐹max=4.0.

If the sum of two Dirac delta functions is invoked, then 𝑠^𝐿1/2 drops at certain 𝐹₁<1 and 𝐹₂>1 (see bold solid line), because both 𝑓₁ and 𝑓₂ are increased by 𝑓 (this trend is shown for 𝜉=𝑐 in Figure 11) and 𝑓₁ (or 𝑓₂) reaches the lean (or rich) flammability limit when 𝐹=𝐹₁ (or 𝐹=𝐹₂). Accordingly, if 𝐹<𝐹₁ (or 𝐹>𝐹₂), then 𝑓₁ (or 𝑓₂) is outside the lean (or rich) flammability limit, hence, 𝑠_𝐿(𝑓₁)=0 [or 𝑠_𝐿(𝑓₂)=0], and, consequently, 𝑠^𝐿1/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 𝑠^𝐿1/2 on 𝐹, calculated by invoking the beta-function and R-2𝛿 PDFs, differ substantially from one another if, for example, 𝐹_min=0.1 and 𝐹_max=4.0 (cf. thin and bold dotted curves).

Figure 17 shows that the probability 𝑝_𝑓=1−𝑝_𝑞 of finding flammable mixture is very sensitive to the shape of presumed PDF 𝑃(𝑓). If the beta-function PDF is invoked, 𝑝_𝑓 depends smoothly on 𝐹 (thin curves), whereas the probability calculated using the sum of two Dirac delta functions drops sharply (bold curves) as 𝑓₁ or/and 𝑓₂ moves beyond the flammable range 𝑓_𝑙<𝑓<𝑓_𝑟.

Figure 17 

Dependencies of the probability of finding flammable mixture composition on the mean equivalence ratio 𝐹, calculated invoking the beta-function (thin curves) and the R-2𝛿 (bold curves) PDFs. 𝑔=0.2. (1)𝐹min=0.4 and 𝐹max=2.2, (2)𝐹min=0.3 and 𝐹max=2,5, (3)𝐹min=0.1 and 𝐹max=4.0.

5. Conclusions

The use of a presumed combustion-progress-variable PDF 𝑃(𝑐) for averaging heat release rate in a premixed turbulent flame characterized by sufficiently large 𝑐^2/[𝑐(1−𝑐)]>0.4 does not seem to be justified, because the mean rate is sensitive to the shape of 𝑃(𝑐) even if first 𝑐, second 𝑐^2, and third 𝑐^3 moments are kept constant when changing the PDF shape.

Because 𝑞(𝑓) differs substantially from 𝑞(𝑓) and 𝑞(𝑓) depends sufficiently weak on the shape of a presumed mixture-fraction PDF 𝑃(𝑓,𝑓,𝑓^2), with 𝑓 and 𝑓^2 being kept constant and 𝑓^2/[𝑓(1−𝑓)]<0.4, the use of 𝑃(𝑓) 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.

If 𝑔>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.

Acknowledgment

This work was supported by the Combustion Engine Research Center (CERC).