A New Analytic pdf for Simulations of Premixed Turbulent Combustion

A new reaction rate source term ωm(c)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega _m(c)$$\end{document} for modelling of premixed combustion with a single progress variable c is proposed. ωm(c)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega _m(c)$$\end{document} mimics closely the Arrhenius source term ωA(c)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega _A(c)$$\end{document} for a large range of activation energies and density ratios while admitting analytic evaluation of many quantities of interest. The analytic flame profile cm(ξ)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_m(\xi )$$\end{document} very closely approximates the numerically integrated Arrhenius flame profiles cA(ξ)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_A(\xi )$$\end{document}. An important feature of cm(ξ)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_m(\xi )$$\end{document} is that it is analytically invertible into a ξm(c)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\xi _m(c)$$\end{document}. Analytic estimates of the laminar flame Eigenvalue Λ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Lambda$$\end{document} and of the Le dependence of the laminar flame speed sL\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_L$$\end{document} are derived, which are more accurate than classic results based on asymptotic analyses. The flamelet pdf p(c)=1/(Δ∗c∗(1-cm))\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p(c)=1/(\Delta *c*(1-c^m))$$\end{document} for a flat laminar flame front in a LES cell of width Δ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta$$\end{document} is derived. The exact mean of the reaction rate ω(c)¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{\omega (c)}$$\end{document} is compared to a beta pdf result, which is shown to be inaccurate for large ratios of filter width to flame thickness Δ/δf\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta /\delta _f$$\end{document} and particularly for high activation energy flames. For multidimensional flame wrinkling we derive the exact relationship p(c)=p1D(c)I(c)Ξ(c)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p(c)=p_{1D}(c)I(c)\Xi (c)$$\end{document} between the 3D pdf p(c), the 1D flat flame pdf p1D(c)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_{1D}(c)$$\end{document}, a correction factor I(c) for change of inner flame structure and a geometrical wrinkling factor Ξ(c)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Xi (c)$$\end{document}. We show that the c dependence of these quantities cannot be neglected for small Δ/δf\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta /\delta _f$$\end{document}. A simple model of a sinusoidally wrinkled flame front qualitatively demonstrates the effect of flame wrinkling on p(c). We show that for large Δ/δf\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta /\delta _f$$\end{document}, a wrinkling of the reaction zone almost constantly increases p(c) in the reaction zone by a wrinkling factor Ξ∗\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Xi ^*$$\end{document} (defined for the surface of the isosurface of maximum heat release) while reducing it near c=0,1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c=0,1$$\end{document} as required for normalisation of p(c). The 1D p(c) evaluated using a reduced filter width Δ′=Δ/Ξ∗\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta '=\Delta /\Xi ^*$$\end{document} may be a good approximation of the wrinkled flame pdf for evaluation of ω(c)¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{\omega (c)}$$\end{document} for such cases.

1D flat flame f Flame

