An extended flamelet-based presumed probability density function for predicting mean concentrations of various species in premixed turbulent flames

Abstract Direct Numerical Simulation (DNS) data obtained by Dave and Chaudhuri (2020) from a lean, complex-chemistry, hydrogen-air flame associated with the thin-reaction-zone regime of premixed turbulent burning are analyzed to perform a priori assessment of predictive capabilities of the flamelet approach for evaluating mean species concentrations. For this purpose, dependencies of mole fractions and rates of production of various species on a combustion progress variable c , obtained from the laminar flame, are averaged adopting either the actual Probability Density Function (PDF) P ( c ) extracted from the DNS data or a common presumed β -function PDF. On the one hand, the results quantitatively validate the flamelet approach for the mean mole fractions of all species, including radicals, but only if the actual PDF P ( c ) is adopted. The use of the β -function PDF yields substantially worse results for the radicals’ concentrations. These findings put modeling the PDF P ( c ) on the forefront of the research agenda. On the other hand, the mean rate of product creation and turbulent burning velocity are poorly predicted even adopting the actual PDF. These results imply that, in order to evaluate the mean species concentrations, the flamelet approach could be coupled with another model that predicts the mean rate and turbulent burning velocity better. Accordingly, the flamelet approach could be implemented as post-processing of numerical data yielded by that model. Based on the aforementioned findings and implications, a new approach to building a presumed PDF is developed. The key features of the approach consist in (i) adopting a re-normalized flamelet PDF for intermediate values of c and (ii) directly using the mean rate of product creation to calibrate the presumed PDF. Capabilities of the newly developed PDF for predicting mean species concentrations are quantitively validated for all species, including radicals.

Presumed probability density function Modeling Direct numerical simulation a b s t r a c t Direct Numerical Simulation (DNS) data obtained by Dave and Chaudhuri (2020) from a lean, complex-chemistry, hydrogen-air flame associated with the thin-reaction-zone regime of premixed turbulent burning are analyzed to perform a priori assessment of predictive capabilities of the flamelet approach for evaluating mean species concentrations. For this purpose, dependencies of mole fractions and rates of production of various species on a combustion progress variable c, obtained from the laminar flame, are averaged adopting either the actual Probability Density Function (PDF) PðcÞ extracted from the DNS data or a common presumed b-function PDF. On the one hand, the results quantitatively validate the flamelet approach for the mean mole fractions of all species, including radicals, but only if the actual PDF PðcÞ is adopted. The use of the b-function PDF yields substantially worse results for the radicals' concentrations. These findings put modeling the PDF PðcÞ on the forefront of the research agenda. On the other hand, the mean rate of product creation and turbulent burning velocity are poorly predicted even adopting the actual PDF. These results imply that, in order to evaluate the mean species concentrations, the flamelet approach could be coupled with another model that predicts the mean rate and turbulent burning velocity better. Accordingly, the flamelet approach could be implemented as post-processing of numerical data yielded by that model. Based on the aforementioned findings and implications, a new approach to building a presumed PDF is developed. The key features of the approach consist in (i) adopting a re-normalized flamelet PDF for intermediate values of c and (ii) directly using the mean rate of product creation to calibrate the presumed PDF. Capabilities of the newly developed PDF for predicting mean species concentrations are quantitively validated for all species, including radicals.

