Subgrid-Scale Modeling of Reaction-Diffusion and Scalar Transport in Turbulent Premixed Flames

ABSTRACT A numerical study of premixed flame-turbulence interaction is performed to investigate the effects of turbulence on the structural features of the flame and the subgrid-scale (SGS) effects on vorticity dynamics, energy transfer mechanism, and turbulent transport across the flame. We consider a freely propagating methane-air turbulent premixed flame interacting with a decaying isotropic turbulence under three different initial conditions corresponding to the corrugated flamelet (CF), the thin reaction zone (TRZ), and the broken/distributed reaction zone (B/DRZ) regimes. We employ the well-established linear eddy mixing (LEM) model in large-eddy simulation (LEMLES), a new subgrid closure for reaction-diffusion occurring in the small-scales based on LEM (RRLES), and a quasi-laminar chemistry based closure in large-eddy simulation (QLLES) to simulate flame-turbulence interactions and to compare predictions with direct numerical simulation (DNS). We assess the accuracy and robustness of the closures by comparing statistical features to highlight their abilities and limitations. The newly proposed RRLES subgrid closure uses a dual-resolution grid for solving the species transport equations. Such an approach is shown to improve the existing LEMLES subgrid model especially at low Reynolds numbers. All SGS closures reveals good agreement, although there are some differences due to the closure used for convective transport of the scalar field and the reaction rate. Further analysis of the DNS dataset shows that there is a significant contribution by dilatation and baroclinic torque terms across the flame. In particular, at higher Karlovitz number, there is an abrupt change in the sign of the dilatation term, which is related to the competing effects of thermal expansion due to heat release and enhanced molecular mixing by the turbulence across the flame brush region. The enhanced mixing leads to localized pockets of cold reactants surrounded by hot products, which is only partly captured by the employed closures. The analysis of SGS kinetic energy and scalar dissipation rates indicates the presence of back-scatter of turbulent kinetic energy, and we also observe counter-gradient transport across the flame. The results suggest that further improvement of the traditional closures is needed to accurately capture the dynamics of flame-turbulence interaction.