Introduction
In the CFD simulation of turbulent premixed combustion using RANS and LES, the size of computational cells is usually too big to fully resolve the laminar flame structure embedded in the turbulent flow field. The ratio of flame thickness to cell size decreases further at elevated pressures due to the drop of diffusivities and heat conductivity with pressure. Thus combustion (subgrid) models are required for the simulation of premixed combustion processes in most technical applications. For better understanding of basic phenomena and for model development, it is often useful in premixed combustion to study the case of a single reaction progress variable c = (T − T u )∕(T b − T u ) assuming single-step irreversible chemistry and adiabatic combustion. T u , T b are the unburnt and fully burnt temperatures, respectively. The c transport equation is given by where , u, c are density, velocity and progress variable, and , c p are the heat conductivity and specific heat at constant pressure. For Arrhenius chemistry and Lewis number Le = 1 , the chemical source term can be written as Poinsot and Veynante (2005): represents the normalized temperature raise and = T a1 ∕T b is a measure of the activation temperature T a1 . The temperature exponent 1 in Eq.
(2) is usually taken as 1 = 0 or 1 = 1 . The continuity equation requires u = const. = u s L at steady state and we have ∼ 1∕T ∼ 1∕(1 − (1 − c)) for constant pressure combustion. One can rescale the spatial coordinate (Poinsot and Veynante 2005) according to = ∫ x 0 u s L c p ∕ dx , yielding a simpler differential equation for c: The prefactor Λ in Eq. (4) represents the eigenvalue of the transport equation which needs to be chosen such that the boundary conditions c = 0 for → −∞ and c = 1 for → +∞ are fulfilled. Asymptotic analysis of Eqs. (3,4) in the limit of large by Zeldovic, Frank Kamenetski and von Karman yielded Λ ∼ 2 ∕2 . An improved result was derived by Williams (2018) for 1 = 1: The turbulent flow field will wrinkle and stretch the flame front, modifying the volume specific fuel consumption rate. For high activation energies, the reaction front will remain thin even for quite high Karlowitz numbers (Nilsson et al. 2019). A variety of subgrid models have been developed to represent the effect of subgrid flame folding/stretching on the fuel consumption rate. The artificially thickened flame model (Colin et al. 2000) makes the flame front resolvable on the LES grid by increasing the diffusion and heat conduction coefficients and reducing the strength of the reaction term, leaving the flame speed unchanged. The effect of non resolved subgrid flame wrinkling needs to be taken into account by an efficiency function. Other approaches are flame surface density (transport) models, where the RHS of Eq. (1) is replaced by u ⟨s c ⟩Σ f with flame surface density Σ f and ⟨s c ⟩ being a surface averaged flame speed. Σ f is either determined by a transport equation or approximated as Σ f = Ξ | ∇c | and evoking algebraic models for the wrinkling factor Ξ (Poinsot and Veynante 2005;Ma et al. 2013). Other researchers have proposed to pretabulate the filtered source term (Fiorina et al. 2010) for given LES cell size Δ or extract more accurate expressions for (c) by 1-D approximate deconvolution (Domingo and Vervisch 2015).
Probability density functions p(c) with normalization condition ∫ 1 0 p(c)dc = 1 allow to evaluate cell averages of any quantity z(c) as z(c) = ∫ 1 0 z(c)p(c)dc . In classic BML theory the pdf is assumed to take the form with (c) ≪ 1 . For (c) → 0 one obtains B ∼ c and A ∼ (1 − c) . Unfortunately, in the limit (c) = 0 accurate mean values of quantities vanishing at both ends c = 0, 1 (like the chemical source term (c) ) cannot be evaluated by p BML (c) . Consequently, the calculation of the mean fuel consumption rate requires knowledge of (c) . A 1-D flamelet pdf can be defined through (Bray et al. 2006): where c − , c + represent the c values on the edges of the computational cell. All ingredients of p f (c) can be evaluated only numerically for Arrhenius chemistry. Bray et al. (2006), Salehi and Bushe (2010) use ad hoc lower and upper bounds c − = > 0 , c + = 1 − to determine N, claiming that N depends only weakly on . Another popular presumed pdf is the beta pdf where the coefficients a, b > 0 can be chosen to reproduce given values of mean and variance of c. In case of a beta pdf, c − = 0 and c + = 1 . The beta pdf is mostly used to model 1 3 diffusion processes and does a good job in representing typical pdf shapes occurring in this case. As a representation of the flamelet pdf, it has been shown (Bray et al. 2006) to deliver inaccurate results for the mean reaction term particularly for large c variance. Knudsen et al. (2010) find a good correlation between the DNS source term filtered to a LES grid and the beta pdf value (with mean and variance in the beta pdf evaluated from the DNS) for small ratios of Δ LES ∕Δ DNS < 5 , while large errors result for larger ratios Δ LES ∕Δ DNS . Domingo et al. (2005) report that their flamelet pdf with ad hoc choice of the integration limits c − , c + delivered unphysical parameters A, B for small values of c variance. They switched to a beta pdf in those cases. In the current contribution we hope to clarify some of these issues.
The paper is structured as follows: We first present the simplified reaction source term introduced by Ferziger and Echekki (1993), which yields analytical results for many quantities of interest. We then introduce a new source term m (c) with similar analytic capability, which however approximates the Arrhenius source term and the flamelet profile much closer. As first applications, we derive expressions for the Λ eigenvalue of the Arrhenius source term and for the dependence of the laminar flame speed s L on Lewis number and Arrhenius parameters, which are more accurate than classic expressions from the literature. We compare flamelet and beta pdf's for different filter widths Δ∕ f . We discuss the relation between pdf, wrinkling factor and correction factor for multidimensionally wrinkled flames and we demonstrate the modification to the pdf by flame wrinkling using a simple model of a sinusoidally wrinkled flame. Finally, we give some conclusions and an outlook to further work. Ferziger and Echekki (1993) proposed to set the source term to zero in the pure diffusion region 0 ≤ c < 1 − 1∕ and to use a linear term (c) = Λ E (1 − c) in the reaction region 1 − 1∕ ≤ c ≤ 1 . Due to the simple form of the source term, the solution of the differential equation

