NUMERICAL SIMULATION OF A BACKWARD-FACING STEP COMBUSTOR USING REYNOLDS-AVERAGED NAVIERSTOKES / EXTENDED PARTIALLY STIRRED REACTOR MODEL

The authors adapt recently developed a large eddy simulation / extended partially stirred reactor (LES/EPaSR) model by Sabelnikov and Fureby for simulation of turbulent combustion to Reynolds-averaged Navier Stokes (RANS) equations. The proposed RANS/EPaSR model is validated against experimental database created at ONERA for an air methane premixed §ame stabilized by a backward-facing step combustor. The RANS/EPaSR model is compared also with the following RANS-based combustion models: ( i ) quasi-laminar model with reduced chemical mechanism (QL RCM); ( ii ) premixed §amelet tabulated chemistry (PFTC) without taking into account the turbulencechemistry interaction (TCI); and ( iii ) a PFTC with a presumed β probability density function (PDF) for a progress combustion variable.

widespread approach for calculating turbulent reactive §ows is based on the RANS equations [1], despite the well-known limitations of such an approach. The RANS equations are not closed and numerous terms have to be modeled: the Reynolds stress tensor, turbulent scalar §uxes, and mean chemical source terms.
The modeling of the mean chemical source terms needs a careful consideration of the TCI and usually is considered as the main problem of moment methods in turbulent combustion. State-of-the-art in the modeling of turbulent combustion is presented, e. g., in [1,2]. We limit ourselves mentioning the following models: (i) laminar §amelet models (for premixed and nonpremixed combustion): turbulent §ame brush is presented as an ensemble of thin laminar §ames; (ii) RANS/PDF models employing presumed PDF or transported PDF [3]; and (iii) phenomenological ¦nite rate chemistry (FRC) models such as a partially stirred reactor (RANS/PaSR) [46] and an eddy dissipation concept (RANS/EDC) [710] models.
Each of these models has certain features that limit its usability for engineering applications. Therefore, it is desirable to develop more versatile and cost-e¨ective RANS combustion models. In this paper, with that aim, the LES/EPaSR model recently developed by Sabelnikov and Fureby [11,12] has been adapted to the RANS simulation. It is worth to mention that PaSR model was improved as well in [5,13] using the additional transport equations for the reactor species. The concept of the model is based on the experimental results of Batchelor and Townsend [14] showing that at high Reynolds numbers, turbulent ¦ne structure is not uniformly distributed in space but concentrated in small isolated regions, whose volume is a small fraction of the total volume. Chomiak [15] has used these results for qualitative description of the turbulent §ame, assuming that the most intensive molecular mixing, chemical reactions, and heat release take place mainly in the ¦ne structures. Later, this model has been re¦ned in [710]. The spatial no uniformity of the turbulent ¦ne structures has been veri¦ed by a direct numerical simulation in nonreactive [16] and reactive [17,18] §ows. The organization of the paper is as follows: the RANS/EPaSR model is presented in section 2. In section 3, the simpli¦ed assumptions introduced for that the model to be adapted to the structure of the industrial code CEDRE developed in ONERA [19]. The validation of the RANS/EPaSR model on the experimental database created at ONERA for an airmethane premixed §ame stabilized by a backward-facing step combustor is presented in sections 4 and 5.