Introduction
Turbulent premixed combustion occurs in several practical devices such as internal combustion engines, gas turbines, and swirl combustors. Premixed flames under a turbulent environment experience a modification to their structure due to stretching and wrinkling induced by the turbulent eddies. In particular, at moderate to high turbulence intensity, the preheat zone of the flame front shows broadening, whereas at higher level of turbulence intensity, the reaction zone may get disrupted leading to a local/global extinction (Dunn, et al. 2009;. Various types of burning have been categorized in the past using velocity-(u 0 =S L ) and length-scale (l=) ratios, and also using non-dimensional numbers such as Reynolds number (Re), Karlovitz number (Ka), and Damköhler number (Da). Here, u 0 is the turbulence intensity, S L is the unstretched laminar flame speed, l is the integral length scale, and is the Zeldovich flame thickness. Often called the Borghi diagram (Borghi, 1985), the exact boundaries between corrugated flamelet (CF), thin-reaction-zone (TRZ), and broken/distributed reaction zone (B/DRZ) regimes are not well defined (Driscoll, 2008) but characterization in terms of these parameters provides some guidance for defining modeling requirements. For example, a higher value of Ka is typically associated with intense turbulence, where smaller turbulent eddies can penetrate within the flame region and disrupt its structure from the classical flamelet structure. These eddies enhance heat and mass transport away from the reaction layer toward the reactants through the process of turbulent micro-mixing (Sankaran et al., 2007). However, some recent experimental studies (Skiba et al., 2015) point to the possibility of local flamelet-like structure within the TRZ or even in the B/DRZ regime, pointing to the inherent multi-scale nature of these turbulent flames. Subgrid models using large-eddy simulation (LES) have to deal with these different physics in a consistent manner.
Direct numerical simulation (DNS) (Aspden et al., 2011;Chakraborty et al., 2009;Nikolaou and Swaminathan, 2015;Poludnenko and Oran, 2011;Sankaran et al., 2007;Shim et al., 2011;Tanahashi et al., 2000;Trouve and Poinsot, 1994;Veynante et al., 1997) of premixed flame-turbulence interactions is typically relegated to canonical configurations where the focus is to provide insight to the physics of such interactions. Often, such studies consider further simplifications such as low-to-moderate Re flows, lower values of Ka, simplified kinetics, two-dimensional geometry, etc. LES, on the other hand, offers potential for application to much higher Re flows (and to practical devices), but validated subgrid scale models that work on all regimes without ad hoc tuning still remain to be established. SGS models are needed not only for scalar processes (reaction, diffusion, mixing) but also for momentum and energy subgrid fluxes. Typically, classical eddy viscosity models used for momentum and energy subgrid flux closures (Schumann, 1975;Smagorinsky, 1963) have been extended using an eddy diffusivity concept for the scalar subgrid flux (Rose, 1977) assuming that a gradient closure is still applicable to subgrid scalar flux. However, since reaction kinetics, molecular diffusion, and subgrid turbulent mixing interactions occur at the subgrid scales (that are not resolved on a typical LES grid), there is significant anisotropy at the small-scales; therefore, a gradient transport closure may not be always appropriate. Refining the LES grid to resolve the flame thickness or the reaction zone would lead to a DNS grid requirement, which is not practical. Although flame front-based geometrical models for flamelet regimes have proven quite successful, they tend to be inadequate when small-scale turbulence can penetrate the flame structure (Dunn et al., 2007). Sometimes, finite-rate kinetics have to be simulated (e.g., to investigate extinction-ignition physics, lean blow out, and other extreme transient events) and then SGS closure for the reaction rate is explicitly needed. Many closures that address reaction kinetics have been developed, for example (citations are representative and not included to be comprehensive): partially stirred reactor (PaSR) (Fureby, 2009), thickened flame model (TFM) (Colin et al., 2000), flame surface density (FSD) (Boger etal., 1998), conditional moment closure (CMC) (Klimenko and and Bilger, 1999), conditional source estimation (CSE) (Steiner and Bushe, 2001), transported probability density function (TPDF) (Colucci et al., 1998;Gao and O'Brien, 1993), multienvironment PDF (MEPDF) (Fox and Raman, 2004), one-dimensional turbulence (ODT) (Echekki et al., 2001), and linear-eddy mixing (LEM) model (Menon and Kerstein, 2011;Menon et al., 1993), etc. The present article focuses primarily on the LES strategy that allows direct inclusion of finite-rate kinetics and in particular on the LEM based subgrid models.
LES subgrid models also have to capture the key features of scale interactions seen in DNS. For example, a priori analysis of the SGS transport and vorticity dynamics using DNS of non-premixed mixing layer (O'Brien et al., 2014a) and planar premixed flame (O'Brien et al., 2014b) showed that negative "dissipation" and counter-gradient transport are present in the flame zone. Analysis of enstrophy balance across the flame has been used to assess SGS modeling (Bobbitt and Blanquart, 2015;Lipatnikov et al., 2014), since enstrophy represents the kinetic energy of small scale turbulence and these scales have to be modeled in LES. Particularly, in a planar premixed flame configuration, DNS (Bobbitt and Blanquart, 2015;Louch and Bray, 2001) suggests that (in a statistical sense) pressure and density gradient induced by the heat-release generates baroclinic torque, which leads to a large contribution in the total balance of the enstrophy. SGS models that are currently employed do not account for these effects explicitly and further study is necessary.
To address some of these observations and also to assess some of the existing SGS modeling strategies for such flows, we perform DNS and LES [using two different multiscale implementation of LEM as a SGS model and one based on quasi-laminar (QL) kinetics within a conventional LES (QLLES) (Grinstein and Kailasanath, 1995)] of a freely propagating planar flame in decaying isotropic turbulence. The two multi-scale implementations of LES as a SGS model include the well established LEM at the subgrid in LES (LEMLES) (Menon and Kerstein, 2011) and a novel reaction-rate closure in LES (RRLES). The strengths and limitations of LEMLES are well known (Menon and Kerstein, 2011) and RRLES is devised to address some of those limitations. For example, the large-scale molecular diffusion across the LES cell faces is ignored in the LEMLES formulation, which leads to inaccurate results when the turbulence is resolved at the LES level as in a DNS or if the turbulence is negligible (laminar conditions). In the RRLES formulation, we use subgrid LEM equations to obtain the filtered reaction rate; however, the large-scale molecular diffusion is included in a manner similar to the conventional LES where typically a closure for the filtered reaction-rate term is used. Therefore, the formulation can asymptote to DNS or laminar conditions. Another key advantage of the RRLES formulation comes from the use of a dual-grid approach (discussed below), where LEM subgrid fields evolve on a coarser LES grid, which has more unresolved scales leading to an improved prediction by solving the subgrid LEM equations. A major limitation of the RRLES formulation is that it can not capture counter-gradient turbulent transport-a unique feature of the LEMLES formulation due to the employed large-scale Lagrangian transport (discussed below in the section on "LES SGS Closure: LEMLES").
In this study we consider three different turbulent premixed flame conditions corresponding to the CF, TRZ, and D/BRZ regimes based on the initial conditions. This is a well known and established canonical setup used for both DNS and LES in the past and also for earlier LEMLES studies Echekki and Chen, 1996;Rutland and Trouvé, 1993;Sankaran and Menon, 2005). Configurationally, it is a very simple constant pressure problem with inflow of reactants ahead of the flame and outflow of products behind it. Initially, the turbulence ahead of the flame is well in the CF, TRZ, and D/BRZ regimes, although as flame-turbulence interaction evolves, the turbulence intensity decreases leading to conditions in the wrinkled flamelet regime for the case originally in the CF regime and TRZ for the other two cases. Note that we do not force the background turbulence, as the flame-turbulence interaction is a multiscale and non-linearly coupled phenomena, which may become affected due to the artificial (unphysical) forcing of the background turbulence. We explore both structural and statistical properties of the flames and assess the role of SGS terms on the enstrophy transport, back-scatter of turbulent kinetic energy, and counter-gradient transport of the scalars.
This article is arranged as follows. The next section describes the governing equations and the subgrid closure formulations considered in this study. The numerical methodology and problem description are presented after that. The following section describes computational assessment of the employed SGS formulations. The discussion of results in terms of flame structure and its statistics, flame-turbulence interaction, and the role of SGS closures are then presented. Finally, the outcomes of this study are summarized.

Mathematical formulation
The compressible LES equations are obtained by spatially filtering the Navier-Stokes equations using a top-hat Favre filter, appropriate for finite volume schemes (Vreman et al., 1994), such that " f is the filtered quantity for a field variable f , and e f is the Favrefiltered quantity defined by e f ¼ ρf =ρ, where ρ represents the local fluid density. The final set of LES equations can be written as: Here, ρ is the density, ðu i Þ i¼1;2;3 is the velocity vector in Cartesian coordinates, P is the pressure, T represents the temperature, and Y k is the mass fraction for the k th species.
Also, N s is the total number of species in the flow. The total energy of the system is the sum of the internal energy (e) and the kinetic energy. As a consequence, the filtered total energy is given as the sum of the filtered internal energy e e, the resolved kinetic energy, 1 2 ½ e u i e u i , and the subgrid kinetic energy k sgs ¼ 1 2 ½g u i u i À e u i e u i . The thermally perfect gas equation of state is used to close the equations as P ¼ ρ e R e T þ T sgs and the filtered internal energy per unit mass is obtained as where e R, the mixture gas constant, is obtained from the species molecular weights ðW k Þ k¼1;...;N S , and R u the universal gas constant, as e enthalpy e h is defined as " ρ e h ¼ " ρe e þ P and for thermally perfect gas is expressed as: The specific heat at constant pressure for the k th species, c p;k , is obtained from classical temperature curve-fits. The filtered viscous stress tensor, ij , and the filtered heat-flux vector, q i are approximated as: Here, e S ij is the resolved rate of strain, given as e S ij ¼ 1 2 @e u i @x j þ @e u j @x i . Finally, e V i;k is the filtered diffusion velocity for the k th species, and is modeled as where W is the mixture molecular weight, and D k and e X k are diffusion coefficient and mole fraction of the k th species, respectively. The diffusion coefficient for a species is obtained through the well known mixture-averaged formulation (Poinsot and Veynante, 2005).
All the subgrid-scale terms, denoted with a "sgs" superscript, are unclosed, and therefore, require specific modeling. These terms are: It should be noted that, in the expressions for θ sgs i;k , q sgs i;k ; and E sgs k , the repeated index k does not imply summation. The strategy to model the terms denoted as subgrid-scale is presented next.

Closure for momentum and energy subgrid flux
The subgrid stress and enthalpy flux and the subgrid viscous work terms are present even in non-reacting flows, and their modeling follows past effort (Menon and Kim, 1996). The unclosed terms in the momentum equation, the subgrid stresses sgs ij , is closed as: Here, we employ a one-equation model for the subgrid kinetic energy (Kim and Menon, 1999;Menon and Kim, 1996) to determine the subgrid eddy viscosity as t ¼ C ffiffiffiffiffiffiffi k sgs p Δ, where Δ is the grid filter width and C is a coefficient calculated using a localized dynamic approach (LDKM) (Génin and Menon, 2010;Kim and Menon, 1999). The subgrid kinetic energy is obtained from a transport model for k sgs as: The closure approximations described above correspond to low Mach conditions and are appropriate for the present study. Further details of this approach and the fully compressible formulation are given elsewhere (Génin and Menon, 2010). Note that the last term in Eq. (20) is the subgrid dissipation of k sgs , which is modeled as sgs ¼ ρC ðk sgs Þ 3=2 =Δ, where C is another model coefficient obtained using the LDKM approach.
There are a few advantages of the LDKM approach using the k sgs approach: (a) the model does not assume equilibrium between production and dissipation at the SGS scale [as implicit in the Smagorinsky model (Smagorinsky, 1963)]; (b) the dynamic approach is localized and stable without requiring any ad hoc averaging; (c) as the grid is refined or the flow becomes laminarized, k sgs becomes naturally negligible since realizability conditions are enforced (Génin and Menon, 2010); and (d) the dynamic coefficients are always positive, resulting in a positive eddy viscosity and hence produces only forward-scatter [unlike the sometimes negative eddy viscosity reported in the dynamic Smagorinsky model (Germano et al., 1991)].

Scalar modeling
All the other SGS terms, namely, Y sgs i;k , θ sgs i;k , q sgs i;k , T sgs , and E sgs k , require closure as they specifically contain the species contribution. Typically, q sgs i;k , T sgs , and E sgs k , are neglected in LES (Fureby and Möller, 1995) and we also do it here as well. The other two terms, Y sgs i;k and θ sgs i;k , are important for closure and are discussed here. Consider the exact (i.e., unfiltered) species equation (written in a slightly different form): In the above equation, we have decomposed the instantaneous velocity field as where e u i is the LES-resolved velocity field, ðu 0 i Þ R is the LESresolved subgrid fluctuation (can be obtained from k sgs ), and ðu 0 i Þ S is the unresolved subgrid fluctuation. Now consider the filtered form of Eq. (21): The second term on the left-hand-side of Eq. (22) is the convective scalar flux that is resolved by the LES grid, while the first term on the right-hand-side is the resolved contribution to the subgrid scalar flux occurring inside the sub-grid (or inside the LES cell). The last two terms represent the filtered diffusion flux and reaction rate for the k th species. We consider SGS closures for each of these terms below. A conventional LES Favre-filtering of the above equation and rewriting by combining with the filtered conservation of mass results in: where S k represents all the SGS effects: There are some assumptions implicit in this form: (a) the spatial filtering and derivatives are assumed to commute (which is problematic on non-uniform grid), (b) the resolved diffusion flux velocity e V k;i is determined in terms of the filtered variables and the filtered scalar gradient, and (c) all Favre filtered subgrid velocity fluctuations are zero. Other terms that appear to still require modeling are the subgrid scalar flux Y sgs k;i , the subgrid diffusion flux θ sgs k;i and the filtered reaction rate _ ω k . Note that the commutation error associated with the filtering and derivative operators occurs when a non-uniform grid is used for spatial discretization, for example, while simulating the flow in a complex geometry. Typically, the commutation error term is lumped together with the SGS flux term prior to modeling (Grinstein and Fureby 2004). We discuss the closure for Eq. (21) or Eq. (23) in the following subsections, but before doing this, we briefly summarize the LEM model and its salient capabilities.

LEM subgrid model
The original LEM model (Kerstein, 1989) was a stand-alone model to account for concurrent interactions between turbulence, molecular diffusion, and reaction kinetics. The model as originally envisioned by Kerstein was one-dimensional (and hence, computationally very efficient), and contained three fundamental physics of scalar mixing and reaction kinetics: (a) molecular diffusion, (b) reaction kinetics, and (c) stirring by turbulent eddies. These processes are implemented by solving the reaction-diffusion model along a notional 1D domain: where s is coordinate along the 1D domain, and then concurrently modifying the scalar distribution by the stirring events that punctuate the deterministic evolution in the reaction-diffusion process. This turbulent stirring is implemented as stochastic events [based on the so-called triplet maps (Kerstein, 1989)] that attempts to capture the effect of turbulent stirring (based on Kolmogorov scaling) on the scalar field and mimics the effect of vortices on the scalar field. Successive folding and compressive motions are modeled so that an initially monotonic profile of, say, mixture fraction can develop multiple extrema with stoichiometric points in between, each corresponding to a flame location. At high turbulence intensity, the time scale of folding, compression, and diffusive mixing might become short enough relative to chemical time scales so that broad reaction zones, stirred internally by small eddies, can be formed. This ability of LEM is unique in the mixing model, but it is also apparent that this feature of LEM limits its applicability to turbulent flows where inertial range scaling is applicable. In 1D LEM, the instantaneous maps represent the outcome of an individual eddy motion, although such a literal connection between maps and eddies is not required. An LEM simulation time advances the 1D unsteady diffusion-reaction equations, including associated dilatation along the 1D domain, and this advancement is punctuated by instantaneous rearrangements of property profiles by mapping operations of a specified form. In effect, the outcome of each map constitutes a new initial condition for further time advancement.
The details of the 1D LEM as a standalone model (Menon and Kerstein, 2011;Menon et al., 1993) have been reported in details elsewhere and therefore are not repeated here. It has also been implemented as a subgrid model for LES (called LEMLES) (Menon and Kerstein, 2011). There are advantages and disadvantages to this implementation. Since the model is in 1D, the resolution can be fine enough to resolve the smallest scales and therefore, reaction kinetics and molecular diffusion can be simulated without any closure. Turbulent stirring in high Re flows recovers turbulent diffusivity (Sankaran and Menon, 2000;Smith and Menon, 1997), and the model has shown the ability to capture both high and low Schmidt number dependencies (Chakravarthy and Menon, 2001) without requiring any model adjustments. On the other hand, reduction to a 1D notional dimension (typically considered aligned normal to the scalar gradient) limits its ability in cases where the flame has to propagate in 3D as opposed to fluctuate around a statistically mean location. In particular, the tangential molecular diffusion is ignored. However, at high Re the turbulent diffusion usually dominates the molecular diffusion, which is captured by the 1D LEM model. Furthermore, many combustion problems of practical interest, as in swirl combustors, premixed burner, etc., involve relatively stationary flames, and therefore this closure has shown ability in many applications in the past. In this study we demonstrate a much simpler SGS closure using the 1D LEM model and assess its potential by comparing it to the more well established LEMLES. These closure strategies are discussed next.

LES SGS closure: LEMLES
The LEMLES is an approach that has been previously discussed (Menon and Kerstein, 2011) and therefore is only described briefly here. To summarize LEMLES (see Figure 1 as well), we can rewrite Eq. (21) as a two-step problem: where Δt LES is the LES time step and the superscript "n" indicates the current time level. It can be seen that the term in the brackets in the right-hand side of Step 1 is a subgrid Figure 1. A typical workflow in LEMLES and RRLES strategies. Note that in the present study, LEMLES implementation uses a single grid approach. evolution process that includes turbulent mixing, molecular diffusion, and finite-rate kinetics all occurring at the subgrid level (i.e., at scales smaller than the LES grid filter Δ). Within each LES cell, the SGS field evolves up to the CFL time step at which time the subgrid fields are convected across the LES cell faces using a Lagrangian transport model in Step 2. The SGS model in Step 1 is the same LEM model described in Eq. (25), while the second step is analogous to the advection step in a conventional finite-difference or finitevolume scheme. However, the implementation of Step 2 in LEMLES is distinctly different to the conventional finite-volume approach, as is reported elsewhere (e.g., Menon and Kerstein, 2011;Menon et al., 1993) since it employs a Lagrangian transport approach. This transport, implemented using a "splicing" algorithm (Menon et al., 1993), is designed primarily for high Re statistically stationary flames where turbulent mixing is much more intense than the molecular mixing at the large (LES-resolved) scales. As a result, this approach does not account for molecular diffusion across LES cell faces (even though Step 1 accounts for SGS molecular diffusion processes). In highly turbulent flows with significant turbulent mixing and transport, this limitation is not apparent, nor does it cause significant statistical errors. However, when most of the turbulence is resolved at the LES level or the flow is locally laminar, then Step 2 cannot account for the molecular diffusion at the resolved scales by the term @ @x i " ρ e Y k e V k;i Â Ã , and moreover, since subgrid stirring by u 0 i À Á S is negligible, the 1D LEM representation of subgrid molecular diffusion is not sufficient. These limitations of LEMLES are well known and are of concern as they do not asymptotically achieve laminar or DNS limits as originally implemented only as a true LES SGS closure. The effect of heat release is also included within the SGS model by carrying out an explicit expansion of the 1D domain, ensuring proper mass conservation (Menon and Kerstein, 2011). For computational and parallel efficiency, the expanded domain is regridded to the original LEM resolution, and this step is known to result in some numerical diffusion that must be carefully assessed. Past studies have shown that this error is tolerable when the flame is statistically stationary, as in the present case.
Finally, in LEMLES, after the completion of the two-step process, the filtered species field e Y k (necessary in the LES equations) is obtained from the subgrid 1D LEM fields as: where Y k;m is the mass fraction of the k th species for the LEM cell m and ρ m is the density of the LEM cell. This filtered species field e Y k is used in the remaining LES equations to complete the closure. Thus, Eq. (4) is not solved in LEMLES and is replaced by the twostep procedure described above. There is no need for explicit closure for Y sgs i;k or θ sgs k;i in this approach.
There are advantages and disadvantages to this type of closure.
Step 1 closure allows for direct resolution of small scale mixing, molecular diffusion in the subgrid and reaction rate without requiring explicit closure (albeit in 1D), and Step 2 allows transport of subgrid scalar gradients across LES cell faces without invoking a finite-volume type differencing. Thus, co-gradient and counter-gradient transport of scalars in the smallscales are captured. On the other hand, large-scale molecular transport across LES cell faces is ignored in this closure, although this is not a major issue when the turbulent transport dominates and the flame is statistically stationary. However, when all of the turbulence is in the resolved grid (as in DNS) or negligible (as in laminar conditions), then LEMLES cannot account for these effects. A remedy for this is investigated and described in the next section.

LES SGS closure: RRLES
In the new reaction-rate closure in LES approach, i.e., RRLES, Eq. (23) is used as in most conventional LES, except that the term S k is modeled differently. The subgrid scalar flux Y sgs k;i now requires a closure and a conventional gradient hypothesis is used as (Rose, 1997): where the eddy diffusivity D t is obtained from the eddy viscosity as D t ¼ t =Sc t , with Sc t being the turbulent Schmidt number (that can also be obtained dynamically but assumed unity here). An advantage of this approach is that when k sgs is negligible (which can occur under laminar or DNS limit), t ! 0 and the above term disappears. A major disadvantage is that in turbulent flows where SGS closure is necessary, only gradient transport is allowed and hence counter-gradient effects cannot be accounted for (unlike in LEMLES). With this closure for the subgrid scalar flux, the Eq. (23) can be rearranged as: The left-hand side of Eq. (30) contains all terms resulting from the convective transport and are now closed with the SGS model, while the right-hand side contains the LESresolved molecular diffusion flux gradient @ @x i " ρ e Y k e V k;i Â Ã , the resolved SGS diffusion flux θ sgs k;i , and the filtered reaction rate _ ω k . Since the LES-resolved diffusion flux is closed using properties at the resolved scale (and in the laminar or DNS limit reduces to the classical molecular diffusion term), only the two terms in S k given by Eq. (24) require explicit modeling. In the current implementation we specifically provide a closure for S k using the subgrid LEM closure as described below. The advantages and disadvantages of the proposed approach are also discussed. The RRLES approach also employs the 1D LEM subgrid model to obtain a closure for the subgrid processes (the terms in Step 1 of LEMLES). Once the subgrid reactiondiffusion equation [Eq. (25)] as in the Step 1 in LEMLES is solved on a 1D domain within each LES cell, the term S k is approximated as described below. Then, the species governing equations in Eq. (30) are solved along with the other LES equations using conventional finite-volume scheme.
Since the Lagrangian transport ( Step 2) of the LEMLES is avoided here, implementation of RRLES requires some specific coupling process. For example, in LEMLES the subgrid scalar fields evolve continuously from initialization, but in RRLES, to enforce conservation of mass and species mass fractions ð e Y k Þ that are evolving on the LES grid, the 1D fields need to be reconstructed every time step based on the corresponding LES-resolved quantities. This reconstruction approach is required since the resolved level does not explicitly contain any subgrid information.
Here, we propose a model for reconstruction that builds on the notional idea that the 1D LEM domains are aligned along the direction of maximum scalar gradient (which has a physical basis and is used in LEMLES) and combine with a dual-grid approach (shown in Figure 2) so that the LEM subgrid scalar fields evolve on a coarser LES grid with resolution Δ LEM ¼ 2Δ (similar to the test filter approach used in LES) while the LES governing equations are still solved on the finer grid Δ. This dual-grid strategy (which could be extended to multiple grids) has the advantage that the estimate for the SGS variation in the 1D LEM (at Δ LEM ) can be constructed from the solution on the finer LES grid. The subgrid turbulence used for stirring in the 1D LEM can also be estimated by the test filtering process. Figure 1 shows a flow chart of the RRLES strategy in comparison to the LEMLES strategy. The key differences between these methods are observable primarily as the replacement of Step 2 in LEMLES by the estimate of S k from RRLES. The reconstruction of the subgrid field is the key issue in RRLES. At present, the gradient based reconstruction is summarized by the following equation: where R is the restriction operator (see Figure 2) that does a conservative filtering from the fine grid (Δ) to the coarse grid (2Δ) to initialize the LEM fields (some sensitivity to how this operator is employed is discussed in the next section). The superscripts "c" and "f " denote the coarse and the fine grids, respectively. The filtered reaction rate within the 1D LEM is computed as _ ω k;LEM ¼ where, ΔV m is the volume of the m th LEM cell. Since the 1D LEM equations given by Eq. (25) are solved at the subgrid level, the effects of the subgrid molecular diffusion and turbulent stirring (as in Step 1 of LEMLES) are implicitly included in Eq. (32). Afterward, Eq. (32) is projected to the fine grid (see Figure 2) using a prolongation operation, such that: There are advantages and disadvantages of this approach. For example, the aforementioned limitation of the LEMLES approach in regions of low turbulence is avoided since there are more unresolved scales at the coarser grid where LEM is used. Also, since the "splicing" algorithm of LEMLES ( Step 2) is avoided, all the limitations of this approach in the limiting case of DNS or in laminar regions are avoided since molecular diffusion at the LES level by the term @ @x i ρ e Y k e V k;i Â Ã is explicitly included. On the other hand, countergradient transport is not feasible by this closure.
Another subtle issue observed is that in the asymptotic limit of low turbulence, the filtered reaction rate term does not approach the quasi-laminar reaction rate limit (Grinstein and Kailasanath, 1995). In the quasi-laminar chemistry-based closure in large-eddy simulation, which is sometimes also called a no-model approach for the reaction-rate term, the filteredreaction rate is simply expressed as _ ω k;QL ¼ _ w k P; e T; e Y 1 ; e Y 2 ; . . . ; e Y N s . To achieve the asymptotic limit under both high and low turbulent conditions, a hybrid approach for the filtered reaction rate is used and Eq. (33) is replaced as where ¼ t = is the ratio of the turbulent eddy viscosity and the kinematic viscosity, respectively, and f is a blending function in the range ½0; 1. We have done some preliminary assessment of the effect of the functional form of f ð Þ on the flame structure in the section "Assessment of RRLES Modeling Strategy." The requirement in Eq. (34) is that at low or negligible SGS turbulence level, f ð Þ ! 0 whereas at higher turbulence level, which is characterized by ) 1, f ð Þ ! 1. At moderate turbulence levels when the ratio ¼ Oð1Þ, the reaction rates computed from the LEM field and the fine grid LES are comparable and therefore both can contribute to the LES closure. The function f ð Þ is determined in a dynamic manner (since t is obtained using a localized dynamic approach, LDKM) and therefore, this approach adapts with the solution. Although a more comprehensive assessment of this hybrid approach will be reported in the future, the results reported in this article show the potential of this approach as a multi-scale SGS closure for finite-rate kinetics modeling.