Ferziger/Echekki Source Term
where H(x) represents the Heaviside jump function, are pure exponentials in . Choosing the location = 0 for c = 1 − 1∕ puts the pure diffusion region at negative and the reacting region at positive . Matching function values and derivatives of c( ) at = 0 , the c( ) profile becomes The eigenvalue of Eq. (9) is given by Λ E = ( − 1) . Figure 1 shows a comparison of the Echekki c E ( ) profile and its source term with the Arrhenius source term and the numerically integrated c A ( ) profile for 1 = 0 and a relatively low = 6 , = 9∕11 corresponding to a density jump T b ∕T u = 4.5 . The c( ) profiles are in reasonable agreement.

New Analytical Flamelet Profile c m () and Source Term ! m (c)
Although the Echekki solution admits analytical evaluations of a number of quantities, which are available only numerically with the Arrhenius source term, it has some disadvantages: c E ( ) is only a reasonable approximation to c A ( ) for small to medium -the pdf is not a good approximation to the Arrhenius pdf in the reactive c region -the piecewise definitions of c E ( ) and its inversion E (c) makes calculations quite cumbersome It is therefore our goal to find an improved approximation to the Arrhenius A (c) , which guarantees the correct boundary conditions of c( ) and still admits analytical integration of Eq. (3). This can be achieved by postulating functional forms of c( ) with analytical inverse (c) and calculating (c) according to: and (c) = z( (c)) . Obviously, the resulting (c) should be positive over the whole interval 0 < c < 1 . The Echekki c E ( ) fulfills this requirement but has the above disadvantages. A more elegant and admissible c( ) is given by: Similar to the Echekki c E ( ) with parameter , c m ( ) has a free parameter m which can be used to mimic c A ( ) for different Arrhenius parameters , , 1 . For large m, c m ( ) approaches a step function H( ) centered at = 0 , which is the correct thin flame limit. The profile c m ( ) can be inverted easily: The profile yields a thermal flame thickness of:  and the reaction source term is given by: Note that this source term is identical to that of the KPP-Fisher equation for m = 1 , generalizations exist which correspond to the case m ≠ 1 . Equation (11) represents the travelling wave equation of these equations with suitably rescaled spatial coordinate . Solutions of this equation and their behaviour have been studied in Kyrychko andBlyuss (2009), Fan (2002).

Determination of Parameter m
To determine values of m which "optimally" represent a given Arrhenius source term A (c) with parameters , , 1 , we generated numerical solutions of Eq. (3) for 1 = 0, 1 for a range of parameters ( , ) and determined m by numerical minimization of Note that the numerical shift 0 , which is necessary to align numerically generated c A ( ) with c m ( ) , is irrelevant since the position of the flame front on the axis is arbitrary. The fitted m are found to be well represented by the following expressions valid in the physically relevant ranges of 4 < < 30 and 1∕2 < < 1: Figure 2 shows comparisons of A (c) with m (c) and of c A ( ) with c m (c) for 1 = 0, = 9∕11, = 6, 18 using these m correlations. Results for 1 = 1 are too similar to be shown here. The new source term is always shifted slightly towards c = 1 compared to the Arrhenius one. Note however the excellent agreement between c A ( ) and c m ( ).  Table 1 shows typical m values for small and large activation energy flames as analysed in Poinsot and Veynante (2005).

