#### 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]: 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 , 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, .

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 and 1 [4], does not allow us to average , because . Accordingly, the following discussion is solely restricted to PDFs presumed in the range of , while the BML PDF, which does not address the intermediate range of , 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 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 is used in order to satisfy the normalizing constraint of .

Equation (3) involves two unknown parameters and , which are commonly evaluated using the following expressions: where is a segregation factor. To do so, the first and second moments are determined by solving proper balance equations. If Favre-averaged balance equations are integrated to compute and , then (i) the right hand side (RHS) of (3) is used to approximate the Favre PDF and (ii) and are substituted with and , 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 and , respectively, see Figure 1. For instance, if (or ), then the PDF given by (3) tends to infinity and resembles Dirac delta function at (or ). However, contrary to Dirac delta function, substitution of (3) into the integral with yields vanishing probability of finding , for example, vanishing probability of finding unburned mixture in the case of premixed combustion (). Indeed, if and differs markedly from zero and unity, then, , , and the above integral is approximately equal to . Moreover, if , then, , because and . Therefore, in the considered example, the integral vanishes as , 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 was introduced by Libby and Williams [27] to simulate partially premixed flames. The second term on the RHS is multiplied with factor in order to satisfy the normalizing constraint of .

The use of (7) substantially simplifies averaging reaction rates, for example, if we consider a double-delta-function PDF for a single random variable . However, evaluation of the parameters , , is more difficult, because the first two moments and are not sufficient to determine the three unknown parameters and the system of three nonlinear algebraic equations and, for example, should be solved for this purpose. To avoid such numerical complications, Ribert et al. [22] have proposed a simplified model that yields 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 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 , 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 , where and are presumed limits of fluctuations in the equivalence ratio, is the heat release rate normalized using its maximum value calculated by varying the combustion progress variable, is the laminar flame speed normalized using its maximum value 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)**

**(a)**

**(b)**

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 the normalized densities and the normalized turbulent burning velocities 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 integrals on the RHSs, with being associated with in all the three cases. Here, 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 , and the well-known complementary-error-function profile [28] 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 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 , or to the normalized maximum (for various ) heat release rate , where 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 , , 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 , it was evaluated invoking the beta-function PDF 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 calculated using the beta-function PDF with the same and the same . However, obtained numerical results indicate that the beta-function and R-2 PDFs yield almost equal for the same and the same if (cf. curves 1 and 2 in Figure 4(c) and note that curves 1 and 3 are indistinguishable). When is decreased, the difference in given by the two PDFs is increased (see Figures 4(a) and 4(b)).

(a) |

(b) |

(c) |

#### 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 (see Figure 5(b)) and becomes extremely strong as , see Figure 5(c).

(a) |

(b) |

(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 , 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 .

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 (thin dashed curve does not reach in Figure 7). At slightly lower , the parameter rapidly grows with (see thin dashed curve) and reaches a value close to that is associated with the maximum for , 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 . 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 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 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 associated the peak for the considered mixture (. 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 , which is roughly equal to . 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, , see dotted-dashed curves in Figure 9, then even for low and a further increase in reduces not only , but also on the RHS of (8). Accordingly, a decrease in with increasing becomes much more pronounced, overwhelms the widening of the computed -curves, and results in decreasing .

**(a)**

**(b)**

**(c)**

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 (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 due to (17), but, on the other hand, the Favre-averaged is equal to by definition of a Favre-averaged quantity .

(a) |

(b) |

(c) |

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 (). In this case, (dashed curves) and (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.

Second, using (9) with , one can easily obtain In the limit case of and, hence, , (20) results in Because , the RHS of (21) may be within the limits of if either (i) and , or (ii) and , or (iii) and . For brevity, we consider here. In cases (i) and (ii), and , respectively, and the first equality in (20) does not hold. In case (iii), and (8) yields . 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 and in the case of . Accordingly, even a small error in computing and/or 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 with a balance equation for [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 . 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 and calculated -curve by varying , keeping constant, and using (3), (4), and (14). Second, we specified and and computed -curve by varying , keeping constant and equal to the above , and using the following equations: Subsequently, we transformed the latter -curve to -curve using Results are shown in Figure 12 in thin and bold curves, respectively. While the -curves computed invoking look symmetrical with respect to and reach peak values around , 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 at least at larger . However, such a method is more complicated, because it requires evaluation of and based on computed and .

#### 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 on the mean equivalence ratio even if is as large as . It is worth remembering that, contrary to fluctuations in the combustion progress variable, which may be characterized by , fluctuations in mixture fraction are characterized by significantly lower in a typical stratified turbulent flame.

**(a)**

**(b)**

**(c)**

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

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 , then computed is more sensitive to the shape of presumed PDF. For instance, Figure 16 indicates that a decrease in from 0.5 to 0.4 and an increase in from 1.9 to 2.2 make 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 drops at certain and (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, [or ], and, consequently, 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 on , calculated by invoking the beta-function and R-2 PDFs, differ substantially from one another if, for example, and (cf. thin and bold dotted curves).

Figure 17 shows that the probability 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 .

#### 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 does not seem to be justified, because the mean rate is sensitive to the shape of even if first , second , and third 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 , with and being kept constant and , 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 , 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).