Governing Equations
The EPaSR model, similar to the PaSR [46] and the EDC [710] models, is based on the assumption that at high turbulent Reynolds numbers, the chemical reactions take place essentially in the ¦ne turbulent structures. These structures are not uniformly distributed in space and its volume is a fraction of the total gas volume. This implies that the ¦ne-structure regions are embedded in a surrounding §uid. The most intensive molecular mixing, chemical reactions, and heat release take place mainly in the ¦ne structures. Henceforth, the quantities inside of the ¦ne structures and in a surrounding §uid will be denoted by * and 0 , respectively. The principal di¨erence between the EPaSR and the PasR models is that the former is formulated using di¨erential equations for the evolution of the quantities in the ¦ne structures and surroundings and the last is based on the algebraic equations. In other words, contrary to the PaSR model, the EPaSR model takes into account the history e¨ects in the turbulent §ow.
The derivation of the LES/EPaSR model [11,12] makes use of the similarities with the mathematical treatment of multiphase §ows. The reader is referred to [11,12] for the details. The RANS/EPaSR model is a direct modi¦cation of the LES/EPaSR model [11,12]. The equations of the RANS/EPaSR model follow directly from the LES/EPaSR after change of characteristic turbulent time and length scales. Let us denote the composition space Here, Y = [Y 1 , . . . , Y Nsp ] are the mass fractions of N sp species and h s is the sensible enthalpy: where C p,i is the heat capacity at constant pressure of species i and T is the temperature. The evolution of ψ is governed by the local balance equations of mass fractions and total energy [1]. The Favre averages, ¦ne structures, and surroundings quantities are related by the following expression: where γ * is the ¦ne-structure volume fraction; and ψ * and ψ 0 are the quantities inside of the ¦ne structures and surroundings, respectively. Also, let us denote " f = ρf /ρ to be the Favre (mass-weighted) average of f , overline designates Reynolds average. The transport equations for ψ * and ψ 0 used in LES/EPaSR model are given in [11,12] where it is also noted that more convenient and straightforward approach is, however, to solve the balance equations for ψ * together with the LES balance equations for " ψ. This approach will be used for RANS/EPaSR model as well. The set of balance equations for ψ * together with the RANS balance equation for " ψ reads: where p is the pressure; τ is the viscous stress tensor; b if the Reynolds stress tensor; and vectors k i , i = 1, . . . , N sp , and k Nsp+1 are the molecular species mass and heat §uxes, respectively. Vectors b i and b Nsp+1 are the turbulent species mass and heat §uxes, respectively. The terms b and b i are modeled with an eddy viscosity and di¨usivity following Boussinesq approach. The exchange terms M * i and ' m in Eqs. (5) and (6) describe the exchange between the ¦ne structures and the surroundings (through immaterial surface separating the ¦ne structures and surroundings). The terms M * i and ' m must be modeled. Neglecting chemical reactions in the surroundings, the mean reaction rates for the Favre-averaged species mass fractions and for the source for the Favreaveraged sensible enthalpy read: where h 0 f,i is the mass enthalpy of formation of species i.
Next, let us consider brie §y modeling the exchange terms M * i and ' m, referring the interested reader to [11,12]. As follows from the core physical considerations, M * i should contain two kinds of terms. The former term, here denoted by - * i , is due to the exchange rate of mass between the ¦ne structures and the surroundings. If the mass exchange rate annuls, i. e., ' m = 0 (the dynamic equilibrium state), - * i annuls as well. The second term, here denoted by Ÿ * i , is a micromixing term. It is similar to the micromixing term in the transport equation for species PDF and in Lagrangian stochastic models of turbulent §ows [3]. The term Ÿ * i is due to molecular di¨usion through the interface between the ¦ne structures and surroundings. Indeed, even if the mass exchange rate is absent, there is exchange through the interface only due to molecular di¨usion. If ψ * = ψ 0 , Ÿ * i , however, turns to zero as in this case, there is no molecular di¨usion. Based on the above considerations, the model for the exchange terms reads: where τ * stands for the micromixing time (or the ¦ne structure residence time); γ * eq is the equilibrium ¦ne-structure volume fraction; ω = " ε/ " k is the turbulent frequency with " ε being the turbulent dissipation rate and " k being the kinetic energy of the turbulence; and C ω is the tuned constant, its value will be selected at the end of this subsection.
It can be concluded from Eq. (10) that if ' m > 0 and γ * < γ * eq , the exchange rate of mass is driven from the surroundings to the ¦ne structures. Otherwise, if ' m < 0 and γ * > γ * eq , then the exchange rate of mass is directed from the ¦ne structures to the surroundings. It can be emphasized that 0 < γ * eq < 1 and γ * tends to its equilibrium state with the characteristic time τ * . Equation (6) for the ¦ne-structure volume fraction γ * can be obtained summing ¦rst N sp Eqs. (5) and taking into account the following identities: To close Eqs. (2)(6), it is necessary also to provide closure models the equilibrium ¦ne-structure volume fraction, γ * eq and the ¦ne-structure residence time, τ * . The LES/EPaSR model [11,12] uses the following submodel for τ * : Here, τ K = (v/" ε) 1/2 and τ -= -/u ′ are the Kolmogorov time scale and the subgrid velocity stretch time, respectively, where -and u ′ are the ¦lter width and the subgrid velocity §uctuations, respectively. Similarly to Eq. (11), the micromixing time for the RANS/EPaSR can be de¦ned as the geometric mean of the Kolmogorov time τ K and the integral time τ t = " k/" ε = ω −1 , that is (see also [5]), The equilibrium reacting ¦ne-structure fraction is taken from [11,12] without modi¦cation: where τ ch is the characteristic chemical time scale. For the premixed combustion, τ ch is taken as τ ch ≈ δ L /S L where δ L and S L are the laminar §ame thickness and the §ame speed, respectively. In addition to the de¦nition of the micromixing time τ * by geometric mean, Eq. (12), another choice was tested. Namely, τ * was taken proportional to the Kolmogorov time scale, i. e., where C is the tuned constant. If C = 1, Eq. (14) is just the Kolmogorov time scale and if C = 0.41, one gets the ¦ne structure residence time introduced in [710]. It is worth noting that the micromixing time calculated as the geometric mean, Eq. (12), is larger than calculated from Eq. (14). Therefore, γ * eq , as follows from Eq. (13), is larger with use of Eq. (14).
The RANS/EPaSR model, Eqs. (2)(6), can be simpli¦ed (similar to LES/ EPaSR model [11,12]) signi¦cantly by assuming that the unsteady, convective (left-hand side of Eq. (5)) and §ux terms (¦rst term in the right-hand side of Eq. (5)) are neglected which, in turn, implies that Therefore in such a limit, Eq. (5) of the RANS/EPaSR model reduces to algebraic RANS/PaSR model: ' It should be noted that Eq. (14) coincides as well with EDC model [710]. Based on this observation, the constant C ω in Eq. (9) is set to 10.5, close to the recommendation of [710].

Reynolds-Averaged NavierStokes / Extended Partially Stirred
Reactor Adaption to CEDRE The structure of CEDRE software (CFD package for diphasic reactive §ows) [20] does not allow implementing directly presented in the subsection 2.1 EPaSR model for the Favre-averaged sensible enthalpy and variables for the ¦ne structures and, therefore, requires for the original EPaSR model to be modi¦ed. In the framework of CEDRE, the mean composition space is de¦ned as [21] " i. e., instead of Eq. (5) for the mean sensible enthalpy, the evolution equation for Favre-averaged total energy is used. The composition space in the ¦ne-structure regions is given by t is the total enthalpy (which includes the mass enthalpy of formation of species). The Favre-averaged total enthalpy " h t and the total energy " e t are related by The temporal variation of the mean pressure in Eq. (7) is neglected, i. e., it is assumed that The Reynolds stress tensor b, the turbulent §uxes b i , i = 1, . . . , N sp , and b Nsp+1 in Eqs. (3) and (4) are closed with the Boussinesq approximation (see, e. g., [1]). For simplicity, it is assumed that the molecular di¨usion coe©cients D i for species i are the same and similar assumption is done for turbulent di¨usion coe©cients D ti , i. e., Here, v is the kinematic viscosity; v t is the turbulent viscosity; and Sc i is the Schmidt number of the species i. It is assumed that the turbulent Schmidt and Prandtl numbers are equal, i. e., Using the Boussinesq hypothesis (relations (16)(18)), the balance equation for the Favre-averaged " ψ i and its counterpart inside of the ¦ne structures of species i read: It is worth noting that the source term ' ω * Nsp+1 in the equation for h * t is equal to zero due to assumption (16).

PREMIXED FLAME STABILIZED BY A BACKWARD-FACING STEP
The proposed RANS/EPaSR model is validated against experimental database created at ONERA for an airmethane premixed §ame stabilized by a backwardfacing step combustor data [22,23]. The RANS/EPaSR model is compared with a few contemporary RANS combustion models such as QL RCM, PFTC no TCI, and PFTC β-PDF. It is worth nothing that the EPaSR model presented above, similar to EDC model, oversimpli¦es the §amewall interaction phenomena. In particular, the special treatment of the micromixing term Ÿ * i in Eq. (8) may be required close to the walls. Especially, it is important for water-cooled combustor (as considered in this section) to address to heat §uxes. This important problem will be considered in future.

Domain, grid, and numerical method
The geometry of the §ow problem was chosen in accordance to the experimental setup of ONERA [19,22]     of the reference frame X = Y = 0 is ¦xed at the nose of the step. Figure 2 represents schematically the main features of the §ow. The streamwise position of the reattachment point is denoted by X r . The two-dimensional (2D) computations were carried out. Five computational grids composed of rectangular cells have been generated (Table 1).
These ¦ve grids have the same topology, being re¦ned in the X and Y directions at the corner of the step and in the recirculation region. Figure 3 shows a close-up view of the step extremity of a 2D grid having 3.7 · 10 4 rectangular cells. In this study, only the results from the intermediate grid are presented as this grid has been found to result in predictions virtually identical to that of the ¦ne grids, whereas the coarse grid shows the results that deviate slightly from the results obtained on the ¦ner grids.

Boundary conditions
At the inlet, methaneair mixture is injected at the uniform temperature " T = 525 K with equivalence ratio equal to 0.8. Two types of inlet conditions for the mean streamwise velocity " u x , the turbulent energy k, and the turbulent length scale l were speci¦ed: These pro¦les were obtained by numerical simulation of the channel of height L y = 0.065 m and length L channel x = 30L y = 1.95 m. The grid was uniform consisting of 3 · 10 2 × 3 · 10 2 = 9 · 10 4 cells. The calculation of this channel was performed with inlet velocity " u x = 53 m/s, which yielded maximum outlet velocity of approximately 58 m/s. The standard wall function method was used in the near-wall region. It is assumed that at the inlet, the ¦ne-structure volume fraction is equal to its equilibrium value, i. e., γ * = γ * eq , Eq. (13). The §ux of the ¦ne-structure volume fraction at the combustor walls was set equal to zero. The channel walls are assumed to be adiabatic. At the outlet, the pressure is set to p = 1 bar. The quenching of the §ame at the combustor walls was modeled letting the chemical source terms to be zero inside the ¦rst cell (0.0005 m from the wall). It should be noted that in the experiment, the combustor walls were cooled by water. The wall temperature was not measured. Assumption of adiabatic walls overestimates the gas temperature near the walls. Figure 4 shows vertical pro¦les of " u x , the root-mean-square velocity §uctu-

Turbulence models
The following models are considered for nonreactive §ow case: kl with Boussinesq viscosity model (BVM), kω with BVM, and kω with explicit algebraic Reynolds stress model (EARSM). Turbulent Prandtl and Schmidt numbers are considered to be constant: Pr t = Sc t = 0.9. The kl turbulence model with BVM was employed for the reactive case. The laminar thermal conductivity and mixture viscosity are calculated with the aid of the Eucken£s and the Suther-land£s laws, respectively. The Schmidt numbers for all the species are assumed to be constants and equal to unity. Detailed description of the physical models can be found in CEDRE documentation [20].

NONREACTIVE CASE
The nonreactive backward-facing step §ow calculations were performed to initialize the mean velocity, temperature, and turbulent characteristics ¦elds, and to test the turbulence models.
The backward Euler temporal integration scheme along with the second-order monotonic upstream scheme for conservation laws scheme in space is applied. The time step is ¦xed to -t = 0.1 ms.
In Fig. 5, streamline patterns are displayed, showing the size of the recirculation region behind the step for di¨erent turbulence models. It is seen that the mean reattachment of a separated §ow behind the step numerical results depends on the used turbulence model. Table 2 presents the comparison of the calculated mean reattachment lengths X r /h obtained with di¨erent turbulence models and the experimental data [22,23] (the error of estimation of X r is about ±0.2h). One can conclude that the kl turbulence model gives the best agreement with the experimental data. Mean velocity ¦elds have been compared with the available experimental measurements (not given here). It is found that the kl (BVM) turbulence model yields a more accurate estimation of the velocity ¦elds and the recirculation region length than the kω turbulence model (BVM or EARSM). The nonhomogeneous inlet pro¦les improve agreement between calculations and experimental pro¦les in the upper part of the chamber.

REACTIVE CASE
The nonuniform boundary conditions, depicted in Fig. 4, were used in the reactive case. The explicit second-order RungeKutta temporal integration scheme along with the ¦rst-order ODFI scheme (the Riemann invariants-based §ow decentering operator scheme for hyperbolic equations) in space are applied. The time step is ¦xed to -t = 0.4 µs.

Reaction Mechanism Model
Westbrook and Dryer two-step reaction mechanism of the methane combustion [24], where the oxidation of CO to CO 2 is reversible, was used for the QL RCM and EPaSR combustion models: The reverse reaction for CO 2 decomposition is used in order to reproduce proper heat release and pressure dependence of the [CO]/[CO 2 ] equilibrium. It should be noted that the two-step reaction o¨ers only an approximate description of methane oxidation and, in particular, overestimates the Favre-averaged temperature " T and underpredicts the ignition delay. The temperature overestimation at the equivalence ratio ϕ = 0.8 is relatively small (about 2%) and can be neglected. The underprediction of the ignition delay does not play in the considered simulation important role because of the use of adiabatic walls conditions. This results in the adiabatic §ame temperature, T ad ≈ 2150 K of the combustion products in the recirculation zone behind the step. As it can be seen from Figs. 6 and 7, T ad is much larger than the temperature in the recirculation zone of the combustor with water cooled walls (1500 < " T < 1750 K at X < 0.04 m and 1900 < " T < 2000 K at X = 0.1 m). In the simulation, the §ame appears immediately downstream of the step edge in the mixing layer formed between incoming reactants §ow and hot combustion products in the recirculation zone. Therefore, all turbulent combustion models tested in the present work yield more fast axial temperature increase with respect to experimental data. The implementation of a detailed kinetic mechanism in a CFD code is still prohibitive due to the high associated computational cost. In order to partially overcome this di©culty, the §amelet tabulated chemistry (FTC) approaches were also used in two versions: (i) tabulated chemistry model (PFTC); and (ii) PFTC without TCI (no TCI).
It should be stressed that the tables are created at the base of a free laminar premixed §ame without heat losses. For better predictions of the combustion, in addition to the boundary conditions taking into account the thermal §uxes at the combustor wall, new tables must be generated introducing an additional entry parameter linked with heat losses.
The FTC approaches allow performing a su©ciently accurate description of the main thermochemical phenomena with a relatively low computational cost of methane combustion modeling. The precalculated tables describe approximately the detailed chemistry by a reduced number of the transported mass fractions. The species Y CH 4 , Y O2 , Y N2 , Y CO , Y CO 2 , and Y H2O are used in CEDRE. The methaneair reaction mechanism is taken from GRI-Mech 3.0 [25] (53 species, 325 reactions). The quasi-laminar table is done with two entry parameters: mixture fraction " Z and progress variable " C. In order to generate laminar premixed FTC table, one-dimensional laminar unstretched premixed §ames are ¦rst solved in physical space, for di¨erent values " C, using CANTERA solver [26]:

The quasi-laminar table has 200 points in "
C. The turbulent table has the same size as quasi-laminar but in addition, it possess 25 points in the variance progress variable space C ′′2 uniformly distributed between 0 and " C(1 − " C). Afterwards, β-PDF integration is performed.
Closing this section, it is worth noting that the EPaSR model solves the equations for " ψ and ψ * . The additional computational cost is of about 60% of the EPaSR model compared with the QL model. T at di¨erent locations X downstream of step calculated using the EPaSR (see Fig. 6) and QL RCM and PFTC (no TCI) models (see Fig. 7) with experimental data [22,23]. In the case of the EPaSR model, the calculations were performed with the di¨erent de¦nitions of the micromixing time τ * :  The calculations were performed for three values of the constant C ω in Eq. (9): 5, 10.5, and 15. It is seen that the calculated mean temperatures in the region nearby the bottom wall of the chamber (Y < 0) are overestimated due to using adiabatic boundary conditions. This overestimation does not depend on the combustion model. It is worth noting that for the EPasR model, the temperature overestimation can be partially explained by the underprediction of the ignition delay length.

