A Layered Description of a Premixed Flame Stabilized in Stagnating Turbulence

ABSTRACT In the thin-flame regime of turbulent premixed combustion, instantaneous progress variable gradients are essentially fixed by propagating flamelets. The laminar flamelet internal structure thus imposes, at least to some extent, the progress variable one-point one-time statistics, i.e., the progress variable PDF . In the present study a generalized presumed probability density function (PDF) shape is considered to account for possible departures from the thin-flame limit. The parameters that characterize this PDF shape depend on the ratio of the Kolmogorov length scale to laminar flame width. The corresponding framework may be associated to a layered description of the turbulent flame brush (TFB) with the thicknesses of the different sub-layers determined from the knowledge of both the Reynolds and Karlovitz numbers, ReT and Ka. It incorporates boundary layers on both sides of the TFB, which are associated to the finite thickness of the local flamelets, i.e., small but finite values of the Karlovitz number. This description leads to a modeling proposal for premixed turbulent combustion. A tabulation is constructed based on canonical one-dimensional planar unstrained laminar premixed flame computations and the different quantities are tabulated as functions of ReT and Ka. The final closure is applied to the numerical simulations of premixed flames stabilized in stagnating turbulence.


Introduction
The development of computational fluid dynamics (CFD) tools over the last decades has allowed their current applications to the numerical simulations of various industrial turbulent flows featuring complex geometries and realistic conditions. However, the physical phenomena involved in turbulent reactive flows display such a wide range of characteristic scales that the statistical, i.e., Reynolds Averaged Navier-Stokes (RANS), representation or the filtered, i.e., large eddy simulations (LES), description of these flows are still required to proceed with their numerical simulations. The direct numerical simulation (DNS) of such reactive flows in complex geometry remains indeed unlikely in the foreseeable future. The use of LES may provide a satisfactory description of the large-scale turbulent dynamics at the price of important computational efforts and costs so that its practical use for engineering purposes still remains limited. Accordingly, RANS and U-RANS approaches remain the reference tools to proceed with the industrial design of large-scale energy conversion devices. In practice, the resort to any of these two frameworks (i.e., RANS or LES) does involve the use of some numerical and modeling parameters so that the corresponding computations cannot be considered as fully predictive thus far. Therefore, the models involved in these numerical simulations and their behaviors must be perfectly understood so as to obtain results that recover relevant trends and can be useful to complement experimental data. The objective of the present study is to provide the analysis and detailed sensitivity investigation of a recent modeling proposal through its application to the numerical simulation of an impinging premixed turbulent flame.
One specific difficulty associated with the numerical simulation of turbulent reactive flows lies in the domains of validity of turbulent combustion models, which remain quite restricted. For instance, the main part of the flame brush can be often described as belonging to the flamelet regime, which is characterized by infinitely small chemical time scales compared to characteristic turbulent time scales. However, on the one hand, the local increase of turbulent kinetic energy from one given location to the other may lead to some thickening of the flame associated to finite-rate chemistry effects. On the other hand, towards the fresh reactants side, the characteristic chemical time becomes infinitely large in comparison with the characteristic turbulent time scales, thus making irrelevant the flamelet assumption to describe this part of the flame brush, i.e., the leading edge. Therefore, the turbulent combustion models proposed for the numerical simulation of premixed reactive flowfields must be as simple as possible but must also be realistic enough to handle such a wide range of variations of the turbulent combustion regime.
Here, we choose to analyze a simple model for the mean reaction rate, able to change its functional form inside the computational domain, depending on the turbulent combustion regime encountered. The relevant criterion accounting for these possible changes is based herein on a layered description of the turbulent flame brush structure that depends upon both Reynolds and Karlovitz numbers.

