The role of presumed probability density function in the simulation of non premixed turbulent combustion

Flamelet Progress Variable (FPV) combustion models allow the evaluation of all thermo chemical quantities in a reacting flow by computing only the mixture fraction Z and a progress variable C. When using such a method to predict a turbulent combustion in conjunction with a turbulence model, a probability density function (PDF) is required to evaluate statistical averages (e.g., Favre average) of chemical quantities. The choice of the PDF is a compromise between computational costs and accuracy level. The aim of this paper is to investigate the influence of the PDF choice and its modeling aspects in the simulation of non premixed turbulent combustion. Three different models are considered: the standard one, based on the choice of a beta distribution for Z and a Dirac distribution for C; a model employing a beta distribution for both Z and C; a third model obtained using a beta distribution for Z and the statistical most likely distribution (SMLD) for C. The standard model, although widely used, doesn't take into account the interaction between turbulence and chemical kinetics as well as the dependence of the progress variable not only on its mean but also on its variance. The SMLD approach establishes a systematic framework to incorporate informations from an arbitrary number of moments, thus providing an improvement over conventionally employed presumed PDF closure models. The rational behind the choice of the three PDFs is described in some details and the prediction capability of the corresponding models is tested versus well known test cases, namely, the Sandia flames, and a test case for supersonic combustion provided by Cheng.


Introduction
The industrial and scientific communities are devoting major research efforts to identify and assess innovative technologies for advanced propulsive concepts. Among such technologies, hydro-carbon combustion has been assumed as a key issue to achieve better propulsive performance and lower environmental impact. In order to improve the know-how to build more efficient engines with lower emissions it is necessary to enhance the knowledge of the combustion phenomena. In this context the simulation of turbulent reacting flows is very useful to cut down experimental costs and to achieve a thorough comprehension of the physical mechanisms involved. Turbulent combustion is a multi-scale problem, where the interaction between chemical kinetics, molecular, and turbulent transport occurs over a wide range of length and time scales. The numerical simulation of such phenomena with detailed chemistry is today prohibitive, so that a reduction model is often employed to condensate the reaction mechanisms and cut down the computational costs. Therefore, different approaches have been proposed to address this problem, such as the reduction of the chemical scheme in intrinsic low dimensional manifolds (ILDM) [6]; the flamelet-based approaches such as the flameletprogress variable (FPV) [7] or flame prolongation of ILDM (FPI) [8]; and Flamelet Generated Manifolds approach (FGM) [9].
Our interest is devoted here to diffusive, either partially premixed or non-premixed, flames which constitute a specific class of combustion problems where fuel and oxidizer are not mixed before they enter into the combustion chamber. In this case mixing must bring reactants into the reaction zone so as to activate and maintain the combustion process. Non-premixed flames can be characterized by a local balance between diffusion and reaction [10]. Their structure can be described by a conserved scalar, the so-called mixture fraction. A diffusive flame can be viewed as an ensemble of thin locally one-dimensional structures embedded within the flow field. Each element of the flame front can then be described as a small laminar flame, also called flamelet. In this paper we focus on FPV approach for turbulent non-premixed flames. The FPV approach is based on the use of only two degrees of freedom, namely, the mixture fraction and the progress variable, that are employed to map all of the thermo-chemical quantities involved in the process. For the case of a turbulent flame one needs to define a probability density function (PDF) to compute the Favre average of the thermo-chemical quantities. In particular, in the present work we are interested in the modeling of the PDF function required to evaluate chemical Favre averaged quantities. The definition of such a function is critical due to the poor knowledge of the two independent variables behaviour. The aim of this work is to provide an extension of the standard FPV model for turbulent combustion, applying the statistically most likely distribution (SMLD) [11] approach to the progress variable PDF, maintaining a good compromise between computational costs and accuracy level. In the second section of this paper we present the rational for the definition of such a PDF. Then, three PDF models are considered and their role in the evaluation of non-premixed flames is analysed. In the third section the numerical results obtained in the simulation of the Sandia flames [4] and of a supersonic combustion [5] are discussed. The paper closes with summary and conclusions.

