Large-Eddy Simulation of Premixed Combustion in the Corrugated-Flamelet Regime

ABSTRACT Large-eddy simulation (LES) is applied to a fuel-lean turbulent propane-air Bunsen flame in the corrugated-flamelet regime. The subgrid-scale (SGS) modeling includes a previously developed treatment of the total enthalpy along with three different SGS velocity, , models. In addressing the filtered reaction rate, a presumed probability density function (PDF) approach is employed for the reaction-progress variable, closed by a transport equation for its SGS variance. The statistics obtained using the three models are in good agreement with the measurements and do not differ significantly from each other for first-order moments suggesting that commonly used SGS modeling may be adequate to get the mean velocities and reaction progress variable. However, all three SGS velocity models fail to reflect a measured bimodality of the PDF of the radial component of the velocity in the central portion of the flame. This emphasizes a need for further development of models required at the reaction rate closure level for practical LES of combustion in the corrugated-flamelet regime.


Introduction
Lean premixed combustion is attractive for practical combustion devices because of its potential to deliver high efficiency and low emissions simultaneously. The sensitivity of lean premixed flames to fluctuations in fluid-dynamic and thermochemical conditions of the reactant mixture, however, causes numerical modeling of such combustion processes to be challenging. Many different modeling avenues have been developed; a summary of the approaches to large-eddy simulation (LES) of premixed combustion is given in Poinsot and Veynante (2005), Pitsch (2006), Gicquel et al., (2012), Bray (2011), andCant (2011). These approaches can be broadly categorized into flamelet and nonflamelet approaches, the former being selected here for the present investigation. Flamelet approach can be further subdivided into geometrical and statistical approaches, with the G-equation and flame-surface-density-based methods belonging to the geometrical category and presumed probability density function (PDF) methods to the statistical category. Flamelet models based on artificial thickening or filtering of laminar flame are also possible (Duwig, 2007;Nguyen et al., 2010;Wang et al., 2011). The present contribution falls within PDF methods, working with a presumed PDF for a reaction-progress variable.
Premixed combustion occurs in different regimes, depending on the competing effects of turbulence and finite-rate chemistry. The flame structures in different regimes differ significantly, ranging from wrinkled flamelets at low turbulence intensities to distributed reaction zones at high intensities (Peters, 2000). The corrugated-flamelet regime lies adjacent to the wrinkled-flamelet regime and is separated from regimes at higher turbulence intensities by the thin-reaction-zone regime. The latter is typical of piloted flames, having high Reynolds numbers and experiencing strong turbulence resulting predominantly from shear-production mechanisms, representative of most practical applications. For this reason, many LES studies aimed at predicting experimental observations have been focused on the thin-reaction-zones regime combustion (Colin et al., 2000;De and Acharya, 2009;Gubba et al., 2012;Wang et al., 2012) and a direct numerical simulation (DNS) with tabulated chemistry approach has been also been attempted for an experimental configuration (Moureau et al., 2011). There is also interest, however, in the corrugated-flamelet regime, in which, just as in the thin-reaction-zones regime, experimental Reynolds numbers are too large for DNS to be feasible, and for which the heat-release effects of the flames on the turbulence may be proportionally greater as a consequence of the lower incoming turbulence levels.
As a continuation of a series of experiments by Furukawa and co-workers (Furukawa et al., 2002;Furukawa and Williams, 2003;Furukawa et al., 2013aFurukawa et al., , 2013b, recently, Furukawa et al. (2016) reported measured statistics of gas velocities, with reaction-progress-variable fields deduced therefrom in a non-piloted lean propane-air premixed Bunsen flame in an open environment. This is one of four flames in the corrugated-flamelet regime that they investigated, and they suggested that it would be of interest to pursue numerical modeling of the experimental results. Because of the aforementioned importance of fuel-lean combustion, the present work was performed with the aim of shedding light on the structure of this flame. This provides a good opportunity to gather insight into the capabilities and limitations of statistical flamelet models for subgrid-scale (SGS) combustion in LES of this flame. Since the shear-produced turbulence does not affect this flame, as one shall see in the results section, the turbulence level experienced by the flame results mainly from flame-intermittency and thermal-expansion effects. These are SGS phenomena in a practical LES, and thus they would affect the SGS velocity that needs to be modeled. It may be debatable whether one needs to include dilatation effect in this velocity modeling or not as propositions have been made in the past to exclude this effect in general (Colin et al., 2000). Three SGS velocity models treating the thermal-expansion effects differently are considered to assess their influence on the premixed combustion from both physical and numerical viewpoints.
This article is organized as follows. The experimental flame is described next in the next section. Following that the numerical model including the SGS combustion and velocity closures with the numerical procedures and boundary conditions are discussed. Results are presented and discussed in section the next, and concluding remarks are made in the final section. Figure 1 shows a direct photograph and a schematic diagram of the flame considered here, a propane-air flame having an equivalence ratio of ϕ ¼ 0:85. Calculations with the PREMIX code (Kee et al., 1985) and complex chemical-kinetic schemes (Sung et al., 1998;UC San Diego, 2014) produce a laminar flame speed of s L ¼ 0:34 m=s and a laminar-flame thermal thickness of th ¼ 0:4 mm. The reactant mixture issues from a pipe of diameter D ¼ 26 mm and spreads into an open environment, as shown schematically in Figure 1b. The bulk mean velocity based on the mass flow rate at cold conditions is U b ¼ 4 m=s. The shear layer present at the interface between jet and surrounding air encloses the flame as sketched in the figure, and this flame configuration is quite different from piloted turbulent premixed flames, where the flame resides inside the shear layer between reactant and pilot streams. Since the flame is not in the shear layer for this configuration, the turbulence experienced by the flame comes from the inlet initially produced by the fully developed turbulent pipe flow but modified by thermal-expansion effects and flame-front intermittency. In addition, the turbulence produced by the shear through the Kelvin-Helmoltz instability is not expected to affect the flame, as will be seen later.

Experiment
A 3D LDV (laser Doppler velocimetry) system having a spherical resolution volume of about 0.2 mm diameter was used to measure three components of velocity, from which various characteristics of the flame were deduced. Measurements along the centerline were taken from h ¼ 10 mm to 110 mm and along the radius for r 15 mm at given height, h, ranging from 30 mm to 90 mm. These measurement regions, indicated in Figure 1b, did not include the shear layer or region where the shear layer interacts with the core of the jet (marked as "Transition" in Figure 1b). In other words, the measurements were focused in regions where shear-generated turbulence effects are weak. This experiment offers a good challenge for combustion modeling because the turbulence induced by the flame intermittency and the turbulence-flame interaction is important in these flames, affecting the required physics of submodels.

Governing equations
The Favre-filtered conservation equations for mass, momentum, and total enthalpy are solved along with other equations required for premixed combustion. These additional equations include filtered progress variable and its SGS variance, to be discussed in the next subsection. The instantaneous progress variable, c, is defined as where Y f ;u is the fuel mass fraction in the unburned mixture and c takes a value of 0 and 1 in the fresh and burned mixtures, respectively. The subgrid stresses are closed using a localized dynamic Smagorinsky model (Germano et al., 1991;Lilly, 1992) and the SGS scalar fluxes are computed using a gradient hypothesis with a dynamic Schmidt number approach (Lilly, 1992). The modeling constant in the localized dynamic Smagorinsky model used here decreases to zero in laminar regions making this model suitable also for low-turbulence regions investigated here. However, other choices like Vreman (Vreman, 2004) or static Sigma (Nicoud et al., 2011) models are possible and these models are expected to give improvements only close to walls or in transitional regions (Kemenov et al., 2012;Rieth et al., 2014;Vreman, 2004) (not relevant for this study). Furthermore, these models involve additional model constants (thus their tuning) and so they are not considered for this study. Also, the SGS scalar fluxes are treated using gradient transport hypothesis as in many earlier LES and the influence of counter-gradient transport for these fluxes, which may exist for most part of the filtered flame under the conditions investigated here, is of interest for a future study. However, adequacy of the gradient hypothesis can be investigated using the experimental comparisons shown in later sections of this article. The transport equation for the filtered total enthalpy, e h, is: where is the thermal diffusivity. The Favre-filtered temperature, e T, is computed as Δh 0 f ;mix Þ= e C p;mix , where f Δh 0 f ;mix is the enthalpy of formation for the mixture and e C p;mix is the specific heat at constant pressure for the mixture, which depends on temperature as described by Ruan et al. (2014), where further details may be found. The mixture density is computed from the state equation, ¼ p e W mix =ðR e TÞ, where p is the filtered pressure, e W mix is the Favre-filtered molecular weight of the mixture and R is the universal gas constant.
In the code, in general, the values of e C p;mix , f Δh 0 f ;mix , and e W mix are computed on the basis of the formula: where e Φ mix refers to any one of the three quantities above, Φ air is its value in air, and e Φ reac is its value in the reacting mixture. This mixing rule is used to include the dilution of burned mixture with the entrained air. A transport equation is used to compute the Favrefiltered passive marker, e ψ, and this equation is: where D ψ is the diffusivity of ψ. The values of e Φ reac are computed from laminar-flamelet values using an equation similar to Eq. (4) to be discussed later. The flamelet values are obtained by computing unstrained planar laminar flames using the PREMIX code (Kee et al., 1985) and complex chemical-kinetic schemes for propane-air combustion. Two mechanisms (Sung et al., 1998;UC San Diego, 2014) were considered but no substantial differences in the laminar-flame structures, speed, or thickness were observed.

Combustion closure
The filtered reaction rate and its related terms in the SGS variance equation are closed using an unstrained flamelet approach. The filtered reaction rate is modeled as: where PðζÞ ¼ e PðζÞ and e PðζÞ is the density weighted subgrid PDF with ζ as the sample space variable for c. The flamelet reaction rate is _ ω ¼ _ W, which is obtained from the freely propagating unstrained planar laminar flame for a given thermochemical If this fiugure was previously printed in another publication please obtain permission to reprint from the original publisher and provide documentation.
condition. The subgrid PDF is modeled using the beta function, which is expressed as: e PðζÞ ¼ ζ aÀ1 ð1 À ζÞ bÀ1 Ψða; bÞ with Ψða; bÞ ¼ The parameters a and b are related to the filtered progress variable, e c, and SGS variance, σ 2 c;sgs , through: a ¼ e c e cð1 À e cÞ σ 2 c;sgs À 1 The choice of a beta-PDF is not uncommon for LES of premixed combustion and has been employed with good results in many past studies, see, for example, Landenfeld et al.  Nambulli et al. (2014), and . Fiorina et al. (2015) remarked that the beta-PDF is not relevant when the SGS flame wrinkling is completely resolved. When this wrinkling is fully resolved (which is a DNS) then the SGS PDF would be a delta function. One would need a very fine grid to resolve the SGS wrinkling, which is not the objective in this article, where the aim is to keep the LES practical. Moreover, the beta-PDF automatically captures the bimodal behavior due to flame intermittency when the SGS variance, which is high due to intermittency, is well predicted as one shall see next. However, other choices of PDF are possible. The values of e c and σ 2 c;sgs are obtained by solving their transport equations. The progress variable equation is: where D is the molecular diffusivity of c. The transport equation for σ 2 c;sgs ¼ e c 2 À e c 2 is written as: where ν t and Sc t are the SGS viscosity and Schmidt number. The third term of Eq. (8) related to chemical reaction is closed by using: which is consistent with the closure in Eq. (4). The subgrid dissipation rate (SDR) is defined from the equation e ε c ¼ DðÑc Á ÑcÞ À DðÑe c Á Ñe cÞ h i and is closed using an algebraic model (Dunstan et al., 2013): This model for the SDR e ε c has been used in conjunction with an algebraic reaction-rate closure for LES of piloted turbulent Bunsen flames . A similar model with corrections for nonunity Lewis numbers has also been studied using DNS (Gao et al., 2014a(Gao et al., , 2014b and LES (Butz et al., 2015;Ma et al., 2014). The factor F in Eq. (10), defined as F ¼ 1 À exp Àθ 5 Δ þ ð Þ, with θ 5 ¼ 0:75 and Δ þ ¼ Δ= th , ensures that e ε c ! 0 when Δ þ ! 0. The filter width, Δ, is assumed to be Δ ¼ V 1=3 i following common practice, where V i is the volume of cell i. The subgrid Damköhler number is Da Δ ¼ t sgs =t c , where t c ¼ th =s L is the chemical time scale and t sgs ¼ k sgs =ε sgs is the SGS flow time scale, k sgs and ε sgs being the SGS kinetic energy and dissipation rate, respectively. These two SGS quantities are related to Δ and the SGS velocity u 0 Δ , which requires modeling. The heatrelease parameter is ¼ ðT ad À T u Þ=T u , where T ad is the adiabatic flame temperature and T u is the reactant temperature. The term σ 2 c;sgs = c describes the effects of flame-front curvature induced by turbulence, chemical reaction, and molecular dissipation. The scaledependent parameter c has previously been obtained dynamically , but that approach cannot be used for the flame investigated here because of the low turbulence level. A dynamic procedure usually invokes scale-similarity assumption, which holds if the energy of a fluctuating reacting scalar decreases monotonically around the scale of interest, which is the flame scale. If the turbulence level is low then the energy at SGS level may be significantly larger than that at the smallest resolved scale. Then, the dynamic procedure becomes ill-posed as noted by Charlette et al. (2002) for a dynamic power-law formulation. This issue is different to dynamic formulation for SGS stresses, like the dynamic Smagorinsky model, where the turbulent kinetic energy cascade to the small scales implies that the energy level between two subsequent scales is always similar. The dynamic formulation for c may lead to its singularity [the denominator of Eq. (4) in  may approach zero], which is unphysical. This was verified through preliminary simulations and for these reasons a static value of c is used, as detailed in the first subsection under "Results." The terms (K c s L = th ) and C 3 À C 4 Da Δ ð Þ ð u 0 Δ =ΔÞ arise from fluctuating dilatation and strain rate resulting from competing effects of turbulence and heat release, respectively, and they are all influenced by SGS turbulence. The parameters describing the influences of thermochemical and turbulence processes and their interplay are K c ¼ 0:79, Dunstan et al., 2013). The basis for these functional forms is discussed in detail elsewhere (Chakraborty and Swaminathan, 2007;Dunstan et al., 2013;Kolla et al., 2009).
The SGS velocity, u 0 Δ , is modeled in three different ways: (1) A model based on SGS stress-related closure proposed by Lilly (1967), given by: is the turbulent viscosity modeled using a localized dynamic Smagorinsky model. The model constant is C L ¼ 0:094.
(2) A second model based on scale-similarity for velocity is (Pope, 2000): where e e u is the velocity vector field obtained using a Gaussian test filter of width b Δ ¼ 2Δ. Different estimates for C q have been given in the literature (Bardina et al., 1980;Colucci et al., 1998;Métais and Lesieur, 1992), and its dynamic evaluation is possible (Cook, 1997), but that requires an estimate for the integral length scale, which is not easy. Also, the turbulence level is low for the experimental flame. For this reason the value C q ¼ 1, which is comparable to the previous estimates, is used here.
(3) Colin et al. (2000) derived another model by expanding Eq. (12) in a Taylor series up to second order and then applying the curl operator to exclude dilatation effects. The resulting expression is: where c 2 ¼ 2 is a typical value suggested by Colin et al. (2000).
One could also use a transport equation for SGS kinetic energy, k sgs , to compute both SGS stresses in the momentum equation and to evaluate the SGS velocity through This transport equation involves pressure correlation terms, which are important for combustion in corrugated-flamelet regime. These terms need modeling and no satisfactory model for LES of premixed combustion is available in the literature. For these reasons, this transport equation is not employed for this study.

Numerical method and boundary conditions
The above conservation and modeled equations are solved with a finite-volume code, PRECISE-MB (Anand et al., 1999). The spatial derivatives are discretized by a secondorder central-differencing scheme for all variables except the progress variable and its variance, which are discretized using a nonoscillatory bounded HLPA (hybrid linear/ parabolic approximation) scheme (Zhu, 1991;Zhu and Rodi, 1991). The time advancement is performed using a second-order scheme with a constant time step so that the CFL number is smaller than 0.3 in the entire computational domain. The velocity-pressure coupling is achieved by using the SIMPLEC algorithm (Doormaal and Raithby, 1984).
The cylindrical computational domain with D c Â L c ¼ 1 m Â 1 m is discretized on a nonuniform grid having 278 Â 506 cells in the radial and axial directions, respectively, with 44 cells for the jet diameter, D, and 80 cells in the azimuthal direction. This grid results in 5.4M cells in total for the chosen computational volume. The filter size, based on the numerical cell volume, has a minimum width of Δ % 1:3 th (0.52 mm) close to the jet exit. Another refined grid having about 23M cells, obtained by increasing the number of cells by about 60% in all three directions, is used to test grid dependency of the flame statistics.
The progress variable is assigned to be 0 and 1 for the reactant jet and entrainment inlet boundaries, respectively. A small velocity of 0.1 m/s and a temperature of 298 K are assigned for the entrained air at the respective inlet boundary. Slip-flow and zero-gradient conditions are assigned on the lateral boundaries for velocity and scalars, respectively. A zero-gradient condition is applied at the outlet for all variables but pressure, which is assigned to be atmospheric.
The cold-flow velocity measurements had provided an exit velocity profile with a maximum velocity on the centerline of U 0 ¼ 4:5m=s (Furukawa et al., 2016). This value, however, is inconsistent with the measured and calculated centerline value of 5.5 m/s at about one exit-diameter distance from the jet exit, implying a strong acceleration over this distance if it were 4.5 m/s at the jet exit. This observation prompted further measurements of jet tube-wall temperature since wall heating of the reactant mixture inside the tube when the flame is present could cause the exit velocity to increase. To capture these influences correctly, a preliminary RANS simulation using Fluent code and k À ε model was performed for the computational domain and conditions illustrated in Figure 2. This cylindrical domain is discretized using 31 Â 775 cells in radial and axial directions, respectively. In terms of wall units, the first cell close to the wall at the pipe exit has y þ ¼ yu =v % 6, where y is the radial coordinate, ν is the kinematic viscosity of the flow, and u ¼ ffiffiffiffiffi w p = is the shear velocity, w being the shear stress at the wall, and these quantities are evaluated a posteriori. Thus, the boundary layer is not resolved and the wallenhancement treatment in Fluent has been used to account for the proper resolution of the flow near the wall. The measured mass flow rate of the propane-air mixture at 298 K is specified at the inlet with flat profiles. A linear temperature profile in the vertical direction along the wall, as shown in the figure, was specified near the exit, based on experimental wall-temperature estimates. This linear profile gives a temperature gradient of about 3 K/ mm along the wall.
The radial variations of axial velocity and temperature at the jet exit obtained from this RANS calculation are shown in Figures 3a and 3b, respectively. The temperature at the jet exit is nonuniform and accelerates the flow. The centerline value of axial velocity was calculated to be about 5.4 m/s, which is close to 5.5 m/s observed in the experiment at 1D from the jet exit. The axial velocity variation in the cold condition shown in Figure 3a depicts the flow acceleration when the wall temperature varies axially. These computed profiles for nonisothermal flow were specified as boundary conditions for the LES.
Although the rms value of 0.22 m/s, reported by Furukawa et al. (2016), is low at the jet exit, it is important to include this because the flame experiences only this turbulence and does not see the shear-generated turbulence as noted earlier. A synthetic turbulence having similar characteristics as in the experiment, therefore, was generated using the digital-filter technique (Klein et al., 2002) and was added to the velocity U shown in Figure 3a for the jet exit. The variations of rms velocity in three directions are imposed using the turbulent kinetic energy obtained from the RANS computation. The value in the axial direction is taken to be twice that in the other two directions. This is acceptable because the centerline value in the axial direction, u 0 rms % 0:24m=s, is close to the measured value of 0.22 m/s. This gives u 0 rms =U 0 % 4:4%, which is slightly smaller than the measured value of 4:8% for the cold flow.
The radial variation of u rms at the jet exit is shown in Figure 3c. The measured longitudinal integral length scale of Λ x ¼ 10:9 mm and the transversal length scales of Λ y ¼ Λ z ¼ D=4 are used. The scales are taken to be uniform as required by the digital filter technique (Klein et al., 2002). The assigned turbulence may thus have some uncertainty and so a sensitivity analysis has been performed by varying the inlet rms values. An increase of 30% of the rms produces an almost 13% decrease in flame length, while a zero inlet rms results in a flame length almost 50% longer than the measured value and also an increased flame 'flapping' leading to the PDFs to exhibit skewness when this artificial inlet condition is applied. This analysis enlightens a strong sensitivity of the flame to the inlet turbulence as anticipated. The uncertainty at the inlet will be removed to some extent in this study by carefully selecting the modeling parameter c to yield the measured flame length (see the next section).
The grid size was also increased by a factor of 1.6 in all directions, resulting in 23M cells for the same computational volume as noted earlier to test grid sensitivity. The CFL number was kept below 0.3 everywhere in the domain by using constant time steps of 30μs and 50μs, respectively, for the 23M and 5.4M grids. The LES are run for 1:6s of real time using 208 cores for 23M and 32 cores for 5.4M grids, of which the latter half ð0:8sÞ is used for collecting statistics. These simulations took about 72 h on a wall clock using 2.60 GHz eightcore Sandy Bridge E5-2670 processors of the Darwin cluster at Cambridge University.

Role of SGS velocity modeling
The SGS velocity will be influenced by dilatation resulting from combustion, whence u 0 Δ values obtained from Eqs. (11) to (13) can differ, since these models presume the role of the dilatation differently. Any difference in u 0 Δ gives a differentε c according to Eq. (10), which influences σ 2 c;sgs . This can influence the filtered reaction rate [see Eqs.
(4) to (6)], which, in turn, can affect the u 0 Δ . Thus, the combustion model employed here has two-way coupling between the flame and the turbulence. Figure 4 shows the sensitivity of the SDRε c (normalized by th =s L ) to u 0 Δ (also normalized), while keeping other parameters in Eq. (10) to be constant, with Δ þ ¼ 1:4, a representative value for the flame regions based on the grid used here. The quantity σ 2 c;sgs = c is taken to be unity for the purpose of this estimate. The variation is weakly nonlinear for the range shown, a 100% change in u 0þ Δ producing only about a 20% change inε þ c . This implies that the same variation could be obtained by changing c by about 20%.
Since the incoming turbulence influences u 0 Δ , any difference in its value obtained synthetically can also influenceε þ c . The subgrid velocity can be written as u 0 with a constant C Þ 0 and u 0Ã Δ given by Eqs. (11) to (13) without their respective model constants. Because of the near linearity observed in Figure 4, Eq. (10) can be approximated as: where e ε c;0 ¼ e ε c ðu 0 Δ ¼ 0Þ and a is to be found to minimize the error (for example, using a ¼ ðe ε c;u 0 Δ ¼2s L À e ε c;0 Þ=2s L gives a maximum error below 2% in the range u 0 Δ ¼ ½0; 2s L ). This is not necessary because the objective here is to demonstrate that e ε c and u 0 Δ models can be tuned simultaneously as discussed below, and Eq. (14) is not actually used for computation. From Eq. (14), the variation of dissipation rate can be expressed as Δe ε c ¼ e ε c À e ε c;0 ¼ aCu 0Ã Δ ¼ Cf 1 ðu 0Ã Δ Þ= c . The variation ofε c for C ¼ 1 is shown in Figure 4. From Eq. (14) one can define 0 c % c =C, which can be used instead of C and c because of the near linearity ofε c with u 0 Δ observed in Figure 4 over the range ½0; 2 of u 0þ Δ expected for the confguration investigated in this study. The value of 0 c for the various u 0 Δ models is chosen by matching the computed and measured flame lengths, estimated to be the location of the peak axial velocity from the jet exit along the centerline, to complete a comparative analysis. It is worth noting that such analysis would not be possible without considerable uncertainty for quantities with no experimental data. The above tuning was necessary to overcome the difficulty produced by the modeling constants of SGS velocity and SDR models, and to some extent that of the inlet turbulence, which affects directly the flame length as discussed in the last section. Alternatively, when dynamic modeling is not possible one might use values suggested by experiments or DNS data available in the literature [see Gao et al. (2014a) for analysis of c using DNS data], which are, however, not universal. This was unnecessary in this study where experimental data is available, and some tuning is probably the best choice to minimize uncertainties for a meaningful comparative analysis.
The values of 0 c obtained using the three u 0 Δ models given in Eqs. (11) to (13), respectively, Lilly, Pope, and Colin et al. models, are shown in Table 1. The 0 c for Lilly model is larger than those for the other two models, suggesting that the SGS velocity estimated by it is also different as will be seen later. The computed axial velocity along the centerline is compared to the measured values in Figure 5 for the three u 0 Δ models. Practically there is no difference among the computed values, which also agree well with measurements.
The variation of e ε þ c within the filtered flame is shown in Figure 6 for the three u 0 Δ models. The values obtained from Pope, Eq. (12), and Lilly, Eq. (13), models for u 0 Δ look Table 1. Optimum values of 0 c found using the three SGS velocity models in Eqs. (11) to (13).   similar, with the latter exhibiting more scatter. The greater scatter with larger values is due to the effects of intermittency of the interface between the reacted jet fluid and entrained air at locations closer to the jet exit. The SDR with Lilly model, Eq. (11), is smaller as a result of the higher value of 0 c . This can be also interpreted to mean that this model needs reduced dissipation to achieve the same flame length or that it overestimates u 0 Δ when 0 c ¼ 6 is used.
Further insights can be obtained by studying the spatial variation of hu 0Ã Δ i, which is the time-averaged u 0Ã Δ ; typical variations are shown in Figures 7a and 7b for the mid-plane. Formally, this quantity is u 0Ã Δ [see Eq. (14)], and the plots are normalized by the respective maximum values to force the variation from 0 to 1, facilitating comparisons, needed because of the differences in 0 c . A significantly larger value and wider region is seen for the shear layer in the model of Colin et al., compared with that of Pope, even in nonreacting regions. Although the model of Colin et al. in Eq. (13) is derived from that of Pope in Eq. (12) there is no guarantee that the SGS velocities predicted with these two models will be the same for nonreacting regions. The maximum values of u 0Ã Δ given in the  A second region having nonzero u 0 Δ is seen in Figures 7a and 7b for Popes' model, Eq. (12), resulting from flame-induced effects. This flame induced u 0 Δ inside the flame brush is of a similar magnitude as that in the shear-layer region. The flame-induced SGS velocity is not visible for the model of Colin et al. because it is designed to exclude the dilatation effect. Nevertheless, the dimensional hu 0Ã Δ i inside the flame brush obtained using Pope, Eq. (12), and Colin et al.,(13) are both found to be about 0.2 m/s, masked for Eq. (13) by the very large value in the nonreacting region.
For Lilly model in Eq. (11), the value of hu 0Ã Δ i in the region of the flame brush is about 0.02 m/s, an order of magnitude smaller than those of the other two models. Since the viscosity is computed dynamically, the dilatational part in Lilly model coming from the strain rate is masked by the fact that the Smagorinsky constant goes to zero in quasilaminar regions. For the same reason u 0 Δ approaches zero even in the shear region close to the jet exit, as seen in Figure 4b. It is worth noting that Lilly model was originally derived from the stresses without involving any dynamic constant. One can think of using this model with a sensible value for Smagorinsky constant (c s ¼ 0:1), which is known to work well when the turbulence level is sufficiently large. However, this approach is not followed here because it is less suited for quasi-laminar and inhomogeneous flows, and also would introduce additional uncertainty through the value for c s . The present study is only concerned to show that Lilly model is unsuitable for predicting effects of dilatation on the SGS velocity when dynamic Smagorinsky model is used for the SGS stresses. However, it would be of interest to investigate the influence of SGS stresses models on the SGS velocity, which is a case for a future study.
One can also compare the u 0 Δ values with a turbulent velocity scale, is the resolved turbulent kinetic energy to assess the relative magnitude between SGS and resolved energy. This is shown in Figure 4c for Pope and Colin et al. models, which gives very similar predictions of u 0 turb . Lilly model also gives very similar predictions and thus is not shown. Using the data in Figure 4 it is easy to verify that a significant part of energy is at SGS scales (use Pope model). An estimation of the Pope's criterion for turbulent kinetic energy, k sgs =K 20% for regions of turbulence without combustion, is also possible and one can easily verify that this criterion is satisfied in the shear region for Lilly and Pope models but not Colin et al. model. However, this is irrelevant for the study here as the SGS velocity constants are scaled through the tuning of parameter 0 c as discussed at the beginning of this section. It is also worth noting that by definition k sgs is a residual kinetic energy and not residual turbulent kinetic energy.
The above analysis shows the potential of the various SGS velocity models considered here. The flame length in this analysis is tuned by selecting an appropriate value for 0 c , and the value of this parameter alters the amount of subgrid SDR in Eq. (14), which in turn alters the value of σ 2 c;sgs and thus _ ω c . Since 0 c ¼ c =f ðCÞ, the role of 0 c can be interpreted as that of controlling the magnitude of u 0 Δ to maintain the correct value of e ε c . These comparisons of the u 0 Δ models suggest that either the dilatation effect should be included in Eq. (13), or a different value for c 2 is required in that equation. Lilly model in Eq. (11) is not suitable for the present flame conditions because the dynamic viscosity used in that model approaches zero, causing the magnitude of u 0 Δ to approach zero, becoming less than the fresh-mixture rms velocity fluctuation irrespective of flame-induced effects.
It will be seen later, however, that none of these models address properly a dilatationrelated subgrid velocity-fluctuation effect that is dominant in the center of the flame zone in the corrugated-flamelet regime.
Comparison of statistics Figure 8 shows axial velocity contours, both instantaneous and time averaged, along with a contour of mean progress variable of he ci ¼ 0:5. These contours indicate the positions of the shear layer and the flame, showing clearly that these two regions do not overlap. The jet becomes highly turbulent at x % 200 mm, but the flame height is only about 100 mm. The acceleration through the flame tip is visible at x % 100 mm, which is consistent with the experimental data, as seen in Figure 5.
The computed rms of radial velocity fluctuation along the centerline is compared with the experimental data in Figure 9 for the three SGS velocity models. The values shown here are only the resolved parts, found to be insensitive to the u 0 Δ model selected. The computed values are seen to agree reasonably with the experimental data. Some underestimate close to the jet exit arise from the uncertainties in the inlet turbulence. There is also some underestimate in the flame region x % 100 mm ð Þand the reason for which will become clear in the next section.
The radial variations of computed mean and rms velocities are compared with the experimental results in Figure 10. These statistics also are insensitive to the u 0 Δ model, and hUi, ffiffiffiffiffiffiffiffiffi ffi hu 02 i p , and ffiffiffiffiffiffiffiffiffi hv 02 t i p agree reasonably well with the experimental data. Some underprediction is seen in the mean (the maximum error being less than 10% for axial velocity) and especially in the rms radial velocities close to the jet exit and for 10 r 15 mm, i.e.,  in the flame region. This underprediction is related to the dilatation effects and a bimodal radial velocity PDF, which is discussed in detail in the next section. The mean velocities are predicted reasonably well also in nonreacting regions surrounding the flame, however, the rms values are underpredicted. The reasons for this might be related to the influences of dilatation effects.
Reynolds-and Favre-averaged progress-variable statistics were also extracted from the experiments by deducing the progress variable from bimodal velocity data (Furukawa et al., 2016). This approach could not be tested by the LES because that broadens the flame at least to a size of Δ > th . The transported propane-based progress variable is used here. The calculated radial variations of the Favre-averaged progress variable are compared with those from the experiment in Figure 11, where the solid curve corresponds to Pope model, Eq. (12), and the dashed curves, close to them, to Lilly and Colin et al. models,Eqs. (11) and (13). Error bars are shown for the experimental data at x ¼ 40 mm, and similar errors are expected for other heights. All SGS velocity models predict similar variations here, although there are small differences in the flame region ð10 r 15 mmÞ, lying within the error limits. The Reynolds averaged values are also shown for x ¼ 40 mm, computed using (Kolla and Swaminathan, 2010): where σ 2 ¼ σ 2 c;sgs þ σ 2 c;res is the Favre variance of c, the resolved variance being obtained as σ 2 c;res ¼ e c À he ci ð Þ 2 . Furukawa et al. (2016) were also able to extract the Reynolds average from their experimental data, although it exhibited error bars larger than those for the Favre average. The experimental data for Reynolds averages were reported only for x ¼ 40 mm. The difference between the computed Reynolds and Favre averages for other locations are similar to that shown for x ¼ 40 mm and one can estimate this Reynolds average using the Favre average and its variance (to be discussed next) through Eq. (15). The significant difference between Favre-and Reynolds-averaged values observed in the figure implies that the variance of the progress variable in both the experiment and the LES is large. Since the flamelets are quasi-laminar, this variance is predominantly produced by combustion, which operates at the SGS level for the LES. All of the agreements between LES predictions and the experiment seen in Figure 11 are good in that the differences are less than the experimental uncertainties, which are similar to those shown for x ¼ 40 mm (Furukawa et al., 2016).
The radial variations of resolved and SGS variances are shown in Figure 12. The results appearing here from the 23M grid will be discussed in a later section. The SGS variance is nearly four to five times the resolved part for the location in the flame closest to the jet exit, and the resolved part is smaller than SGS part for all locations shown here. The relative contribution of the resolved part increases with increasing height because of the associated relative increase in the turbulence level and flame thickness. These observations suggest that models for σ 2 c;sgs deduced from a passive-scalar scenario would be inappropriate. This is expected in low turbulence as such models, typically written as σ 2 c;sgs;model % AΔ 2 Ñe c Á Ñe c ð Þ (Pierce and Moin, 1998), require the most of the energy to be at super-grid scale. However, this model is shown to be inappropriate even for flames with high turbulence levels  as it does not consider the effect of the flame, which is typically at subgrid level in LES. Hence, its use for LES of premixed combustion should be discouraged. Furukawa et al. (2016) also extracted from their data an average turbulent-flame central contour based on a Favre-averaged progress variable of 0.5 and an average turbulent-flame thickness based on the radial extent of radial-velocity bimodality. Figure 13  those results with LES predictions, using he ci ¼ 0:5 and t ¼ 1=j@he ci=@rj max in the computations. Again, the dashed curves marked as 23M will be discussed in a later section. The agreements are seen to be good in that none of the differences between experimental and LES results exceed experimental uncertainties. Furukawa et al. (2016) refer to the bulge seen in the range of the experimental results at the top of Figure 13 as a "bubble," indicating that there is some experimental evidence for the associated double inflection point in the contour. They indicate, however, that the experimental uncertainties are large enough that one cannot be sure that this "bubble" exists. The LES does not exhibit this "bubble." The flame-brush thickness predicted by the LES is slightly thinner with the Lilly and Colin models, consistent with the results in Figure 11, but the differences are small, and predictions are in agreement with the experiment, within experimental error.
In summary, the results suggest that the flame statistics are not influenced by the choice of u 0 Δ modeling if the parameter c is chosen carefully to match the measured flame length and the SGS variance of c is obtained using its transport equation. Thus, the unstrained flamelet model can be used with careful selection of SGS modeling to capture these physical aspects. Figure 14 shows PDFs of velocities at five radial positions that span the turbulent flame brush, at a representative mid-flame height of 50 mm. The LES results shown are those for the resolved part only. A good agreement is seen between the LES results and measurements for the tangential velocity PDF. The computed axial velocity PDF also compares well the measurement except for the small shift in the location of the peak. This shift is consistent with the results in Figure 10. For the radial component, however, the experiment exhibits a strong bimodality in the center of the turbulent flame that is not reflected in the simulation results. The LES does show skewness of the radial-velocity PDFs in the flame region, with a high-velocity tail near the fresh mixture and probably a low-velocity tail near of burned gas sides of the flame brush. Since the mean dilatation produces higher radial velocities in the burned gas, the skewness may arise from the occasional penetration of reacted mixture into the position at low average completion of combustion, as is seen in the computed time-dependent behavior. The bimodality at the measurement positions where it is seen, however, cannot come from these resolved-scale components with the practical grids used here. It is a subgrid phenomenon on these grids (indeed small-scale phenomenon for ultra fine grids or DNS) for the corrugated-flamelet regime, that is not present in any of the three subgrid velocity models tested here.

Bimodality of radial velocity PDFs
This radial-component bimodality, coming from flame intermittency, that is, back and forth flamelet movement in the radial direction, as suggested by Furukawa and Williams (2003), is strong at r ¼ 11mm in the experiment. The computational results do not give this PDF because the filtered flame width is similar to the flame-brush width, as shown in Figure 15, and the flame-brush width is about 6 mm as shown in Figure 13b. This implies Further insight can be gleaned by considering the measured bimodal PDF as two Gaussians. The experimental data suggests that the means are μ L % 0:63 and μ R % 1:70, and the variances are σ 2 L % 0:026 and σ 2 R % 0:039, respectively, for the unburned and burned mixtures. The subscripts L and R denote the left (unburned mixture) and right (burned) parts of the PDF. The width of these two Gaussians, represented by the respective variances, are influenced by turbulence, and the separation of the two peaks are controlled by flame intermittency. The turbulence contribution is captured quite well in the simulations and this can be seen by writing the joint PDF of c and v r as Pðφ; ζÞ ¼ PðφjζÞPðζÞ, where φ is the sample space variable for v r . If one presumes the marginal PDF to be BML as has been done by Libby and Bray (1980) and Bray et al. (1985) then: where f ðζÞ is the burning mode part of the marginal PDF, PðζÞ, which can be evaluated by integrating the beta PDF used for reaction-rate closure between the limits of ζ ¼ 0:05 x  Figure 15. Isocontours of Favre-filtered (black lines) and Favre-averaged (green lines) progress variable at values 0.1 and 0.9 for the 5.3M grid. An enlarged view is shown for radial positions marked in Figure 14 in the inset. The Favre-filtered isoline is shown for two different times. and 0.95, but it is not easy to get σ f as it is related to the effects of combustion on SGS velocity fluctuation. The sum ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi σ 2 L þ σ 2 R p compares well with the computed rms value of 0.26 at r ¼ 11mm shown in Figure 16 (solid line). Nearly 78% of the total variance, hv 02 r i, at this radial location comes from the flame intermittency (evaluated using the distance between the peak locations in the experimental PDF), which suggests that the subgrid part must be σ 2 sgs ¼ σ 2 f % 0:23. Even with the gross assumption that the total u 0 Δ contributes to hv 02 r i, the computed variance is increased by only a small amount of 3:6 Â 10 À3 , which is negligible compared to σ 2 sgs . One can artificially include the intermittency effect by increasing the size of data sampling region to 6 mm as shown in Figure 16 (symboled line). Since this is artificial the corresponding PDF is not shown in Figure 14.
The maximum error for ffiffiffiffiffiffiffiffiffiffi hu 0 2i p =hUi is found to be less than 2% although the difference between measured and computed values in Figure 16 for the axial velocity rms seems high. There is some difference for the tangential component, which results from the flame intermittency. This contribution will be highly anisotropic in the corrugated-flamelet regime and thus including u 0 Δ fully for the radial component is not justifiable. Similarly, other corrections such as Deardorff correction (Deardorff, 1971;Ribault et al., 1999) for the rms values cannot be used because they generally require isotropy, which is not met in flames. Therefore, improved SGS velocity modeling strategies are required to capture the strong anisotropic effects of flame intermittency in practical LES.

Influence of numerical grid
Strictly, one should include the contribution of SGS fluctuations when constructing the velocity PDF, and these fluctuations are unavailable in LES framework unless the grid is refined to a level so that they becomes negligible. Indeed, one can refine the grid to have Δ < th , so the filtered flame width becomes comparable to the flamelet thickness to capture the bimodal behavior. The minimum filter width based on the numerical cell volume is about Δ þ ¼ 1:4 for the 5.2M grid and to have a filter width of Δ þ % 0:5, which is comparable to the experimental resolution length, one needs a grid with about 120M cells, which would be a very expensive calculation and thus it is not attempted here. However, another grid with about 23M cells is used to test the sensitivity of the statistics discussed above.
The first-order statistics obtained using the 23M grid are very similar to those obtained from the 5.4M grid (for example, seec h i ¼ 0:5 contour shown in Figure 13), and thus they are not shown here. The second-order statistics show some sensitivity to grid resolution as seen in Figure 16. Since the filter size decreases with an increase in grid resolution, the SGS velocity fluctuations are expected to decrease, which is seen in this figure. However, u 0 Δ values are still high compared to the resolved rms values for the three velocity components, and this is because the flame effects are responsible for the predominant part of the SGS kinetic energy in these regions. The resolved rms of axial velocity fluctuation for the 23M grid compares reasonably well with the experimental data, the tangential component does not show any sensitivity, and there is some improvement for the radial component. There is, however, a substantial difference between the computed and measured ffiffiffiffiffiffiffiffiffiffiffi v 0 r 2 q because of the missing subgrid part, which is related to the flame intermittency effects, as noted above. The improvements in the rms value for the radial component is reflected in the PDFs shown in Figure 14 for locations inside the flame brush. The SGS variance, σ 2 c;sgs , does not seem to be sensitive to the grid resolution and the resolved part increases with grid resolution as one would expect. The difference in the resolved part also becomes less sensitive to the grid further downstream. All of these points can be observed in results for the 23M grid shown in Figure 12. This suggests that the SGS combustion process is captured well by the reaction-rate closure used for this study. However, the influence of the flame intermittency and resulting acceleration effects on the SGS velocity fluctuation is not represented directly by the commonly used SGS velocity models. Thus, the experimentally observed bimodality of the radial velocity is not captured fully.

Summary and conclusions
A non-piloted Bunsen flame of a propane-air mixture with an equivalence ratio of 0.85 and low turbulence level, lying in the corrugated-flamelet regime, has been simulated using an unstrained-flamelet model for the filtered reaction rate. The turbulencechemistry interaction is driven by flame-induced effects, which is different from the shear-produced turbulence effects in standard piloted jet flames. Thus, the modeling of SGS velocity is expected to play a role, and three different models treating the dilation effects differently are considered. The SGS variance of a reaction-progress variable required for the presumed PDF is obtained by using its transport equation. The dissipation rate required for this equation is closed by applying an algebraic model developed by Dunstan et al. (2013). The grid sensitivity of the statistics is also assessed using 5.4M and 23M grids.
The computed statistics compare reasonably well with the experimental results, and the grid sensitivity of the statistics is found to be negligible. This suggests that the unstrained flamelet closure is robust and adequate for practical LES of premixed combustion. There is an underprediction of the variance of the radial component of velocity in the central part of the flame, where the respective PDF was observed to be bimodal in experiments. Detailed analyses of LES and experimental PDFs suggest that the turbulence effects are captured quite well in the LES, whereas the effects of flamelet intermittency is not. These strong and anisotropic effects in the corrugated-flamelet regime of combustion are not represented by the commonly used SGS velocity models, and improved models are required.