Application: Derivation of Improved 3 Eigenvalue
As first applications of this model we present a new approach to derive improved approximations for the laminar flame eigenvalues Λ 0,1 . In contrast to results from classic asymptotic analyses (which are derived in the limit of large activation energy, i.e. large values), the new expressions are also accurate for moderate values of .
We use the fact that by construction m already guarantees the correct limiting behaviour c( → −∞) = 0 and c( → +∞) = 1 . The optimal m values are found to be closely correlated to c max , the location of the maximum of A (c) , which is independent of prefactor Λ . Therefore a choice of the Arrhenius Λ which minimizes the difference between A (c) and m (c) for given m should provide a good approximation to the correct Arrhenius flame eigenvalue Λ . Normally one would look for a minimimum of the integral of the squared difference between A (c) and m (c) . Since this does not admit an analytic solution we solve the following equation for Λ: This equation is acceptable in this case due to the undulary difference of A (c) and m (c) . For 1 = 0 we get with Ei(z) represents the exponential integral function. To simplify the expressions in Eqs. (19,20), we note that exp[ ∕( − 1)] → 0 for large and < 1 and that the function f (z) = 1∕(e z Ei(−z)(z + 2)z + z + 1) with z = ∕ can be approximated accurately by a second order polynomial in z in the interesting range 4 < z < 40 . The following expressions for Λ 0 , Λ 1 approximate Eqs. (19, 20) very closely in this z range: For comparison to classical results, we perform a Taylor series in 1∕ , yielding the following leading order terms: and Note the similarity of the first two terms in Eq. (24) to the Williams result Eq. (5), which was derived by asymptotic expansions in 1∕ and assuming 1 = 1 . The Williams factor (3 − 1.344) is replaced by (3 − 6∕5) with additional terms in only. For 1 = 0 , the term multiplying is replaced by 2 − 5∕4 with different terms. While the Williams Λ works well for > 8 and 1 = 1 , the new expressions are also valid for smaller > 4 and provide a much better approximation in the case of 1 = 0 . Figure 3 shows the numerically integrated c A ( ) � s using the Arrhenius source term A (c) for = 9∕11, = 6, 1 = 0 with Eigenvalues Λ 0 and Λ Williams . The violation of the boundary condition c( ) = 0 for c → −∞ when using Λ Williams is apparent. 1 3

Application: Lewis Number Dependence of Laminar Flame Speed
In the case of Le ≠ 1 , an additional transport equation for the normalized fuel concentration Y(x, t) needs to be solved: (25) can be cast in the following form: emphasising the similarity to Eq. (1) except for the denominator 1/Le. In the chemical source term (c) the factor (1 − c) in front of the exponential in Eq. (4) needs to be replaced by (1 − ) (Poinsot and Veynante 2005;Ferziger and Echekki 1993): Again rescaling the spatial coordinate according to = u s L c p ∕ , we obtain: and where now the chemical source term depends not only on c but also on : The prefactor Λ Le = ∕(c p 2 u s 2 L ) in the source term represents the Lewis-dependent eigenvalue of the transport equations, which needs to chosen such that the boundary conditions c, = 0 for → −∞ and c, = 1 for → +∞ are fulfilled simultaneously. For Le = 1 , it is obvious that ( ) ≡ c( ) solves both Eqs. (28, 29) with the given boundary conditions and Λ chosen as described above.
From numerical integration of Eqs. (28, 29), varying Λ Le to fulfill the boundary conditions we find (maybe not surprisingly) that c( ) � s for Le ≠ 1 are almost identical to c( ) for Le=1, while the ( ) are more diffuse for Le < 1 and less diffuse for Le > 1 . We find also, that (1 − ) is closely proportional to (1 − c) in the region where the source term (c, ) is notably different from zero. Inspection of Eq. (30) shows that if I Le (1 − ) = (1 − c) with constant I Le , the source term in Eq. (28) would be identical to Eq. (4) and the solution for c( ) would indeed be the same as in the case Le = 1 provided that Λ Le ∕I Le = Λ Le=1 ( , ) . Equation (29) can be cast in a formally similar form to Eq. (28) by rescaling the spatial coordinate in Eq. (29) as = Le : which is formally identical to Eq. (28) if we replace by c, by � = I Le , and by � = I Le . The Eigenvalue of the equation to fulfill the boundary condition must therefore be equal to Λ Le=1 ( � , � ) . This yields the following equations for Λ Le and I Le : A similar slightly more complex relation can be obtained for s L in the case of 1 = 1 . The dependence of s L (Le)∕s L (Le = 1) on Le is much stronger than on , . Figure 4 shows s L (Le)∕s L (Le = 1) versus Le for = 0.818 , typical values of and 1 = 0, 1 . We can see that the differences between s L for 1 = 0 and 1 = 1 in this range of Le numbers are minimal (of the order of 1% at Le = 2 , but increasing at higher Le). A Taylor series of s 2 L in inverse powers shows that to leading order s L (Le)∕s L (Le = 1) →

