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 āˆ«š‘Š=10š‘Š(š‘)š‘ƒ(š‘)š‘‘š‘ [6ā€“10], or āˆ«10š‘Š(Ģƒš‘,š‘“)š‘ƒ(š‘“)š‘‘š‘“ and šœŒš‘=āˆ«10šœŒš‘(š‘“)š‘ƒ(š‘“)š‘‘š‘“ [11ā€“13], or š‘†š‘”=āˆ«10š‘†š‘”(š‘“)š‘ƒ(š‘“)š‘‘š‘“ and šœŒš‘ [14ā€“19], or āˆ¬š‘Š=10š‘Š(š‘,š‘“)š‘ƒ(š‘,š‘“)š‘‘š‘š‘‘š‘“ and āˆ¬šœŒ=01šœŒ(š‘,š‘“)š‘ƒ(š‘,š‘“)š‘‘š‘š‘‘š‘“ [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 āˆ«10š‘ƒ(šœ‰)š‘‘šœ‰=1.

Equation (3) involves two unknown parameters š‘Ž(š±,š‘”) and š‘(š±,š‘”), which are commonly evaluated using the following expressions: š‘Ž(š±,š‘”)=šœ‰ī€·š‘”āˆ’1ī€øī‚€āˆ’1,š‘(š±,š‘”)=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, Ī“(š‘Ž)ā‰ˆš‘Žāˆ’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 [29ā€“31] 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)šœ=šœŒš‘āˆ’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š‘=1āˆ’2ī‚€āˆšerfcī‚šœ‹šœ(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 š‘2ā‰¤1 (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<š‘2ā‰¤1 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, š‘1ā‰Ŗ1 (dashed curves) and 1āˆ’š‘2ā‰Ŗ1 (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āˆ’š›¼)š‘2ī€»ī€ŗ1āˆ’š›¼š‘1āˆ’(1āˆ’š›¼)š‘2ī€»,š‘ī…ž2ī€·š‘=š›¼1āˆ’š‘ī€ø2+ī€·š‘(1āˆ’š›¼)2āˆ’š‘ī€ø2=š›¼(1āˆ’š›¼)2ī€·š‘1āˆ’š‘2ī€ø2+š›¼2ī€·š‘(1āˆ’š›¼)1āˆ’š‘2ī€ø2ī€·š‘=š›¼(1āˆ’š›¼)1āˆ’š‘2ī€ø2.(20) In the limit case of š‘”ā†’1 and, hence, š‘ā€²2ā†’š‘(1āˆ’š‘), (20) results in š‘š›¼āŸ¶2ī€·1āˆ’š‘2ī€øš‘2ī€·1āˆ’š‘2ī€øāˆ’š‘1ī€·1āˆ’š‘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 [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.

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.


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