Finite thickness flamelets statistical description
In this study, fully premixed conditions are studied so that a unique progress variable c (x, t), defined to be c ; 0 in the fresh unburned mixture and c ; 1 in combustion products, is considered to characterize the composition of the mixture. Accordingly, Lewis numbers for all species are close to unity, and heat losses as well as composition inhomogeneities are not addressed. Thus, the Favre averaged progress variablec ¼ ρc=ρ is the solution of the following balance equation: @ ρc @t þ @ @x k ρṽ kc þ ρv 00 k c 00 À ρD where v k denotes the velocity field, D the molecular diffusion coefficient, and Ω c ¼ ρω c is the chemical reaction rate expressed in kg · m -3 · s -1 . In the above equation, the quantities c 00 ¼ c Àc and v 00 k ¼ v k Àṽ k denote Favre fluctuations of the progress variable and velocity components, respectively. The quantity ρv 00 k c 00 denotes the turbulent fluxes of the progress variable; the reader is referred to Robin et al. (2011) and references therein for further details concerning its modeling. Under the low Mach number assumption, the thermodynamic contribution p 0 ðx; tÞ to the pressure field pðx; tÞ is considered to be homogeneous, i.e., p 0 ðx; tÞ ¼ p 0 ðtÞ. Therefore, density may be expressed as a function of the composition variable only, i.e., ρðcÞ ¼ ρ u =ð1 þ τcÞ, where τ ¼ ðρ u =ρ b Þ À 1 is the expansion factor.
As it is well known, the closure of the mean chemical rate Ω c appearing in the righthand side (RHS) of Eq. (1) is one of the most important sub-problems in the field of turbulent combustion modeling and various expressions or functions of the mean progress variablec, i.e., Ω c ðcÞ ¼ Ω c ðcÞ ¼ ρðcÞω c ðcÞ, can be found in the literature. Most of the available closures have been developed under the flamelet hypothesis (Bray, 2016). In the corresponding flamelet regime of turbulent combustion (Borghi, 1985;Peters 1986;Williams, 1976), it is standardly assumed that the molecular processes, i.e., molecular diffusion and chemical reactions, are the dominant processes. However, the laminar flame thickness l f does not provide any firm information about the real thickness of flamelets present in the turbulent flows. Therefore, we proposed recently to introduce a characteristic length scale l x , defined as the relevant turbulent flamelet length scale at which molecular processes involved in the flame become the dominant processes (Kha et al., 2016). In the first step of the analysis, the value of l x has been approximated by the Kolmogorov length scale η K . The ratio between the laminar flame thickness and the turbulent flamelets length scale Rf x is thus considered to be of the order of the square root of the Karlovitz number: . Accordingly, the extreme values of the progress variable that delineate the flamelet profile are no longer zero and unity but will be denoted c min and c max . Beyond these limits, the laminar flame structure is indeed expected to be altered by the turbulent flowfield (Anand and Pope, 1985;Mura et al., 2003). These values will depend on the Karlovitz number. For instance, as established in the recent analysis by Kha et al. (2016), we have c min ¼ exp À1=Rf x ð Þand c max À 1 in the limit of infinitely large values of the Zel'dovich number, such that δ c ¼ c max À c min ¼ 1 À exp À1=Rf x ð Þ. Figure 1 displays the evolution of δ c as a function of Ka for two distinct approximations of l x . The first corresponds to the one discussed above, i.e., l x ¼ η K , which leads to δ c ¼ 1 À expðÀ1= ffiffiffiffiffiffi Ka p Þ. The second assumption implies l x ¼ l G with l G being the Gibson scale introduced by Peters (1986Peters ( , 1992. It may be worth recalling here that the Gibson scale is determined by setting the eddy turnover velocity to the laminar flamelet propagation velocity S L , in such a manner that l G ¼ l T : S L =v RMS ð Þ 3 with v RMS the characteristic turbulent velocity fluctuation such that Re T ¼ v RMS : l T =ν. In this case, it can be readily shown that δ c ¼ 1 À expðÀ1=Ka 2 Þ.
Considering that the flamelets characteristic scales l x may differ from the laminar flame thickness, l f , the following expression of the progress variable PDF is introduced: where α c and β c denote the probability that c equals c min and the probability that c equals c max , respectively. This PDF can be thought as a generalization of the well-known BML form (Bray and Moss, 1977), which is recovered in the limit associated to infinitely large values of the Damköhler number ðγ c ! 0Þ and infinitely small values of the Karlovitz number c min ! 0 ð Þand c max ! 1 ð Þ . In the above expressions, the functionf ðcÞ represents the PDF conditioned inside the internal structure of flamelets. This function can be deduced from the structure of a one-dimensional planar unstrained premixed laminar flame:f where C Ã 1 ¼ C 1 l f is a nondimensionalized normalization constant: As emphasized above, the values of c min and c max may differ from zero and unity, i.e., the Karlovitz number is not infinitely small. Thus, the progress variable PDF, as defined by Eq.
(2), as well as the internal PDF and its nondimensionalized normalization constant C Ã 1 , take into account the finite thickness of flamelets.
Accordingly, firstand second-order moments of the progress variable are: with and It seems worth emphasizing that, according to Eq.
(2), we havec ¼ c min forc 2 0; c min ½ andc ¼ c max forc 2 c max ; 1 ½ , while the second-order centered moment f c 002 is zero for mean progress variable valuesc that lie in 0; c min ½ [ c max ; 1 ½ . The probability to find flamelets is proportional to the intermittency factor between the fresh reactants and fully burned gasescð1 ÀcÞ and to the departure from the maximum segregation rate value S max À S [see, for instance, Bray et al. (2006)]: where the segregation rate S reaches its maximum value S max ; which, in the limit of large values of the Reynolds number, can be expressed as follows: The above equality indeed provides a generalized expression of the maximum segregation rate consistent with Eq. (2). Thus, the expression of the mean chemical rate can be written as: where I ω c is the mean chemical rate conditioned inside flamelet structures, a quantity obtained from laminar flame computations: Finally, the modified form of the PDF proposed above evidences the existence of boundary layers associated to the flame edges.

Properties of finite thickness flamelets
The properties of laminar flamelets are now studied in detail. Different levels of description can be retained to conduct the analysis. The limit of infinitely large values of the Zel'dovich number β provides a theoretical framework to get analytical results. For instance, in this limit, the progress variable evolution is given by cðx Ã Þ ¼ expðx Ã Þ; where the normalized coordinate x Ã ¼ x l f has been introduced.
For the purpose of analyzing the influence of finite values of the Zel'dovich number, laminar flame calculations based on simplified chemistry and transport properties may be first considered. The numerical calculations are performed by using the Mathematica TM software, and the conservation equations are solved in a frame of reference attached to the flame so that only the steady-state solution is evaluated. Thus, the transport equation of the progress variable writes: where v u denotes the velocity in the fresh reactants, and the reaction rate is given by: where B is the reduced pre-exponential factor and α ¼ ðT b À T u Þ=T b is a temperature factor (Williams, 1985). The chemical reaction rate can be recast in Ω c ¼ ρ u B F ðc; τ; βÞ with the function F ðc; τ; βÞ defined as follows: Then, the dimensionless form of the transport equation (14) is: where ω c ¼ Ω c D u ρ u v 2 u ¼ EF ðc; τ; βÞ is the dimensionless chemical reaction rate. The quantity D u denotes the molecular diffusion coefficient at the cold edge of the flame and the nondimensional number E is defined by E ¼ B D u v 2 u . The above equation is solved in the domain c 1 À by following the procedure of Bray et al. (2006). More realistic laminar flame structures may be evaluated from a detailed chemical mechanism using the Cantera library. A methane air flame is considered here for a value of the equivalence ratio Φ ¼ 1:0. Chemistry is described with the detailed GRI3.0 mechanism (Bowman et al., 1997). The resulting profiles of the progress variable c(x) and the corresponding scalar gradient dc=dx are presented in Figure 2 with the progress variable evolution deduced from ρ ¼ ρ u ð1 þ τcÞ. The burning velocity value is 38.6 cm/s. The detailed chemical mechanism leads to a slow evolution of the progress variable from c ¼ 0:9 to c ¼ 1:0 on the burned gases side. This is clearly visible on the scalar gradient profiles (see Figure 2). It is expected that turbulence may alter the laminar flamelet profile in this specific zone where the scalar gradients are rather weak and some experimental data bring support to this picture (Chen and Bilger 2002;Chen et al., 1996).
The conditional scalar moments I 1 and I 2 , see Eqs. (7) and (8), can be calculated from the knowledge of the laminar flame structure and associated internal PDF shapef ðcÞ, and the flamelet boundaries c min and c max . As emphasized above, these flamelet boundaries c min and c max are related to the value of the Karlovitz number. Therefore, the quantities defined by Eqs. (7), (8), and (13) are related to both the flamelet inner structure and the value of the Karlovitz number. In the case of infinitely large values of the Zel'dovich number, i.e., β ! 1, the upper boundary is such that c max ! 1 and the lower limit c min is related to the Karlovitz number and may be expressed from c min ¼ expðÀ1 ffiffiffiffiffiffi Ka p Þ. Moreover, in this case, the scalar gradient through the laminar flame structure is given by dc=dx Ã ¼ c, except in the reaction zone ðc ! 1Þ. Thus, the conditional scalar moments I 1 and I 2 can be obtained analytically: where the expression of the PDFf ðcÞ has been used and, according to Eq. (4), the normalization constant C Ã 1 has been approximated by ffiffiffiffiffiffi Ka p . Then, considering the relation between c min and the Karlovitz number, the conditional scalar moments can be expressed as functions of c min :  20) or as functions of the Karlovitz number: Equations (21) and (22) provide the firstand second-order scalar moments conditioned inside flamelets for large values of the Zel'dovich number as functions of the Karlovitz number only. Flamelet structures obtained from the numerical calculations discussed at the beginning of this section can now be used to evaluate these conditional quantities. Figure 3 displays the integral quantities as functions of the Karlovitz number. It also shows the difference between I 1 and I 2 , which is an important quantity for the evaluation of the mean chemical rate [see Eq. (12)]. It is worth noticing that the ratio I 2 /I 1 is also an important quantity since, in the limit of infinitely small values of Ka, one can show that it imposes the segregation rate at the leading edge 1 : where the following two identities have been used: f Kha et al. (2016) for further details].
The first part of Figure 3 (on the top) corresponds to the case of large values of the Zel'dovich number, i.e., eqs. (21) and (22), the middle one to the case of finite Zel'dovich number ðβ ¼ 10Þ but simplified chemistry and the third figure (on the bottom) corresponds to detailed chemistry calculations ðΦ ¼ 1Þ. These results show that, for vanishingly small values of the Karlovitz number, I 1 and I 2 vanish for either large values of the Zel'dovich number or finite Zel'dovich number provided that simplified chemistry is considered. However, when detailed chemistry is considered, I 1 and I 2 tends to unity for vanishingly small values of the Karlovitz number, which is ascribed to the specific behavior of the progress variable gradient observed at the burned gas side when detailed chemistry is considered. It is worth noting that, even in the small Karlovitz number limit, the behaviors of I 1 and I 2 do depend on the level of description retained for the chemistry, i.e., single step or detailed chemistry. We will see below that this dependence, however, disappears as soon as we consider the difference between I 1 and I 2 , and their ratio I 2 =I 1 .
The quantities I 1 and I 2 indeed do not appear separately in the mean chemical rate expression but through the difference I 1 À I 2 as well as the ratio I 2 /I 1 . Whatever the description of the chemistry, i.e., single step or detailed chemistry, it is indeed 0.2 0.4 0.6 0.8 1 I 1 I 2 I 1 -I 2 I 2 /I 1 Ka Figure 3. Evolutions of laminar flamelet properties I 1 , I 2 , I 1 À I 2 , and I 2 =I 1 as functions of the Karlovitz number. The results reported on the top correspond to the infinitely large Zel'dovich number limit ðβ ! 1Þ, those in the middle to one-step chemistry calculations conducted with β ¼ 10, and those presented at the bottom have been obtained from detailed chemistry computations ðΦ ¼ 1:0Þ. remarkable that, for small values of the Karlovitz number Ka, the difference I 1 À I 2 vanishes and the ratio I 2 /I 1 tends to unity (see second and third graphs of Figure 3). The behavior observed for I 2 /I 1 in the infinitely large Zel'dovich limit appears to be more specific since, as previously underlined in the recent study by Kha et al. (2016), it tends to 0.5 when Ka ! 0. Figure 3 also shows that taking into account the reaction zone, i.e., finite values of the Zel'dovich number, slightly modifies the value of the ratio I 2 /I 1 only for small values of the Karlovitz number. Another important result is that the difference I 1 À I 2 is a constant approximately equal to 0.2 for a large range of Karlovitz number ð0:2 < Ka < 1Þ. Figure 4 shows the evolution of the flamelet properties I 1 À c min and c max À I 1 with the Karlovitz number. These two quantities behave in a way very similar to the one observed for I 1 and I 2 , as shown in Figure 3. Moreover, for large values of the Karlovitz number ð0:5 < Ka < 1Þ, they are almost constant whatever the laminar flame structure considered. We shall see later on that these two differences I 1 À c min and c max À I 1 dictate the thickness of the leading and trailing edge layers of the turbulent flame brush. Figure 4. Evolution of two flamelet properties I 1 À c min (top) and c max À I 1 (bottom) as functions of the Karlovitz number for different laminar flame structure.