Analytic Flamelet pdf's
In the following, we derive the subgrid flamelet pdf of a freely propagating flat premixed flame within a CFD cell of width Δ . Figure  The flamelet pdf p(c) is given by with normalisation factor N. The normalisation condition for p(c) yields The new profile c m ( ) yields: and The functional form of the new flamelet pdf p m (c) may look similar to a beta pdf p (c) but there are notable differences. The first factor 1/c could be realized in p (c) by setting a = 0 , which is however not admissible in this context since the pdf would not be integrable at c = 0 . Note that the factor 1/c arises as a consequence of the pure diffusion in the preheat region c → 0 and is independent of the Arrhenius parameters , , 1 . The second factor 1∕(1 − c m ) is responsible for the characteristic asymmetry of the flamelet pdf, its integral also diverges logarithmically c = 1 . In contrast to the beta pdf, the power m is on c, not on (1 − c).

Mean Value of Chemical Source Term
For the new c( ) flame profile and pdf, the mean of the reaction source term m (c) evaluates analytically as: An analytical result can also be derived for the mean of the sum of laminar diffusion and reaction source terms, which is often modelled together e.g. in flame surface density models: This result again is true for all flamelet pdf's. For small Δ , this expression obviously approaches the correct limit dc∕d . In flame surface density models the mean of the sum of diffusion and source terms is modelled as u s L | ∇c | , which in our context is equivalent to dc∕d because u = 1 = s L due to the scaling of .

Evaluation of Means and Variances of c
For evaluation of mean and variance of c from p m (c) the following integrals are defined: Specifically, we have where 2 F 1 (a, b, c;z) is the hypergeometric function, which can be evaluated numerically from its power series definition for z < 1 . Note that its last argument c m ∈ [0, 1] . The last equality in Eq. (48) results from addition theorems of 2 F 1 (a, b, c;z) . For integer m, I 1 (c, m) reduces to expressions containing only powers of c, logarithms and trigonometric functions, e.g.: The means of c and c 2 in a cell of width Δ evaluate as: and from which the variance c 2 m , Δ −(c m , Δ ) 2 can be calculated.

Numerical Evaluation of c
The numerical evaluation of I n (c, m) becomes numerically ill-defined for c values approaching 1. The reason is that the series expansion of the hypergeometric function 2 F 1 (a, b, c;z) does not converge well for z → 1 . c can be evaluated in c or space: where + = − + Δ , c − = c m ( − ), c + = c m ( + ) and − = m (c − ) , + = m (c + ) . The indefinite integral of c m ( ) is given by: We now choose a small such that 1 − c m ( ) < . By replacing c m ( ) in the integrand of Eq. (53) by 1 for > , we can assure that the hypergeometric function is only evaluated for values of its last argument far enough from 1 so that the series still converges. In the results presented here, we have used = 10 −6 . Explicitly, = log (1 − (1 − ) −m )∕m.

Evaluation of Favre Averages
At constant pressure, the ratio of density to the unburnt density u is proportional to 1/T, yielding and for the Favre average of c Those two quantities cannot be evaluated analytically for arbitrary values of m, but integer (and some half-integer) values of m yield analytic results in terms of powers of c, logarithmic and trigonimetric functions. Results are too unwieldy to be shown here. Figure 6 shows c as function of c for various values of Δ together with the BML limit c = c∕(1 − c) for m = 4, = 3∕4 . We can see that the BML limit is approached only for Δ ≫ f ( f = 1.87 for m = 4).

Analytic Results for ˇ-pdf
The mean of the new chemical source term m (c) can also be analytically evaluated if the beta pdf p (c) is used instead of the flamelet pdf p m (c): Since an explicit expression is desirable for fast evaluation in the CFD code and for the study of the behaviour of c − (c, Δ) , we provide the following approximation. With This approximation has an error below 1% and represents c − , c + considerably better at large Δ∕ f . Obviously, since c is available analytically as function of c − and Δ through Eq. (50), the relation betweem c and c − can also be inverted numerically e.g. by 1-D-re-interpolation. Since c is a strictly monotonic function of c − and Δ , this procedure is always well defined.

