Modeling of Progress Variable Variance Transport in Head-On Quenching of Turbulent Premixed Flames: A Direct Numerical Simulation Analysis

ABSTRACT The statistical behavior of the transport of reaction progress variable variance has been analyzed using three-dimensional direct numerical simulation (DNS) data for different values of Damköhler, Karlovitz, and global Lewis numbers in the context of head-on quenching of turbulent premixed flames by an inert isothermal wall. It has been found that reaction rate contribution to the variance transport acts as a leading order source, whereas the molecular dissipation term remains as the leading order sink for all cases considered here. However, all of the terms of the variance transport equation decay significantly in the near-wall region once the quenching starts. The existing models for the turbulent transport, reaction, and dissipation contributions to the variance transport do not adequately capture the near-wall behavior. The wall effects on the unclosed terms of the variance transport equation have been analyzed using explicitly Reynolds averaged DNS data and the existing closures of the unclosed terms have been modified to account for the near-wall effects. A-priori DNS analysis suggests that the proposed modifications to the existing closures for the unclosed terms of the variance transport equation provide satisfactory predictions both away from and near to the wall.


Introduction
The variance of reaction progress variable plays a pivotal role in the modeling of turbulent premixed combustion (Bray, 1980;Swaminathan and Bray, 2011). The magnitude of the variance of reaction progress variable is often necessary for the modeling of mean reaction rate in the context of Reynolds Averaged Navier-Stokes (RANS) simulations (Bray et al., 2006;Linstedt and Vaos, 1999;Mantel and Bilger, 1995;Moler et al., 1996;Mura et al., 2007;Ribert et al., 2005;Robin et al., 2006). The reaction progress variable c can be defined in terms of a suitable reactant mass fraction Y R in the following manner: According to Eq. (1), c increases monotonically from zero in the unburned gas (subscript 0) to unity in fully burned products (subscript 1). The variance of reaction progress variable is given by: f c 00 2 , whereq ¼ ρq=" ρ and q 00 ¼ q Àq represent the Favre average and Favre fluctuation of a general quantity q, respectively, with ρ being the gas density, and the Reynolds averaging is shown by the overbar. The scalar variance f c 00 2 is one of the important quantities for the flamelet (Bray et al., 1985;Linstedt and Vaos, 1999;Swaminathan and Bray, 2011) and conditional moment (Klimenko and Bilger, 1999;Swaminathan and Bilger, 2001) based closures. Consequently, the variance f c 00 2 is often required for the well-known eddy break up (EBU) models (Linstedt and Vaos, 1999;Swaminathan and Bray, 2005). Furthermore, f c 00 2 is an essential gradient of the tabulated chemistry based modeling of turbulent premixed combustion (Domingo et al., 2005;Savre et al., 2008).
Based on a presumed bi-modal probability density function (pdf) of c with impulses at c ¼ 0 and c ¼ 1:0 according to the Bray-Moss-Libby (BML) model (Bray et al., 1985), one obtains: f c 00 2 ¼c 1 Àc ð ÞþO γ c À Á where O γ c À Á is the burning mode contribution. The contribution of O γ c À Á can be neglected and f c 00 2 assumes its maximum possible valuec 1 Àc ð Þwhen P c ð Þ can be approximated by a bi-modal distribution with impulses at c ¼ 0 and c ¼ 1:0, and this condition is realized for high values of Damköhler number (i.e., Da ) 1), where the flame front is thinner than the Kolmogorov length scale, and the turbulent eddies do not affect the flame structure. However, O γ c À Á cannot be neglected for small values of Da (i.e., Da<1) and subsequently f c 00 2 remains smaller thanc 1 Àc ð Þ, and thus, it is necessary to solve variance transport equation along with other modeled conservation equations in the context of RANS simulations.  and Malkeson and Chakraborty (2010) analyzed the statistical behaviors of scalar variance transport in turbulent premixed and stratified flames, respectively. Furthermore,  demonstrated that global Lewis number Le has significant influences on the various terms of the transport equation of f c 00 2 . However, all of the aforementioned analyses have been carried out for flames, which are away from the wall, and the analysis of reaction progress variable variance f c 00 2 transport in the near-wall region during flame-wall interaction is yet to be addressed in existing literature. A number of previous studies Rutland, 1998, 2002;Bruneaux et al., 1996Bruneaux et al., , 1997Dabireau et al., 2003;Gruber et al., 2010Gruber et al., , 2012Lai and Chakraborty, 2015;Poinsot et al., 1993) analyzed flame-wall interaction using three-dimensional (3D) DNS data. Poinsot et al. (1993) discussed about possible wall functions in the context of flame surface density (FSD)-based closure using 2D DNS data of head-on quenching of turbulent premixed flames. Bruneaux et al. (1997) proposed near-wall modifications to the models of the unclosed terms of the FSD transport equation based on channel flow DNS data. Rutland (1998, 2002) addressed the near-wall closure of FSD and turbulent scalar flux for turbulent premixed flame-wall interaction. Dabireau et al. (2003) analyzed the statistical behavior of wall heat flux for both premixed and diffusion flames based on DNS data. Gruber et al. (2010Gruber et al. ( , 2012 carried out 3D detailed chemistry DNS of turbulent premixed flame interaction with an inert isothermal wall and indicated the presence of flame instabilities, which are yet to be understood in detail. Recently, Lai and Chakraborty (2015) have proposed near-wall modifications to a well-known scalar dissipation rate (SDR)-based mean reaction rate closure proposed by Bray (1980) and also extended an algebraic SDR closure for accurate prediction of the near-wall behavior using DNS data of head-on quenching of statistically planar turbulent premixed flames by an isothermal inert wall. Lai and Chakraborty (2015) also analyzed the effects of global Lewis number Le on the statistical behaviors of wall heat flux, flame quenching distance, and the closures of SDR and mean reaction rate in the context of head-on quenching of turbulent premixed flames. It is worth noting that none of the aforementioned analyses on flame-wall interaction concentrated on the statistical behavior of the variance f c 00 2 transport in the near-wall region. In order to address this gap in existing literature, the present analysis concentrates on the analysis of the variance f c 00 2 transport in the near-wall region for head-on quenching of statistically planar turbulent premixed flames with global Lewis number Le ¼ α T =D ¼ λ= ρC P D ð Þranging from 0.8-1.2 for different normalized values of root-mean-square turbulent velocity u 0 =S L and integral length scale l=δ th , where S L , δ th , α T , λ, C P , and D are the unstrained laminar burning velocity, thermal flame thickness, thermal diffusivity, thermal conductivity, specific heat at constant pressure, and mass diffusivity, respectively. Thus, the main objectives of the current analysis are: (1) To identify the near-wall effects on the statistical behavior of the unclosed terms of the scalar variance f c 00 2 transport equation.
(2) To propose modifications to the existing models of the unclosed terms of the variance f c 00 2 transport equation in order to account for the near-wall behavior.
The rest of the article will be organized as follows. The information related to mathematical background and numerical implementation will be provided in the next two sections. Following this, results will be presented and subsequently discussed. The main findings will be summarized and main conclusions will be drawn in the final section of this article.

Mathematical background
The present analysis considers a single-step Arrhenius-type irreversible chemical reaction (i.e., Reactants ! Products) for the purpose of computational economy as 3D DNS simulations with detailed chemistry remain prohibitively expensive for a detailed parametric analysis (Chen et al., 2009). An actual combustion process includes a number of species with different values of Le, but often the Lewis number of the deficient reactant is considered to be the characteristic Lewis number (Im and Chen, 2002;Mizomoto et al., 1984). Law and Kwon (2004) proposed a method for estimating the effective Lewis number based on heat release rate measurements, whereas Dinkelacker et al. (2011) proposed a methodology of estimating effective Lewis number as a linear combination the mole fractions of the mixture constituents. It is worth noting that several previous analyses (Clavin and Williams, 1982;Han and Huh, 2008;Haworth and Poinsot, 1992;Libby et al., 1983;Rutland and Trouvé, 1993;Sivashinsky, 1977;Trouvé and Poinsot, 1994) analyzed the effects of Le in isolation based on simple chemistry and the same approach has been adopted here. The instantaneous transport of reaction progress variable c is governed by: where _ ω is the reaction rate of reaction progress variable c. Reynolds averaging of Eq. (3) yields: The transport equation of the reaction progress variable variance f c 00 2 can be obtained using relation ρc 00 2 ¼ ρc 2 À " ρc 2 : À2ρu 00 j c 00 @c @x j |fflfflfflfflfflfflffl ffl{zfflfflfflfflfflfflffl ffl} whereε c ¼ ρDÑc 00 Á Ñc 00 =" ρ is the SDR of the reaction progress variable. In Eq. (5), D 1c is a closed term which denotes the molecular diffusion of f c 00 2 . The term T 1c is the turbulent transport term, whereas T 2c denotes the generation/destruction of f c 00 2 by the mean scalar gradient. The term T 3c is the reaction rate contribution and ÀD 2c ð Þ is the molecular dissipation term. The terms T 1c , T 2c ; T 3c and ÀD 2c ð Þare the unclosed term in the context of f c 00 2 transport. Equation (5) indicates that ÀD 2c ð Þ closure translates to the modeling of SDRε c . The modeling of T 1c , T 2c , T 3c , andε c in head-on quenching of turbulent premixed flames has been addressed here using explicitly Reynolds averaged 3D DNS data.

Numerical implementation
The conservation equations of mass, momentum, internal energy, and reaction progress variable for compressible reacting flows are solved in nondimensional form in the present analysis using a well-known DNS code SENGA (Jenkins and Cant, 1999). The simulation domain is taken to be a rectangular box of size 70:6δ Z Â 35:2δ Z Â 35:2δ Z , where δ Z ¼ α T0 =S L is Zel'dovich flame thickness with α T0 and S L being the thermal diffusivity of the unburned gas and unstrained laminar burning velocity, respectively. The simulation domain is discretized using a Cartesian grid of 512 Â 256 Â 256 ensuring 10 grid points across the thermal flame thickness δ th ¼ T ad À T 0 ð Þ =Max ÑT L , where T ad , T 0 , andT are the adiabatic flame temperature, unburned gas temperature, and the instantaneous dimensional temperature, respectively, and the subscript 'L' is used to refer to unstrained planar laminar flame quantities. A 10th-order central difference scheme has been used for the spatial discretisation for the internal grid points and the order of differentiation gradually drops to a one-sided 2ndorder scheme at the nonperiodic boundaries (Jenkins and Cant, 1999). A 3rd-order Runge-Kutta scheme (Wray, 1990) has been used for explicit time advancement. The reactive flow field has been initialized by a steady unstrained planar laminar premixed flame solution, which in turn has been superimposed on top of an initially homogeneous isotropic field of the turbulent velocity fluctuations away from the wall. The turbulent velocity field away from the wall has been initialized by a homogeneous isotropic incompressible field of turbulence, which was generated using a standard pseudo-spectral method (Rogallo, 1981) following the Batchelor-Townsend Spectrum (Batchelor and Townsend, 1948). The left-hand-side boundary in the x 1 -direction (i.e., x 1 ¼ 0) is taken to be the chemically inert isothermal wall with temperature T w ¼ T 0 , where no-slip boundary conditions are imposed and zero mass flux is specified in the wall normal direction. A partially nonreflecting outlet boundary condition is specified in the right-hand-side boundary in the x 1 -direction. Transverse directions are considered to be periodic. The nonperiodic boundaries have been specified using the Navier-Stokes Characteristic Boundary Conditions (NSCBC; Poinsot and Lele, 1992) technique. Three different global Lewis numbers (Le ¼ 0.8, 1.0, and 1.2) have been considered for this analysis and standard values are chosen for Zel'dovich number (i.e., β ¼ T ac T ad À T 0 ð Þ =T 2 ad ¼ 6), Prandtl number (i.e., Pr ¼ 0:7), and the ratio of specific heats (i.e., γ ¼ 1:4), where T ac is the activation temperature. The heat release parameter ¼ T ad À T 0 ð Þ =T 0 is taken to be 6.0. The simulations have been carried out for different initial values of normalized root mean square value of turbulent velocity u 0 =S L , Damköhler number , and turbulent Reynolds number Re t ¼ ρ 0 u 0 l=μ 0 (where ρ 0 and μ 0 are the unburned gas density and viscosity, respectively), which are listed in Table 1. It can be seen from Table 1 that the cases A, C, and E (B, C, and D) have the same values of Da (Ka), and Ka (Da) is modified to bring about the changes in Re t . The simulations have been carried out for a time when the maximum, mean, and minimum values of wall heat flux assume identical values following the flame quenching. The simulation time remains different for different cases but the simulations for all cases have been continued for t ! 12δ Z =S L , where 12δ Z =S L corresponds to 21, 30, 21, 15, and 21 initial eddy turnover times (i.e., t e ¼ l=u 0 ) for cases A-E, respectively. The nondimensional grid spacing next to the wall y þ ¼ u τ Δx=ν remains smaller than unity for all turbulent cases, where u τ ¼ ffiffiffiffiffiffiffiffiffi ffi τ w =ρ p , τ w , and ν are the friction velocity, wall shear stress, and kinematic viscosity, respectively.
The Reynolds/Favre averaging has been carried out by ensemble averaging the quantities over statistically homogeneous x 2 À x 3 directions at a given x 1 -location. The statistical convergence of the Reynolds/Favre averaged values have been assessed by comparing the values obtained using the full sample size and half of the sample size in the x 2 À x 3 directions, and a good agreement has been obtained for all cases. The results obtained using the full sample size will only be shown in the next section for the sake of conciseness.

Lewis number effects on flame-wall interaction
The distributions of reaction progress variable c, nondimensional temperature T ¼T À T 0 À Á = T ad À T 0 ð Þ , and the nondimensional reaction rate _ ω Â δ Z =ρ 0 S L in the central x 1 À x 3 plane are shown in Figure 1 for turbulent case D for Le ¼ 0:8, 1.0, and 1.2. For the unity Lewis number case, c and T are identical when the flame is away from the wall (e.g., t ¼ δ Z =S L ), which is not valid in the Le ¼ 0:8 and 1.2 cases even when the flame remains away from the wall. The quantities c and T become significantly different from each other in the near-wall region once quenching starts even for Le ¼ 1:0. It can be seen from Figure 1 that local super-adiabatic temperatures (i.e., T > 1:0) are obtained for the nonunity Lewis number flames (i.e., Le ¼ 0:8 and Le ¼ 1:2). The high values of temperature in the Le ¼ 0:8 flames are associated with the flame-wrinkles, which are convex towards the reactants, whereas high temperature values in the Le ¼ 1:2 flames are obtained for flame-wrinkles that are concave towards the reactants. This behavior is consistent with several previous findings Cant, 2005, 2009;Rutland and Trouvé, 1993). A combination of strong focusing of reactants and weak defocusing of heat is responsible for high values of temperature and reaction rate magnitude at the wrinkles, which are convexly curved towards the reactants in the Le ¼ 0:8 cases. Just the opposite mechanism is responsible for high values of temperature in the regions which are concave towards the reactants in the Le ¼ 1:2 cases. It can be seen from the reaction rate distributions in Figure 1 that the chemical reaction rate _ ω drops significantly once the flame reaches in the vicinity of the wall due to a large amount of heat loss through the wall. It has been found that chemical reaction ceases to exist in a region in the x 1 =δ Z < Pe min , where Pe min is the minimum Peclet number, where Pe ¼ X=δ Z is the wall Peclet number with X being the wall normal distance of T ¼ 0:9 isosurafce as defined by Poinsot et al. (1993). The minimum Peclet number for head-on quenching of laminar premixed flames ðPe min Þ L has been found to increase with decreasing Le (e.g., ðPe min Þ L ¼ 3:09, 2.83, and 2.75 for Le ¼ 0:8, 1.0, and 1.2, respectively). The rate of mass diffusion is greater than the thermal diffusion rate in the Le ¼ 0:8 case, and thus the reactants from the vicinity of the wall diffuses faster into the approaching flame than the rate of propagation of isotherms towards the wall. As a result, the minimum Peclet number for the laminar flame ðPe min Þ L in the Le ¼ 0:8 case has been found to be greater than in the unity Lewis number case. The rate of thermal diffusion is greater than the rate of mass diffusion of fresh reactants from the vicinity of the wall in the Le ¼ 1:2 case, and as a result the isotherms can reach closer to the wall before quenching than in the Le ¼ 1:0 case. This, in turn, leads to smaller value of ðPe min Þ L in the Le ¼ 1:2 case than in the Le ¼ 1:0 case. The values obtained for ðPe min Þ L are consistent with several previous experimental (Huang et al., 1986;Jarosinsky, 1986;Vosen et al., 1984) and computational (Poinsot et al., 1993) findings. The minimum value of wall Peclet number Pe min in turbulent flames remains comparable to the corresponding laminar flame value Pe min ð Þ L for Le ¼ 1:0 and 1.2 cases but Pe min in the turbulent Le ¼ 0:8 cases assumes smaller magnitude than the corresponding Pe min ð Þ L . The flame quenching is initiated in the turbulent Le ¼ 0:8 cases when the flame-wrinkles, which are convex towards the wall and associated with super-adiabatic temperatures, approach the wall. The high rate of chemical reaction due to super-adiabatic temperature and relatively weak thermal diffusion rate in turbulent Le ¼ 0:8 cases enable the aforementioned flame-wrinkles to come closer to the wall than in the corresponding laminar flame before quenching. Although super-adiabatic values of temperature are obtained in the turbulent Le ¼ 1:2 cases, these high temperature zones are associated with flame-wrinkles, which are concave towards the reactants and as a result quenching of other regions of the flame (which acts to reduce the temperature) starts to occur before these zones interact with the wall. Thus, the minimum Peclet number for the turbulent Le ¼ 1:2 cases remains comparable to the corresponding laminar values. Interested readers are referred to Lai and Chakraborty (2015) for a more detailed discussion on the influences of Le on the minimum Peclet number Pe min , which is not repeated here for the sake of conciseness.
The temporal evolutions of f c 00 2 andc 1 Àc ð Þ in the direction normal to the wall are shown in Figure 2. According to Eq.
(2), f c 00 2 becomes equal toc 1 Àc ð Þfor a presumed bimodal distribution of c with impulse functions at c = 0 and c = 1.0, which is strictly valid for Da ) 1 flames (Bray et al., 1985). The difference between f c 00 2 andc 1 Àc ð Þ provides a measure of the extent of the deviation of the reaction progress variable pdf P c ð Þ from the presumed bi-modal pdf of c with impulse functions at c = 0 and c = 1.0. Figure 2 shows that f c 00 2 remains smaller thanc 1 Àc ð Þ for all cases even when the flames are away from wall (e.g., t ¼ 2δ Z =S L ) for all Lewis number cases due to small values of Damköhler number (i.e., Da < 1:0). However, later on f c 00 2 drops significantly during flame quenching and eventually vanishes in the regions close to the wall, wherec 1 Àc ð Þ continues to assume nonzero values (i.e.,c 1 Àc ð ÞÞ 0). It has been shown by Lai and Chakraborty (2015) that P c ð Þ does not resemble a bi-modal distribution in the near-wall region x 1 =δ Z Pe min ð Þ L (not shown here for this reason) and thus it is not possible to obtain f c 00 2 by an algebraic closure (e.g., Eq. (2)) and under this situation one has to solve a modeled transport equation in order to evaluate f c 00 2 . Statistical behavior of the variance f c 00 2 transport The variations of T 1c ; T 2c ; T 3c , and ÀD 2c ð Þwith normalized wall normal distance x 1 =δ Z are shown in Figure 3 for all cases considered here. The following observations can be made from the variations of T 1c ; T 2c ; T 3c , and ÀD 2c ð Þ with x 1 =δ Z for all cases considered here: • For all cases the reaction rate term T 3c and the molecular dissipation term ÀD 2c ð Þ remain a leading order source and sink terms, respectively, in the variance f c 00 2 transport equation when the flame is away from the wall. The magnitudes of both the terms decrease with time as flame starts to quench. The terms T 3c and ÀD 2c ð Þ remain the same order of magnitude away from the wall but T 3c vanishes in the region given by x 1 =δ Z <Pe min due to flame quenching, whereas ÀD 2c ð Þ continues to act as a dominant sink term even when T 3c disappears. However, ÀD 2c ð Þ eventually vanishes when the flame is completely quenched. • The mean scalar gradients term T 2c acts as a sink for all cases considered here. It can be seen from Figure 3 that the turbulent scalar flux ρu 00 1 c 00 shows counter-gradient transport (i.e., ρu 00 1 c 00 Â @c=@x 1 >0 and -ρu 00 1 c 00 Â @c=@x 1 <0) throughout the flame brush for all cases considered here.
• The turbulent transport term T 1c shows negative values close to the wall but assumes positive values away from the wall during early stages of flame quenching. However, T 1c assumes positive values in the near-wall region and negative values away from the wall as a result of the reversal of flow direction (after quenching the flow is directed towards the wall in contrast to the flow away from the wall before quenching) at later stages of flame quenching. • One obtains the following scaling estimates of T 1c ; T 2c ; T 3c and ÀD 2c ð Þ according to the scaling arguments of Swaminathan and Bray (2005): where the gas density is scaled using the unburned gas density ρ 0 , the turbulent velocity fluctuations associated with scalar fluctuations are scaled using the unstrained laminar burning velocity S L , the mean gradients are scaled using the turbulence integral length scale l, and the length scale associated with gradient of fluctuating quantities is scaled using the flame thickness δ th . In Eq. (6), _ ω is scaled as _ ω,ρ 0 S L =δ th . It can be seen from Figure 3 that the magnitudes of the turbulent transport and mean scalar gradient terms T 1c and T 2c remain smaller than those of T 3c and ÀD 2c ð Þ, especially when the flame is away from the wall before flame, which is consistent with the scaling estimates presented in Eq. (6). Furthermore, it can be seen from Figure 3 that the magnitudes of T 3c and ÀD 2c ð Þ increase with decreasing Le, which is consistent with previous findings by .
The modeling of the terms T 1c ; T 2c ; T 3c , and ÀD 2c ð Þ will be discussed next in this article.