Combustion model 2.1 The flamelet approach
The FPV model proposed by Pierce [7] is used in this work to evaluate all of the thermo-chemical quantities involved in the combustion process. This approach is based on the parametrization of the generic thermochemical quantities, φ, in terms of two variables: the mixture fraction Z and the progress variable C: Equation (1) is taken as the solution of the steady laminar flamelet equation: where χ is the scalar dissipation rate modeled in terms of molecular diffusivity of mixture fraction D Z , χ = 2D Z (∇Z) 2 ; ρ is the density;ω φ is the source term related to φ [7]. Each solution of equation (2) is a flamelet and the solution variety over χ = χ st , called S-curve, is shown in figure 1. From equation (1) one can obtain the Favre-averages of φ using the definitions: where P (Z, C) is the density-weighted PDF, P (Z, C) is the PDF and ρ is the Reynolds-averaged density. As usual, φ is decomposed as: and, where φ and ρ are the fluctuations. This ensures that the filtering process does not alter the form of the conservation laws.
The choice of such a PDF plays a crucial role in the definition of the model, being a compromise between computational costs and accuracy level. In this respect, this paper provides an extension of the standard FPV turbulent combustion model combined with a RANS equation solver [12], where different fundamental hypotheses are used to define the PDF function for the progress variable, C. The influence of the different PDFs in the simulation of non-premixed turbulent combustion is the final aim of this research.

Presumed probability density function modeling
In order to investigate the role of the presumed PDF one can, first of all, use the Bayes theorem and take the PDF as the product between the marginal PDF of Z and the conditional PDF of C|Z: Therefore, one has to presume, or evaluate, the functional shape of such PDFs. Let us consider the marginal PDF, P (Z). It has been shown, by several authors, that the mixture fraction is best described by a passive scalar and that the PDF of a passive scalar can be approximated by a β distribution function [13,14,15]. The two parameter family of β-distribution in the interval x ∈ [0, 1] is given by: where Γ(x) is the Euler function and a and b are two parameters related to x and x 2 respectively For all the three models presented here, the β-distribution is employed for P (Z).
To presume the functional shape of the distribution of a reacting scalar, one needs to make some constitutive hypotheses. To simplify the problem, in this work we assume the statistical independence of Z and C, so that, for all of the considered models, P (Z, C) = P (Z) P (C), namely C = C|Z. The most widely used hypothesis (model A), implying great simplification in the theoretical framework, consists in assuming that the conditional PDF, P (C), can be modeled by a Dirac distribution. It can be shown that there is only one solution of equation (2) for each chemical state. With this criterion the Favre-average of a generic thermo-chemical quantity is given by: Therefore, one has only three additional transport equations (for Z, Z 2 and C) to evaluate all thermochemical quantities in the flow thus avoiding the expensive solution of a transport equation for each chemical species. However, it is well known that a reactive scalar [3], such as C, depends on a combination of solutions of equation (2) for each chemical state and therefore its PDF cannot be accurately approximated by a Dirac distribution. Therefore, the second model (model B) is given by assuming that Z and C are distributed in the same way, namely, using a β -distribution, thus giving the following joint PDF: This avoids the simplification seen before and, consequently, the model requires the evaluation of an additional transport equation for C 2 . Moreover, the probability distribution of a reacting scalar is often multi-modal, unlike the β function, and its functional form depends on the turbulence-chemistry interaction. Therefore, one can think about a distribution built considering, as constraints, the only available informations, namely the value of Z, Z 2 , C and C 2 . The third model (model C) is obtained evaluating the conditional PDF as the statistically most likely distribution (SMLD) [3]. It can be shown that if one knows only its first three moments, the PDF can be evaluated using "Laplace's principle of insufficient reason" [11]. The technique is developed following the statistical mechanics arguments presented by Heinz [2]. Relying on the knowledge of the first three moments of P (C), a unique measure, S, of the predictability of a thermodynamic state can be defined. S is an entropy function depending on P (C), S = S( P (C)) [16] that can be thought of as the Boltzmann's entropy: where Q(C) is a bias density function to integrate information when no moments are known. In this paper the form of Q(C) proposed by Pope [3] is assumed. The goal is to construct a PDF that maximizes the entropy S. Following the Lagrangian optimization approach, the functional S * is defined by involving the constraints on the moments: In the above equation µ n are the Lagrange's multipliers while the last fraction term is introduced to normalize P (C). The expression for P (C), obtained evaluating the maximum of S * , reads: where: The role of PDF in non-premixed flames At this point the model still needs an additional assumption to be closed. Here we assume that the first and the last point of P (C) are equal to the first and last points of β(C) evaluated with the given values of the mean and variance: This assumption does not affect the multi-modal nature of the distribution, but simplifies the model implementation (there is no need to evaluate the roots of a non-linear system). The major advantage of the SMLD approach over conventionally employed presumed PDF closure models is that it provides a systematic framework to incorporate an arbitrary number of moment information. It is noteworthy that, since C is used instead of C|Z as argument of P , also this model assumes statistical independence of Z and C.
For the case of a turbulent flame, equation (1) must be written in terms of the Favre averages of Z and C and in terms of their variance. Using the model A one can tabulate all chemical quantities in terms of Z, Z 2 and C because of the properties of the δ-distribution. On the other hand, models B and C express φ in terms of C 2 too and therefore they need to evolve a transport equation also for C 2 . Therefore, in order to evaluate a diffusive flame with the three models described above, one has to define a transport equation for each of the flamelet variables, namely, Z, Z 2 , C and C 2 . Only in the case of model A, employing the δ-distribution for P (C), the transport equation for C 2 is not needed and only three equations are solved. The transport equations read: where D is the diffusion coefficient for all of the species, given as D = ν/P r evaluated assuming a unity Lewis number; ν is the cinematic viscosity and P r the Prandtl number; D t Z = D t C 2 = D t C = D t C 2 = ν/Sc t are the turbulent mass diffusion coefficients and Sc t the Shmidt turbulent number;ω C is the source for the progress variable. The gradient transport assumption for turbulent fluxes is used and the mean scalar dissipation rate, χ t appear as a sink term in the equation (21). At every iteration, the values of the flamelet variables of the model are updated and the Favre-averaged thermo-chemical quantities are defined, using equation (3). Such solutions provide the mean-mass-fractions which are used to evaluate the density, the enthalpy and all of the transport properties of the fluid.