Numerical approach and problem description
The LES equations are solved using a finite-volume solver called LESLIE. The numerical method to solve these equations is well established and has been used in the past to study both canonical and complex flows (Génin and Menon, 2010;Masquelet and Menon, 2010;Patel and Menon, 2008;Sankaran and Menon, 2005;Undapalli et al., 2009). The scheme is nominally fourth-order accurate in space and time, and implemented in parallel for computational efficiency. Although the spatial discretization at the DNS/LES level is fourth-order accurate, in the LEMLES and the RRLES formulations, due to the use of splicing and dual-grid strategy, respectively, formally the method is second-order accurate. Thermally perfect gas is considered in this study, and all the transport properties such as viscosity, specific heat, etc., are computed by the mixture-averaged formulation. We use a moderately complex four-step and eight-species methane-air mechanism (Smith and Menon, 1997) to include effects of finite-rate chemical kinetics. Figure 3a shows a schematic of the premixed planar flame configuration. We consider three different turbulent premixed flames corresponding to the CF, TRZ, and B/DRZ regimes (based on the initial turbulence and flame conditions) to analyze the resolved and SGS dynamics of flame-turbulence interaction. These cases are indicated on the premixed regime diagram  in Figure 3b at both the initial time and the time at which the flame statistics are analyzed. The initial flame front is obtained from a laminar premixed flame solution, and is specified near the center of the domain with reactants and products on its left and right sides, respectively. The extent of the computational domain is L Â L Â L in the streamwise (x), transverse (y), and spanwise (z) directions, where L ¼ 0:0055 m. The flow field is initialized using an isotropic turbulent flow field obtained using the Kraichnan spectrum (Kraichnan, 1970). It is further superimposed with the one-dimensional planar flame solution obtained at ϕ ¼ 0:8, T ref ¼ 570 K, and P ref ¼ 1 atm. Here, ϕ denotes equivalence ratio of the methane-air mixture, T ref is the temperature on the reactants side, and P ref is the reference pressure. The flame conditions, particularly the preheated conditions and the equivalence ratio chosen here, are nominally based on past studies, which are typical of gas turbines, spark-ignition engines, and combustors (Sankaran et al., 2007;Sankaran and Menon, 2005). A characteristic based inflow-outflow boundary condition is used in the streamwise (x) direction and periodic boundary condition is used along the spanwise (z) and the transverse (y) directions. Table 1 summarizes the simulation parameters of all the cases in terms of turbulence and flame parameters. The three premixed flames are estimated to be in the CF, TRZ, and B/DRZ regimes based on initial turbulence intensity, and are respectively labeled as Case A, Case B, and Case C. We conduct four simulations employing DNS, LEMLES, RRLES, and QLLES. These cases are labeled as Case A m , Case B m , and Case C m , where subscript "m" corresponds to different closures. Here DNS and LES cases are conducted using the same code and the same scheme. Note that the LES code is capable of simulating more complex flows (Masquelet and Menon, 2010;Patel and Menon, 2008;and Undapalli et al, 2009). The LES cases are performed using the various closures described above. In Table 1, l is the integral length scale, ¼ =S L is the Zeldovich flame thickness, u 0 is the turbulence intensity, S L is the laminar flame speed, and Re, Ka, and Da are respectively, the integral Reynolds number, the Karlovitz number, and the Damköhler number defined as Re , and Da ¼ S L l u 0 . These values are estimated based on the initial conditions of the simulation.
The grid resolution, N x , N y , and N z along x-, y-, and z-directions, respectively, is chosen based on the past studies, and for the conditions reported here is sufficient to reach k max η ! 1 for DNS, where k max is the largest wave number and η is the Kolmogorov length scale. In particular, k max η ¼ 7:23, 2.6, and 2.1 for Case A 1 , Case B 1 , and Case C 1 , respectively. Note that the grid resolution of a well resolved LES of the flame-turbulence interaction is dictated by two requirements. First, the turbulence on the LES grid should be resolved in a manner so that it satisfies the Pope criterion (Pope, 2004), which requires that the resolved turbulent kinetic energy should be approximately 80% of the total kinetic energy. Second, the flame structure should be adequately resolved, which requires approximately 10 points across the the thermal flame thickness, i.e., L ¼ T b À T u ð Þ = ÑT j j max , when no closure is used for turbulence-chemistry interaction (Colin et al., 2000). Here, subscripts "b" and "u" denote burned and unburned regions, respectively. The LES cases for Case A, Case B, and Case C approximately resolve about 90%, 80%, and 70% of the total turbulent kinetic energy at t=t 0 ¼ 2, 2, and 3, respectively. Here t 0 ¼ u 0 =l is the initial eddy turnover time. In addition, with the grids employed in the present study, it is estimated that the thermal flame thickness is resolved by around by 20 points in DNS and 5 points in LES in Case A. Clearly, the LES grid does not resolve the thermal flame thickness and can be considered as a borderline, particularly for the QLLES cases. Also, 12 LEM cells are used in each LES cell for both LEMLES and RRLES. Estimates suggest that, in the SGS scales, eddies of size of the order of η=2 are resolved. Simulation using an increased resolution with 24 LEM cells does not show any significant effects on the flame statistics (not shown here), and, therefore, all the results reported here are based on simulations employing 12 LEM cells.
The simulations are carried out long enough to allow flame-turbulence interaction to evolve and all the results are compared after two (three) initial eddy turnover times for Case A/B (C), i.e., t=t 0 ¼ 2ð3Þ. Even though turbulence decays in time in the present study, there is a period during which flame-turbulence interaction attains a quasi-stationary state, and therefore, the dynamics and statistics associated with the flame-turbulence interaction can be analyzed during such period. Typically, the initialization with the Kraichnan spectrum evolves to a physical state in about one eddy turnover time and additional one or two eddy turnover time is sufficient for flameturbulence interaction to evolve (Savre et al., 2013). The time evolution of the turbulent flame speed shown in Figure 6 (discussed later) is used to assess the evolution of flame-turbulence interaction to a quasi-stationary state. The DNS data are filtered to the LES grid scale using the top-hat filter to have a consistent comparison.