Introduction
To protect the environment and mitigate the threat of global warming, there is a strong need for replacement of engines that utilize chemical energy bound in fossil fuels with highly efficient, flexible, and ultra clean engines capable for utilizing chemical energy bound in renewable carbon-free fuels. Hydrogen is considered to be one of the most attractive carbon-free fuels due to a variety of H 2 production technologies [1e5] and unique characteristics of hydrogen flames, such as a low ignition energy, a wide range of flammability limits, a high laminar burning velocity, etc. [6]. Accordingly, burning of hydrogen in various engines has been attracting significant amount of attention [6e14], thus, making combustion of hydrogen (including emissions from hydrogen flames [15,16]) an important and rapidly growing area of hydrogen energy. These recent developments have been motivating fundamental research into basic characteristics of laminar combustion of hydrogen [17e20]  To increase an impact of the fundamental science on applied research into future ultra-clean and highly efficient engines that burn hydrogen, there is a strong need for development of Computational Fluid Dynamics (CFD) models that (i) allow for complex (and fuel-specific) chemistry of turbulent combustion, and (ii) are capable for predicting mean concentrations of various (not only major reactants and products, but also radicals such as O, OH, H, etc.) species in a turbulent premixed flame under a wide range of conditions. Accordingly, over the past two decades, several approaches, e.g. a transport equation for a Probability Density Function (PDF) [47,48], thickened flame model [49], or flamelet concept [50] coupled with a presumed PDF [51e54], were advanced for predicting concentrations of various species in turbulent flames. In particular, the flamelet concept coupled with a presumed PDF is implemented into major commercial CFD codes and is widely used in applied research, e.g. see Table 4 in Ref. [55]. Such an approach is also in the focus of fundamental studies [56e82] and is supported by recent experimental and Direct Numerical Simulation (DNS) data that, as reviewed elsewhere [83,84], indicate that the domain of validity of the flamelet concept is significantly wider than earlier expected.
The flamelet concept consists in adapting results (the socalled flamelet library) of numerical simulations of a set of laminar premixed flames, performed by invoking a sufficiently detailed model of molecular transport and a sufficiently detailed chemical mechanism. By adopting Flamelet Prolongation of Intrinsic Low-Dimensional Manifolds (FPI) [85] or Flamelet Generation Manifolds (FGM) [86] techniques, a flamelet library is often stored in a form of dependencies of temperature T L ðxÞ, density r L ðxÞ, mole (or mass) fractions X n;L ðxÞ and rates W n;L ðxÞ of consumption/production of n ¼ 1; …; N species on a set x of independent variables. The set x may consist of a single combustion progress variable c, which varies from zero in fresh reactants to unity in equilibrium combustion products, but may also involve mixture fraction in the case of partially premixed combustion, mixture enthalpy in the case of non-adiabatic flame, pressure in the case of non-isobaric burning (e.g., in a piston engine), stretch rate in highly turbulent flames, etc. Henceforth, we restrict ourselves to the basic case of a single independent scalar variable c.
To adapt a flamelet library in a Reynolds Averaged Navier-Stokes (RANS) or Large Eddy Simulation (LES) study of turbulent burning, the library should be coupled with a model of the influence of turbulence on combustion. Such a coupling is commonly implemented (i) using the following 1 RANS equations (1) Tðx; tÞ ¼ (or basically similar LES equations, which are not written for brevity) for mean production/consumption rates, mole fractions, temperature, and density, respectively, and (ii) invoking a presumed PDF Pðc; x; tÞ. Here, W and X are N-dimensional vector-functions that encompass reaction rates W n and mole fractions X n , respectively, for 1 n N species. As reviewed elsewhere [87,88], the presumed PDF is commonly modeled by (i) assuming its general shape, which depends on a few unknown parameters, and (ii) evaluating these parameters by comparing the first two moments of the cðx; tÞ-field, calculated using the PDF, with the moments obtained by solving their transport equations. The PDF shape is 1 In applications, such equations are commonly adapted for the mass-weighted quantities and PDF. Here, conventional mean quantities and PDF are addressed for simplicity, and because the issues and solutions discussed in the following are basically similar for conventional and mass-weighted quantities. i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 5 ( 2 0 2 0 ) 3 1 1 6 2 e3 1 1 7 8 presumed adopting a sum of Dirac delta functions [52,56,58,63] or combinations of Dirac delta functions and a flamelet PDF [53,54,61,66,73,78], but the b-function PDF P b ðc; x; tÞ [51] is most widely accepted [53,57,59,60,62,64,65,67e72,74e77,79e82], because its shape is very flexible and, depending on the magnitudes of the aforementioned moments, P b ðc; x; tÞ can vary from a quasi-bimodal PDF associated with the flamelet regime of premixed turbulent combustion [89] to a quasi-Gaussian PDF associated with extreme turbulence.
Despite the wide use of the flamelet concept jointly with a presumed PDF, such an approach definitely requires better validation and further development. On the one hand, recent simulations aimed at a posteriori assessment of the approach have showed its limited capabilities for predicting mean concentrations of intermediate species [59,62,65,70,71,75,76,80,81]. On the other hand, these results are not sufficient to draw the negative conclusion regarding all constituents of the discussed approach. For instance, predictive capabilities of Eq. (1), Eqs. (2)e(4), and a presumed PDF can be significantly different, as will be discussed later. Moreover, substantial disagreement between simulated and measured (or DNS) data, observed in the papers cited above, could stem not only from eventual limitations of the flamelet concept and/or the invoked PDF, but also from limitations of other models invoked in a posteriori study. For instance, as reviewed elsewhere [90,91], capabilities of available models for predicting thermal expansion effects in turbulent flames are still very limited and such limitations could contribute to the discussed disagreement.
Therefore, there is still need for a priory study that will allow researchers to assess predictive capabilities of Eq. (1), Eqs. (2)e(4), and/or presumed PDFs not only jointly, but also separately. However, a priory quantitative assessment of Eqs.
(1)e(4) and a presumed PDF have yet been very limited: the presumed b-PDF was quantitatively compared with DNS data by Donini et al. [72], whereas Lapointe and Blanquart [79] explored Eq. (1) for a single rate W c , i.e. the source term in the transport equation for the mean combustion progress variable.
The present paper reports results of a study originally conceived with a view to fill the discussed knowledge gap by analyzing recent DNS data [92] in order to perform a priory quantitative assessment of Eq. (1), Eqs. (2)e(4), and the presumed b-PDF for various species. It is worth stressing that the work is not limited to exploring these three submodels altogether but aims also at testing each submodel separately. The point is that the flamelet concept yielding Eqs. (1)e(4) and the presumed PDF approach are based on different reasoning. Therefore, Eqs. (1)e(4) could perform better than the presumed b-PDF or vice versa. Moreover, despite their apparent similarity, Eqs. (1) and (2) aim at solving basically different problems, i.e., prediction of the mean rate of product creation W c and evaluation of the mean mole fractions of various species. Accordingly, hypotheses and models developed to solve the former problem, which was also attacked in many studies that did not invoke Eq. (1), may differ significantly from hypotheses and models developed to solve the latter problem. Therefore, predictive capabilities of Eqs. (1)e(4) may be significantly different, see Section A priori assessment of conventional presumed PDF approach and its three constituents.
As will be discussed later, results of the present study did quantitatively validate Eqs. (2)e(4) but supported neither Eq.
(1) nor the presumed b-PDF. The latter finding called for advancing the presumed PDF approach by developing a better model of PðcÞ. Therefore, the focus of the study was extended, and another major goal was pursued. Thus, the present work aims not only at assessing Eqs. (1)e(4) and the presumed b-PDF, but also at developing and quantitatively validating a new presumed PDF model. Moreover, the present study compares different choices of the combustion progress variables.
The structure of the paper is as follows. In the next section, the DNS attributes [92] are briefly reported. Results of a priory quantitative assessment of Eq. (1), Eqs. (2)e(4), and the presumed b-PDF are discussed in the third section. A new flamelet-based presumed PDF approach is developed in the fourth section, followed by conclusions.