Comparison of Different pdf's
A comparison of p A (c) (evaluated numerically from c A ( ) ) with the Echekki p E (c) and the new pdf p m (c) is shown in Fig. 8 together with the respective source terms A (c) , E (c) and m (c) . It is evident that p E (c) , featuring a slope discontinuity at c = 1∕ is not a particularly good approximation to the real p A (c) in contrast to the new p m (c).  Figure 9 compares the new pdf p m (c) and the beta pdf p (c) for two values of Δ . The exact mean and variance of the c distribution within the cell was used to calculate the parameters a, b in p (c) . Similar to Bray et al. (2006), we find that p (c) does not fit the flamelet pdf well in the region where (c) > 0 , causing an over prediction of (c) for large c variance, i.e. large Δ∕ f . Note that in contrast to p (c) , the shape of the flamelet pdf p m (c) is actually independent of c and Δ ; those parameters only influence the limits c − and c + and the normalisation factor N = Δ.
Both pdf's approach the correct limit of p(c) → (c − c) in the DNS limit ( Δ << f ). For small ratios Δ∕ f , only a small part of p m (c) is cut out near c while for larger Δ∕ f , c − and c + move nearer towards c = 0, 1 and reveal a larger portion of the complete 1∕(dc m ( )∕d ) shape. In contrast, the beta pdf resembles a Gaussian near c for small variance (i.e. small Δ∕ f ) and switches to a double-delta type behaviour for larger Δ∕ f . Figure 10 (left) shows a comparison of (c) evaluated using the beta pdf (with exact mean and variance of c to evaluate a, b) with the exact (c) (evaluated using p m (c) ) vs. Δ for two different values of m. On the right we plot the ratio (c)∕ (c) . (c) over predicts the exact value of (c) for Δ > 0.75 . In general (c)∕ (c) increases with Δ and saturates for Δ → ∞ to a value of m∕3 + 0.4 . The error introduced when using the beta pdf thus increases with increasing activation energy (i.e. , i.e. m).

Effect of 3D Flame Wrinkling
In the thin reaction zone regime, large turbulent eddies of size l e ≫ f will wrinkle the reaction zone (Fig. 11) while leaving its internal structure largely intact. This will increase the mean reaction rate in the cell approximately by the wrinkling factor Ξ = Σ∕Δ 2 , where Σ is the area of the wrinkled reaction zone within the LES cell and Δ 2 is the area of the flat flame propagating in direction.
For a given c( ) field within an LES cell volume Ω of size Δ 3 , the area of an isosurface Σ of a certain (c * ) value is defined by (Osher and Fedkiw 2005): Note that (Boger et al. 1998) We can define a generalized wrinkling factor Ξ(c * ) through: Note that Ξ(c * ) can be smaller than 1 for c * isosurfaces which are not or only partially contained within the LES cell. The fine grained pdf of a given c( ) field is given by Poinsot and Veynante (2005), Gao and OBrien (1993)

Evaluation of !(c) in the Thin Reaction Zone Regime
Using the definition of the fine-grained pdf, we can evaluate (c) as We can multiply the delta function by 1 = |∇c| |∇c| : In the thin reaction zone limit, c isosurfaces within the reaction zone are largely parallel and | ∇c( ) | c * ≈| dc∕d | 1D,c * . Replacing the c gradient in the denominator by the 1D gradient and moving it outside the Ω integral (note that | dc∕d | 1D,c * is constant during the Ω integration since it depends on c * only) yields: The second integral is just Σ(c * ) , see Eq. (63) and the first term is equal to the flat flame pdf: We obtain the following approximation in the thin reaction zone limit: using Eq. (65). For large LES cells ( Δ > f ), Σ(c * ) and Ξ(c * ) might be approximately constant in the c region where (c) > 0 as we show later using a simple flame wrinkling model. Equation (70) can then be further simplified as proposed in Boger et al. (1998): The assumption of a constant wrinkling factor (independent of c) is the basis of many algebraic flame surface density models.
Obviously, Ξ * in Eq. (71) should represent the wrinkling of the c * isosurface representing the maximum heat release. Note that the wrinkling of c isosurfaces outside the reaction region ( (c) > 0 ), e.g. within the preheat zone, will not affect (c).

Relation Between pdf, Wrinkling Factor and Correction Factor
We can also multiply the delta function in Eq. (67)  This decomposition nicely separates the contributions of the freely propagating flat flame pdf p 1D (c * ) , correction factor I(c * ) (providing the modification of the flame inner structure due to e.g. strain) and wrinkling factor Ξ(c * ) (contributing the geometrical effect of 3D wrinkling of c * isosurfaces). This exact equation is valid independent of the number of spatial dimensions and of the form and size of the LES cells. Such relationships make pdf formulations attractive and conceptually more general than e.g. 1-D flamelet filtered tabulations or 1-D approximate deconvolution methods. The general applicability also makes pdf methods suitable candidates for LES simulations of high pressure premixed flames, where a large amount of subgrid flame wrinkling is to be expected.