Flow equations and numerical solution
The numerical method developed [12] has been employed to solve the steady-state RANS equations with the k-ω turbulence closure. For an axisymmetric multi-component reacting compressible flow of n species the system of the governing equations can be written as: where t is the time variable; x and r are the axial and the radial coordinate, respectively; Q=(ρ, ρ u x , ρ u r , ρ H, ρk, ρω, ρ R n ) is the vector of the conserved variables; ρ, ( u x , u r ), H indicate the Favre-averaged values of density, velocity components and specific total enthalpy, respectively; k and ω are the turbulence kinetic energy and its specific dissipation rate; finally, R n is a generic set of conserved variables related to the combustion model; as described above, R n is the set of independent variables of the flamelet model, namely, Z, Z 2 , C, C 2 ; E, F , and E v , F v are the inviscid and viscous flux vectors [17], respectively; S is the vector of the source term. A cell-centered finite volume space discretization is used on a multi-block structured mesh. The convective and viscous terms are discretize by the third-order-accurate Steger and Warming [18] flux vector-splittingscheme and by second-order-accurate central differences, respectively. An implicit time marching procedure is used with a factorization based on the diagonalization procedure of Pulliamm and Chaussee [19], so as to allow a standard scalar alternating direction implicit (ADI) solution procedure [20]. Only steady flows are dealt with in this paper. Therefore, the unsteady terms are eliminated from the governing equations and the ADI scheme is iterated in the pseudo-time until a residual drop of five orders of magnitude for all of the conservation-law equations (24) has been achieved. Additional details, when needed, are provided in the following for each single application.

Numerical results
This section provides the comparison among the results obtained using the three combustion models so as to assess the influence that the PDF choice may have in the prediction of turbulent non-premixed flames. The well-knonw subsonic Sandia flames [4] are considered at first, then a supersonic test case is analysed, whose experimental data are available in the literature [5]. The steady flamelet evaluations, for both test cases, have been performed using the FlameMaster code [21].    The temperature field obtained by the three combustion models for the case of the Sandia Flame F is presented in figure 2. Moreover, figure 3 provides the temperature distributions along the axial direction. It appears that in the near-burner region model C is in better agreement with the experimental data than the other two models. Moving away from the burner (x > 20 d ref ) the agreement deteriorates; this is probably be due to the accuracy limits of the RANS approach in the prediction of the mixing process that greatly affects combustion. From this two set of figures one can see that there is an improvement provided by model C in the evaluation of the flame core and of the flame shapes. The radial distributions of temperature, mixture fraction and CO 2 mass fraction at the axial coordinate x/d ref equal to 1 and 15 are shown in figures 4 and 5 respectively. One can see that the results of model B (corresponding to presume that P (Z, C) = β(Z)β(C)) and model C (corresponding to the choice P (Z, C) = β(Z)P SM L,2 (C)) are in better agreement with the experimental data [4] than the results of model A (corresponding to the standard choice P (Z, C) = β(Z)δ(C − C|Z)). In particular, a sligth improvement in the evaluation of temperature is observed. In fact, considering figure 4, model A predicts a maximum temperature of about 2250 K that is not present in the experimental data and is much less evident in the results of the other two models, providing a spike of about 2050 K. The improvement in the prediction of the flame core is evident in the figure 3, where the differences in the results by the three models are more relevant. Moreover in figure 5 one can see the shift of the location of temperature peaks that is considerable in Flame E and Flame F because of the higher inlet velocity (with respect of flame D case). It is interesting to observe that for the case of Flame D the lower probability of local extinction is registered; consequently, the statistical independence of Z and C is not negligible [23] and none of the models is able to give results in good agreement with the experimental data even in the first section, see figure 4. Envisaging a model that does not require the hypothesis of statistically independence of Z and C should considerably enhance the prediction accuracy also in the regions close to the burner.