Assessment of RRLES modeling strategy
The RRLES formulation requires a reconstruction strategy for the subgrid field and a blending function as discussed in the section "LES SGS Closure: RRLES." We first perform some limited numerical experiments with Case C to assess sensitivity of the prediction to the choice of the reconstruction strategy and the blending function parameters. Then we freeze the RRLES parameters for all the simulations discussed later in the "Results and Discussion" section.
The reconstruction method used in RRLES relies on the notional idea that the LEM fields are aligned along the direction of the maximum scalar gradient. The idea has a strong physical basis (Kerstein, 1989), but we investigate how sensitive is the filtered reaction rate to the initial form of reconstruction. Although more studies are warranted and are being carried out, here we discuss three reconstruction strategies. The first is the baseline approach in which the two-grid approach is used to reconstruct the subgrid scalar fields at the 2Δ resolution. In the second approach, a Gaussian distribution of the subgrid field is specified using the mean and variance of the distribution determined from the resolved level values. For the subgrid variance we use a model proposed (Pierce and Moin, 1998) of the following form: where ϕ is a scalar and f ϕ 002 is its subgrid variance. The constant C is taken to be 0.5 based on numerical experiments. Often this model for the scalar variance is not adequate enough, and typically an additional transport equation for the scalar variance is solved Poinsot and Veynante, 2005). However, in the present study, we have only used this model as one of the approaches to initialize the subgrid fields at the LEM level so that sensitivity of initialization on the flame statistics can be assessed. In the third approach, no subgrid variation is considered. The entire subgrid field is set to the resolved level values. The different forms of the reconstruction and the subgrid structure after one LEM evolution step are shown in Figure 4 for Case C. The LEM evolution step involves the stirring events accompanied with reaction and diffusion. Figure 5a compares the prediction of the filtered reaction rate with the reference DNS results. It is clear that all forms of initial reconstruction match reasonably well with DNS, implying that the final value of the filtered reaction rates might not be very sensitive to the initial form of reconstruction for this particular case, although there are some noticeable differences. In particular, the uniform reconstruction marginally over-   predicts in the post-flame region, and in the preheat zone 0:2 e c 0:5 ð Þwhere the turbulence-chemistry interaction is significant. On the other hand the initial reconstruction based on Gaussian distribution consistently under-predicts the reaction rate. More studies are needed, but for now we use the gradient-based reconstruction for all the subsequent studies.
The second input for the RRLES is the blending function, f ð Þ, is used in Eq. (34) to ensure the model asymptotically achieves both turbulent and laminar limits. Again, this is an area that requires further assessment but for now, we use a smooth blending function given by where a, b, and c are blending parameters that determine the position and the rate at which the blending transitions from the LEM to the QL model predicted reaction rates. We currently fix a ¼ 100 and c ¼ À50, and investigate the sensitivity to b in Figure 5b. For b ¼ À0:475, the blending achieves a reasonable match with the DNS results in the preheat region and this value is used for all the simulations described in the next section.