Illustration of Wrinkled Flame pdf
In the appendix, we present a simple model for the 3-D c( ) field of a sinusoidally wrinkled flame, which allows analytical calculations of isocontours. Implict is the assumption that the c gradient on all c * isosurfaces is equal to that of the flat flame at c = c * . A wrinkling of the reaction zone c * isosurface can push other isosurfaces partly or totally out of the LES cell as illustrated in Fig. 12. This reduces the corresponding Σ(c) 's and translates into a reduction of p(c) as shown in Fig. 13. Parameters used to produce these plots are given in the "Appendix". We see that for small enough wrinkling amplitude A∕Δ , p(c) is enhanced almost uniformly by the wrinkling factor Ξ * in the reactive c region, where (c) > 0 , while p(c) is decreased below p 1D (c) already for c > c − and c < c + . Also shown in Fig. 13 are 1D pdf's where we replaced Δ by a smaller filter size Δ∕Ξ * (which automatically moves cut-off values c − , c + closer towards c ). These surrogate pdf's appear to be a good approximation of the model pdf in the reactive c region for not too large wrinkling amplitudes and large Δ∕ f .
Shapes of premixed flame pdf's derived from DNS data (Salehi and Bushe 2010;Salehi et al. 2013;Jin et al. 2008) and from application of the (stochastic) linear eddy model to numerically derived 1D flamelet pdf's (Tsui and Bushe 2014) look similar to the pdf's derived from this simple analytical model. It appears that the latter can reproduce the effects of turbulent flame wrinkling on p(c) at least qualitatively.

Effect of Flame Stretch and Curvature on p(c)
In a locally stretched flame, the reaction layer will become thinner, increasing | ∇c | over (dc∕d ) 1D and thus decreasing I(c) and p(c). In contrast, a locally thickened reaction zone (e.g through turbulent mixing by very small vortices or through differential diffusion) should lead to an increase of I(c) and p(c). Such an increase of I(c) with increasing Karlowitz number has been observed in highly turbulent flame DNS's (Luca et al. 2019).
Stronger modifications of local p(c)'s are expected in large LES cells containing regions of large flame curvature (e.g. cusps, tips of flame fingers) and regions of subgrid flameflame interaction. Here the filtered flame front might not be represented well by the 1-D flat flame structure.