Mean temperature ¦eld
The transverse gradient of the mean temperature ∂ " T /∂Y is very large for X < 0.25 m, independently on models. The calculated temperatures in the region X < 0.25 m exceed the experimental data. As was already discussed in subsection 5.1, it is explained by the use of adiabatic wall condition. For larger distances 0.25 ≤ X < 0.71, large temperature gradients are preserved for the QL RCM, the EPaSR with micromixing time scale de¦ned by Kolmogorov, and Magnussen time scales, Eq. (14). The comparison of the EPaSR and the §amelets models shows that the temperature obtained with the former model is closer to the experimental data at X = 0.46 m.
The EPaSR model with the geometric micromixing time scale τ * = √ τ K τ t and C ω = 10.5 yields the best agreement with experimental data. The variation of the constant C ω for the geometric mixing time shows the following impact on the calculated temperature. Decreasing C ω from 10.5 to 5 results in the temperature decrease owing to the reduction of micromixing between fresh mixture and combustion products. On the contrary, increasing C ω from 10.5 to 15 results in the temperature increase because of augmentation of the micromixing between fresh mixture and combustion products. Figure 8 illustrates the dependence on the de¦nition of the micromixing time τ * on the pro¦les of the volume fraction of the ¦ne structures γ * and its equilibrium value γ * eq . the lowest γ * and γ * eq , in correspondence with the remark given at the end of subsection 2.1, are obtained for the geometric micromixing time, Eq. (12). The Magnussen de¦nition of the micromixing time (C = 0.41 in Eq. (14)) yields the values of γ * and γ * eq approximately 1.3 times larger with respect to γ * and γ * eq calculated with the Kolmogorov scale time (C = 1 in Eq. (14)). In the bottom region just downstream of the step, the values of γ * and γ * eq calculated with three di¨erent de¦nitions of the micromixing time are lower than 0.5. Then, γ * and γ * eq increase downstream, at X = 0.71 m, nearby the bottom wall γ * eq = γ * = 0.82 for the geometric mean micromixing time and close to γ * eq = γ * = 0.92 for the Magnussen and the Kolmogorov micromixing times. In the middle of the channel, γ * is slightly larger than its equilibrium Figure 8 Pro¦les of volume fraction γ * of the ¦ne structure (curves) and its equilibrium state γ * eq in reactive backward-facing step §ow (signs); RANS/EPaSR; Cω = 10.5: 1 ¡ Kolmogorov micromixing time; 2 ¡ Magnussen mixromixing time; and 3 ¡ geometric micromixing time value γ * eq for the geometric mean micromixing time √ τ K τ t . The better correspondence of the numerical predictions with experimental data obtained using the geometric mean micromixing time √ τ K τ t , Eq. (12), can be explained by the lower values of γ * with respect to the de¦nition by Eq. (14) within the mixing layer where the chemical reactions mainly take place. It permits to smooth the temperature pro¦le, i. e., to decrease the transverse temperature gradient ∂ " T /∂Y . This e¨ect is illustrated in Fig. 9. Figure 9 compares the pro¦les of the mean temperature with the temperature in the ¦ne structures. It is seen that the pro¦les of T * are much steeper with respect to " T , because T * > " T in accordance with the relation (see Eq. (1)): Outside of the mixing layer, the impact of submodel for micromixing time weakens because the mean chemical source terms are relatively small. In addition to the calculations with the EPaSR, PFTC no TCI, and QL RCM models, the calculations with the PFTC β-PDF models were performed as well. Di¨erent closure relations were applied to dissipation and source terms in the transport equation of the combustion progress variable variance C ′′2 . The best results were obtained with PFTC BML (BrayMossLibby) model. In this model, the BML closure for the dissipation term in which a wrinkling factor is approximated by constant 5 and the source in the transport equation of the variance of progress variable is found directly from the table. In general, the PFTC BML overestimates the values of temperature in the combustor.
The comparison of EPaSR, QL RCM, and PFTC BML models is presented in Fig. 10. This ¦gure shows the isolines of temperature " T = 1500 K. The PFTC BML and QL RCM models overestimate the §ame angle in the entire domain except the region where the §ame impacts the upper wall. The §ame angle calculated with the EPaSR model is close to the experimental data in the interval 0 < X < 0.4 m. It is overestimated in the interval 0.4 < X < 0.7 m and underestimated in the interval 0.7 < X < 0.8 m. The divergence between the EPaSR model and experimental data nearby the combustor top wall can be explained by the use of the simple §amewall interaction model ( §ame quenching) and by the behavior of γ * nearby the wall. It is seen from Fig. 8 that the di¨erence 1 − γ * is small (because the ratio of two time scales τ * /τ ch in Eq. (13) decreases approaching to the walls) and, as a result, the EPaSR model reduces to the QL approach. The correspondence with the experimental data can be improved by elaborating micromixing model close to the wall. Figure 11 shows the temperature " T ¦eld in the combustion chamber obtained with the EPaSR model with the geometric mean micromixing time, Eq. (12). The thickness and position of the §ame brush are in correspondence with Figs. 6 and 10. The possible reason of more gradual §ame development with respect to the experimental data nearby the combustor top wall in the region 0.7 < X < 0.8 m is, as was already noted, the oversimpli¦ed §ame quenching model.