DNS attributes
The DNS data [92] were computed adapting the Pencil code [93], with the simulation attributes being similar to the attributes of earlier simulations performed by the same group and discussed in detail elsewhere [94]. Accordingly, we will restrict ourselves to a brief summary of these attributes.
Both studies [92,94] dealt with statistically 1D, planar premixed flames that propagated along the x-axis in a parallelepiped (19:18 Â 4:8 Â 4:8 mm) discretized using a uniform mesh of 960 Â 240 Â 240 nodes. The periodic boundary conditions were set on the transverse sides and Navier-Stokes Characteristic Boundary Conditions [95] were set at the inlet and outlet.
At t ¼ 0, a planar laminar flame (H 2 -air mixture with the equivalence ratio F ¼ 0:81 at 310 K, the laminar flame speed S L ¼ 1:84 m/s and thickness d L ¼ ðT b ÀT u Þ =maxjVTj ¼ 0:36 mm) was embedded into the computational domain at x ¼ x 0 . Subsequently, the 3D continuity, Navier-Stokes, species and energy transport equations were numerically solved using a detailed reaction mechanism (9 species, 21 reactions) by Li et al. [96].
Homogeneous isotropic turbulence was pre-generated using large-scale forcing in a cube with the periodic boundary conditions and was evolved until a statistically stationary state was reached. During combustion simulations, the turbulence entered the computational domain through the inlet boundary and decayed in the direction of the mean flow. The sole difference between the earlier [94] and recent [92] DNSs consisted of the inlet turbulence characteristics. In the present case [92], the rms velocity u dissipation of turbulent kinetic energy, averaged over the cube; time scale; subscript u or b designates unburned or burned mixture, respectively; and the summation convention applies to repeated indexes. Four different combustion progress variables were defined as follows: The mean profiles qðc k Þ of various quantities q were evaluated as follows. First, the qðx; tÞ and c k ðx; tÞ-fields were averaged over each transverse plane x ¼ const at each instant t. Second, obtained dependencies of qðx; tÞ on x were transformed to dependencies of < qjx ¼ c k ðtÞ > on a sample variable x using the averaged profiles of c k ðx; tÞ divided in 101 intervals, i.e. a transverse plane x ¼ const contributed to the value of < qjx j ¼ c k ðtÞ > provided that c k ðx; tÞ Àx j < 0:005 withx j ¼ 0:01j and j ¼ 0; …; 100: Third, the profiles of < qjx j ¼ c k ðtÞ > were averaged over various instants (54 snapshots stored each 5 ms, 1:401 ms t 1:566 ms), thus, yielding qðc k Þ reported in the following.
To examine Eqs. (1)e(4), PDFs P k ðx; x; tÞ were sampled from grid points characterized by c k ðx; tÞ Àx j < 0:005 for each J is a sample variable for the instantaneous c k ðx; tÞ-fields and J ¼ 100 if the opposite is not stated. Subsequently, the instantaneous PDFs were transformed to P k ðx; x ; tÞ, as discussed earlier, and were averaged over all instants. Finally, the time-averaged PDFs P k ðx; c k Þ were substituted into Eqs.
The presumed b-function PDFs P b;k ðx; c k Þ were modeled as follows [51,53].
by extracting the profiles of c 2 k ðc k Þ from the DNS data, as dis- z aÀ1 e Àz dz is the gamma function, and g k ¼ ðc 2 k Àc 2 k Þ =½c k ð1 Àc k Þ is the segregation factor. Subsequently, the PDFs P b;k ðx; c k Þ were applied to average W L ðc k Þ; X L ðc k Þ; T L ðc k Þ, and r L ðc k Þ using Eqs. (1)e(4).

Results and discussion
A priori assessment of conventional presumed PDF approach and its three constituents  (4), respectively, for the fuel-based c F ≡c 1 , product-based c 3 , temperaturebased c 4 and both for the actual PDF and for the b-function PDFs. The use of the oxygen-based c 2 jointly with the b-function PDF P b;2 ðx; c 2 Þ results in significantly underestimating the mean density, see dotted line in Fig. 2b. This underestimation is associated with non-monotonous variation of the mole fraction of O 2 in the laminar flame, as will be discussed later.
The left column in Fig. 3 quantitatively validates Eq. (2) for the major reactants, H 2 and O 2 , and product H 2 O for all four c k and for both PDFs. However, the mean mole fractions of radicals, see the middle and right columns in Fig. 3, are well predicted only when (i) the fuel-based combustion progress variable is adopted, see Fig. 3b, and (ii) the actual PDF P 1 ðx; c 1 Þ is used, cf. dashed and solid lines. The use of other c k yields worse results for radicals, especially for H. For HO 2 , OH, and O, the worst results are obtained for the oxygen-based c 2 , see the second row. Recent analysis [97] of DNS data obtained from highly turbulent lean H 2 -air flames also indicated that defining c based on H 2 (O 2 ) was the best (worst, respectively) choice.
Thus, Figs. 1e3 indicate that, while the analyzed case is characterized by Ka > 10 and is associated with strong preferential diffusion effects [26], the flamelet library is useful for predicting mean species concentrations provided that the library is properly adopted. The present results, as well as recent experimental and DNS studies reviewed elsewhere [83,84], imply that the domain of the flamelet concept validity is larger than earlier expected.
On the contrary, Fig. 4 indicates that Eq. (1) performs poorly for the mean rate W c for any of the four combustion progress variables. It is of interest to note that, for the temperaturebased c 4 , the use of P b;4 ðx; c 4 Þ yields better results, cf. dotted and solid lines in Fig. 4d, than the actual PDF does, see dashed line. This example illustrates that, if Eqs. (1) and (5) are jointly tested, a wrong conclusion regarding their validity could be drawn due to mutual cancellation of two types of errors (limitations of the flamelet library for the reaction rates and a wrong PDF). Therefore, to assess Eqs. (1)e(5) in a solid manner, each of the equation should also be tested separately.
Limitations of the flamelet library for the reaction rates are also indicated in Fig. 5, which reports evolution of turbulent burning velocities evaluated as follows for the species-based combustion progress variables c 1 (H 2 ), c 2 (O 2 ), and c 3 (H 2 O), or for the temperature-based combustion progress variable c 4 .
Here, M k is the molecular weight of species k. For each c k , the actual U T;k ðtÞ, see solid lines, differs substantially from U T;k ðtÞ yielded by Eq. (7) or (8) and Eq. (1) with the actual or b-function PDF, see dashed and dotted lines, respectively. Thus, while both Eq. (1) and Eqs. (2)e(4) stem from the same flamelet concept, Figs. 1e5 clearly show that the former equation performs significantly worse under conditions of the present study. The point is that variations in the mole fractions X, temperature T, or density r in a flame are substantially smoother than variations in the rates W that depend on the temperature in a highly non-linear manner. For instance, Fig. 6 reports dependencies of mole fractions X L ðc F Þ of various species and rates W L ðc F Þ of their production/consumption on the fuel-based combustion progress variable c F . These results have been obtained from the unperturbed stationary, planar, 1D laminar flame and have been re-normalized using the peak values of X L ðc F Þ or W L ðc F Þ, respectively. Comparison of Fig. 6a and b indicates significantly smoother variations in X n;L ðc F Þ when compared to W n;L ðc F Þ for each species n and a similar trend appears to hold in the studied turbulent flame. Accordingly, eventual errors associated with the flamelet concept, i.e. reduction of Xðx; tÞ and Wðx; tÞ to X L ½cðx; tÞ and W L ½cðx; tÞ, respectively, and eventual errors in modeling the PDF PðcÞ could result in significantly larger errors when averaging the non-linear rates W when compared to averaging the mole fractions. X:

Implications for modeling
Since Figs. 1e5 indicate that predictive capabilities of Eqs.
(2)e(4) are significantly better when compared to Eq. (1), the mean mole fractions of various species, the mean temperature and the mean density can be evaluated by adapting Eqs.
(2)e(4) independently of Eq. (1), e.g., by invoking a single-stepchemistry model, which yields a closure relation for the mean rate W c and performs better than Eq. (1). In such a case, X; T; and r can be calculated at a post-processing stage after computations of the mean flame speed, thickness, and structure.
Discussion of such a closure relation for W c is beyond the scope of the present work and the interested reader is referred to Refs. [88,89,98e102]. Here, we restrict ourselves to noting that predictive capabilities of simple single-step-chemistry models should not be underestimated. Indeed, first, as hypothesized by Prudnikov [103], reviewed elsewhere [88,98], and supported by more recent experimental data [104,105], various premixed turbulent flames have the same mean structure, i.e., variations of the Reynolds-averaged c along the normal to a mean flame surface are well approximated by the complementary error function if the distance is normalized with the mean flame brush thickness. Single-step-chemistry models can predict this fundamental feature of premixed turbulent combustion [88,98,106e108], whereas any effect of combustion chemistry on the mean flame structure has not yet been revealed.  Second, as hypothesized by Karlovitz et al. [109], reviewed elsewhere [88,98,103] and supported by more recent data [105,110], the growth of mean turbulent flame brush thickness follows the turbulent diffusion law in various experiments. Single-step-chemistry models can predict this fundamental feature of premixed turbulent combustion [88,98,108], whereas any effect of combustion chemistry on the growth of mean turbulent flame brush thickness has not yet been documented.
Third, while some influence of combustion chemistry on turbulent burning velocity is expected due to local quenching of combustion in extreme turbulence, the flamelet concept does not seem to hold under such conditions. As far as less, but still intense turbulence is concerned, single-stepchemistry models can predict turbulent burning velocity and flame speed reasonably well, as reviewed elsewhere [88,98,108]. It also worth stressing that, in recent DNS studies [111,112], a weak influence of combustion chemistry on U T has been documented at Karlovitz numbers significantly higher than unity.
To apply Eq. (2) at the post-processing stage, not only a closure relation for the mean rate W c , but also a PDF PðcÞ are required and modeling the PDF still challenges the combustion community. The issue could be addressed by developing approaches that deal with a transport equation for the PDF [47,48]. This research direction appears to be prioritized from the fundamental perspective and the present results provide additional motivation for such investigations. Alternatively, from the application perspective, the presumed PDF approach could be improved, as the use of the b-function PDF has not allowed us to predict the mean concentrations of radicals, see Fig. 3. This request is responded in the rest of the present paper by taking the following opportunity, which has yet been overlooked. The opportunity stems from the proposed combination of the flamelet Eqs. (2)e (4) with an independent closure relation for the mean rate W c .
Indeed, if Eq. (1) is not used to close W c , but another model of W c is invoked, the presumed PDF could involve one more unknown parameter, because the following three equations cðx; tÞ ¼ which (or their mass-weighted counterparts) are typically adopted to evaluate parameters of a common presumed PDF Pðc; x; tÞ, could be supplemented with one more constraint of where the mean rate W c ðx; tÞ is considered to be a known quantity yielded by an advanced model of the influence of turbulence on premixed combustion, as discussed earlier. In the rest of the present paper, W c is extracted from the DNS data. The use of the extra constraint given by Eq. (12), which is similar to Eq. (1), but aims at modeling the PDF Pðc;x;tÞ, rather than the mean rate W c ðx;tÞ, offers an opportunity to consider a wider set of presumed PDFs, which may involve four parameters. Alternatively, Eq. (11) could be substituted with Eq. (12). In such a case, the number of the PDF parameters is not increased, but a transport equation for the second moment c 2 ðx; tÞ could be skipped, thus, allowing us to circumvent modeling [102] of the scalar dissipation rate in that transport equation.
Moreover, since consumption of a fuel is mainly localized to reaction zones, substitution of Eq. (11) with Eq. (12) offers an opportunity to build a PDF that better predicts the probability of finding that zone and, hence, better predicts mean concentrations of various species. The mean rate W c ðx; tÞ is directly linked with the aforementioned probability, whereas such a link appears to be doubtful for the variance c 002 ðx; tÞ [87,88]. The latter claim is obvious in the Bray-Moss-Libby (BML) limit [113], where the variance is solely controlled by the probabilities of finding unburned reactants and equilibrium combustion products.