Ka Ka
Layered description of the turbulent flame and tabulation of finite thickness flamelets characteristics The turbulent flame brush structure as a whole must be now analyzed in the light of the previous discussion. Figure 5 shows the typical sketch of a turbulent flame brush consistent with the finite thickness flamelets presumed PDF of Eq. (2). Accordingly, the flame brush is bounded by c min <c < c max and we use the following layered description of the turbulent flame brush (Kha et al., 2016): • Upstream of the turbulent flame brush (c < c min ), a turbulent preheat layer must be considered. The chemical reactions are exponentially small in this zone and can therefore be discarded from the analysis; the PDF support reduces toPðc; xÞ ¼ δðc ÀcÞ. • The appearance of the first flamelet structure corresponds to c min <c < c le ; this is the leading edge of the turbulent flame brush. In this region, the PDF is approximated byPðcÞ ¼ ð1 À γ c Þδðc À c min Þ þ γ cf ðcÞ. • In the region c le <c < c te , the composition changes abruptly from c ¼ c min to c ¼ c max . The molecular processes involved in flames occur at the smallest scales of the turbulent flow and the probability γ c to lie inside flamelets is small enough to consider that second-order averaged quantities (variance, scalar turbulent fluxes, etc.) are proportional to ðc max ÀcÞðc À c min Þ. Then the mean chemical rate may be evaluated from Eq. (12). • The trailing edge of the flame brush ðc te <c < c max Þ is associated to the disappearance of the last flamelet structure. In this region, the PDF is approximated bỹ According to this description of the structure of the turbulent flame brush, the mean chemical rate is expressed in five different sub-layers, respectively: the turbulent preheat layer, the leading edge, the main (intermittent) layer, the trailing edge, and the weak reactive layer. It is necessary to evaluate the boundaries of each of these layers, namely c min , c le , c te , and c max . It has been shown earlier that the value of the Karlovitz number determines the flame brush boundaries, i.e., c min and c max . Equation (5), once applied to the leading and trailing edges, yields: and In these expressions the approximation introduced previously by Kha et al. (2016) is used to evaluate γ le and γ te ; we assume that the probability to be inside these smallest structures is proportional to the ratio of their thickness l x ¼ l ηκ to the largest scale of the turbulence, i.e., the integral length scale l T , which is responsible for the largest wrinkles of the flame front. The sensitivity of the computational results to the parameter κ will be evaluated afterwards. From these expressions, one can deduce the thickness in the composition space of the leading and trailing edge layers as δ le ¼ c le À c min ¼ γ le ðI 1 À c min Þ and δ te ¼ c max À c te ¼ γ te ðc max À I 1 Þ, respectively. Considering the dependency of the different boundaries on the Reynolds and Karlovitz numbers, the mean quantities can be tabulated as functions of these two parameters. From a general point of view, tabulations may either be completed before the numerical simulations or generated on-the-fly during the simulation. The tabulation methodology we use here follows the first of these two strategies and is based on the canonical 1D unstrained laminar premixed flame computations that have been presented in the previous section. However, in the present study it is parameterized by the values of the Reynolds and Karlovitz number, an original feature. From the couple (Ka, Re T ), the values c min (Ka), c le (Ka, Re T ), c te (Ka, Re T ), and c max (Ka) are determined as well as the different sub-layers of the turbulent flame brush. These evaluations are based on the Cantera computations of premixed laminar flamelets incorporating detailed chemistry. The corresponding profiles of chemical species and temperature are themselves parameterized as functions of a progress variable c, defined as a marker of the heat released by chemical reaction through ρ ¼ ρ u =ð1 þ τcÞ and weighted by the progress variable PDFPðcÞ introduced above. The tabulations includes not only the boundaries of the sub-layers but incorporates also mean species mass fractions as well as any other quantity that is required to evaluate the mean chemical rate, especially I 1 , I 2 , and I ω c .
From the above layered description, the mean chemical rate expressions are given by: ω c ðcÞ ¼ ω c ðcÞ % 0c 2 ½0; c min ð1 À γ c Þω c ðc min Þ þ γ c I Ã ω cc 2 ½c min ; c le ðS max À SÞcð1 ÀcÞI ω c =ðI 1 À I 2 Þc 2 ½c le ; c te ð1 À γ c Þω c ðc max Þ þ γ c I ÃÃ ω cc 2 ½c te ; c max ω c ðcÞc 2 ½c max ; 1 8 > > > > < > > > > : where I Ã ω c and I ÃÃ ω c are defined by I c ÃÃ ðcÞ ω c ðcÞf ðcÞdc. However, it has been experienced in the numerical simulations that the leading and trailing edges sub-layers are quite difficult to discriminate from the turbulent preheat layer ðc < c min Þ and weak reactive layer ðc > c max Þ. Due to their small thicknesses it is indeed difficult to capture these thin layers on conventional computational meshes. Therefore, a reduced three-layer model of mean chemical rate is used to conduct the numerical simulations:ω c ðcÞ ¼ ω c ðcÞ % 0c 2 ½0; c le γ c I ω cc 2 ½c le ; c te ω c ðcÞc 2 ½c te ; 1 8 < : ( Thus, in an attempt to generalize available fast chemistry closures, the present proposal consists in modifying the closure of the mean chemical at the leading and trailing edges of the turbulent flame brush. As emphasized above, in the present application the leading edge sub-layer (LEL) is merged with the turbulent preheat layer (TPL) and it is acknowledged that in the latter the mean progress variable evolution is driven by mean convection and turbulent transport. It is noteworthy that some specific features of the turbulent flame exhibit some sensitivity to the quenching criterion c le . For instance, previous one-dimensional traveling wave (TW) analyses have established that, due to the quenching influence, it is quite difficult to recover the classical Kolmogorov-Petrovskii-Piskunov (KPP) solution (Kolmogorov et al., 1937) from a computational estimation of the turbulent flame speed, see Catlin and Lindstedt (1991) and Sabelnikov et al. (1998). The detailed investigation conducted by Corvellec (1998) establishes that the convergence to the theoretical KPP solution, as the quenching value is reduced, is quite slow. Different propagation characteristics of the mean reacting flow may be obtained depending on the accuracy of the numerical integration and on the treatment of the boundary conditions. This study also emphasized that, in practical CFD computations of turbulent premixed flames, i.e., 2D and 3D computations performed with standard discretization schemes, it is almost impossible to get rid of the computational procedure influence and the numerical solution will remain altered by the mesh resolution and boundary conditions (Sabelnikov et al., 1998). The present procedure is quite a contrast to the above discussion. Indeed, the present study is aimed at avoiding the systematic resort to the infinitely fast chemistry approximation (i.e., the telegraphic functional formcð1 ÀcÞ of the mean chemical rate) in the whole turbulent flame brush thickness (TFB). As emphasized above, such an hypothesis does not hold on the leading and trailing edges of the TFB but is implied in the KPP TW-analyses. A quenching value is finally introduced but it is based on physical grounds and consideration of small but finite values of the Karlovitz number. It involves the thicknesses in the progress variable space of both the turbulent preheat layer, from 0 to c min , and leading edge layer, from c min to c le . The closure of the mean reaction rate provided by Eq. (28) is now assessed through its application to a practical geometry.

Description of the numerical setup
The experimental setup that is considered to proceed with the numerical simulations corresponds to a premixed flame stabilized in stagnating turbulence. This geometry has been studied inter alia by Cho et al. (1988Cho et al. ( , 1989. In these experiments a turbulent stream of reactants at room temperature impinges on a refractive surface so as to stabilize a planar premixed turbulent flame. For many years it has been recognized that such impinging flames provide a reference test case for the validation of premixed turbulent combustion models (Bray et al., 1991(Bray et al., , 1994(Bray et al., , 2000Dong et al., 2013;Lee and Huh, 2004;Sabelnikov and Lipatnikov, 2011). In the reference papers by Cho et al. (1988Cho et al. ( , 1989) data on the velocity and temperature characteristics of these flames are provided. The corresponding set of experimental data is considered as a reference database to assess turbulent combustion closures and it has been used to this purpose in many other studies [see, for instance, the recent analyses conducted by Sabelnikov and Lipatnikov (2011) or Dong et al. (2013)]. A uniform axisymmetric flow of methane and air is provided by a 50-mm diameter nozzle. This jet of reactants is surrounded by a co-flow of air with an outer diameter of 100 mm. The turbulent flames are stabilized in the stagnating flow thanks to a refractive plate placed at a distance of 75 mm above the jet exit plane. The bulk flow velocity of the jet is held constant at 5.0 m/s and the incoming turbulence is generated by a grid. The characteristic integral scale of turbulence l T can be estimated 2 to lie between one and two millimeters (Dong et al., 2013). The fluctuating velocity RMS v RMS is approximately equal to 0.3 m/s. The equivalence ratio of the methane-air mixture is unity ðΦ ¼ 1Þ.
From the computational point of view, a 2D mesh, made of approximately 15,900 computational cells, is used to represent a small part of the physical domain (see Figure 6). The inlet condition at the bottom of the computational domain is imposed by using the experimental data while an outlet boundary condition is applied on the top right side. Boundary conditions on the left side of the computational domain are given by the symmetry of the flow and on the right side it is approximated by a special condition, which can be either an inlet or an outlet. The numerical simulations are performed with the Code Saturne, a parallel general-purpose 3D low-Mach-number computational fluid dynamics (CFD) solver developed by EDF (Archambeau et al., 2004). It relies on a finitevolume approximation of the Navier-Stokes equations. The set of equations to be solved consists of the averaged Navier-Stokes equations completed with equations for turbulence modeling and for additional scalars (enthalpy, concentration of species, mixture fraction, as well as any other user-defined scalars). The time marching scheme is based on a velocity field prediction step followed by a pressure-correction step. Equations for turbulence and scalars are solved separately afterwards. The discretization in space is based on the fully conservative, unstructured finite-volume framework, with a fully colocated arrangement for all variables. For the present application balance equations for the mean momentum are solved in conjunction with a classical k-ε closure and a transport equation for the mean progress variable [see Eq.
(1)]. The boundary conditions for turbulent quantities k and ε at the inlet jet of methane and air are, respectively, k 0 ¼ 3ðu RMS Þ 2 =2 and ε 0 ¼ C μ k 3=2 0 =l T with C μ = 0.09. A passive scalar, i.e., mixture fraction , is also considered to follow the mixing between the jet of reactants and the coflowing jet of air.

Reference case
To proceed with the numerical simulation, the mean chemical rate in Eq. (1) is represented by the closure discussed in the previous section, see Eq. (28). It is worth noting that, in a first step of the analysis, the progress variable scalar flux ρv 00 k c 00 is closed by using the standard, i.e., gradient-diffusion law: ρv 00 k c 00 ¼ À ρD T ð@c=@x k Þ with D T ¼ ν T =Sc T . The value of the turbulent Schmidt number is set to its reference value Sc T ¼ 0:7 (Launder, 1976). The present set of results has been obtained with l T = 1 mm and a value of the constant κ equal to 0.05, see Eq. (26). The sensitivity of the computational result to the choice of this value will be evaluated subsequently.
It is noteworthy that the model used to perform the present set of numerical simulations has been set forth for fully premixed conditions only. Further work is clearly required to extend it to partially premixed conditions and we therefore did not attempt to describe turbulent combustion in regions where the flow of premixed reactants is subject to dilution by ambiant air. The mean chemical rate is thus set to zero as soon as the mixture fraction significantly differs from its value at the nozzle exit. This can be seen in Figures 7 and 8 reporting the field of the mixture fraction superimposed with progress variable isolines and the mean progress variable fieldc superimposed with mean-chemical rate isolines, respectively. In Figure 8 we also compare the mean field of the progress variable obtained from the transport equation (1) to another expression of the progress variable that can evaluated from the tabulated detailed chemistry: where Y eq CO and Y eq CO2 denote the CO and CO 2 mass fractions at equilibrium. The profiles of c min , c le , c te , and c max , evaluated along the symmetry axis, are reported in Figure 9. This figure shows that the values of c min are too small to resolve the turbulent preheat sublayer [0, c min ]. It also firmly confirms that it would have been quite difficult to discriminate between c te and c max , except at the price of prohibitively high computational costs. The profiles of the mean axial velocity component and mean  progress variable are presented together with experimental data in Figure 10. A quite satisfactory quantitative agreement is obtained between the numerical simulation results and the experimental data. The sensivity of the corresponding results to changes in the turbulent Schmidt number value Sc T is illustrated in Figure 11. The sensivity of the computational results to a modification of the equivalence ratio Φ of the incoming stream of reactants has been also examined and, in Figure 12, results obtained for Φ = 0.8 are compared to those obtained for stoichiometric conditions, i.e., Φ = 1, which correspond to experimental operating conditions. It is interesting to see that, for a value of κ held constant, the results confirm that the model is sensitive to mixture composition variations. Finally, Figure 13 provides a comparison with another set of results issued from a computation performed with c le ¼ 0 and c te ¼ 1:0. It evidences the influence of the consideration of the leading and trailing edges sublayers onto the computational results.

Impact of the CGD modeling on the first moment quantities
Since it is well known that the experimental geometry under consideration is subject to counter-gradient diffusion (Bray et al., 1998), we examine the possible influence of these phenomena on the computational results. Therefore the next set of computational results has been obtained by using an algebraic closure that accounts for the countergradient diffusion (CGD) influence. The corresponding closure may be written as follows, ρv 00 k c 00 ¼ ðτS L ρc 00 2 λ À ρD T jjÑcjjÞ e M k ; where e M k denotes the component of the mean orientationM of the flux, which is approximated by e M ¼ Ñc=jjÑcjj and λ is a parameter related to the fluctuations of the local flame front. Such a closure was previously used by Dong et al. (2013) to conduct the numerical simulation of the same experimental test case. Figure 14 reports the scalar flux as obtained (i) from the classical gradient-diffusion law (GD), (ii) from the above generalized closure (CGD), and (iii) experimental data. From this figure representing the scalar flux as a function of the mean progress variable, it cannot be disputed that, in the present conditions, CGD processes are  Cho et al. (1988). Bottom: Comparison between the axial profile of the mean axial velocity component and experimental data documented by Cho et al. (1988). essential to represent satisfactorily the scalar fluxes. A comparison between the two sets of computational results obtained with a closure that incorporates or not the CGD effects is reported in Figures 15. These figures report the profiles of the mean progress variable and the mean axial component of the velocity respectively. It is quite interesting to see that, compared to the results obtained from the generalized CGD closure, the standard GD closure leads to a satisfactory level of agreement with experimental  Figure 11. Top: Comparison between the axial profile of the mean progress variable and experimental data documented by Cho et al. (1988). Bottom: Comparison between the axial profile of the mean axial velocity component and experimental data documented by Cho et al. (1988). Different values of the turbulent Schmidt number Sc T are considered.
data. Therefore, for the purpose of the present investigation, which is not focused on the representation of counter-gradient diffusion (CGD) and associated flame-generated turbulence (FGT) effects, the standard GD approximation will be kept throughout the rest of this study. x Figure 12. Top: Comparison between the axial profile of the mean progress variable and experimental data documented by Cho et al. (1988). Bottom: Comparison between the axial profile of the mean axial velocity component and experimental data documented by Cho et al. (1988). Calculations were performed for two distinct values of the equivalence ratio: Φ ¼ 0:8 and Φ ¼ 1:0; the latter corresponds to experimental conditions.

Leading edge layer thickness: A sensitivity analysis
In the model for the mean chemical rate introduced above, an important and critical parameter of the final closure is the value of the constant κ, which relates through Eq. (26) the probability to lie in the flamelet structure to the turbulent Reynolds number. Equation (25)  x Figure 13. Top: Comparison between the axial profile of the mean progress variable and experimental data documented by Cho et al. (1988). Bottom: Comparison between the axial profile of the mean axial velocity component and experimental data documented by Cho et al. (1988). Calculations were performed by considering either the multilayered flame brush representation (continuous line) or by forcing c le and c te to zero and unity, respectively (dotted line).
the value of this parameter. It is therefore decided to analyze the dependency of the computational results on the specification of the value of κ. Figure 16 reports the profiles of the progress variable fluxes and mean axial component of the velocity respectively. These results have been obtained by using three distinct values of κ, which are κ ¼ 0:05, 0:5, and 5. This offers quite a large range of variations for this parameter since the maximal value is 100 times larger than the smallest one. As expected, the obtained results display some dependency on the value used to perform the numerical simulation but it must be acknowledged that the level of sensitivity does not seem to be critical and we will see below that the sensitivity of the computational results to the specification of the inlet boundary conditions is much larger.
Sensitivity of the computational results to the specification of the inlet flow integral length scale Figure 17 reports the profiles of the progress variable fluxes and mean axial component of the velocity respectively. These results have been obtained for three distinct values of the integral length scale imposed at the inlet boundary (jet nozzle exit). The corresponding values are l T ¼ 0:5 mm, l T ¼ 1 mm, and finally l T ¼ 2 mm. Unsurprisingly, the obtained results display an important sensitivity to the specification of l T , which further emphasizes the crucial need for well-documented experimental databases featuring also well-defined boundary conditions.

Conclusions
In this work, a presumed PDF of the progress variable has been used to account for the influence of finite thickness flamelets and associated progress variable boundaries. This leads to the introduction of a layered description of the turbulent flame brush (TFB) based on the consideration of both Reynolds and Karlovitz numbers effects. This description is subsequently used to propose a first-order closure of the mean x Figure 15. Top: Comparison between the axial profile of the mean progress variable and experimental data documented by Cho et al. (1988). Bottom: Comparison between the axial profile of the mean axial velocity component and experimental data documented by Cho et al. (1988). GD denotes computational results obtained with the standard diffusion-like closure of the progress variable turbulent fluxes. CGD denotes computational results obtained with Eq. (30).
progress variable. The closure is derived under the assumption of sufficiently large values of the turbulence Reynolds number, with the molecular fluxes neglected with respect to its turbulent counterpart, but accounts for finite values of the Karlovitz number. The corresponding formalism has been used in conjunction with tabulated chemistry to compute Laboratory flames and the preliminary results obtained are rather encouraging.
A sensitivity analysis has also been conducted with respect to the main parameter of the modeling proposal and, as expected, the obtained results display some dependency Comparison between the axial profile of the mean progress variable and experimental data documented by Cho et al. (1988). Bottom: Comparison between the axial profile of the mean axial velocity component and experimental data documented by Cho et al. (1988).
on the value used to perform the numerical simulation. However, it has been undoubtedly shown that there exists an equivalent sensitivity to the specification of the inlet boundary condition, i.e., the turbulence integral length scale at the nozzle exit. To conclude, it is quite clear that it would be useful to conduct additional validations of the proposed closure against other well-documented experimental data. It would be especially interesting to analyze the response of the computational model in conditions featuring stronger departures from the thin laminar flamelet combustion regime. In  Figure 17. Sensitivity of the computational results to the specification of the integral length scale l T at the nozzle exit. Top: Comparison between the axial profile of the mean progress variable and experimental data documented by Cho et al. (1988). Bottom: Comparison between the axial profile of the mean axial velocity component and experimental data documented by Cho et al. (1988). this respect, the highly turbulent premixed V-shaped flames studied by Renou and coworkers (Mazellier et al., 2010) provide an interesting test case so as to proceed with this additional assessment of the proposed model.