Mean velocity ¦elds
The comparison of the calculated pro¦les of the mean velocity " u x with the experimental data at di¨erent X cross sections is presented in Fig. 12 for the EPaSR model with di¨erent speci¦cations of the micromixing time τ * : (i) geometric mean, Eq. (12), for three values of C ω : 5, 10.5, and 15 in Eq. (9) for the micromixing term; and (ii) Kolmogorov and Magnussen τ * , Eq. (14), for one value of C ω = 10.5. The best correspondence with the experimental data is observed with the geometric mean micromixing time. It was found that the variation of the C ω constant has a minor importance on the mean velocity ¦eld. Figure 13 presents the comparison of the calculated pro¦les of " u x using QL RCM and PFTC no TCI with the experimental data. From Figs. 12 and 13, one may conclude that the di¨erence between the EPaSR and PFTC no TCI models appears for X ≥ 0.25 m. The EPaSR model with the geometric mean micromixing time yields the better results. The pro¦les of " u x for the QL RCM are very similar to the EPaSR model with Kolmogorov and Magnussen micromixing times τ * . All these models overestimate the Favre-averaged streamwise velocity. The comparison of the calculated pro¦les of the transverse mean velocity " u y , with the experimental data at di¨erent X cross sections is presented in Figs. 14 and 15.
One may observe essential di¨erence between calculated and experimental pro¦les close to the step: the calculated velocity is negative at Y > 0 while in the experiment, it is positive. In the middle part of the combustion chamber, the calculated velocity " u y with the geometric mean micromixing time τ * has a satisfactory agreement with the experimental data. The peaks with minimal values of the transverse velocity are shifted to the top wall. For X ≥ 0.25 m, the absolute values of " u y are overestimated. The model with the Kolmogorov micromixing time better approximates the experimental results than the model with Magnussen micromixing time. The variation of C ω has a minor importance on the mean transverse velocity.
It should be noted that in the experiments, the ration between the mean reattachment length X r and the size of the step h, noted as X r /h, is located within 2.85 and 3.42. It follows from Figs. 1215 that X r /h is shorter than the experimental results. In particular, for all considered models, X r /h < 2.7. In the case of the EPaSR model, the relative error of the numerical results with the experimental value is about 8%23% that is the smallest value among all the considered models. In addition, on the considered grid, none of the models captures the secondary recirculation region. This fact can explain the di¨erence Figure 15 Pro¦les of Favre-averaged transverse velocity " uy in reactive backwardfacing step §ow with homogeneous inlet pro¦les; RANS/quasi-laminar approaches: 1 ¡ experiments; 2 ¡ PFTC no TCI; and 3 ¡ QL RCM (to be continued) uy in reactive backward-facing step §ow with homogeneous inlet pro¦les; RANS/quasi-laminar approaches: 1 ¡ experiments; 2 ¡ PFTC no TCI; and 3 ¡ QL RCM between the numerical and experimental transverse mean velocity " u y near the step.

CONCLUDING REMARKS
The present study describes and evaluates a new RANS ¦nite rate chemistry combustion model ¡ the RANS/EPaSR model. This model is the adaptation of the LES/EPaSR model, recently developed in [11,12]. The proposed RANS/EPaSR model is validated against experimental database created at ONERA for an air methane premixed §ame stabilized by a backward-facing step combustor. The RANS/EPaSR model is compared also with the following RANS-based combustion models: (i) QL RCM; (ii) PFTC without taking into account TCI; and (iii) a PFTC with a presumed (β-function) PDF for progress combustion variable.
For this combustor, the overall best performing model is the RANS/EPaSR model.