Summary and Conclusions
In this paper we present a chemical source term m (c) to be used in the transport equation of a single reaction progress variable c, which closely mimics the Arrhenius A (c) but yields an analytical solution c m ( ) with a simple analytical inverse m (c) . The parameter m can be adapted so that c m ( ) closely mimics the Arrhenius c A ( ) for the relevant range of its parameters , , 1 .
We derive the flat flame premixed flamelet pdf p m (c) = 1 Δ 1 c(1−c m ) and provide an analytic correlation to evaluate the limits of the c integration c − , c + as function of c and Δ for efficient implementation of p m (c) into CFD codes. The integral of p m (c) diverges logarithmically at c = 0, 1 , but these limits are never accessed. Due to the closed form of the expressions no piecewise definition of the pdf as in Salehi and Bushe (2010), Domingo et al. (2005) or functions at c = 0, 1 like in BML theory are required. The explicit form of the flamelet p(c)'s show that the 1/c-behaviour near c = 0 is universal, independent of Arrhenius parameters, since it is caused by the pure heat diffusion within the preheat zone.
Many quantities of interest can be evaluated analytically with this approach. We derive estimates of the laminar flame Eigenvalue Λ and of the Le dependence of the laminar flame speed s L which are more accurate than classical results from the literature based on asymptotic analysis.
The cell averages of source term m (c) and of the sum of source and diffusion terms (equivalent to dc∕d ) are evaluated analytically with the new flamelet pdf and compared to beta pdf results, which are also available in closed form when using m (c) as source term. We show that for Δ∕ f < 1 (i.e. small c variance) the beta pdf results are reasonable, while for large Δ∕ f the beta pdf strongly over predicts (c) , particularly for flames of large activation energy. Due to the simplicity of the new pdf, we can expect that p m (c) could replace the beta pdf in presumed pdf modelling of the progress variable in premixed combustion.
The effect of flame wrinkling on the pdf is investigated using the fine-grained pdf formulation. We derive the exact relationship p(c) = p 1D (c) * I(c) * Ξ(c) , separating the effects of the 1-D flamelet pdf from those of flame thickening and geometrical flame wrinkling. The c dependence of I(c) and Ξ(c) needs to be taken into account for large wrinkling factors and if isosurfaces for different values of c are wrinkled differently, e.g. in the case of large Karlowitz number flames.
For illustration of the effect of flame wrinkling on the 3D pdf we derive a simple model of sinusoidally wrinkled thin reaction fronts where for a given c * value | ∇c | is identical the 1 3 flat flame gradient. We show that for large Δ∕ f and moderate flame wrinkling the main effect on the pdf is an increase of p(c) in the reactive region by a wrinkling factor Ξ * , representing the wrinkling of the isosurface at maximum heat release. p(c) decreases near c = 0, 1 as required by its normalisation condition. In this limit the wrinkled flame pdf's resemble flat flame pdf's evaluated at a smaller filter Δ � = Δ∕Ξ * in the reaction region c ≈ c * .
As the wrinkling factor Ξ(c) is mostly determined by the turbulent flow field, it may be simpler to model Ξ than other quantities like the c variance or the c scalar dissipation rate which depend on the turbulent field and the inner flame structure. Many models for (constant) Ξ * are available in the framework of flame surface density theory (for an overview, see e.g. Ma et al. 2013). Ultimately, the choice of modelling style may however be a matter of taste or experience.
In the future we plan to compare the derived pdf's with data extracted from DNS's of premixed flames at ambient and elevated pressures (e.g. Klein et al. 2018). We will investigate modelling strategies for Ξ(c) and I(c) and we will attempt to extend the theory to the case of stratified / partially premixed flames.
1 3 repeated parabola sections. Since the results are lengthy and the model is only approximate we omit the details here. Figure 14 shows the wrinkling of a c * front for (n = 1) together with isolines at different distances d. The cell size is chosen as Δ = 4 and located between −2 < < 2 and −2 < < 2 as shown by the dashed frame. The unwrinkled c * contour would be located at ≡ 0 . It is evident that curves at larger distances d are folding back unto themselves (plotted here as thin lines). They need to be cut at the edges of the cell (here = ±2 ) to guarantee a unique c field. Isosurfaces sufficiently near to the c * isosurface (plotted in full gray) will not fold onto themselves and will not have to be cut in this way. Isosurfaces at larger distances d will be also cut at the top/bottom boundaries (here = ±2 ) of the cell, see dashed curves in Fig. 12.
Due to the assumption of | ∇c |=| dc∕d | 1D,c * on all c * isosurfaces, the correction factor I(c) = 1 and the ratio of the wrinkled flame pdf to the flat flame one becomes p(c)∕p 1D (c) = Ξ(c) = Σ(c)∕Δ 2 . In the present model, Σ(c) is equal to the length of the sinusoidal isoline contained within the cell times the width of the cell Δ . Therefore p(c)∕p 1D (c) is equal to the length of the wrinkled isolines within the 2D cell divided by Δ . Note that this is a purely geometrical quantity dependent only on A∕Δ and n in the present simple model. The c value corresponding to an isoline at distance d can be calculated as c(d) = c m ( = d).
For illustration of the behaviour of the generalized wrinkling factor Ξ(c) , we place a low activation energy flame (m = 4.45) in an LES cell of width Δ = 4 as shown in Fig. 15. This arrangement yields c = 0.7, c − = 0.135, c + = 0.99997 . The generalized wrinkling factor Ξ(c) for wrinkling with n = 1, 2 cosine waves inside the cell at amplitudes A = 1 and A = 2 is shown in Fig. 16. Also shown are results of the same flame but in a cell of double size Δ = 8. Figure 17 shows a comparison of the resulting pdf's. We see that for Δ = 4 and the smaller wrinkling amplitude A = 1 , p(c) is increased by a nearly constant factor (equal to the ratio of the length of the sine function to a flat isoline) in the reaction region (c > 0) . The pdf is reduced towards c → c − and c → c + (compared to the flat flame pdf) as required by p(c) normalisation. For larger wrinkling amplitude A = 2 , the maximally possible enhancement of p(c) due to isocontour wrinkling is not realized over the full c range between c − and c + since some isocontours within the reaction front are cut by the upper/lower cell boundaries. For A = 2 in a larger cell with Δ = 8 , a smaller portion of  Fig. 17. The over prediction of (c) by the beta pdf reduces with larger flame wrinkling, provided that the  exact mean and variance of c are used. This can be understood through the observation that the 1D pdf resembles one at a reduced cell size Δ � = Δ∕Ξ * . Note that the over prediction of (c) by the beta pdf is less pronounced at small Δ∕ f .