Results and discussion
We first compare the flame structure and its statistics predicted by different closures and then assess the ability of various closures to capture the SGS effects on the flameturbulence interactions in terms of enstrophy transport, back-scatter of turbulent kinetic energy, and type of turbulent transport.

Structural features of flames
The flame structure and its characteristics are assessed using a progress variable c, which is defined based on the mass fraction of the fuel (methane) as where Y CH 4 denotes the mass fraction of methane, and subscripts "u" and "b" denote its value in the fresh reactants and the burned products sides, respectively. The value of c varies from 0 in the fresh reactants to 1 in the burned products. As in past studies (Sankaran and Menon, 2005;Savre et al., 2013), we identify a notional flame surface by an iso-level of c ¼ 0:8 and the bounds of the flame brush correspond to c 2 ½0:01; 0:99. While comparing LEMLES and RRLES predictions, the filtered progress variable e c is defined in manner similar to Eq. (37) but using the filtered mass fraction of the fuel. Direct comparison is difficult since e c is obtained from LEMLES with the mass fraction filtered from the SGS LEM field, whereas in RRLES the filtered mass fraction is directly obtained from the LES. The statistics of flame-turbulence interactions are analyzed after such interactions reach a quasi-stationary state. This is identified in terms of time evolution of the turbulent flame speed S T , which is defined as where V is the computational volume and the subscript "0" denotes conditions at the inlet (reactants side). The turbulent flame speed is a global quantity, which specifies the mean speed at which reactants must be provided to keep the position of the flame front statistically stationary. There are alternate definitions of the turbulent flame speed (Poinsot and Veynante, 2005), and each of these definitions predicts a similar qualitative behavior. Figure 6 shows the time evolution of the normalized turbulent flame speed S T =S L for Case A 1 , Case B 1 , and Case C 1 . We can observe that the flame speed initially increases, and then saturates to a quasi-stationary value. In Case A 1 , B 1 , and C 1 , the value of S T =S L is 1.2, 1.4, and 2.9, respectively, during such stage, which are reached at approximately t=t 0 ¼ 2 for Case A 1 and B 1 and 3 for Case C 1 . Therefore, all the discussion of the flame statistics in the present study for Case A, B, and C are presented at these times. Figure 7 shows instantaneous structure of the flame brush obtained from DNS for the three flames. We can observe intense wrinkling and stretching of the initially planar flame surface by the turbulent eddies. The effect of heat release and associated thermal expansion across the flame is apparent in terms of increase in the overall length scales associated with the protruding structures on the reactants and the products sides. A typical role of turbulent eddy is to transport pockets of cold reactants toward the reaction zone, which sustains the flame. As expected, the flame structure is continuous, thus precluding any local extinction. The width of the flame brush varies spatially, due to the straining effect of the large-scale eddies, and the increased flame surface area results in an increased fuel consumption rate (Poinsot and Veynante, 2005;Sankaran and Menon, 2005;Yuen and Gülder, 2009). Another key feature to observe is a broadening of the flame brush as the Karlovitz number increases from Case A to Case C. This occurs due to enhanced transport of heat and mass away from the reaction zone toward the fresh reactants by the turbulent eddies (Lipatnikov and Chomiak, 2002). Finally, an important feature for Case C in Figures 7e and 7f is that due to penetration of eddies within the flame brush region and enhanced mixing, localized cold pockets of reactants are surrounded by hot products, which in turn affect the dynamics of flame-turbulence interaction through modification of dilatation and baroclinic torque (discussed further in the section "Enstrophy Transport").  The wrinkling and stretching of the flame is predicted by all closures in a consistent manner, and the increase in the wrinkling along with enhanced transport of heat from the reaction zone from Case A to Case C is qualitatively captured by all closures. All simulations show a continuous flame structure and no sign of local extinction for these conditions. Qualitatively, there are some noticeable differences between different closures. In particular, DNS, RRLES, and QLLES show more fine-grained features compared to the LEMLES formulations, which is particularly apparent in Case C. Note that the interpretation of the LEMLES results for the scalar field differs from the other formulations due to its multi-scale nature, as the LEMLES-based resolved scalar field corresponds to the resolved value of the corresponding LEM field on the LES grid obtained through Eq. (28), whereas the resolved scalar field is explicitly solved in the other formulations. A key feature to notice is that, in Case A, LEMLES does not capture the flame structure in a manner similar to the other approaches, where we can observe significant differences. This occurs due to a very low or negligible subgrid turbulence in Case A, which affects the prediction of the flame-turbulence interactions by the LEMLES formulation. Note that LEMLES neglects the laminar diffusion of the scalar fields across the LES cell faces, which becomes important under low turbulence conditions as in Case A. However, the RRLES formulation is able to capture the flame structure as it includes the laminar diffusion, and additionally, due to the blending approach employed by the formulation through Eq. (34), it asymptotes toward the quasi-laminar value of the filtered reaction-rate term when turbulence is either low or resolved.
The wrinkling of the flame structure by turbulent eddies is further evident in Figure 9, where instantaneous streamwise profile of the filtered progress variable e c is shown for all cases. The averaged profile of e c is also shown for reference (discussed further in the next section). In each case, we observe a significant spreading about the averaged profile showing the effect of turbulent eddies that are penetrating and broadening the preheat zone. The flame broadening and disruption of the preheat zone with an increase in Ka is consistent with the past studies (Lipatnikov and Chomiak, 2002). The instantaneous flame structure, which is predominantly monotonic in Case A and Case B, becomes non-monotonic (multiple flame crossings) in Case C due to increased turbulent mixing, and this effect is captured by all the LES closures. However, in Case A, the predictions by LEMLES show a very sharp flame structure compared to the other cases, which demonstrates the limitations of the LEMLES formulation when the subgrid turbulence level is either small or resolved. Note that the multiple flame crossing in Case C occurs at the subgrid level in both LEMLES and RRLES, but their effects are captured differently due to the nature of closure used here. A key feature to observe is the flame excursion on the reactants side, which is marginally larger in QLLES in comparison to the other cases. This occurs due to a higher turbulent flame speed (not  shown here) and the filtered reaction-rate (discussed below) predicted by the QLLES formulation, particularly in Case C. Table 2 summarizes some key structural characteristics of the flames, such as wrinkling factor Ξ ð Þ, mean flame curvature (" κ), and mean reaction zone thickness " 0:5 À Á . The wrinkling factor ΞðtÞ is defined as ΞðtÞ ¼ AðtÞ=Aðt ¼ 0Þ, where AðtÞ denotes area of the flame front at time t identified by e c ¼ 0:8. The flame curvature κ is defined as κ ¼ @n i @x i c¼c 0 , where n i ¼ À @c @x i = @c @x i is the unit-normal vector pointing toward the reactants and c 0 ¼ 0:8 is the value of the progress variable used to mark the flame surface. A negative (positive) value of the curvature implies a concave (convex) flame with respect to the reactants. The reaction zone thickness 0:5 is defined in terms of the thermal flame thickness corresponding to a specific value of the progress variable (Yuen and Gülder, 2009) as 0:5 ðtÞ ¼ T b À T u ð Þ = ÑT j j c Ã ¼0:5 , where T b and T u denote temperature in the burned products and unburned reactants sides, respectively, and c Ã is another progress variable defined as c Ã ¼ T À T u ð Þ= T b À T u ð Þ . Similar to the progress variable c given by Eq. (37), c Ã varies from 0 to 1 from fresh reactants to burned products sides, respectively.
In all cases, the wrinkling factor Ξ>1, indicating an increased flame surface area due to wrinkling effects of the turbulent eddies, which is apparent from Figures 7-9. Overall, there is a very good agreement between different closures for all the cases except LEMLES, which shows a higher (lower) value of Ξ indicating increased (decreased) wrinkling of the resolved flame surface in Case B (C). The mean value of the curvature predicted by all the closures (except LEMLES) is close to zero in Case A and Case B, implying the presence of both concave and convex surface elements. Clearly, LEMLES is not able to capture the curvature effects at the smaller level of turbulence. This can be attributed to the role of laminar diffusion in Case A and Case B of the scalar field across the LES cells, which is ignored in the LEMLES formulation. With an increase in Ka, the mean value of the curvature increases. This trend is consistently captured by all the closures, but there are noticeable quantitative differences. In particular, QLLES over-predicts this quantity in Case C. RRLES is able to capture the curvature in good agreement with the DNS data in all the cases. Finally, the reaction zone thickness does not show much variation in all the cases and RRLES is able to consistently capture the trend as Ka is increased. These results clearly demonstrate that RRLES can asymptote toward laminar or turbulent conditions in a continuous manner and can accurately capture the flame statistics. While QLLES performs reasonably, some of the statistics, particularly curvature, indicate the need to have a closure model for turbulence-chemistry interaction as Ka is increased.

Flame statistics
The global flame structure is examined in terms of the spatially averaged quantities, where the spatial averaging of a quantity q is performed along the homogeneous y-and z-directions through qðx; y; z; tÞdydz: (39) Figure 10 shows the streamwise profile of the spatially averaged quantities for all cases. In all the cases, we observe a good agreement of the results with the DNS data, and all the methods are able to capture the global flame structure, except in Case A where LEMLES under-predicts the turbulent flame speed, thus affecting the flame statistics. Additionally, there are some noticeable differences in results from different closures. In Case A with lower value of Ka, LEMLES under-predicts the turbulent flame speed, which leads to a different spatial structure of the flame structure. In particular, the filtered reaction rate term is under-predicted and much narrower in extent compared to other cases. In Case B, LEMLES under-predicts the filtered reaction rate term, but overall the flame structure compares well with the other cases. However, at higher value of Ka (Case C), we observe that the LEMLES agrees with the DNS data. This clearly demonstrates the need to have a (c) Case C Figure 10. Streamwise profile of the spatially averaged scaled temperature, fuel mass fraction, and fuel reaction rate. closure that can asymptote to lower turbulence or laminar conditions. QLLES tend to over-predict the reaction rate, which manifests into the flame excursions towards the reactants side as observed in Figure 9l. RRLES shows an overall good agreement with the DNS results for the cases, implying that the blending employed by the formulation allows it to include contribution from LEM at the subgrid level as the turbulence level increases and with a decrease in turbulence level it takes contributions from the quasi-laminar chemistry through Eq. (34). With an increase in Ka, we observe the flame-broadening homogenization of the preheat zone and a drop in the peak value of the heat release. With enhanced transport of mass and heat, a reduction in fuel consumption is evident from a wider distribution of the progress variable compared to a sharp front observed in the laminar solution. All of the SGS closures demonstrate these features albeit with noticeable differences.
The coupling of differential diffusion (non-unity Lewis number) effects with the flame curvature κ significantly impacts the local burning rate. For example, in case of Le < 1, the diffusion of reactants is higher than the thermal diffusion. Therefore, regions of flame having convex/concave shape toward reactants leads to a focusing/de-focusing effect, which increases/decreases the local reaction rate and the flame propagation speed. The PDF of the normalized flame curvature shown in Figure 11 suggests that the mean value of κ is close to zero, implying the presence of both convex and concave elements in Case A. With an increase in Ka, we can notice an increase in the mean and skewness of the curvature. Overall, the features of the PDF of curvature predicted by RRLES compare reasonably well with the DNS and are consistent with past studies (Sankaran et al., 2007;Yuen and Gülder, 2009). However, the PDF predicted by QLLES differs in both cases, suggesting that SGS modeling does play a role. The RRLES result is closer to QLLES in Case A and Case B and is closer to LEMLES in Case C, demonstrating the asymptotic nature of the formulation. Figure 12 shows the conditional variation of the filtered reaction rate for methane. The peak location is shifted toward higher value of e c for all cases. The improved predictions by the RRLES formulation are noticeable in all the cases, where we can observe a good agreement with the DNS data. The asymptotic behavior of the RRLES closure is again evident here, as we can observe that its prediction approaches DNS and QLLES in Case A when the level of turbulence is small. These predictions by the various closures highlight their abilities and limitations. The unphysical results predicted by LEMLES in Case A demonstrate the limitations of the formulation when the local flow conditions tend to have a low subgrid-scale turbulence.

Enstrophy transport
The effect of thermal expansion due to heat release on the small scales is shown in Figure 13 in the form of instantaneous contours of vorticity magnitude in the central x À y plane for the DNS cases. The vorticity field is associated with the smallest turbulent scales where heat release occurs (Bobbitt and Blanquart, 2015;andSwaminathan, 2007, Steinberg et al., 2012) and plays an important role in the flame stretching, as the vorticity vector preferentially lies in the plane of the flame leading to an extensional strain of the flame (Rutland and Trouvé, 1993). It is apparent in Figure 13 that heat release dissipates small scale vorticity across the flame leaving only larger scales. Further in the post-flame region, the flow tends to recover toward an isotropic state but with relatively larger scales.
Since small scales get dissipated by the heat release, we investigate dynamics near these scales to assess model performance. Enstrophy, which is defined as Ω ¼ ω i ω i , where ω i is the vorticity field, can be used to assess interactions at the small scales since it represents turbulent kinetic energy at the dissipation scales. An enstrophy transport equation has been proposed earlier in the context of RANS modeling (Gorski and Bernard, 1996) and was used recently to post-process DNS data (Bobbitt and Blanquart, 2015) to investigate scaling of various terms in the transport equation. Here, we use a similar approach to assess the current results. The unfiltered enstrophy equation can be written as: The terms on the right-hand side of Eq. (40) are, respectively, vortex stretching, dilatation, baroclinic torque, and viscous dissipation. Note that vortex stretching and viscous dissipation terms are present in both non-reacting and reacting flows. However, dilatation and baroclinic torque terms play a major role in reacting flows and, therefore, are investigated here. Additionally, the dilatation term in Eq. (40) is a scaled value of the flow field dilatation ðÑ Á uÞ by the enstrophy. Figure 14 shows normalized dilatation and baroclinic torque variation across the flame region for the DNS cases. The normalization of these terms is performed by using the where denotes a particular term in the enstrophy equation and c k is a parameter with a value of 10 3 , 10 5 , and 10 7 , respectively, for Case A, Case B, and Case C. The dilatation and baroclinic torque terms are small in the post-flame (or product) region due to the heat-release-induced thermal expansion. Baroclinic torque term is significant only directly ahead of the flame surface in the preheat zone, particularly in Case A and Case B, whereas the dilatation term varies significantly ahead of the flame due to turbulent transport of heat and mass away from the flame surface toward reactants, particularly in Case B and Case C. In the post-flame region, the small scales get dissipated leading to a spatial homogeneity of the density, and therefore, a smaller value of the dilatation term. Note that the baroclinic torque term occurs due to a non-alignment of the gradient of density and pressure fields. Compared to Case A, these terms are larger in Case B and Case C, which is associated with the increased stretching of the flame surface and the associated spatial non-homogeneity of the density field. Figure 15 shows profiles of the four terms in the right-hand-side of Eq. (40) conditioned on the progress variable. In Case B and Case C, we observe that the vortex stretching term monotonically decreases with respect to c and as expected the viscous dissipation term always remains negative in all the cases, which is consistent with past observations (Bobbitt and Blanquart, 2015). However, the dilatation term changes sign across the flame region as Ka is increased (see Figure 16b, and Figure 18c). In Case C, the dilatation term in Eq. (40), which is the scaled value of the flow field dilatation ðÑ Á uÞ by the enstrophy (a positivesemidefinite quantity by definition), shows that there is expansion ðÑ Á u > 0Þ ahead of the flame region (0 < c % < 0:5), which becomes compression ðÑ Á u < 0Þ in the later part of the flame region (0:5 % < c < 1), before becoming negligible in the post-flame region. This change in behavior is consistent with the results shown in Figure 7-9, where the presence of localized pockets of cold reactants surrounded by expanding hot products is evident. Such a disruption of the preheat zone by the turbulent eddies manifests into a negative value of the dilatation term in the later half of the flame brush, which is discussed below further in terms of the flow field dilatation. This behavior is not evident in Case A and Case B, where eddies are not strong enough to disrupt the structure of the flame brush. This is further evident from the PDF of the dilatation term at different values of the progress variable (c) for Case C 1 , which is shown in Figure 16a. We can clearly observe a consistent transition in the shape of the PDF, which demonstrates the appearance of sign change of the dilatation term within the flame brush region as observed in Figure 15c. Table 3  dilatation term in Eq. (40) at different values of c, which shows that the mean values become negative, implying the presence of compression effects for larger values of c. In addition, the skewness becomes increasingly negative, which further demonstrates the transitional behavior of the dilatation term across the flame brush region at higher values of Ka. To further understand the role of turbulence intensity on this behavior, we consider the time evolution of the dilatation term conditioned with respect to the progress variable for Case C 1 at three different times, namely t=t 0 ¼ 2, 3, and 5, which is shown in Figure 16b. It is evident that as the turbulence decays, the magnitude of the negative value of the dilatation term reduces and a predominantly positive value of the dilatation term occurs within the flame brush region as also observed in Case A 1 and Case B 1 . This clearly indicates a transient nature of the competing effects of mixing induced due to enhanced turbulence and thermal expansion due to heat release, which leads to multi-scale and anisotropic flame-turbulence interactions. Furthermore, it is apparent that the transition of the dilatation term in Eq. (40) from a positive to a negative value occurs due to a transitional behavior of the flow field dilatation ðÑ Á uÞ across the flame brush, which is examined next.
The flow field dilatation can be expressed in terms of the effects of thermal expansion due to heat release and molecular diffusion (mixing), and is given by the dilatation equation (Swaminathan et al., , 1997 where ¼ ðT b À T u Þ=T u is the heat release parameter, W c Ã is the reaction rate, D ; Ñ Á ραÑc Ã is the molecular mixing, and c Ã ¼ ðT À T u Þ=ðT b À T u Þ is a progress variable based on the temperature field. Note that the dilatation equation is based on the assumption that the molecular weight is uniform across the flame brush and is equal to the unburned value. This assumption is reasonable to analyze the competing effects of thermal expansion and molecular diffusion in the present study. Figure 17 shows the conditional variation of the three terms in Eq. (41) normalized by the laminar flame speed S L and thermal flame thickness L . We can observe that the effect of the thermal expansion due to heat release is positive throughout the flame brush, whereas the molecular mixing term transitions from a positive value on the unburned side (0 < c % < 0:5) to a negative value on the later part of the flame brush. The net effect of these two terms indicates a change in sign in the flow field dilatation. Such a transition from expansion (Ñ Á u > 0) to compression (Ñ Á u < 0) is clearly an effect of an increased level of mixing, a characteristic feature of high Ka premixed flames (Aspden et al., 2011;Savre et al., 2013), which brings pockets of unburned reactants close to the region where heat-release is occurring. The change in sign of the flow field dilatation is also a transient phenomenon as shown in Figure 16b through the dilatation term of the enstrophy transport equation. Since the competing effect of thermal expansion and molecular mixing is transient in nature, it may not be observed in statistically quasi-stationary premixed flame turbulence interaction simulations, where the background turbulence is maintained through an artificial forcing. Figure 18 shows that the conditional variation of dilation shows quantitative differences between DNS and the LES closures. In Case A, all the closures under-predict the magnitude of dilatation and baroclinic torque terms. LEMLES tends to under-predict the magnitude of dilatation and baroclinic torque terms in Case B and Case C. RRLES and QLLES under-predicts in the reactants side of the flame brush, i.e., 0 e c 0:5, whereas toward later part of the flame brush these closures approach the DNS predictions in Case B and Case C. In Case C, the abrupt change in the behavior of dilatation term is captured in all the simulations although there are noticeable differences in prediction by all the closures in comparison to the DNS data. So, qualitatively RRLES and QLLES capture the trends in both the cases, but the differences are apparent in the preheat zone, particularly in Case C where the multi-scale turbulence-chemistry interaction occurs and there is a presence of anisotropic scales of motion.
These results clearly highlight limitations of these models and warrant further investigation. Since the assessment of SGS closures is done in a posteriori manner, these are considered predictions by models that do not explicitly include these effects. Future studies will address these results to determine if SGS models can be devised for both fluid and scalar subgrid fluxes to capture these effects.

Energy transfer mechanism
Heat release and thermal expansion at the small scales impact energy transfer mechanisms in turbulence. Although forward-scatter or energy cascade from large to small scales always occur, in reacting flows there is also some feedback. To characterize these effects, we consider the transport equation for the resolved turbulent kinetic energy, which is given by  Figure 17. Conditional variation of the normalized dilatation ( L Ñ Á u=ðS L Þ), molecular diffusion ( L D=ðρ u S L Þ), and reaction-rate ( L ρW c Ã =ðρ u S L )) with respect to the progress variable (c Ã ) at t=t 0 ¼ 3 for Case C 1 . where Here, α Å and Å are the pressure-velocity correlation terms, α is the work done by the viscous stress, α sgs is the work done by the subgrid stress, is the kinetic energy dissipation rate due to the molecular viscosity, and sgs is the subgrid dissipation rate. In general, sgs > 0 implies a forward-cascade of the turbulent kinetic energy from large to small scales. However, it is possible for the SGS energy dissipation rate to become a locally negative value implying a back-scatter of turbulent kinetic energy from small to large large implying presence of a backward-scatter. These features are not included in conventional SGS closures and remain a modeling challenge. The behavior of sgs and sgs;c is further examined in terms of the single point statistics by extracting PDF of these quantities at different locations within the flame brush, which is shown in Figure 20. In Case A, we can clearly observe the presence of backwardscatter across the flame brush region. In Case B and Case C and at all values of e c, we can clearly observe the presence of both forward and backward-scatter. However, the mean value remains positive as shown in Figure 19. A key feature to note is that for 0:2 e c < 0:8, i.e., in the preheat zone, the amount of backward-scatter increases with e c. At the conditions considered in this study, the turbulent eddies can enter the preheat zone leading to homogenization, and their residence time within the flame brush region  Figure 19. Conditional variation of sgs (a, c, e) and sgs;c (b, d, f) for Case A 1 , Case B 1 , and Case C 1 . Solid and dashed curves denote the conditional averaged value and the standard deviation about the averaged value, respectively. is adequate to allow for cross-scale exchange of energy between the exothermic smallscales and the large-scales. However, for e c ! 0:8, i.e., in the post-flame region, the energy transfer predominantly shows the forward-scatter due to a predominant presence of the large-scales and relaxation of turbulence toward isotropy.
To assess the impact of eddy diffusivity closures employed in the current LES, we first investigate the degree of correlation between sgs and sgs;c in Figure 21 in the DNS data. We can observe a significant correlation between both these quantities (even when both are negative). The variation of sgs with respect to sgs;c shows a presence of substantial back-scatter in Case A and marginal back-scatter in Case B and Case C. This implies that, statistically, forward-scatter of the turbulent kinetic energy and scalar variance is predominant at high Ka, which is apparent in Figure 19, but back-scatter may have to be explicitly included. A model for back-scatter within this closure strategy is being explored and will be reported in the near future.
In the current LES, both sgs and sgs;c are positive by virtue of their closure, but analysis of the DNS data shows that there can be negative values for these quantities. So to assess this further, we compare the modeled and the computed quantities in Figure 22. In addition to the filtered DNS data, we use the DNS data to compute the modeled quantities. It is apparent that in Case A, none of the closures can capture the back-scatter phenomenon as evident from the DNS results. In Case B and Case C, both sgs and sgs;c remain positive but at lower (higher) Ka, i.e., Case B (Case C), LES over-predicts (under-predicts) in comparison to the DNS data. A need for modified SGS dissipation closure is apparent from these results.

Turbulent scalar transport
Another important aspect of flame-turbulence interaction is that both co-and countergradient (CG) turbulent transport phenomena are typically observed in turbulent premixed flames . This CG transport is usually associated with the gas expansion, whereas turbulent mixing is considered to be responsible for the co-gradient transport (Wenzel and Peters, 2000). We analyze turbulent transport by comparing SGS scalar flux Y sgs i;c and the resolved scalar gradient @e c=@x i . Without loss of generality, we examine only the terms along the x-direction (x;x 1 ). Figure 23 shows a scatter plot of the scaled Y sgs 1;c with respect to the scaled @e c=@x 1 obtained from the filtered DNS data. The scaling is based on the maximum magnitude of the corresponding terms. The scatter data in the first (both positive) and the third quadrants (both negative) indicate CG transport, which is more prominent for Case A and Case C in comparison to Case B. In Case A, the CG transport is evident throughout the flame brush region due to a lower level of turbulence and the effect of thermal expansion. In Case B, the CG transport is associated mainly with the preheat zone of the flame, whereas in Case C, it is again observed all across the flame region. This is notably different from expectations as with increase in turbulence level, role of turbulent mixing should be dominant. However, this can be attributed to the stretching effects of turbulence, which lead to creation of sharp gradients, which coupled with gas expansion lead to CG transport. To assess if the original LEMLES model handles such effects albeit indirectly, we filter the LEMLES resolved fields with a filter of size 2Δ and recompute the SGS scalar flux. We . Figure 24 shows variation of hY sgs 1;c i with respect to the resolved gradient of the filtered progress variable in the x-direction for only Case B and C. Since LEMLES leads to unphysical results for the flame-turbulence interaction in Case A, therefore, the analysis of counter-gradient transport has not been considered for this case. We observe that LEMLES does predict both co-and counter-gradient transport, and as the turbulence level is increased from Case B to Case C, the presence of countergradient transport is much more apparent. These results demonstrate one unique capability of the LEMLES that needs to be exploited further in RRLES since the latter model does not have this capability.

Conclusions
Numerical studies of subgrid-scale dynamics of premixed flame interaction with decaying isotropic turbulence are carried out using two variants (LEMLES and RRLES) of a subgrid mixing model based on LEM and quasi-laminar chemistry approximation within LES. The two variants based on LEM as a subgrid reaction-diffusion-mixing model differ in their coupling to the resolved LES scales. We consider three premixed flames at different levels of initial turbulence intensity and analyze flame-turbulence interactions as the Karlovitz number is increased. The RRLES formulation uses a dual-grid framework that alleviates some of the limitations associated with the original LEMLES formulation. The strengths and limitations of both models are assessed by comparing with DNS data obtained using the same code. Overall, both closures capture the flame structure in reasonable agreement with the DNS data but with some noticeable differences in quantitative predictions, which occur primarily due to differences in the type of scalar transport employed by the two models. In particular, the RRLES formulation includes large-scale molecular diffusion, which is ignored in the LEMLES formulation, and therefore, the model can be used in cases where the turbulence is either resolved at the LES level as in a DNS or if the turbulence is negligible. Additionally, due to use of a dual-grid strategy in the RRLES formulation, we have more unresolved scales available at the subgrid LEM level, thus leading to improved predictions. However, a major disadvantage of the RRLES formulation is that it cannot capture counter-gradient turbulent transport, which is a unique feature of the LEMLES formulation. This aspect of the model will be investigated further in future studies.
Further analysis suggests that both these SGS closures still cannot capture all the critical physics associated with the dynamics of flame-turbulence interaction. Using the enstrophy balance equation dilatation and baroclinic torque contributions are identified as key properties strongly affected by thermal expansion with a significant SGS contribution that current closures do not account for quantitatively. Forward/backward-scatter of energy is assessed in terms of the SGS dissipation rate of turbulent kinetic energy and scalar variance, and it is shown that when the effect of turbulent diffusion reduces in comparison to the thermal expansion due to combustion, back-scatter was evident. The back-scatter implies an increased dissipation of the small-scales in the post-flame region and enhanced large-scales, which are ignored in the SGS closures considered in this study. The presence of countergradient transport is also observed in both the cases, suggesting co-gradient SGS scalar transport may not be adequate. Future studies will explicitly focus on these issues.