Modeling of turbulent transport of scalar variances
According to BML (Bray et al., 1985) the joint pdf between velocity vectorũ and reaction progress variable c can be expressed as: where α c ; β c , and γ c are the weights associated with the pdf contributions, P Rũ ; 0 ð Þ and P Pũ ; 1 ð Þ are the conditional velocity pdfs in reactants and products, respectively, and fũ; c ð Þ originates from the interior of the flame. For high Damköhler number flames the third contribution can be ignored and in the case of unity Lewis number flames one gets: Bray et al., 1985). Based on Eq. (7) one gets the following expressions for high Damköhler number (i.e., Da ) 1) flames (Bray et al., 1985): where u i ð Þ R and u i ð Þ P are the ith components of mean velocity conditional on reactants and products, respectively. The last terms on the right-hand side of Eqs. (8a)-(8c) can be ignored for Da ) 1.  demonstrated that ρu 00 i c 00 1 À 2c ð Þ does not adequately predict ρu 00 i c 00 2 obtained from DNS data for low Damköhler number (i.e., Da < 1) combustion and proposed an alternative model as: where m ¼ 0:3 is a model parameter. It is worth to noting that O γ c À Á contribution for low Damköhler number (i.e., Da < 1) combustion is represented by 2 f c 00 2 = f c 00 2 þc Á 1 Àc ð Þ h i in Eq.  Figure 4 shows the variations of ρu 00 1 c 00 2 with normalized wall normal distance x 1 =δ Z as obtained from DNS data along with the predictions of Eq. (9) for all cases considered here. Equation (9) mostly provides satisfactory performance away from the wall but this model underpredicts the magnitude of the negative contribution of ρu 00 1 c 00 2 in the near-wall region when the flame starts to interact with the wall (see Figure 4). Based on this observation, Eq. (9) has been modified in the following manner: þ 2 is the model parameter, which remains active close to the wall wherecÞT but the magnitude of A w increases with increasing wall normal distance and asymptotically approach 1.0 away from the wall wherec %T. It can be seen from Figure 4 that the model given by Eq. (9) starts to underpredict the magnitude of ρu 00 1 c 00 2 at an early stage of flame quenching (e.g., t = 8δ Z =S L for Le ¼ 1:0 and t = 6δ Z =S L for Le ¼ 0:8). Furthermore, Eq.
(9) starts to predict the wrong sign of ρu 00 1 c 00 2 at later stages of flame quenching in Le ¼ 0:8 cases (e.g., t = 6δ Z =S L for Le ¼ 0:8). The sign of ρu 00 1 c 00 2 is incorrectly predicted when ρu 00 1 c 00 2 = ρu 00 1 c 00 1 À 2g mc ð Þ h i becomes negative. In order to avoid this discrepancy A 3 w À 2g mc À Á is introduced in Eq. (10), which assumes a negative value in the near-wall region where 1 À 2g mc ð Þremains positive. The term A w remains active in the near-wall region wherec and T are different from each other as a result of flame quenching. The nonzero value ofc ÀT À Á arises due to different boundary conditions used for the reaction progress variable and nondimensional temperature at the isothermal inert wall (i.e., Dirichlet boundary condition for nondimensional temperature and Neumann boundary condition for reaction progress variable). The ðc ÀT) dependence of A w ensures that the effects of enthalpy loss due to wall heat transfer are reflected on both the qualitative and quantitative variations of ρu 00 1 c 00 2 depending on the distance of the flame from the wall. The quantities,c andT approach each other away from the wall (i.e., x 1 =δ Z ) Pe min ), butcÞT andc >T in the near-wall region during flame quenching. A model parameter similar to A w was previously used in the context of FSD-based closure for flame-wall interaction (Bruneaux et al., 1997). The quantitiesc andT approach each other away from the wall (i.e., x 1 =δ Z ) Pe min ), which leads to A w ¼ 1:0 and thus Eq. (10) reduces to Eq. (9) away from the wall. The wall normal distance at whichc andT approach each other, and the discrepancy between the prediction of Eq. (9) and DNS data depend on Le (e.g., the discrepancy is greater in extent in the Le ¼ 0:8 case than in the Le ¼ 1:0 and 1.2 cases), and thus, A w is taken to be Lewis number dependent. It can be seen from Figure 4 that the underprediction of ρu 00 1 c 00 2 by Eq. (9) in the near-wall region can be eliminated by the modification proposed in Eq. (10).
It is worth noting that the success of the model given by Eqs. (8c), (9), and (10) depend on appropriate modeling of turbulent scalar flux ρu 00 i c 00 . Furthermore, the modeling of ρu 00 i c 00 plays a pivotal role in the evaluation of T 2c . The near-wall modeling of turbulent scalar flux ρu 00 i c 00 will be addressed next in this article.
Algebraic closure of turbulent scalar flux ρu 00 i c 00 Using Eq. (8a) one obtains Cant, 2009, 2015;Malkeson and Chakraborty, 2012): The slip velocity u i ð Þ P À u i ð Þ R n o can be expressed as (Chakraborty and Cant, 2009): where M i ¼ À @c=@x i ð Þ= Ñc j jis the ith component of the flame normal vector based on the Favre averaged reaction progress variable, Δu turb is the contribution to the slip velocity arising from turbulence, and Δu hr is the contribution to the slip velocity arising from heat release. Using Eqs. (11) and (12) one obtains: Using Δu turb ¼ Àα 2 ffiffiffiffiffiffiffiffiffi ffi 2k=3 q (where α 2 is a model parameter andk ¼ ρu 00 j u 00 j =2" ρ is the turbulent kinetic energy) (Chakraborty and Cant, 2009;Veynante et al., 1997) one obtains: The quantity Ñc j j can be scaled as Ñc j j,1=δ b , where δ b is the flame brush thickness. Accordingly, the velocity jump due to heat release over a distance equal to the flame thickness based on reaction progress variable gradient for a corresponding laminar flame (i.e., δ L ,1= Ñc j j L ) can be estimated as: where Ñc j j L is estimated as Ñc j j L ,c 1 Àc ð ÞLe=δ th , in which δ th =Le provides an estimate for the laminar flame thickness based on the reaction progress variable gradient.
According to Veynante et al. (1997), can be expressed as M i , which upon using in Eq. (8b) yields Cant, 2009, 2015): where Chakraborty and Cant, 2015) are the model parameters with Re L ¼ ρ 0k 2 =μ 0ε andε ¼ μ @u  Figure 5 shows the variations of ρu 00 1 c 00 with normalized wall normal distance x 1 =δ Z as obtained from DNS data along with the predictions of Eq. (14) for all cases considered here. It can be seen that ρu 00 1 c 00 is positive throughout the flame brush and gradually reduces zero at the wall. The positive value of ρu 00 1 c 00 is indicative of counter-gradient transport as @c=@x 1 remains positive in the positive x 1 -direction. It can be seen from Figure 5 that Eq. (14) satisfactorily predicts the qualitative behavior of ρu 00 1 c 00 when the flame is away from the wall but this model significantly overpredicts ρu 00 1 c 00 once the flame approaches the wall, and Eq. (14) predicts nonzero values of ρu 00 1 c 00 at the wall. This starts to happen at an earlier time for higher values of u 0 =S L ,Re 1=4 t Ka 1=2 ,Re 1=2 t Da À1=2 because the flame starts to interact with the wall at an earlier time instant due to greater extent of flame wrinkling. In order to eliminate the inadequacies of Eq. (14) in the near-wall region the following modification has been suggested: is the model parameter. Figure 5 shows that Eq. (14) overpredicts the magnitude of ρu 00 1 c 00 in the near-wall region. The presence of a wall leads to a decay in velocity fluctuation, which gives rise to a reduction of the magnitude of scalar flux ρu 00 1 c 00 . However, this behavior is not sufficiently captured by Eq. (14) and it overpredicts the magnitude of ρu 00 1 c 00 . For this reason, Eq. (14) is revised to propose a new model (i.e., Eq. 15) where the parameter A 1 accounts for the reduction of scalar flux magnitude due to the presence of a wall. The model parameter A 1 is responsible for eliminating the overprediction of ρu 00 1 c 00 in the near-wall region. The functional dependence of A 1 on Leðc ÀT) and x 1 =δ Z ensures that this parameter remains active close to the wall wherec ÞT. The turbulent scalar flux components ρu 00 i c 00 vanish at the wall (i.e., x 1 ¼ 0) because the velocity component fluctuations u 00 i vanish at the wall due to no-slip condition. The model parameter A 1 contains an error function that depends on x 1 =δ Z , which ensures that both ρu 00 i c 00 ¼ 0 and A 1 ¼ 0 at x 1 =δ Z ¼ 0. Furthermore, the error function in A 1 ensures that it increases from 0 at x 1 ¼ 0 with increasing x 1 =δ Z and asymptotically approaches 1.0 away from the wall (i.e., Modeling of reaction rate term T 3c According to Bray et al. (1985), the reaction rate contribution T 3c can be expressed as: where c m is given by c m ¼ ò Þ being the burning mode pdf. This parameter c m has been found to be equal to 0:87, 0:85, and 0:83for Le ¼ 0:8, 1.0, and 1.2, respectively. Bray (1980) proposed the following closure for the mean reaction rate of reaction progress variable _ ω in terms of scalar dissipation rate e ε c for Da ) 1 flames based on a presumed bi-modal pdf of c with impulses at c ¼ 0 and c ¼ 1:0: It was shown by Chakraborty and Cant (2011) and , based on scaling arguments and DNS data, that Eq. (17) also remains valid for Da<1 as long as the flamelet assumption remains valid. Thus, the reaction rate term T 3c can be expressed as: (18) Figure 6 shows the variations of T 3c with normalized wall normal distance x 1 =δ Z as obtained from DNS data along with the predictions of Eq. (18) for all cases considered here. It can be seen from Figure 6 that Eq. (18) satisfactorily predicts T 3c when the flame is away from the wall, but once the quenching starts, Eq. (18) predicts nonzero values at the wall and in the near-wall region, where T 3c either vanishes or assumes negligible values. This behavior originates due to nonzero value of 2" ρe ε c = 2c m À 1 ð Þin the near-wall region, where _ ω vanishes due to flame quenching (Lai and Chakraborty, 2015). Recently, Lai and Chakraborty (2015) extended the expression given by Eq. (18) to predict _ ω accurately in the near-wall region in the following manner: The parameters A 2 ; A 3 , A 4 , and B in Eq. (19a) are given by: .
In Eq. (19b),c w andT w are the Favre averaged values of reaction progress variable and nondimensional temperature at the wall, and Å ¼ Pe min ð Þ L erf 8Le À 6 ð Þ=2 is the parameterization for the minimum Peclet number Pe min for turbulent flames (see Lai and Chakraborty (2015) for further discussion on this). Equation (19) combines the advantages of both FSD (i.e., _ ω ¼ ρ 0 S L AE gen , where AE gen ¼ Ñc j j is the generalized FSD (Boger et al., 1998)) and SDR-based (Eq. (17)) mean reaction rate closures. The FSD-based closure is known to overestimate _ ω in the near-wall region during flame quenching (Bruneaux et al., 1997;Lai and Chakraborty, 2015). Moreover, ρ 0 S L AE gen accurately predicts _ ω away from the wall for unity Lewis number flames but ρ 0 S L AE gen underpredicts (overpredicts) _ ω for the Le ¼ 0:8 (Le ¼ 1:2Þ cases when the flame is away from the wall (not shown here but interested readers are referred to Lai and Chakraborty (2015) for further information). By contrast, 2" ρe ε c = 2c m À 1 ð Þpredicts _ ω satisfactorily for all cases irrespective of Le when the flame is away from the wall. In Eq. (19a), the generalized FSD is estimated as AE gen , ffiffiffiffiffiffiffiffiffi ffi e ε c =D q andc w ÀT w À Á dependence of ψ in A 3 ensures that the prediction of Eq. (19a) captures the correct spatial distribution of mean reaction rate _ ω at different stages of flame quenching depending on the values ofc w andT w . The involvement of Å in the model parameters A 2 and A 3 ensures that _ ω vanishes in the region given by x 1 =δ Z ( Pe min , whereas these model parameters asymptotically approach 1.0 away from the wall (i.e., x 1 =δ Z ) Pe min ). Furthermore, the involvement of π implicitly includes reacting boundary layer information into Eq. (19a). For the Le ¼ 1:0 cases,c andT are identical to each other for x 1 =δ Z ) Pe min ð Þ L , which leads to A 2 e LecÀT ð Þ ¼ 0 and A 3 A 4 ¼ 0, and thus, Eq. (19a) reduces to Eq. (17). For Le Þ 1:0 cases,c andT are not equal to each other even when the flame is away from wall, and the involvement of Le on the first term on the right-hand side of Eq. (19a) accounts for this effect. The involvement of 1=Le B in the second term on the right-hand side of Eq. (19a) compensates the underprediction (overprediction) of _ ω by the FSD-based closure for the turbulent Le < 1 (Le > 1) cases. It can be seen from Figure 7 that Eq. (19a) satisfactorily predicts the mean reaction rate _ ω at different stages of flame quenching for all cases considered here. Substituting _ ω from Eq. (19a) in Eq. (16) yields the following model for the reaction rate term T 3c : ffiffiffi ffi e ε c D r e À0:5 It can be seen from Figure 6 that Eq. (20) satisfactorily predicts T 3c both away and close to the wall for all cases considered here. Equations (19a) and (20) indicate that the satisfactory prediction of _ ωand T 3c depends on accurate evaluation of SDR e ε c . Moreover, the closure of D 2c ¼ 2ρe ε c depends on the modeling of e ε c . Thus, the modeling of SDR e ε c will be discussed next in this article.
In Eq. (21), K Ã c is a thermochemical parameter, which is defined as (Kolla et al., 2009): The parameter K Ã c is equal to 0:74τ, 0:78τ, and 0:80τfor Le ¼ 0:8, 1.0, and 1.2, respectively. It is instructive to present the transport equation of SDR e ε c in order to understand the origin of Eq. (21) Chakraborty and Swaminathan, 2010;Gao et al., 2015;Mantel and Borghi, 1994;Mura and Borghi, 2003;Swaminathan and Bray, 2005): for strengthening of the density variation term T 2 and the contribution of turbulence-scalar interaction T 3 arising from the strain rate due to flame normal acceleration as a result of augmented chemical heat release for small values of Le Swaminathan, 2010, 2011). Figure 8 shows the variations of e ε c with normalized wall normal distance x 1 =δ Z as obtained from DNS data along with the predictions of Eq. (21) for all cases considered here. It can be seen from Figure 8 that Eq. (21) significantly overpredicts e ε c in the nearwall region where the equilibrium between T 2 þ T 3 þ T 4 þ f D ð Þ ½ and ÀD 2 ð Þ is not maintained (Lai and Chakraborty, 2015). Thus, Eq. (21) yields an erroneous value of dissipation rate term ( À D 2c Þ and SDR e ε c in the near-wall region (x 1 =δ Z Pe min ð Þ L ). Lai and Chakraborty (2015) modified Eq. (21) in the following manner in order to account for the near-wall behavior: where the model parameters A ¼ 0:5 erf x 1 =δ Z À π ð Þþ1 ½ and exp À1:2Lec w ÀT w À Á 3 only remain active close to the wall to account for the flame-wall interaction and they asymptotyically approach 1.0 away from the wall. The involvement of π includes reacting boundary layer information into the model given by Eq. (24). Furthermore, ðc w ÀT w ) dependence of Eq. (24) accounts for the effects of nonadiabaticity due to wall heat transfer, which influences both the qualitative and quantitative variations of e ε c depending on the distance of the flame from the wall. Figure 8 shows that Eq. (24) predicts e ε c accurately for both near to and away from the wall.

