Large Eddy Simulation of Premixed Combustion: Sensitivity to Subgrid Scale Velocity Modeling

ABSTRACT An algebraic reaction rate closure involving filtered scalar dissipation rate of reaction progress variable is studied. The filtered scalar dissipation rate closure requires a model for sub-grid scale velocity, , which is estimated using four algebraic models and transported sub-grid scale kinetic energy. A priori analyses using direct numerical simulation (DNS) data show that the filtered dissipation rate, and thus the reaction rate closure, has some sensitivity to the model. The sensitivity of various statistics obtained from large eddy simulation (LES) of three piloted Bunsen flames of stoichiometric methane-air mixture to the modeling of is observed to be weaker compared to that for the DNS analysis. Moreover, analysis using transported sub-grid scale kinetic energy does not indicate a necessity to include flame-generated turbulence in the modeling of for the Bunsen flames in the thin reaction zones regime. The measured and computed flame brush structures are compared and studied and the algebraic closure for the filtered reaction rate is found to be quite good.


Introduction
The dynamics of large-scale turbulent eddies down to a cut-off scale are solved with models to mimic the influences of remaining sub-grid scales in large eddy simulation (LES). This methodology and its modeling are described in a number of earlier studies for non-reacting and reacting flows, see for example the books by Pope (2000) and Poinsot and Veynante (2005). The LES studies on turbulent combustion are reviewed by Pitsch (2006) and Gicquel et al. (2012). The combustion is usually a sub-grid scale (SGS) phenomenon requiring modeling and various modeling approaches used for premixed combustion are reviewed and summarized in earlier studies (Cant, 2011;Gicquel et al., 2012;Poinsot and Veynante, 2005;Swaminathan and Bray, 2011). These approaches can be broadly categorized into two classes, namely, flamelets and non-flamelets or geometrical and statistical (Gicquel et al., 2012). The geometrical category of flamelets includes thickened flame (Colin et al., 2000;De and Acharya, 2009), flame surface density or flame-wrinkling (see, for example, Boger et al. 1998;Chakraborty and Cant 2007;Chakraborty and Klein 2008;Gubba et al., 2012;Hawkes and Cant 2001;Knikker et al., 2004;Wang et al., 2012), and level-set or G equation (Moureau et al., 2008;Pitsch, 2005). The statistical category of flamelets includes approaches, such as EBU (eddy-break-up model), algebraic closure involving scalar dissipation rate (Butz et al., 2015;Dunstan et al., 2013;Gao et al., 2014;Langella et al., 2015;Ma et al., 2014), and presumed probability density function (PDF) methodology with laminar flamelets. The nonflamelets category includes transported PDF and conditional moment closure (CMC) methodologies. The CMC method (Klimenko and Bilger, 1999) was tested rigorously for RANS simulations of premixed flames (Amzin and Swaminathan, 2013;Amzin et al., 2012;Martin et al., 2003) but yet to be applied for LES of premixed combustion. Different PDF approaches have been used for LES of premixed combustion (Bulat et al., 2014;Dodoulas and Navarro-Martinez, 2013;Rowinski and Pope, 2013). Each of these methods has its merits and drawbacks, and the interest here is on the algebraic closure involving scalar dissipation rate discussed briefly in the next section [see Eq. (1)].
The SGS velocity, u 0 Δ , is an input to the algebraic closure as one shall see later and this velocity is typically modeled using Smagorinsky model (Smagorinsky, 1963) for the SGS eddy viscosity in earlier studies (Butz et al., 2015;Gao et al., 2014;Ma et al., 2014) or using a scalesimilarity model of Pope (2000). Since the combustion is a SGS phenomenon, its interaction with turbulence and other physical processes at the SGS level must be represented well in combustion modeling. This is specifically important for typical LES where the majority of the energy containing dynamic scales are resolved. In other words, the typical LES is seen here as a higher fidelity approach compared to URANS as expressed in classical views of LES (Deardorff, 1974;Lilly, 1966Lilly, , 1967Schumann, 1975;Smagorinsky, 1963) rather than as a degenerate from direct numerical simulation (DNS). The flamelet is usually thinner than the smallest scale resolved in such LES and SGS combustion model must represent relevant physical processes and their interactions at the SGS level quite accurately and robustly. It is well known that these processes are small-scale turbulence, chemical reactions or heat release, molecular diffusion, and their close interactions with one another in premixed combustion. Thus, the modeling of u 0 Δ is expected to play an important role and the sensitivity of the algebraic reaction rate closure to this modeling is yet to be explored. Specifically, the adequacy of non-reacting flow modeling of u 0 Δ for reacting flows is unclear. Hence, the objective of this study is to conduct this sensitivity analysis using DNS data and by performing LES of piloted Bunsen premixed flames of Chen et al. (1996), which was investigated in many earlier studies using RANS (Amzin et al., 2012;Herrmann, 2006;Kolla and Swaminathan, 2010b;Lindstedt and Vaos, 2006;Prasad and Gore, 1999;Salehi and Bushe, 2010;Salehi et al., 2012;Heinz, 2008, 2010) and LES (De and Acharya, 2009;Dodoulas and Navarro-Martinez, 2013;Pitsch and de Lagneste, 2002;Wang et al., 2011;Yilmaz et al., 2010) paradigms employing various combustion modeling approaches.
This article is organized as follows. The algebraic closure is discussed in Section 2 and the u 0 Δ models investigated here are presented in Section 3. The experimental flames and their computational modeling are discussed in Sections 4 and 5, respectively. A priori analyses using DNS data is discussed in Section 6.1 and LES results are presented in Section 6.2. The main findings of this study are summarized in the final section.

