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 [13]. 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 [13] 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 𝑊=10𝑊(𝑐)𝑃(𝑐)𝑑𝑐 [610], or 10𝑊(̃𝑐,𝑓)𝑃(𝑓)𝑑𝑓 and 𝜌𝑏=10𝜌𝑏(𝑓)𝑃(𝑓)𝑑𝑓 [1113], or 𝑆𝑡=10𝑆𝑡(𝑓)𝑃(𝑓)𝑑𝑓 and 𝜌𝑏 [1419], or 𝑊=10𝑊(𝑐,𝑓)𝑃(𝑐,𝑓)𝑑𝑐𝑑𝑓 and 𝜌=01𝜌(𝑐,𝑓)𝑃(𝑐,𝑓)𝑑𝑐𝑑𝑓 [2025]. Note that the use of the 𝐺-equation to simulate partially premixed burning [1419] 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 [613, 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 10𝑃(𝜉)𝑑𝜉=1.

Equation (3) involves two unknown parameters 𝑎(𝐱,𝑡) and 𝑏(𝐱,𝑡), which are commonly evaluated using the following expressions: 𝑎(𝐱,𝑡)=𝜉𝑔11,𝑏(𝐱,𝑡)=1𝜉𝑔11,(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, Γ(𝑎)𝑎1, 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.

Second, another presumed PDF, that is, a combination of Dirac delta functions 𝛿(𝜉)𝑃(𝑐,𝑓,𝐱,𝑡)=𝛼(𝐱,𝑡)𝛿𝑐𝑐1𝛿(𝐱,𝑡)𝑓𝑓1+[]𝛿(𝐱,𝑡)1𝛼(𝐱,𝑡)𝑐𝑐2(𝛿𝐱,𝑡)𝑓𝑓2(𝐱,𝑡)(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 10𝑃(𝑐,𝑓)𝑑𝑐𝑑𝑓=1.

The use of (7) substantially simplifies averaging reaction rates, for example, 𝜉𝑊(𝐱,𝑡)=𝛼(𝐱,𝑡)𝑊1(+[]𝑊𝜉𝐱,𝑡)1𝛼(𝐱,𝑡)2(𝐱,𝑡)(8) if we consider a double-delta-function PDF 𝑃𝜉(𝜉,𝐱,𝑡)=𝛼(𝐱,𝑡)𝛿1+[]𝛿𝜉(𝐱,𝑡)𝜉1𝛼(𝐱,𝑡)2(𝐱,𝑡)𝜉(9) for a single random variable 𝜉. However, evaluation of the parameters 𝛼(𝐱,𝑡), 𝜉1(𝐱,𝑡), 𝜉2(𝐱,𝑡) 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[]𝜉(𝐱,𝑡)+1𝛼(𝐱,𝑡)2(𝐱,𝑡),𝜉2(𝐱,𝑡)=𝜉2(𝐱,𝑡)𝜉2(𝐱,𝑡)(10)=𝛼(𝐱,𝑡)𝜉21([]𝜉𝐱,𝑡)+1𝛼(𝐱,𝑡)22(𝐱,𝑡)𝜉2(𝐱,𝑡),(11) and, for example, 𝜉3(𝐱,𝑡)=𝜉3(𝐱,𝑡)3𝜉(𝐱,𝑡)𝜉2(𝐱,𝑡)𝜉3(𝐱,𝑡)=𝛼(𝐱,𝑡)𝜉31[]𝜉(𝐱,𝑡)+1𝛼(𝐱,𝑡)32(𝐱,𝑡)3𝜉(𝐱,𝑡)𝜉2(𝐱,𝑡)𝜉3(𝐱,𝑡)(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𝑔1/2,𝜉(𝐱,𝑡)2(𝐱,𝑡)=𝜉(𝐱,𝑡)+1𝑔𝜉(𝐱,𝑡)1/2(𝐱,𝑡),𝛼(𝐱,𝑡)=1𝜉(𝐱,𝑡)(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 [2931] in Figure 3, while results of validating the mechanism against published data on ignition delay times are reported elsewhere [35].

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 𝑓0, for example, 𝑓1=𝑓2=𝑓0 in (7). Accordingly, we will compare the normalized mean heat release rates 𝑤𝑐,𝑐2,𝑓0=10𝑤𝑐,𝑓0𝑃𝑐,𝑐,𝑐2𝑑𝑐,(14) the normalized densities 𝜌𝑐,𝑐2,𝑓0=10𝜌𝑐,𝑓0𝑃𝑐,𝑐,𝑐2𝑑𝑐,(15) and the normalized turbulent burning velocities 𝑢𝑡=𝛿𝑡𝑤𝑑𝜁=𝛿𝑡10𝑤𝑑𝑐𝑑𝜁1𝑑𝑐(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, 1𝜌=,1+𝜏𝑐(17)𝜏=𝜌𝑏11 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=10𝑞(𝑓)𝑃𝑓,𝑓,𝑓2𝑑𝑓(19) 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 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 𝜉3 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)).

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).

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 𝑝0+𝑝1𝜉+𝑝2𝜉2+𝑝3𝜉3, 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 𝑐21 (thin dashed curve does not reach 𝑐=1 in Figure 7). At slightly lower 𝑐, the parameter 𝑐1 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 𝑐1𝑐𝑚 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 𝑐1 is substantially less than 𝑐𝑚 and 𝑤(𝑐1)max{𝑤(𝑐)}. Accordingly, an increase in 𝑤(𝑐1) 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<𝑐21 in all studied cases, contrary to the 2𝛿-PDF.

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 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 𝑐2 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 𝑐2 is lower (higher) than the value of 𝑐𝑚=0.868 associated the peak 𝑤(𝑐) for the considered mixture (𝐹=1.6). Therefore, there is a 𝑐 such that 𝑐2(𝑐)=𝑐𝑚, 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 𝑐2>𝑐𝑚 even for low 𝑐 and a further increase in 𝑔 reduces not only (1𝛼), but also 𝑤(𝑐2) 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 𝑢𝑡.

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 𝜌1 is equal to 1+𝜏̃𝑐 due to (17), but, on the other hand, the Favre-averaged 𝜌1 is equal to 𝜌1 by definition of a Favre-averaged quantity ̃𝑞𝜌𝑞/𝜌.

As far as each of the two presumed PDFs investigated here is concerned separately, the sum of Dirac delta functions at 𝑐1 and 𝑐2 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, 𝑐11 (dashed curves) and 1𝑐21 (dotted-dashed curves). Therefore, both 𝑤(𝑐1) and 𝑤(𝑐2) 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.

Second, using (9) with 𝜉=𝑐, one can easily obtain 𝑐=𝛼𝑐1+(1𝛼)𝑐2,𝑐1𝑐=𝛼𝑐1+(1𝛼)𝑐21𝛼𝑐1(1𝛼)𝑐2,𝑐2𝑐=𝛼1𝑐2+𝑐(1𝛼)2𝑐2=𝛼(1𝛼)2𝑐1𝑐22+𝛼2𝑐(1𝛼)1𝑐22𝑐=𝛼(1𝛼)1𝑐22.(20) In the limit case of 𝑔1 and, hence, 𝑐2𝑐(1𝑐), (20) results in 𝑐𝛼21𝑐2𝑐21𝑐2𝑐11𝑐1.(21) Because 𝑐1(1𝑐1)0, the RHS of (21) may be within the limits of [0,1] if either (i) 𝑐1=0 and 𝑐2<1, or (ii) 𝑐1>0 and 𝑐2=1, or (iii) 𝑐1=0 and 𝑐2=1. For brevity, we consider 𝑐1𝑐2 here. In cases (i) and (ii), 𝛼=1 and 𝛼=0, respectively, and the first equality in (20) does not hold. In case (iii), 𝑤(𝑐1)=𝑤(𝑐2)=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: 𝑤=10=𝑤(𝑐)𝑃(𝑐)𝑑𝑐10𝑤(𝑐)𝜌𝑃(𝑐)𝜌=𝑑𝑐Γ(𝑎+𝑏)Γ(𝑎)Γ(𝑏)10𝑤(𝑐)𝜌𝜌𝑐𝑎1(1𝑐)𝑏1𝑔𝑑𝑐,𝑎=̃𝑐𝐹1𝑔1,𝑏=(1̃𝑐)𝐹1.1(22) Subsequently, we transformed the latter 𝑤(̃𝑐)-curve to 𝑤(𝑐)-curve using 𝑐̃𝑐,𝑔𝐹=Γ(𝑎+𝑏)Γ(𝑎)Γ(𝑏)10𝑐𝜌𝜌𝑐𝑎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 [3843] 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.

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.

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.

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.

If the sum of two Dirac delta functions is invoked, then 𝑠𝐿1/2 drops at certain 𝐹1<1 and 𝐹2>1 (see bold solid line), because both 𝑓1 and 𝑓2 are increased by 𝑓 (this trend is shown for 𝜉=𝑐 in Figure 11) and 𝑓1 (or 𝑓2) reaches the lean (or rich) flammability limit when 𝐹=𝐹1 (or 𝐹=𝐹2). Accordingly, if 𝐹<𝐹1 (or 𝐹>𝐹2), then 𝑓1 (or 𝑓2) is outside the lean (or rich) flammability limit, hence, 𝑠𝐿(𝑓1)=0 [or 𝑠𝐿(𝑓2)=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 𝑓1 or/and 𝑓2 moves beyond the flammable range 𝑓𝑙<𝑓<𝑓𝑟.

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).