New flamelet-based presumed PDFs
The encouraging results obtained in the present work by testing the flamelet Eqs. (2)e(4) suggest that the presumed PDF approach could be substantially advanced using the flamelet PDF P L ðcÞf1=ðd L jVcj L Þ [53,54], but also invoking Eq. (12), as proposed earlier [87,88,114]. Development and assessment of such a new presumed flamelet-based PDF is the subject of the following discussion. Henceforth, we restrict ourselves to results computed by adopting the fuel-based combustion progress variable c F , because the mean concentrations of the radicals were best predicted using c F , as discussed earlier. Accordingly, subscript F will often be skipped for brevity in the following. Results obtained for the mean density, temperature, and mole fractions of H 2 , O 2 , and H 2 O will not be reported, because, for these quantities, all presumed PDFs discussed in the following yield comparably good agreement with the raw DNS data.
The choice of the flamelet PDF as a starting point is supported in Fig. 7. Indeed, in a large interval of x 1 < c < x 2 , appropriately re-normalized flamelet PDFs P L ðcÞ ¼ b =ðd L jVcj L Þ are close to the actual PDFs PðcÞ extracted directly from the DNS data at various cðx; tÞ, cf. dotted and solid lines, respectively. The re-normalization means that Eq. (9) does not hold, but the factor b is tuned in order to minimize the difference between the flamelet and actual PDFs in an interval of c characterized by P L ðcÞ < 2minfP L ðcÞg. Accordingly, the tuned value of b depends on c. The reasons for the re-normalization will be discussed later. It is worth remembering that the actual PDF has been obtained as follows. First, instantaneous PDFs P½c; cðx; tÞ have independently been extracted from each snapshot, with the PDF variations along the normal to the mean flame brush being characterized by cðx; tÞ. Second, those PDFs have been averaged over time to yield the actual PDF Pðc; cÞ whose spatial variations are characterized with c. Fig. 7 shows that the boundaries x 1 and x 2 of an interval where the two PDFs match well depend on c. In the largest part of the mean flame brush, x 1 is small and x 2 z0:9, see Fig. 7b. Since the fuel consumption in the laminar flame occurs at c 0:9, see solid line in Fig. 6b, the flamelet PDF P L ðcÞ plotted in Fig. 7b appears to fit the actual PDF sufficiently well, as far as averaging W L;c ðcÞ is concerned. Fig. 6a implies that the flamelet PDF P L ðcÞ may also be applied to evaluate mean mole fractions of HO 2 and H 2 O 2 , because these mole fractions are small at c > 0:9 in the laminar flame. However, P L ðcÞ seems to be less suitable for averaging the mole fractions of HO 2 and H 2 O 2 at the trailing edge of the mean flame brush, because x 1 is increased by c, see Fig. 7c. Moreover, the left boundary of the validity domain for the flamelet PDF, i.e. x 1 , may be increased by Re t [115], as will be discussed later. Furthermore, (i) the mole fractions of the radicals H, O, and, especially, OH peak at large c in the laminar flame, see Fig. 6a, and (ii) differences between the actual and flamelet PDFs are well pronounced at c > 0:9, especially if the mean c is low, cf. Fig. 7a and 7c. Accordingly, the flamelet PDF appears to perform worse for the radicals H, O, and OH. Note that the differences between the two PDFs could in part stem from insufficient statistics sampled at cðx; tÞ/1, because the probability of finding cðx; tÞ ¼ 1 is low under the DNS conditions. Fig. 7 elucidates one more issue associated with the flamelet PDF. It grows rapidly at c/0 or c/1 and is much larger than the actual PDF in these edge zones. Accordingly, such zones can significantly contribute to the integrals on the RHSs of Eqs. 9e12. Consequently, at least one of these constraints is expected to be violated. For instance, as will be discussed in detail later, if substitution of P L ðcÞ into Eq. (12) allowed us to satisfy this constraint, Eq. (9) would be violated, i.e., its RHS would be too large. Indeed, the RHS of Eq. (9) is controlled by a large P L ðc /1Þ, whereas the RHS of Eq. (12) is controlled by much smaller P L ð0:1 < c < 0:9Þ, because the rate W c ðcÞ vanishes at large c.
From the physical perspective, the use of the flamelet PDF P L ðcÞ for small and large c appears to be flawed also. As emphasized by Kuznetsov and Sabelnikov [115], at high Reynolds numbers, the local structure of premixed flames at small c is similar to the structure of an inert scalar field, i.e. the gradient Vc and, hence, the PDF PðcÞ are controlled by turbulent straining if c is sufficiently small and Re t is sufficiently high. Indeed, significant perturbations in the local structure and thickness of preheat zones characterized by small c are well documented in highly turbulent flames [83,84]. Similar perturbations are expected to be of substantial importance in the radical-recombination zone [116,117], whose thickness is relatively large, because reactions are relatively slow therein. On the contrary, the local reaction zone appears to be less affected by turbulent eddies, e.g., the zone thicknesses in laminar and turbulent flows are comparable, as reviewed elsewhere [83,84]. Accordingly, adaptation of the flamelet PDF P L ðcÞ to a bounded interval of 0 < x 1 < c < x 2 < 1 seems to be justified better than the use of P L ðcÞ in the entire interval of 0 c 1. For these reasons, the PDFs P L ðcÞ plotted in Fig. 7 are renormalized and Eq. (9) does not hold for them (if intervals of c < x 1 and c > x 2 are not considered, then, integration of the flamelet PDF from zero to unity is not required). We may also note that the flamelet PDFs shown in Fig. 7 vary weakly in wide intervals of c, with the exception of the edge zones. This is consistent with recent experimental and DNS data [83,84] that indicate that reaction zones are thin even in highly turbulent flows. Indeed, if a zone is thin, variations in the actual PDF PðcÞ are expected to be weak in the zone, see solid lines for 0:2 < c < 0:8 in Fig. 7, similar to weak variations in P L ðcÞ in the same interval of c.
Based on the above results and reasoning, the flamelet PDF should solely be applied to a bounded interval of 0 < x 1 < c < x 2 < 1, with the PDF being re-normalized to satisfy Eq. (12). To satisfy other constraints, e.g. Eqs. 9 and 10, and, apparently, Eq. (11) if the PDF invokes four parameters, the behavior of PðcÞ at c < x 1 and c > x 2 should also be modeled in a turbulent flame. The simplest presumed PDF that could satisfy all these requirements appears to be as follows PðcÞ ¼ adðc À x min Þ þ bdðx max À cÞ þ gHðc À x 1 ÞHðx 2 À cÞP L ðcÞ; (13) where dðxÞ and HðxÞ are Dirac delta function and Heaviside function, respectively, 0 x min x 1 , x 2 x max 1, and P L ðcÞ is normalized so that A similar presumed PDF was earlier put forward by Bushe et al. [54,61,66,73], but these authors assumed that x 1 ¼ x min , x 2 ¼ x max and did not adapt Eq. (12) to calibrate the PDF parameters. On the contrary, the use of Eq. (12) is the most important peculiarity of the present approach. Within its framework, the parameters a, b, and g in Eq. (13) and other presumed flamelet-based PDFs discussed in the following are evaluated by adopting Eqs. (9), (10) and (12).
For the PDF presumed by invoking Eq. (13), these three constraints read Results reported in the following were obtained by setting x min ¼ 0 and x 2 > 0:9. Accordingly, the first term on the Left Hand Side (LHS) of Eq. (16) and the first two terms on the LHS of Eq. (17) vanish or are negligibly small. Fig. 8 shows that, if c < 0:6 (or c > 0:6), Eqs. (13e17) with x min ¼ 0 and x max ¼ 1, see curves 3, predict the mean mole fractions of H and O, see curves 2 1, better (worse, respectively) than the conventional b-function PDF does, see curves 2. For OH, the two presumed PDFs yield close results if c < 0:4, but the conventional b-function PDF performs better if c > 0:4. Thus, substitution of the latter PDF with the considered flameletbased PDF does not substantially improve predictions of the mean mole fractions of the radicals H, O, and OH. This result is associated with the fact that Eq. (13) does not allow us to model processes localized to the radical-recombination zones characterized by large c, whereas the radical concentrations change drastically in such zones [116,117].
Comparison of curves 3 and 4 in Fig. 8 also emphasizes importance of appropriately modeling contributions from the radical-recombination zones to X H , X O , and X OH . The sole difference between the models applied to compute the two curves consists of substitution of x max ¼ 1, see curves 3, with x max ¼ ðx 2 þ1Þ=2, see curves 4. In the latter case, agreement with the mean mole fractions of the radicals, evaluated using the actual PDF, see curves 1, is substantially better. Moreover, Eqs. (13e17) with x min ¼ 0 and x max ¼ ðx 2 þ1Þ=2 perform significantly better than the conventional b-function PDF, see curves 2, for H, O, and OH in the largest part of the mean flame brush with the exception of its trailing edge (c > 0:95). The improved predictions yielded by Eqs. (13e17) with x min ¼ 0 and x max ¼ ðx 2 þ1Þ =2 imply that a layer around the surface of cðx; tÞ ¼ ðx 2 þ1Þ =2 is sufficiently representative for averaging 2 Curves 1 in Fig. 8 look more scattered than the curves shown in dashed lines in the middle column in Fig. 3, because the PDFs adopted to compute results reported in Fig. 8 have been calculated by dividing the interval of 0 x 1 into a larger number (J ¼ 250) of bins. ðx 2 þ1Þ =2, to X H , X O , and X OH is partially counterbalanced by overestimation of the contribution from the trailing half of the zone, i.e., ðx 2 þ1Þ =2 < c < 1. Here, the terms "underestimation" and "overestimation" are used, because the actual PDF is larger (smaller) than the flamelet-based PDF in the leading (trailing, respectively) halves of the zone, see Fig. 7.
To allow for processes localized to the radicalrecombination zones, the following simple flamelet-based presumed PDF PðcÞ ¼ adðcÞ þ bHðc À x 2 ÞHð1 À cÞ þ gHðc À x 1 ÞHðx 2 À cÞP L ðcÞ; could also be used. When compared to Eq. (13) with x min ¼ 0, Eq. (18) adopts a constant PDF in the interval of x 2 < c < 1 instead of a single Dirac delta function. Eqs. (9), (10), (12), (14) and (18) read Here, b 0 ≡bð1 Àx 2 Þ and contribution from the second term on the RHS of Eq. (18) to the mean rate W c is neglected, because W c;L ðc > x 2 Þ is small for the values of x 2 used in the present study, see solid curve in Fig. 6b. Curves 5 in Fig. 8a and b shows that Eqs. (18e21) predict X H and X O substantially better than the three other PDFs discussed earlier do. For X OH , Eqs. (13e17) with x min ¼ 0 and x max ¼ ðx 2 þ1Þ =2 and Eqs. (18e21) perform similarly, but much better than the two other PDFs discussed earlier do.
A common drawback of all the flamelet-based PDFs considered above consists of the lack of a solid criterion for selecting x 2 . This drawback not only offers a room for tuning, but could also reduce predictive capabilities of the approach, because selection of x 2 is not a trivial task. To demonstrate  importance of the appropriate choice of x 2 , Fig. 9 shows dependencies of the following two integrals on x 2 (note that Fig. 9b reports results I W;F ðx 1 ; x 2 Þ of integration of the negative rate of production of the fuel). These results have been calculated by setting x 1 so that W c;L ðc ¼ x 1 Þ =maxfW c;L ðcÞg ¼ 0:01, but the integrals I c ðx 1 ; x 2 Þ and, especially, I W ðx 1 ; x 2 Þ depend weakly on x 1 if it is sufficiently small. A rapid increase in I c ðx 1 ; x 2 Þ when x 2 /1, see Fig. 9a, results from a sharp increase in P L ðc /1Þ. A decrease in the magnitude I W;F ðx 1 ; x 2 Þ at large x 2 is explained as follows. If x 2 is increased by retaining the same P L ðcÞ, an integral of is significantly increased due to integration of a large P L ðcÞ over the added interval of c. However, this integral should be equal to unity due to the normalization constraint given by Eq. (14). Consequently, the PDF P L ðcÞ ¼ b=ðd L jVcj L Þ should be renormalized, i.e. b should be decreased with increasing x 2 . The decrease in b results in decreasing P L ðcÞ and, hence, W c;L ðcÞP L ðcÞ in the reaction zone, where the rate W c;L ðcÞ is high. On the contrary, the extension of the integration domain to larger c has a minor direct effect on I W ðx 1 ;x 2 Þ, because the rate W c;L ðcÞ is too low in the added domain. Accordingly, the indirect effect, i.e. the decrease in b, required to satisfy Eq. (14) when increasing x 2 , dominates and the integral I W ðx 1 ; x 2 Þ is decreased at large x 2 . Fig. 9 indicates that the use of Eqs. (13) and (14) requires an appropriate pre-selection of x 2 . For instance, if c is too small, substitution of g given by Eq. (17) or (21) into Eq. (16) or (20), respectively, could yield gI c ðx 1 ; x 2 Þ > c for certain x 2 , thus, requiring negative a or b. Moreover, since g 1 due to Eq. (15) or (19), Eq. (17) or (21) with x 2 considered to be independent of c can be satisfied only if jI W ðx 1 ; x 2 Þj ! maxfW c ðcÞg and this constraint has been adopted when calculating results plotted in curves 3e5 in Fig. 8. Otherwise, jI W ðx 1 ; x 2 Þj can be too small if x 2 is too large, see Fig. 9b. In such a case, Eqs. (15) and (17) or (19) and (21) require g 1 and g > 1, respectively. However, for c associated with W c ðcÞ < maxfW c ðcÞg, lower values of I W ðx 1 ; x 2 Þ can be consistent with Eqs. (12) and (17), or (21). Accordingly, a larger x 2 may be set to extend an interval that the flamelet PDF is applied to. This extension is beneficial, because (i) the flamelet and actual PDFs match in a wide range of c in Fig. 7 and (ii) the former PDF is based on physical reasoning, whereas the second term on the RHS of Eq. (13) or (18) is sufficiently arbitrary. An increase in x 2 results in increasing (decreasing) an interval of c that the third (second, respectively) term on the RHS of Eq. (13) or (18) is applied to. Accordingly, the third and physically sound (second and arbitrary) term plays a more (less, respectively) important role when x 2 is increased. Thus, x 2 can and should vary with c.
Accordingly, the following unified algorithm is proposed to build the flamelet-based presumed PDF by adopting Eqs. (13e17) in the simplest case of x min ¼ 0 and x max ¼ 1. First, the integrals I c ðx 1 ; x 2 Þ and I W ðx 1 ; x 2 Þ are pre-calculated for various x 2 by processing data obtained in simulations of the unperturbed complex-chemistry laminar flame. Second, for each c, the product of gI W ðx 1 ; x 2 Þ is compared with W c ðcÞ and the largest value of x 2;i consistent with Eq. (17) is selected. If this value is also consistent with Eqs. (15) and (16), where 0 a 1, 0 b 1, and 0 g 1, the mean mole fractions are evaluated using Eq. (2). Otherwise, the next lower value x 2;iÀ1 is used and the procedure is repeated until all three constraints, i.e. Eqs. (15e17), are satisfied for certain "optimal" x 2;iÀj , where 0 j < i. While there is no guarantee that all these constraints can simultaneously be satisfied, this is so for c ! 0:06 in the present study. For lower c, the algorithm does not seem to be applicable, because the PDF Pðc; cÞ tends to Dirac delta function dðcÞ as c/0, the rate W c ðc /0Þ/0, and the mean concentrations of radicals are very low at small c. Accordingly, at low c, Fig. 10 e Dependence of the optimal x 2 on the mean value c F of the fuel-based combustion progress variable. i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 5 ( 2 0 2 0 ) 3 1 1 6 2 e3 1 1 7 8 one can assume that b ¼ 0 and find two other PDF parameters adopting Eqs. (15) and (16). Then, the algorithm is straightforward and free from tuning (the choice of x 1 will be discussed later) provided that the PDF shape is presumed, e.g., Eq. (13) with x min ¼ 0 and x max ¼ 1 is invoked.
Results computed by adopting Eqs. (13e17) and the suggested algorithm are shown in curves 6 in Fig. 8. Agreement with the mean mole fractions computed by substituting the actual DNS into Eq. (2) is encouraging for all three radicals. Some scatter of the former results stems from (i) the use of a discrete set of allowed x 2 and (ii) a step growth of the optimal x 2;iÀj with c, see Fig. 10. Note that x 1 has been selected based on the constraint of W c;L ðc ¼ x 1 Þ=maxfW c;L ðcÞg ¼ ε, with almost the same results being obtained for ε ¼ 0:1 and 0.01 under conditions of the present study. It is worth remembering, however, that the use of a larger ε and, hence, a larger x 1 may be necessary at higher turbulent Reynolds numbers [115].
While, under conditions of the present study, the value of ε≪1 weakly affects X H ðcÞ, X O ðcÞ, and X OH ðcÞ calculated by adapting various extended flamelet-based presumed PDFs discussed above, this value should be sufficiently small in order for these PDFs to predict X HO2 ðcÞ and X H2O2 ðcÞ, cf. curves 3 and 4 in Fig. 11. As suggested in Fig. 6a and confirmed by the present analysis of the DNS data, these results are weakly sensitive to the behavior of the PDFs at large c. Accordingly, solely profiles calculated using the simplest Eq. (13) with x min ¼ 0 and x max ¼ 1 are reported in Fig. 11. Comparison of curves 1 and 4 shows that the extended flamelet-based presumed PDF well predicts the mean mole fractions of HO 2 and H 2 O 2 provided that x 1 is sufficiently small. This PDF performs better than the conventional b-function PDF does, see curves 2.
Finally, predictive capabilities of the developed flameletbased presumed PDF approach are quantitatively validated in Fig. 12. The dashed lines show results computed by adopting Eqs. (18e21) with x min ¼ 0, x 2 depending on c, and x 1 resulting from the constraint of W c;L ðc ¼ x 1 Þ =maxfW c;L ðcÞg ¼ 0:01. For all intermediate species with the exception of H 2 O 2 , the approach is clearly superior to the use of the conventional b-function PDF. There is a room for further improvement of the model, e.g., by smoothing variations in x 2 with c in order to smooth curves shown in dashed lines in Fig. 12a. However, such polishing of the developed approach does not seem to be necessary at this stage of research, because the achieved  . agreement between the model and DNS results is already impressive, cf. dashed and solid lines in Fig. 12. Therefore, if the developed approach is adopted in a CFD study, precision of numerical simulations of mean concentrations of various species in a premixed turbulent flame appears to be controlled by other invoked models, e.g., by models of turbulence in premixed flames [90,91].
In other words, Fig. 12 implies that, from the purely application perspective, modeling the PDF PðcÞ is no longer a bottleneck, at least under conditions of the present study. Assessment of the developed approach at higher Reynolds numbers, higher Karlovitz numbers, for other fuels and species, or/and within the LES framework appear to be more important and interesting tasks for future research. Nevertheless, it is worth remembering that the presumed PDF given by Eq. (13) or (18) invokes a very simple approximation of PðcÞ at large c. While such simplifications appear to be justified for applications, at least at the present stage of CFD research into complex-chemistry premixed turbulent combustion, these simplifications may be disputed from the fundamental perspective. Therefore, the behavior of PðcÞ at large c definitely requires further study, e.g., by advancing the transport equation for the PDF. In particular, comparison of results obtained by adapting different PDFs, see Fig. 8, motivates targetdirected investigation of PðcÞ and its transport equation in the radical-recombination zone of complex-chemistry flames. This task is of crucial importance for predicting the mean concentrations of H, O, and OH.