Algebraic closure for filtered reaction rate
The direct relationship between averaged reaction rate and scalar dissipation rate in high Damköhler number flames was shown by Bray (1979Bray ( , 1980 for RANS methodology, which was also shown to hold for low Damköhler number combustion (Chakraborty and Cant, 2011;Chakraborty and Swaminathan, 2011). Dunstan et al. (2013) demonstrated its applicability for LES, which is supported by subsequent analyses using DNS  and LES (Butz et al., 2015;Ma et al., 2014). This statistical model is: where the overbar denotes LES filtering and the Favre filtered scalar dissipation rate of a reaction progress variable c is e N c ¼ ρα Ñc Á Ñc ð Þ=ρ. The model parameter is where f(c) is the burning mode PDF of c and the subscript 'L' refers to planar laminar flame (Bray, 1979(Bray, , 1980. These integrals can be evaluated using laminar flame results and thus C m ' 0:83, obtained using DNS data, is a thermochemical parameter.
The filtered dissipation rate is the sum of resolved and unresolved contributions: " ρ e N c ¼ ρ αðÑe c Á Ñe cÞ þ ρe ε c ; with the unresolved part modeled as: where F ¼ 1 À exp Àθ 5 Δ þ ð Þ with θ 5 ¼ 0:75 and Δ þ ¼ Δ=δ th is the normalized filter width and the factor F ensures that e N c ! N c when Δ þ ! 0 (Dunstan et al., 2013). The planar laminar flame speed and its thermal thickness are denoted as s L and δ th , respectively. The SGS velocity scale, u 0 Δ , is to be modeled and the sub-grid 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 is the SGS flow time scale, which is related to u 0 Δ and Δ. The SGS kinetic energy is K and its dissipation rate is ε sgs . The heat release parameter is defined as τ ¼ ðT ad À T u Þ=T u using the adiabatic flame temperature T ad and reactant temperature T u . The other parameters signifying the influences of thermochemical and turbulence processes and their interplay are Kolla et al., 2009;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 past studies Swaminathan, 2013, 2014;Kolla et al., 2009;Kolla and Swaminathan, 2010b;Ruan et al., 2014;Swaminathan et al., 2012) showed that the above parameterization is not arbitrary and various terms in Eq. (2) are closely related to certain physical processes involved in the transport of scalar dissipation rate Kolla et al., 2009). The terms involving (K c s L =δ th ) and C 3 À τC 4 Da Δ ð Þ ð u 0 Δ =ΔÞ respectively arise due to fluctuating dilatation and strain rate resulting from competing effects of turbulence and heat release. Hence, these parameters are not tuneable for given reactant mixture and turbulence conditions. The termcð1 ÀcÞ=β c comes from combined influences of flame front curvature effects induced by turbulence, chemical reaction, and molecular dissipation processes Chakraborty and Swaminathan, 2007). Although a reasonably robust value of β c ¼ 6:7 was established in past RANS calculations of various premixed flames (Ahmed and Swaminathan, 2013;Amzin and Swaminathan, 2013;Amzin et al., 2012;Kolla et al., 2009;Ruan et al., 2014;Swaminathan et al., 2012), this value cannot be used for LES because the processes influencing β c are effected by SGS turbulence. Thus, it can be evaluated dynamically as suggested by Langella et al. (2015) and Gao et al. (2014). However, a representative value guided by DNS analysis as discussed in Section 6.1 is used here to investigate the influences of u 0 Δ modeling on the reaction rate closure and LES statistics.
3. Modeling of u 0 Δ The u 0 Δ influences the SGS SDR,ε c , in two ways-a direct influence is through the linear dependence in Eq. (2) and additional nonlinear effects come through Ka Δ dependence of C 3 and C 4 . The values ofε c obtained using Eq. (2) is accurate if the DNS values of u 0 Δ are used as one shall see later. The SGS SDR directly influences the reaction rate through Eq. (1) and thus the turbulent flame speed. Hence, one can expect that the later quantity will be predicted well when u 0 Δ is modeled correctly. Furthermore, the role of u 0 Δ and its interaction with flame may differ depending on the combustion regime. For example, the thin flame front in the corrugated-flamelets regime can produce velocity fluctuations at SGS level through flame intermittencies and baroclinic torque mechanisms, and these mechanisms may become weaker compared to the turbulence in thin reaction zones regime combustion and beyond. Also, the role and nature of u 0 Δ in practical combustors with large turbulent Reynolds number are open questions. Thus, the modeling of u 0 Δ needs to be considered carefully and the following models are tested here.
(1) A first model based on a shear stress related closure proposed by Lilly (1967) is: where ν t ¼ C s Δ 2 ffiffiffiffiffiffiffiffiffiffiffi 2 e S ij e S ij q is the turbulent viscosity obtained using localized dynamic model (Lilly, 1992). The model constant is C L ¼ 0:094.
(2) A second model based on scale-similarity for velocity is written as (Pope, 2000): where C q is of order unity and it is 1 for this study.ũ i is the velocity field obtained using a Gaussian test filter.
(3) A third model is based on scale-similarity for kinetic energy proposed by Bardina et al. (1980). It is written in Galilean invariant form as (Ferziger, 1997): with C 0 q ' 0:4 deduced from Liu et al. (1994). (4) A variant of the first model can be deduced by considering a constant energy transfer down to Δ scale as shown in Appendix A. This model is: where C s is the Smagorinsky model constant obtained dynamically. (5) The u 0 Δ can also be obtained from SGS kinetic energy, p by assuming isotropy for the SGS velocity components. This assumption is less likely to hold for corrugated-flamelets combustion because of stronger acceleration across the flame front compared to those along the front. However, this approximation may hold reasonably well for combustion at large Reynolds number. A transport equation using the above definition for K can be written as: The first term is the substantial derivative of K, the second term is the molecular diffusion, the third term is the work done by the sub-grid stress tensor, τ R ij ¼ τ r ij þ 2Kδ ij =3 with τ r ij as defined in Section 5, the fourth term is the turbulent transport of the fifth term is the dissipation of K, which is modeled as k % C " ρ K 1:5 =Δ, and the last term is the pressure-work related term. This modeled equation is similar to that given by Chai and Mahesh (2012) for compressible flow except for the dilatational dissipation resulting from compressibility. The pressure related term is given by: which is modeled as: following the practice in RANS (Kolla and Swaminathan, 2010a) to include dilatation resulting from heat release effects and Pr t % 1 is used here. The modeling of Π and ε k is included systematically as listed in Table 1 to study their influences on LES statistics. The case K1 excludes the pressure related term, Π, and a static value of 0.916 is used for C , whereas the case K2 includes Π. The value for C is determined dynamically in the case K3. The above modeled transport equation for K is ad hoc at this time and improvements can be made by analyzing its individual terms using DNS data and revising the above models. Although such an analysis would be worthwhile, our interest here is on algebraic closures for u 0 Δ and its sensitivity on LES statistics for the algebraic reaction rate closure.
Common u 0 Δ models, as those introduced earlier in this section, are derived originally for nonreacting flows and thus it is of interest to assess their adequacy for combusting flows in which flame-generated turbulence may affect the SGS turbulence. Colin et al. (2000) suggested that the u 0 Δ models involving v t (Smagorinsky model) are dominated by thermal expansion in the absence of turbulence and a scale similarity model as in Eq. (4) can be used by taking the rotational part of this u 0 Δ . A revised form of this model as written by Wang et al. (2011) is u 0 Δ ¼ 2Δ 3 x jÑ 2 ðÑ Â e uÞjð0:1Δ=Δ x Þ 1=3 , where the mesh size estimated from the local numerical cell volume is Δ x . This model was derived for use in conjunction with thickened flame model and its efficiency function requiring turbulence level of unburned mixture that would have existed locally in the absence of heat release rate. It may not be entirely appropriate to use this model with Eq. (2) for the following reason. The local gradients of c will be influenced by combustion even in the absence of turbulence and the local values of Ñc result from a fine balance among turbulent transport, heat release, advection, and molecular diffusion. All of these processes will be influenced by the heat release effects and the local velocity will be affected by the acceleration induced by dilatation resulting from the flame. Thus, it would be inappropriate to exclude the dilatational effects on u 0 Δ required for Eq. (2) and so the model of Colin et al. (2000) is not considered further in this study. On the other hand, it is important to evaluate if the common u 0 Δ models can capture the correct SGS turbulence-flame interaction and its effect on the filtered reaction rate. This is attempted in Section 6.2 for flames in thin reaction zones regime and further insights are in  for the corrugated flamelets regime combustion.
The closure in Eq. (1) is statistical and so one must not insist that _ ω ! _ ω in the limit of Δ ! 0. Also, it would not be physically meaningful to impose that the propagation speed of a resolved flame to be s L when Δ ! 0 for a statistical closure and such expectations are acceptable for geometrical category. However, if one likes to impose the above condition then Eq. (1) must be written as r ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffi where a is a parameter of order unity and _ ω is the flamelet reaction rate. Gao et al. (2014) showed that Eq. (1) recovered DNS results for Δ + > 1 and, hence, Eq. (1) is used here.

Experimental cases
Three piloted stoichiometric methane-air Bunsen flames of Chen et al. (1996) considered in many past RANS (Amzin et al., 2012;Herrmann, 2006;Kolla and Swaminathan, 2010b;Lindstedt and Vaos, 2006;Prasad and Gore, 1999;Salehi and Bushe, 2010;Salehi et al., 2012;Heinz, 2008, 2010) and LES (De and Acharya, 2009;Dodoulas and Navarro-Martinez, 2013;Langella et al., 2015;Pitsch and de Lagneste, 2002;Wang et al., 2011;Yilmaz et al., 2010) studies are chosen for this investigation. The reactant jet with a diameter of D = 12 mm issues into quiescent air and the flame is stabilized by laminar pilot flames of the same mixture. The pilot ring of a diameter D p = 68 mm is water cooled and thus its burned mixture is sub-adiabatic. The turbulence in the reacting region is shear driven as there is no turbulence generating device in the reactant flow path. Three flames (F1, F2, and F3) with bulk mean velocity of U b = 65, 50, and 35 m/s have a Reynolds number based on D and U b of 52,000, 40,000, and 24,000, respectively, and this variation is achieved by changing only U b . The combustion conditions in these flames, as shown in Figure 1, are in the thin reaction zones regime and these conditions are calculated using the centerline values of turbulence RMS velocity, u 0 , and its integral length scale Λ = 2.4 mm near the nozzle exit reported by Chen et al. (1996). The Zeldovich thickness is δ ¼ ν u =s L , where v u is the kinematic viscosity of the reactant mixture. The turbulent Reynolds number is defined as Re t ¼ u 0 Λ=ν u . The symbols D1 to D7 in Figure 1 denote the conditions of DNS flames used for a priori testing discussed in Section 6.1.

Simulation detail
The Favre filtered transport equations for conservation of mass and momentum are solved along with additional filtered equations required for combustion. The sub-grid stresses are modeled using the dynamic Smagorinsky model (Germano et al., 1991;Lilly, 1992) when u 0 Δ is modeled using one of Eqs.
(3) to (6) for the combustion part. A test-filter of size b Δ % 2Δ is employed for all dynamic procedures used here following the common practice, and the filter width is computed using the computational cell volume, is modeled using K then the sub-grid stresses are modeled as τ r ij ¼ À2C k Δρ ffiffiffi ffi K p ð e S ij À e S kk =3Þ, with C k computed dynamically (Germano et al., 1991). The additional equations are for Favre filtered progress variable,c, and total enthalpy, e h. The instantaneous c can be defined using absolute temperature, T, or a species mass fraction. Here, it is defined using the fuel mass fraction, Y f , as c ¼ 1 À Y f =Y f ;u so that it takes a value of 0 and 1 in the unburned and burned mixtures, respectively. This choice avoids a spurious flame that will appear numerically in the mixing layer between the pilot and coflow if c is defined using T or any other species mass fraction. The transport equation forc is: and the filtered reaction rate is closed using the models described in Section 2. One can also use a transport equation for e Y f instead ofc as there is no particular advantage in using either e Y f orc for premixed combustion. The governing equation for e h is: Including this equation allows one to handle the sub-adiabatic pilot stream discussed in the previous section. The sub-grid scalar fluxes appearing in Eqs. (10) and (11) are computed using dynamic Schmidt number approach (Lilly, 1992). The filtered temperature e T is obtained from e h using: where T 0 ¼ 298:13 K and e C p;mix is the specific heat capacity at constant pressure for the mixture, which depends on temperature as described by Ruan et al. (2014). The fluid mixture density is computed using the state equation ρ ¼ p e W mix =ðR e TÞ; where " p is the filtered pressure, e W mix is the Favre-filtered molecular mass of the mixture and R is the universal gas constant. The quantities e C p;mix , f Δh 0 f ;mix , and e W mix are specified through a lookup table constructed using planar unstrained laminar premixed flame of stoichiometric CH4-air mixture using CHEMKIN and GRI 3.0 mechanism with mixture averaged diffusivity formulation. The laminar solution is then convoluted using a Gaussian filter kernel with a width equal to the cube root of the smallest computational cell volume in LES. This table contains e C p;reac , f Δh 0 f ;reac , e W reac ; and filtered mass fractions of various species as a function ofc and is used during the LES. The LES statistics obtained using three times larger width for this filter kernel showed a negligible sensitivity of the statistics to this kernel width. Also, the use of an unstrained flamelet excludes the possible effects of fluid dynamic strain on filtered flame in LES, which will be explored in a future study.
A transport equation, similar to Eq. (11), for a passive fluid marker, e Ψ, is included to account for the mixing or dilution of burned mixture with the entrained air. These effects were assessed to be important (Kolla and Swaminathan, 2010b) to capture the averaged values of various species mass fractions for x ! 4:5D, specifically for the F1 flame having the highest Reynolds number and thus this is included in this study. Briefly, the influence of entrained air is included using a mixing rule: e Φ mix ¼ e Ψ e Φ reac þ ð1 À e ΨÞΦ air (Kolla and Swaminathan, 2010b) where e Φ is a generic value representing e C p;mix or f Δh 0 f ;mix or e W mix or filtered mass fractions. The subscript 'reac' denotes that these values are taken from the lookup table for localc and the subscript 'air' denotes their values in the air stream.

Numerical grid, boundary, and initial conditions
The filtered conservation and combustion modeling equations are solved using PRECISE-MB, which is based on a finite volume methodology (Anand et al., 1999). The spatial derivatives are discretized using a second-order central differencing scheme and the velocity-pressure coupling is achieved using the SIMPLEC algorithm (Doormaal and Raithby, 1984). The discretized equations are time advanced using a second-order accurate scheme with a constant time step and the CFL number is kept below 0.3 everywhere inside the computational domain spanning 40D in the axial and radial directions as shown in Figure 2, giving a total computational volume of π(40D) 3 /4. This volume is discretized using a structured multi-block grid having nonuniform numerical cells. These cells are fine near the burner exit and grow gradually in the downstream and radial directions. A coarse grid having 22 cells for D and about 4 cells within Λ measured at the jet exit gives a total of about 1.5 million cells for the computational volume considered. This grid has 404 cells in the streamwise direction along the centerline. Increasing the cell count to 32 for D and 6 for Λ keeping other grid parameters to be almost the same yields about 4.2 million cells in total. These two grids are used to assess the grid sensitivity of the LES statistics computed here. The coarse grid satisfies the 80% turbulent energy criterion of Pope (2000) and has the smallest normalized filter width Δ þ min ¼ minðΔÞ=s L % 1:3 and the fine grid has Δ þ min % 0:8. Figure 2. A schematic of the experimental (Chen et al., 1996) and computational setup of piloted stoichiometric methane-air Bunsen flames.

Boundary and initial conditions
The measured values at jet exit are used to specify the jet velocity and there is no velocity fluctuation specified at the exit as the turbulence is predominantly shear driven. A uniform velocity of 1.5m/s is specified for the pilot stream on the basis of the volumetric flow rate in the experiments. A small velocity of 0.2 m/s is assigned to the coflowing air to mimic the air entrainment and a 25% change to this velocity is observed to influence the LES statistics negligibly. The pilot temperature is unspecified in the experimental study and values ranging from 1785 to 2248 K were used in past studies (Amzin et al., 2012;De and Acharya, 2009;Dodoulas and Navarro-Martinez, 2013;Herrmann, 2006;Kolla and Swaminathan, 2010b;Lindstedt and Vaos, 2006;Pitsch and de Lagneste, 2002;Prasad and Gore, 1999;Salehi and Bushe, 2010;Salehi et al., 2012;Heinz, 2008, 2010;Wang et al., 2011;Yilmaz et al., 2010) suggesting that the heat loss to the pilot burner varies from 0 (De and Acharya, 2009) to 34% (Dodoulas and Navarro-Martinez, 2013;Lindstedt and Vaos, 2006). A value of 17% was used by Herrmann (2006) and 20% heat loss was assumed by Pitsch and de Lagneste (2002) and Wang et al. (2011). For this study, it is taken to be 16% following Amzin et al. (2012) and Kolla and Swaminathan (2010b), which gives 1950 K for the pilot temperature. These past studies showed that this uncertainty influences the temperature only close to the jet exit (x ≤ 3D), which is also confirmed here by changing this temperature over a range of about 200 K. The filtered progress variable is 0 in the jet exit and 1 for the pilot and coflowing streams. The passive fluid marker, e Ψ, is 1 in the jet and pilot fluids and 0 in the coflowing air. The lateral boundaries are specified to be slip walls and the outlet has zero gradient in the streamwise direction for all of the variables.
The simulations are run using 96 cores running at 2.60 GHz, Intel Sandy Bridge E5-2670 processors on Darwin cluster at Cambridge University for about 0.15s to 0.25s of real flow time in total. The data are collected for about 32, 26, and 22 flow-through times for F1, F2, and F3 flames, respectively. The flow-through time is defined using the computational domain length and the respective U b at the jet exit. These sampling times are substantially larger than what has been used typically in past studies for these flames and these simulations in the above cluster took about a day for 1.5M and two days for 4.2M grids on the wall clock.

Results and discussion
The details of DNS data used for a priori analyses are elaborated elsewhere ) and thus only essential information required here is presented briefly. The DNS considered freely propagating turbulent premixed flames for a range of turbulence and thermochemical conditions given in Table 2 along with other relevant parameters defined as u 0þ ¼ u 0 =s L , Λ þ ¼ Λ=δ th . The Damköhler and Karlovitz numbers are Da ¼ Λ þ =u 0þ and Ka ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi u 0þ 3=Λ þ p , respectively. The initial values of these parameters are given in Table 2 and the combustion conditions of DNS flames shown in Figure 1 are representative of those in the experimental flames. A single irreversible reaction with kinetic parameters representative of preheated methane-air combustion was used.
The instantaneous quantities of interest are filtered using a Gaussian kernel with filter widths in the range of 0:4 Δ þ 2:8 for a priori validation to be discussed next. The test filter width used for this analysis is b Δ þ ¼ 2Δ þ as for the LES. More detail of the DNS data and their post-processing techniques can be found in Gao et al. (2014).
6.1. A priori analysis The various quantities required to get u 0 Δ through Eqs.
(3) to (6) can also be extracted from the DNS data and thus these models can be verified by comparing their estimates to u 0 Δ;DNS . These estimates are denoted as u 0 Δ;model , or u 0 Δ;Eq:4 if Eq. (4) is used for modeling in the following discussion. Figure 3 compares the DNS results with modeled values across the filtered flame of the D5 case for two filter widths, Δ þ ¼ 0:8 and 2.8. These values of u 0 Δ;DNS , shown in this figure, are obtained by averaging u 0 Δ conditional onc value. The D5 flame is shown because its combustion condition is akin to those in the experimental flames. The variation of u 0 Δ;DNS shows that the SGS kinetic energy decreases across the flame and this trend does not depend on the filter widths used. As one would expect, the SGS kinetic energy level decreases as the filter width is decreased. All of these variations are captured by the standard models. The model in Eq. (6) predicts an increase in u 0 Δ , which is expected for flames in the corrugated or wrinkled flamelets regime, but D5 is in the thin reaction zones regime as shown in Figure 1. Although the other models underpredict u 0 Δ , with Eq. (5) giving the lowest value, their estimates are of the same order as the DNS result. The values estimated using Eq. (4) with C q = 2.5 agree very well with the DNS results as shown in this figure, which suggests that the model parameter does not depend on the filter width. Furthermore, this  analysis suggests that the modeling constants derived for nonreacting flows need to be revised for reacting flows, specifically for flames in the corrugated-flamelets regime . However, no such revision may be required for flames in the thin reaction zones regime because of stronger influence from turbulence compared to flame-induced effects on u 0 Δ .
6.1.2. Effect of u 0 Δ onε c and _ ω closures The filtered dissipation can be rewritten as: after averaging it over the DNS computational volume. This averaging operation is denoted by hÁ Á Ái v in the above equation. The typical variation of N with Δ + obtained from the DNS data is shown in Figure 4a. As one would expect N approaches 1 when Δ + goes towards zero implying that the SGS scalar dissipation rate,ε c , is negligible. As Δ + increases, the SGS contribution becomes substantial compared to the resolved component as has been observed in LES (Butz et al., 2015). This SGS contribution becomes about 50% and 2.5 times of the resolved part for Δ + = 1.2 and Δ + = 2.8, respectively. Thus,ε c closure plays an important role for the filtered reaction rate calculation for large filter widths. The behavior of N seen in Figure 4a is also suggestive of a power-law behavior as noted by Dunstan et al. (2013) and Gao et al. (2014), and further details can be found in these references. Figure 4a also compares the DNS values with estimates obtained using various u 0 Δ models discussed in Section 2. As one would expect, N values obtained using u 0 Δ;DNS agree well with its values obtained directly from the DNS, especially for large filter widths,  demonstrating that Eq.
(2) and its parameterization are very good as has been observed by Dunstan et al. (2013) and Gao et al. (2014). The error is observed to be about 24% for Δ þ % 2:4 when Eqs.
(3) or (4) are used and the error increases to about 27% for the scale similarity model of Bardina et al. (1980) in Eq. (5) when Δ þ % 2:4. The maximum error for the model based on the constant energy transfer hypothesis in Eq. (6) is about 10% suggesting that this model works well in an overall sense. A similar behavior is observed for other DNS flames listed in Table 2 also. The influence of u 0 Δ modeling on the local behavior ofε þ c ¼ε c δ th =s L is shown in Figure 4b along with the values extracted directly from the DNS data. Although Figure 4a showed a good match between the DNS and modeled values of N , there are some differences inε þ c when u 0 Δ;DNS is used. The local variations obtained using Eq. (6) for u 0 Δ follows closely the variation obtained using u 0 Δ;DNS . The other three models for SGS velocity underestimateε c significantly, as seen in Figure 4b. This is consistent with the results shown in Figure 3. It is worth noting that C q = 1 is used for Figure 4b and employing C q = 2.5 would improve the model comparison, which can be verified quite easily using Eq. (2) with u 0 Δ from Figure 3. However, the magnitude of this constant is observed not to influence LES results significantly.
The equation (2) can be written as: where B ¼ Fe cð1 À e cÞ=β c , after using K c ¼ 0:79 τ. A direct influence of u 0þ Δ is through the linear term, which is also compounded by nonlinear effects coming through ka Δ involved in C 3 and C 4 . When Ka Δ ( 1 the contributions from the second and third terms in Eq. (14) becomes very small. This behavior is expected because u 0þ Δ becomes very small under this condition and the unresolved dissipation rate is controlled by the thermochemical process. The SGS turbulence related term dominates when Ka Δ ) 1 because C 3 % 1:5 and C 4 % 0. The contributions from the terms involving C 3 and C 4 in Eq. (14) are of order unity for finite values of Ka Δ . The values of these two terms become about 0.38 and 0.13 for u 0þ Δ ¼ 12 and they change to 0.10 and 0.23, respectively, for u 0þ Δ ¼ 4 when Δ þ ¼ 2:8. These estimates can change in a nonlinear manner with Δ þ . However, the SGS turbulence level is linked to the filter width and so one cannot change these two parameters independently. Thus, the SGS velocity scale is expected to play a role in typical LES irrespective of the combustion regime. Furthermore, u 0 Δ model will also influence the fluid dynamic field through the effects of heat release. Thus, one must be cautious in drawing conclusions from DNS analysis and extending them to LES without further testing.
Equation (1) can be used to show that the difference between the DNS and modeled filtered reaction rate is _ ω DNS À _ ω model À Á , e ε c;DNS À e ε c;model À Á , which also holds for the volume averaged values. Thus, the differences between modeled and DNS values observed in Figure 4 apply for the filtered reaction rate also and a direct comparison is shown by Gao et al. (2014).
The above a priori assessments are conducted using β c = 4, see Eqs.
(2) and (14), for the case D5. This parameter takes a similar value [see Table 2 of Gao et al. (2014)] for other DNS cases and SGS dissipation rate obtained using this β c obeys the realizability condition,ε c ! 0. The β c was shown to be influenced by τ  and thus its value for the experimental flames having τ = 6.4 can be quite different. Indeed, an extensive parametric testing using LES of F1, F2, and F3 flames suggests that β c % 7:5 for these flames, which is used for the LES results discussed next.

A posteriori analysis
The computational model used is a ssessed first using nonreacting flow results. The measurements of mean velocity and turbulent kinetic energy for a nonreacting jet are reported in Chen et al. (1996) for the conditions of the F2 flame. The LES results are time averaged first over the sampling period and then averaged in azimuthal direction in order to get radial variation of mean quantities at various axial positions, as has been done by Langella et al. (2015). This averaging procedure, which is also density weighted when required (Poinsot and Veynante, 2005), is denoted by using hÁ Á Ái in the following discussion. Figure 5 compares the computed and measured radial variations of averaged streamwise velocity and turbulent kinetic energy (TKE), hki, for five axial positions. The velocity is normalized using U b ¼ 50m=s and the TKE is normalized using k 0 ¼ 10:8m 2 =s 2 . The TKE is obtained as hki ¼ 0: are obtained using Eq. (3) for the cold flow. The model in Eq. (6) was also used in the analysis but no significant difference was observed in the results.
The computed averaged streamwise velocity agrees very well with measurements as one sees in Figure 5 and the differences between 1.5M and 4.2M grid results are negligibly small. There are some differences between the computed and measured values of hki for r ≤ 0.55D, which decrease as one moves in the downstream direction (i.e., x > 6.5D). As seen in Figure 5, the turbulence levels represented by hki in regions of combustion (one shall see later that this would be r > 0.6D) are captured quite well. There is no turbulence specified at the jet exit and feeding some inlet turbulence improved the comparison of hki for x = 2.5D, but no significant changes were observed for other downstream locations because the turbulence is dominated by the shear generation mechanism in this flow configuration. It is clear from Figure 5 that the 1.5M grid is adequate to capture the flow dynamics and, hence, it is used for further analysis. Nevertheless, 4.2M grid is also used when required.

Sensitivity of filtered reaction rate
The normalized minimum filter sizes, Δ þ min , implied by the minimum cell volume are 1.3 and 0.8 for the 1.5M and 4.2M grids, respectively. These sizes become 2.8 and 1.7, respectively, if one uses the diagonal distance of the smallest grid. Since a spatially nonuniform grid is used, this width varies and typical distribution obtained using the local cell volume is shown in Figure 6 as a histogram for both the grids. This histogram is constructed using the grid information collected from the entire and two subparts of the computational domain having 0 x 20D and 0 r 10D, and 0 x 11D and 0 r 2D, where r ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffi . The most likely width of the implied filter is about 2δ th and 0.9δ th for the 1.5M and 4.2M grids, respectively, for the smallest subdomain. When the subdomain size is increased the most likely and averaged widths are increased substantially, especially for the 1.5M grid. The most likely value of Δ + is about 4 and the average value is about 5 to 6. There is a further increase in these values when the whole domain is considered for this analysis. The larger grid sizes near the lateral and outflow boundaries increase the count around Δ þ % 14 and about 10, respectively, in Figures 6a and 6b. These Δ + are substantially larger than what was typically used in past LES studies of these flames.
The influence of u 0 Δ modeling on the filtered reaction rate computed in the LES is shown in Figure 7 for five models used in this study, as it is important to understand this sensitivity because the averaged reaction rate, thus the turbulent flame speed, will depend on this modeling for the algebraic reaction rate closure, Eq. (1), used in this study. The reaction rate is normalized using ρ u s L =δ th and the color map is shown for _ ω þ > 0:01 for an arbitrarily chosen time. The algebraic models in Eqs.
(3) to (5) give relatively thinner fronts compared to that obtained using Eq. (6) or the K equation. The flame length is more or less the same for all of the models and there is a slightly longer flame for the scale-similarity model in Eq. (4). Lilly's model in Eq. (3) gives the shortest flame. Also, the peak reaction rate near the burner exit is seen only for Eq. (4) and K transport equation. The peak reaction rate predicted by the other models is relatively smaller. The averaged normalized reaction rate shown in Figure 8 is almost similar for the various SGS velocity models. The difference is seen only in the region close to the jet exit (x ≤ 2.5D). The differences observed in Figures 7 and 8 suggest some sensitivity to u 0 Δ modeling. It is not possible to do a quantitative analysis using reaction rate because this quantity is not reported in Chen et al. (1996) as it is not easy to measure it. However, one can make this assessment for overall performance of these models using statistics gathered from LES and measurements as has been done next.

Sensitivity of statistics
A priori analyses presented in Section 6.1 suggest that the sensitivity of LES statistics to u 0 Δ modeling must be investigated. The flame F2 simulated using the 4.2M grid is chosen for this assessment because the results for the 1.5M grid is better than what is shown here and the results for the F1 and F3 flames are very similar. Figures 9 and 10 respectively show the sensitivities of fluid dynamic and scalar quantities to u 0 Δ modeling. The radial variation of hUi=U b , is shown in Figure 9 for five streamwise locations. The averaged normalized TKE is also shown in this figure. The contribution of SGS kinetic energy to hki is obtained using the same u 0 Δ model consistently for the four algebraic cases discussed in Section 2. The sensitivity of hUi to the SGS velocity modeling is negligible whereas hki shows a weak sensitivity. All of the models overestimate hki for r ≤ 0.8D when x ≥ 8.5D and their estimates for the outer region, r > D, compare reasonably well with the measured values. The sensitivities of flame-related quantities, such as the fuel mass fraction and normalized temperature, are shown in Figure 10. The estimates of hY CH4 i by the scale-similarity model in Eq. (4) are closer to the measured values. Although a similar behavior is seen for the normalized temperature there are some large underestimates for r > 0.7D, which is related to the performance of filtered reaction rate closure in Eq. (1) on 4.2M grid having Δ + < 1. The DNS analysis discussed in Section 6.1 suggests that this reaction rate closure underestimates the reaction rate for Δ < δ th, which is also observed in LES (Langella et al., 2015). The results for the 1.5M grid discussed in Appendix B compares well with experimental data.
Typical results obtained using the modeled K equation, see Eq. (7), is shown in Figure 11 for the flame F2 along with the results obtained using Eq. (4). The influence of K-based u 0 Δ modeling on hUi=U b is negligible. The detail of the models K1, K2, and K3 are tabulated in Table 1. A similar behavior is observed for hki=k 0 and the averaged fuel mass fraction shown respectively in the second and third rows of Figure 11. Although there are some small differences in the near field for the normalized averaged temperature, the differences become negligible for downstream locations as seen in the bottom row of Figure 11. From these results, one concludes that the results of K-based modeling of u 0 Δ are very similar to those of the scale-similarity model in Eq. (4). Thus, it does not seem to justify the additional computational effort required for K transport equation. Furthermore, a negligible difference observed for the models K2 and K1 or K3 and K1 suggests that flame-generated turbulence included through Π in Eqs. (7) and (9) does not significantly influence the filtered reaction rate, thus the flame speed. Hence, the influence of flame on u 0 Δ modeling can be neglected for the flames studied here. This is because these flames are in the thin reaction zones regime where the relatively high level of turbulence may mask the effects of flame-generated turbulence and this may not be so for flames in the corrugated flamelets regime as shown in . However, one must acknowledge that the various sub-modeling for κ equation used here are ad hoc and further investigation is required to improve these models or to develop new ones.
Based on the above assessments, the scale-similarity model in Eq. (4) seems appropriate for the conditions of the flames investigated here. The models in Eqs. (3) and (5) are less satisfactory compared to those in Eqs. (4) and (6). Furthermore, it is to be noted that SGS stress-based models consider only the deviatoric part of this stress tensor and not the trace. x/D = 2.5 Figure 9. Influence of SGS velocity modeling on averaged streamwise velocity, hUi, turbulent kinetic energy, hki for F2. Symbols () are experimental results (Chen et al., 1996), lines are from simulations using Eq. (3) ( ), Eq. (4) ( ), Eq. (5) ( ), and Eq. (6) ( ).
Hence, it is debatable if one can use the stress-related models to get u 0 Δ because the SGS kinetic energy per unit volume is related to the trace and so Eq. (4) is used for this study. It is possible to envisage that this model can be used for LES of flames over a wide range of combustion conditions. It is useful to understand the ability of the filtered reaction rate closure to capture the flame brush structure. Since the focus here is on the u 0 Δ modeling sensitivity, the flame brush structure comparison is discussed in Appendix B for the sake of completeness of this investigation. These results show that the flame brush structure is captured quite well using the 1.5M grid, which is smaller than what was used in past studies.

Summary and conclusion
Large eddy simulations of piloted Bunsen flames of stoichiometric methane-air mixture are conducted. The filtered reaction rate is modeled using an algebraic closure involving scalar dissipation rate of reaction progress variable requiring a model for SGS velocity, u 0 Δ . This x/D = 2.5 10.5 Figure 10. Influence of SGS velocity modeling on averaged fuel mass fraction and normalized temperature for F2. The legend is as for Figure 9.
closure belongs to the statistical category of flamelet modeling. Four algebraic models for u 0 Δ are investigated along with a modeled transport equation for the SGS kinetic energy. Although some sensitivity of the SGS scalar dissipation rate,ε c , and thus the filtered reaction rate, closure to the u 0 Δ model is observed in a priori testing discussed in Section 6.1, the LES statistics show a weak sensitivity. The good agreement of flame brush structure computed using the algebraic closure for the filtered reaction rate, _ ω, along with the scale-similarity model, u 0 Δ ' C q je u i À e e u i j, suggests that these models work well for the flame conditions investigated in this study and also that ad hoc models accounting for flame-generated turbulence may not be required when the level of turbulence is relatively high as for the thin reaction zones regimes, which can explain the good comparisons shown in this study. However, the algebraic closure for the reaction rate does not include the finite rate chemistry and flame stretch effects. This closure also has one adjustable parameter, β c , which can be evaluated dynamically as has been done by Langella et al. (2015) but a static value is used for this study to investigate the sensitivity to u 0 Δ model without further complexity. Also, the β c value used here is very close to the average of the values in the flame region obtained using the dynamic approach. The influences of flame stretch and finite rate chemistry can be included using either unstrained or strained flamelets as discussed in .

Acknowledgment
The authors gratefully acknowledge the support from EPSRC, Siemens, and Rolls-Royce.

Funding
This work is funded by the grant numbered EP/I027556/1.

Nomenclature
c progress variable C m combustion model parameter C p specific heat capacity at constant pressure (kJ/kg-K) C L , C q , C 0 q SGS velocity modeling constants kinematic viscosity (m 2 s −1 ) Å pressure-dilatation term in Eq. (7) ρ density (kg/m 3 ) σ 2 T variance of normalized T τ heat release parameter τ R residual stress tensor (kg/m-s 2 ) τ r anisotropic part of τ R Ψ passive fluid marker _ ω reaction rate (kg/m 3 -s)  Dodoulas and Navarro-Martinez (2013) used a computational domain of size 15D × 5D employing about 510k and 225k grids in the physical space and their results on their largest grid, and with 16 stochastic fields, is used here for comparison. Further detail on the combustion modeling and numerical methods used in these studies can be found in these references and the results available there are used for comparison. The radial variation of U h i=U b for various streamwise locations are shown in Figure B1 along with experimental measurements. Only two representative results from past studies are shown here because those results are very similar to one another. The current results are shown for the 1.5M grid and the results of the 4.2M grid are shown only for the F2 flame because the relative behavior is very similar for the other two flames. The difference in U h i=U b obtained using the 1.5M and 4.2M grids is negligibly small and the values obtained using the 1.5M grid compare well with the measurements also. These results agree well with those from past studies; however, a close examination of Figure B1 suggests that the agreement with the experimental data is better for the results of Wang et al. (2011), which used the 8.5M unstructured grid. Nevertheless, the difference among the numerical results and the data is observed to be small.

Superscripts
The radial variation of k h i=k 0 at various streamwise locations is shown in Figure B2 for the three flames. The LES results compare better with experimental data compared to previous RANS results as one would expect. The agreement between the measured and current computational values improves as the Reynolds number increases (going from F3 to F1), which is consistent with the general expectations for LES. The current results are in good agreement with the previous results as seen in Figure B2. However, the results of Wang et al. (2011) agree better with the data because they have used the 8.5M grid. The results obtained using the 4.2M  (Chen et al., 1996) () and simulations ( ) using the 1.5M grid. Results obtained using the 4.2M grid ( ) are shown only for F2. Previous RANS (Kolla and Swaminathan, 2010b) ( ) and LES (Wang et al., 2011) ( ) results are also shown for comparison.
grid agree quite well with the 1.5M results and the measurements for r ≥ 0.8D. The 4.2M grid results are good in general but an overprediction is seen for the inner region of the jet at downstream locations. Generally, second-order statistics such as TKE will have increased sensitivity to mesh resolution compared to first-order statistics as shown in Figure 19 of Wang et al. (2011). Thus, one would expect to see a better agreement for the 4.2M grid and the large difference seen here is because the combustion submodel in Eq. (1) is at its limit for this grid as noted earlier (see the discussion on Figures 4,7,8,9,and 10). This SGS combustion model underestimates the reaction rate, heat release rate effects, and their interactions with turbulence as noted earlier for the 4.2M grid case. For these reasons, the 4.2M grid results will not be considered for the discussion below unless it is required otherwise. The radial variation of averaged fuel mass fraction is compared to the measured values in Figure B3. The results of previous numerical studies are also included in this figure. The thickened flame-LES results of De and Acharya (2009) agree very well with the measurements for x 6:5D for the high Da flame F3 and these computational results gradually move away from the experimental data as seen in this figure for the flame F1. However, the results of Wang et al. (2011) using a dynamic thickened flame approach show a good agreement uniformly. The agreement between the measured and computed values using the LES-PDF method (Dodoulas and Navarro-Martinez, 2013) is not uniformly good. The agreement is relatively better for x ≥ 8.5 for the flames F1 and F2, and for the flame F3 the agreement seems to improve from x = 6.5D as seen in Figure B3. It had been suggested that including the differential diffusion effects could improve these estimates for the LES-PDF method. The mean fuel mass fraction computed in this study using the flamelet-based algebraic model agrees very well with the measurements as shown in Figure B3. A similar comparison is shown in Figure B4 for the normalized temperature, hT þ i. The values shown in Figures 10 and B4 are calculated as follows. First, the normalized filtered temperature is obtained using an expression proposed by Kolla and Swaminathan (2010b): where the RMS of normalized temperature fluctuation is obtained using T 00 ¼ e T þ À h e T þ i and (Poinsot and Veynante, 2005). Strictly, one must also include the sub-grid variance, σ 2 T þ ;sgs , while calculating σ T þ and this SGS variance is unavailable in the modeling framework used here. A transport equation for σ 2 T;sgs needs to be solved to include the influences of reaction rate, dissipation rate, and turbulence on the evolution of this SGS variance and this will be considered in a future study.
Although one should use hT þ i for comparison with measurements, no significant difference between hT þ i and h e T þ i was observed for the flame conditions considered in this study. The flame F3 has the highest Damköhler number in the set and the computed values agree well with the data as shown in Figure B4. A very good agreement among various numerical studies is also observed for this flame. A large overestimate seen for the flame F1 for x 4:5D is because of the uncertainty in the pilot temperature, as it is not reported by Chen et al. (1996). Furthermore, the flame in the near field of the jet is expected to have quasi-laminar unsteady structure since the turbulence in these regions is low. The results in Figure B4 suggest that the compounded effects of these two factors (uncertainty and quasi-laminar unsteady flames) become more apparent as the turbulent Reynolds number increases (compare the difference between the computational and experimental values at x = 2.5D for all the three flames). The values of hT þ i computed in this study agree well with the measurements at down stream locations as shown in the Figure B4. In general, a good agreement is observed for the statistics computed in this study with measurements. The radial variations of averaged mass fractions of major, O 2 , H 2 O, and CO 2 , and minor, OH, species and intermediates, H 2 and CO, are shown in Figures B5-B10 for various axial locations. The current statistics are also compared to past results available in the literature. A 2-step chemistry with empirical rate expression was used by De and Acharya (2009) and 1-step chemistry was used by Wang et al. (2011). The LES-PDF study by Dodoulas and Navarro-Martinez (2013) used a 4-step, 15-step, and GRI3.0 mechanisms and the results of the 15-step mechanism are used here. Because of these differences in the chemistry, a uniform comparison cannot be made for the various species.
The mass fractions shown in these figures are obtained after correcting for the dilution effect as noted in Section 5, which is elucidated using the oxygen mass fraction at axial positions of x/ D = 2.5, 8.5, and 10.5 for the flame F1. In the near field (x/D = 2.5), the influence of this effect is almost negligible and including this effect captures the oxygen mass fraction increase in the outer region, r > 0.8D, due to dilution with the air. The oxygen mass fraction computed using a 2-step chemistry and thickened flame model by De and Acharya (2009) is relatively larger compared to those obtained in the current study and those from the LES-PDF calculation by Dodoulas and Navarro-Martinez (2013). This difference increases with downstream distance. Wang et al. (2011) did not report O 2 and CO 2 mass fractions and thus they are not shown in Figures B5 and B7. The current results are improved over the previous RANS results as one would expect for the flames F1 and F2. This improvement is observed for H 2 O, H 2 , and OH. The averaged mass fraction of water vapor computed in this study agrees well with measured values and those reported in earlier LES studies. The underestimate seen for CO 2 mass fractions in Figure B7 is consistent with many earlier studies (De and Acharya, 2009;Dodoulas and Navarro-Martinez, 2013;Kolla and Swaminathan, 2010b) employing different combustion modeling approach. However, the values obtained using the LES-PDF method are relatively closer to the measurements for downstream locations. Although CO 2 is transported in the calculations reported by De and Acharya (2009) there is a large underestimate inside the flame brush for x ≥ 6.5D for the flame F1. For high Da flame F3, the computed CO 2 mass fractions agree well with the measured values. The OH mass fractions obtained using the LES-PDF (Dodoulas and Navarro-Martinez, 2013) and the current LES-flamelets methods agree quite well with one another and also with measurements for the flame F1. These agreements are relatively poor for upstream locations and for the largest Da flame F3 as shown in Figure B8. Nevertheless, the general trend is captured well by both of these modeling approaches.
The mean mass fraction of H 2 computed in this study agrees well with the measured values as shown in Figure B9. None of the previous LES studies reported this result and thus they are not included in this figure for comparison. The slight overestimate at x/D = 2.5 is because of the uncertainty in the pilot boundary condition and this uncertainty does not seem to be of any importance as one moves downstream.
The agreement for hY CO i shown in Figure B10 is not good and there is a large overestimate for all three flames, which is consistent with the underestimate of CO 2 mass fraction. This could be related to the slow oxidation of CO, which may not be good for flamelet approach. However, a similar behavior is seen in Figure B10 for many previous studies using other combustion modeling techniques. This suggests that the reason may be somewhere else, perhaps in the measurement itself. Nevertheless, further analysis of the combustion model used in this study using other flame configurations and conditions is required. It is worth noting that, although the CO values computed by De and Acharya (2009) seem to agree quite well with the measurements, this seems to be inconsistent with their CO 2 values see Figure B7.