Supersonic combustion test case
The second test case presented is the case of hydrogen-air supersonic combustion proposed by Cheng et al. [5]. The supersonic burner provides an annular, axisymmetric jet of hot, vitiated wet air at Mach number equal to 2, average axial velocity of 1417 m/s, temperature of 1250 K and pressure of 107 kP a. The vitiated air is composed of the following set of mass fractions: Y O2 = 0.245, Y H2O = 0.175 and Y N2 = 0.58. The exit conditions of the air stream are given by a pre-combustion at low temperature [5]. The hydrogen exit is estimated as a chocked flow, average axial velocity of 1780 m/s, temperature of 545 K and pressure of 112 kP a. The diameter of the fuel stream is d ref = 2.362 mm, taken as the reference distance. At the inlet section, the inner diameter of the vitiated air stream is 3.812 mm and its outer diameter is 17.78 mm [24], respectively. The computational domain is axisymmetric and includes the divergent part of the air nozzle; it extends 150 d ref and 50 d ref along the axial and radial directions, respectively, and has been discretized using about 70000 cells. The evaluation of the flamelet library has been performed using the Jachimowski supersonic combustion scheme [25]: 34 sub-reactions upon 11 species. The flamelet library has been computed over a grid with 250, uniformly distributed points, in the Z and C directions, and 50 uniformly distributed points in the Z 2 and C 2 directions. A first consideration is in order. From the temperature contours provided in figure 6 one can see that none of the developed models is able to ignite the mixture. It is argued that this happens because of the absence of OH in the inlet air pilot composition. In fact, the OH fraction has significant impact on the ignition delay, so that, if OH is absent at the inlet ignition cannot occur [26]. Cheng et al. [5] observed OH traces in the inlet composition of the air stream, coming from the pre-combustion; but the OH mass fraction was negligible with respect to that of H 2 O, O 2 and N 2 . A quantitative analysis is provided in figure 7, showing that the standard model A can evaluate with a better agreement than the other two models the data in the first section. On the other hand, in the section at x = 10 d ref there is an evaluated temperature lower than that from the data, because of the non-burning mixture. Moreover, it is interesting to observe that model B and model C produce almost the same results, see figures 6, 7 and 8, in fact the dashed-dotted line is almost overlapped to the solid line. In figure 8 the distributions of Z and Z 2 for two different sections are presented. One can see that the evaluation of the mixture fraction variance is better in the last two models than in the standard one, especially in the first section. The problems encountered in the reproduction of the experiment of Cheng et al. comes from the use of the steady flamelet equation (2). Given the complexity of the flow the implicit assumption of a stoichiometric scalar dissipation, χ = χ st , is a nonoptimal choice. The very high level of local extinction probability gives in turn the statistical independence of Z and C [23], so the best way to evaluate such a flow is, hopefully, to use the unsteady-FPV in conjunction with the model C to obtain a good estimation of the progress variable distribution without the hypothesis of a stoichiometric scalar dissipation rate.

Conclusions
This paper provides an extension of standard FPV model combined with a RANS solver introducing the SMLD approach to describe the progress variable distribution. The work is composed of four sections. The first one is an introduction to the problem of the presumed PDFs in the non-premixed combustion. The second section describes the combustion model developed with a new closure method for the SMLD technique. The third section provides the flow governing equations, the additional transport equations of the combustion models, and numerical solution. The numerical results are discussed in the last section for the case of a subsonic CH 4 -air flame and a H 2 -air Mach 2 flame. The analysis is performed in order to validate the applicability of the developed models. In the first case, the flow is very simple and the steady flamelet equation well represents the phenomenon; in the second case the assumption of a stoichiometric scalar dissipation rate fails because of the complexity of the turbulence effects.