Conclusions
A quantitative a priori assessment of the flamelet approach to evaluating the mean density r, the mean temperature T, the mean mole fractions X n of various species, and the mean source term W c in the transport equation for the mean combustion progress variable c has been performed by analysing complex-chemistry DNS data obtained recently by Dave and Chaudhuri [92] from a lean hydrogen-air flame characterized by u 0 =S L ¼ 3:6 and the Karlovitz number Ka ¼ 13. Four different choices of the combustion progress variable have been probed, with the PDF PðcÞ being either (i) extracted directly from the DNS data or (ii) presumed by invoking the widely-used b-function and adopting the first two moments of the cðx;tÞ-field extracted directly from the DNS data.
The results show that rðcÞ; TðcÞ, and X n ðcÞ for all species (major reactants H 2 and O 2 , product H 2 O, and radicals H, O, OH, HO 2 and H 2 O 2 ) are well predicted by the flamelet Eqs.
(2)e(4) provided that (i) the combustion progress variable is defined using the fuel mass fraction and (ii) the actual PDF extracted directly from the DNS data is substituted into these equations. In line with other recent experimental and DNS data reviewed elsewhere [83,84], this result indicates that the domain of validity of the flamelet concept is substantially wider than earlier expected. The use of the b-function PDF yields substantially worse results for the radicals. The mean rate W c and turbulent burning velocity based on it are poorly predicted for all probed c even if the actual PDF is substituted into Eq. (1).
These findings imply that, in order to evaluate the mean temperature, density and species concentrations, Eqs. (2)e(4) could be coupled with another closure relation for the mean rate W c , whose predictive capabilities are better documented when compared to Eq. (1). In such a case, Eqs. (2)e(4) could be implemented as post-processing of a mean c-field computed by numerically integrating a single transport equation for the mean combustion progress variable.
Therefore, modeling the PDF PðcÞ in Eqs.
(2)e(4) is put on the forefront of the research agenda. Accordingly, an extended flamelet-based presumed PDF approach has been developed by using the mean rate W c for the PDF calibration. If the rate is extracted from the DNS data, the newly developed approach well predicts mean concentrations of all species (H 2 , O 2 , H 2 O, H, O, OH, HO 2 , H 2 O 2 ) without using any tuning parameter.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.