Conclusions
The reaction progress variable variance f c 00 2 transport and its modeling in the context of RANS have been analyzed for head-on quenching of turbulent premixed flame due to an inert isothermal wall using 3D simple chemistry DNS data for global Lewis numbers Le ranging from 0.8 to 1.2. The statistical behaviors of the unclosed terms in the transport equation of f c 00 2 have been analyzed in detail and their relative magnitudes have been explained based on scaling arguments. It has been found that the reaction rate contribution T 3c and the molecular dissipation term ÀD 2c ð Þare the leading order source and sink terms, respectively, in the f c 00 2 transport equation. However, the reaction rate contribution T 3c vanishes in the near-wall region due to flame quenching, whereas ÀD 2c ð Þcontinues to act as a dominant sink. The mean scalar gradient term T 2c acts as the sink term for all cases considered here, since the turbulent scalar flux ρu 00 1 c 00 shows counter-gradient transport in these cases. The turbulent flux of scalar variance ρu 00 1 c 00 2 assumes positive values in the near-wall region but becomes negative away from the wall at early stages of flame quenching but an opposite behavior is observed at the final stage of quenching.
The performances of previously proposed models for turbulent fluxes ρu 00 i c 00 2 and ρu 00 i c 00 , reaction rate contribution T 3c and scalar dissipation rate e ε c have been assessed with respect to the corresponding quantities extracted from DNS data. It has been found that the aforementioned models do not adequately predict the near-wall behavior of the unclosed terms of the variance f c 00 2 transport equation. The existing models for the unclosed terms of the variance f c 00 2 transport Figure 8. Variations ofε þ c ¼ε c Â δ Z =S L obtained from DNS data and the predictions of Eqs. (21) and (24) with x 1 =δ Z at t = 4δ Z =S L , 6δ Z =S L , 8δ Z =S L , 10δ Z =S L for turbulent cases A-E with Le ¼ 0:8, 1.0, and 1.2. Please refer to the table in Figure 2 for the color scheme. equation have been modified to account for the near-wall behavior in such a manner that the modified models asymptotically approach the existing model expressions away from the wall. The functional forms of the modeling parameters have been proposed in such a manner that they follow the asymptotic behavior in terms of normalized wall normal distance x 1 =δ Z . It is, however, likely that they need to be validated further based on both experimental and DNS data for high values of Re t . Furthermore, the proposed models need to be implemented in actual RANS simulations to assess their predictive capabilities. Some of the aforementioned issues will form the basis of future analyses.