Modelling of Turbulent Premixed and Partially Premixed Combustion

A significant increase in our energy consumption, from 495 quadrillion Btu in 2007 to 739 quadrillion Btu in 2035 with about 1.4% annual increase, is predicted (US Energy Information Administraion, 2010). This increase is to be met in environmentally friendly means in order to protect our planet. Despite the renewable energy sources are identified to be the fastest growing in the near future, they are expected to meet only one third or less of this energy demand. Also, the renewable generation methods face significant barriers such as economical risks, high capital costs, cost for infrastructure development, low energy conversion efficiency, and low acceptance level from public (US Energy Information Administraion, 2010) at this time. It is likely that improvements will be made on all of these factors in due course. The role of nuclear technology in the energy market will vary from time to time for many cultural and political reasons, and the perceptions of the general public. In the current climate, however, it is clear that the fossil fuels will remain as the dominant source to meet the demand in energy consumption. Hence, optimised design of combustion and power generating systems for improved efficiency and emissions performance are crucial. The emissions of oxides of nitrogen and sulphur, and poly-aromatic hydrocarbons are known sources of atmospheric pollution from combustion. Their detrimental effects on environment and human health is well known (Sawyer, 2009) and green house gases such as oxides of carbon and some hydrocarbons are also included as pollutants in recent years. The emission of carbon dioxide (CO2) from fossil, liquid and solid, fuel combustion accounts for nearly 76% of the total emissions from fossil fuel burning and cement production in 2007 (Carbon Dioxide Information Analysis Center, 2007). The global mean CO2 level in the atmosphere increases each year by about 0.5% suggesting a global mean level of about 420 ppm by 2025 (Anastasi et al., 1990; US Department of Commerce, 2011) Such a forecasted increase has led to stringent emission regulations for combustion systems compelling us to find avenues to improve the environmental friendliness of these systems. Lean premixed combustion is known (Heywood, 1976) to have potentials for effective reduction in emissions and to increase efficiency simultaneously. Significant technological advances are yet to be made for developing fuel lean combustion systems operating over wide range of conditions with


Introduction
A significant increase in our energy consumption, from 495 quadrillion Btu in 2007 to 739 quadrillion Btu in 2035 with about 1.4% annual increase, is predicted (US Energy Information Administraion, 2010). This increase is to be met in environmentally friendly means in order to protect our planet. Despite the renewable energy sources are identified to be the fastest growing in the near future, they are expected to meet only one third or less of this energy demand. Also, the renewable generation methods face significant barriers such as economical risks, high capital costs, cost for infrastructure development, low energy conversion efficiency, and low acceptance level from public (US Energy Information Administraion, 2010) at this time. It is likely that improvements will be made on all of these factors in due course. The role of nuclear technology in the energy market will vary from time to time for many cultural and political reasons, and the perceptions of the general public. In the current climate, however, it is clear that the fossil fuels will remain as the dominant source to meet the demand in energy consumption. Hence, optimised design of combustion and power generating systems for improved efficiency and emissions performance are crucial. The emissions of oxides of nitrogen and sulphur, and poly-aromatic hydrocarbons are known sources of atmospheric pollution from combustion. Their detrimental effects on environment and human health is well known (Sawyer, 2009) and green house gases such as oxides of carbon and some hydrocarbons are also included as pollutants in recent years. The emission of carbon dioxide (CO 2 ) from fossil, liquid and solid, fuel combustion accounts for nearly 76% of the total emissions from fossil fuel burning and cement production in 2007 (Carbon Dioxide Information Analysis Center, 2007). The global mean CO 2 level in the atmosphere increases each year by about 0.5% suggesting a global mean level of about 420 ppm by 2025 (Anastasi et al., 1990; US Department of Commerce, 2011) Such a forecasted increase has led to stringent emission regulations for combustion systems compelling us to find avenues to improve the environmental friendliness of these systems. Lean premixed combustion is known (Heywood, 1976) to have potentials for effective reduction in emissions and to increase efficiency simultaneously. Significant technological advances are yet to be made for developing fuel lean combustion systems operating over wide range of conditions with desirable characteristics. This is because, ignition and stability of this combustion, controlled strongly by turbulence-combustion interaction, are not fully understood yet. Many of the recent studies are focused to improve this understanding as it has been noted in the books edited by Echekki & Mastorakos (2011) and Swaminathan & Bray (2011). The OEMs (original equipment manufacturers) of gas turbines and internal combustion engines are embracing computational fluid dynamics (CFD) into their design practice to find answers to what if? type questions arising at the design stage. This is because CFD provides quicker and more economical solutions compared to "cut metal and try" approach. Thus, having an accurate, reliable and robust combustion modelling becomes indispensable while developing modern combustors or engines for fuel lean operation. In this chapter, we discuss one such modelling method developed recently for lean premixed flames along with its extension to partially premixed combustion. Partial premixing is inevitable in practical systems and introduced deliberately under many circumstances to improve the flame ignitability, stability and safety. Before embarking on this modelling discussion, challenges in using the standard moment methods for reacting flows, which are routinely used for non-reacting flows, are discussed in the next section along with a brief discussion on three major computational paradigms used to study turbulent flames. Section 1.2 identifies important scales of turbulent flame and discusses a combustion regime diagram. The governing equations for Reynolds averaged Navier Stokes (RANS) simulation of turbulent combustion are discussed in section 2 along with turbulence modelling used in this study. The various modelling approaches for lean premixed combustion are briefly discussed in section 3. The detail of strained flamelet model and its extension to partially premixed flames are presented in section 4. Its implementation in a commercial CFD code is discussed in section 4.2 and the results are discussed in section 5. The final section concludes this chapter with a summary and identifies a couple of topics for further model development.

Challenges
In the RANS approach, the instantaneous quantities are decomposed into their means and fluctuations. The mean values of density, velocity, etc., in a flow are computed by solving their transport equations along with appropriate modelling hypothesis for correlations of fluctuating quantities. These modelling are discussed briefly in section 2. The simulations of non-reacting flows have become relatively easier task now a days. The presence of combustion however, significantly complicates matters and alternative approaches are to be sought. Combustion of hydrocarbon, even the simplest one methane, with air includes several hundreds of elementary reactions involving several tens of reactive species. If one follows the traditional moment approach by decomposing each scalar concentration into its mean and fluctuation then it is clear that several tens of partial differential equations for the conservation of mean scalar concentration need to be solved. These equations will involve many correlations of fluctuations requiring closure models with a large set of model parameters. More importantly, the highly non-linear reaction rate is difficult to close accurately. To put this issue in a clear perspective, let us consider an elementary chemical reaction R 1 + R 2 → P, involving two reactants and a product. The law of mass action would give the instantaneous reaction rate for R 1 asω 1 = −AT b ρ 2 Y 1 Y 2 exp (−T a /T),w h e r eA is a pre-exponential factor and T a is the activation temperature. The temperature is denoted as T and the mass fractions of two reactants are respectively denoted by Y 1 and Y 2 . Let us take b = 0 for the sake of simplicity. In variable density flows, it is normal to use density-weighted means, defined, for example for mass fraction of reactant R 1 ,as Y 1 = ρ Y 1 /ρ and its fluctuation y ′′ 1 . Substituting this decomposition into the above reaction rate expression, one obtainṡ The exponential term can be shown (Williams, 1985) to be Since T a / T is generally large, at least about 20 terms in the above expansion are required to have a convergent series. This is impractical and also T " / T is seldom smaller than 0.01 in turbulent combustion to neglect higher order terms in the above expansion. There are already some approximations made while writing Eq. (1) and furthermore, while averaging this equation to getω one must not forget that ρ exp(−T a /T) = ρ exp(−T a / T).
It is clear that one needs to solve a large set of coupled partial differential equations with numerous model parameters, which poses a serious question on the accuracy and validity of computed solution using the classical RANS approach which usually tracks the first two moments for each of the quantities involved. This is a well-known problem in turbulent combustion and alternative approaches have been developed in the past (Echekki & Mastorakos, 2011;Libby & Williams, 1994;Swaminathan & Bray, 2011). The first two statistical moments of one or two key scalars are computed instead of solving hundreds of partial differential equations for the statistical moments and correlations of all the reactive scalars. The statistics of the two key scalars, typically a mixture fraction, Z, and reaction progress variable, c, are then used to estimate the thermo-chemical state of the chemically reacting mixture using modelling hypothesis. Many modelling approaches have been proposed in the past and the readers are referred to the books by Libby & Williams (1994), Peters (2000), Mastorakos (2011), andBray (2011). Approaches relevant for lean premixed and partially premixed flames are briefly reviewed in later parts of this chapter. The three computational paradigms generally used to study turbulent combustion are (i) direct numerical simulation (DNS), (ii) large eddy simulation (LES) and (iii) RANS. These approaches have their own advantages and disadvantages. For example, the detail and level of information available for analysis decreases from (i) to (iii). Thus DNS is usually used for model testing and validation and it incurs a heavy computational cost because it resolves all the length and time scales (in the continuum sense) involved in the reacting flow. The general background of DNS is discussed elsewhere (Chapter ??) in this book. With the advent of Tera-and Peta-scale computing, it is becoming possible to directly simulate laboratory scale flames with hundreds of chemical reactions. However, direct simulations of practical flames in industry are not to be expected in the near future.
On the other hand, in RANS all the scales of flow and flame are modelled and thus it provides only statistical information and it is possible to include different level of chemical kinetics detail as will be discussed later in this chapter. If the RANS simulations are performed carefully, they provide solutions with sufficient accuracy to guide the design of combustors and engines. These simulations are cost effective and quick, and thus they are attractive for use in industries. The LES is in between these two extremes as it explicitly computes large scales in the flow and it is well suited for certain class of flows. Still models are required for quantities related to small scales. The general background of LES is discussed in Chapter ?? of this book. It should be noted that combustion occurs in scales much smaller than those usually captured in LES and thus they need to be modelled. Most LES models are based on RANS type modelling and this chapter presents combustion modelling for the RANS framework. The details of turbulent combustion modelling depend on the combustion regime, which is determined by the relativity of characteristic scales of turbulence and flame chemistry.

Regimes in turbulent premixed combustion
The characteristic flame scales are defined using the unstrained planar laminar flame speed, s 0 L , and Zeldovich flame thickness given by δ = D/s 0 L ,w h e r eD is a molecular diffusion coefficient. Using these velocity and length scales, one can define the flame time scale as t F = δ/s 0 L . The turbulence scales are defined using turbulent kinetic energy, k = 0.5ρ u ′′ i u ′′ i /ρ, and its dissipation rate, ε = 2ρνs ij s ij /ρ (strictly 2ρν(s ij s ij − s ii s ii /3/ρ for combusting flows), where s ij is the symmetric part of the turbulent strain tensor and ν is the kinematic viscosity of the fluid. The characteristic turbulence length and time scales are respectively Λ = k 1.5 / ε and t T = k/ ε. The RMS of turbulent velocity fluctuation is u ′ = 2 k/3. The viscous dissipation scales are the Kolmogorov scales given by l η =( ν 3 / ε) 1/4 and t η =( ν/ ε) 1/2 . For the sake of simplicity, let us take equal molecular diffusivities for all reactive species and the Schmidt number to be unity. One can combine these characteristic scales to form three non-dimensional parameters, turbulence Reynolds, Damköhler and Karlovitz numbers, given respectively by ( 3 ) Figure 1 illustrates the relationship between these three parameters in (u ′ /s 0 L )-(Λ/δ) space. This diagram is commonly known as combustion regime diagram (Peters, 2000). The left For u ′ /s 0 L < 1, turbulent fluctuation due to large eddy cannot compete with the flame advancement by the laminar flame propagation mechanism, thus laminar flame propagation dominates over flame front wrinkling by turbulence. For u ′ /s 0 L > 1a n dK a < 1, the flame scales are smaller than all the relevant turbulent scales. As a result, turbulent eddies cannot disturb the inner reactive-diffusive structure but only wrinkle it and the flame also remains quasi-steady. The Ka = 1 line represents that the flame thickness is equal to the Kolmogorov length scale, which is known as the Klimov-Williams line. Although the Kolmogorov scales are smaller than the flame thickness, the inner reaction zone is expected to be undisturbed by the small scale turbulence when Ka < 100. However, the preheat zones are disturbed by the small-scale turbulence. This regime is known as thin reaction zones regime. When Ka > 100, the Kolmogorov eddies can penetrate into the inner reaction zones causing local extinction leading to distributed reaction zones. However, evidence for these reaction zones are sparse (Driscoll, 2008). Also shown in the diagram are the combustion regimes for practical engines (Swaminathan & Bray, 2011). The spark ignited internal combustion engines operate in the border between corrugated flamelets and thin reaction zones regimes, whereas the power gas turbines operate in the border between corrugated and wrinkled flamelets. Aero engines do not operate in premixed mode for safety reasons and if one presumes a premixed mode for them then they are mostly in thin reaction zones regime.

Modelling framework for turbulent flames
In the RANS modelling framework, the conservation equations for mass, momentum, energy and the key scalar values are averaged appropriately. These equations are solved along with models and averaged form of state equation. As noted earlier, density weighted average is used for turbulent combustion and the Favre averaged mass and momentum equations are given by ∂ρ ∂t in the usual nomenclature. In adiabatic and low speed (negligible compressibility effects) combustion problems, the Favre averaged total enthalpy equation given by for the reacting mixture is solved. The total enthalpy h is defined as h = ∑ i Y i h i ,w h e r e Y i is the Favre averaged mass fraction of species i. The specific enthalpy of species i is where the standard specific enthalpy of formation for species i and its specific heat capacity at constant pressure are respectively h 0 i and c p,i . The equation (7) is used to calculate the temperature, T, from the computed values of h. The state equation is given by p = ρ R T, where R is the gas constant. One must also know the mean mass fraction fields, which are obtained using combustion modelling. The two key scalars used in turbulent combustion are the mixture fraction, Z, and reaction progress variable, c. These two scalars are defined later and the transport equation for their Favre averaged values are respectively given by The Reynolds stress, ρ u ′′ i u ′′ j , and fluxes in Eqs. (5) to (9) and the mean reaction rate,ω c ,o f c need closure models. In principle, all of the above equations are applicable to both RANS and LES, and the correlations must be interpreted appropriately. The Reynolds fluxes are typically closed using gradient hypothesis, which gives, for example, u ′′ j Z ′′ = −D T ∂ Z/∂x j , where D T is the turbulent diffusivity. This quantity and the Reynolds stress are obtained using turbulence modelling.

Turbulence modelling
The modelling of turbulence is an essential part of turbulent combustion calculation using CFD. A variety of turbulence models are available in the literature (Libby & Williams, 1994;Swaminathan & Bray, 2011, see for example) and an appropriate choice should be guided by a physical understanding of the flow. A standard two-equation model like k-ε model can be used to model simple free shear flows. In this model, the Reynolds stress is related to the eddy viscosity, μ T ,by However, complex geometries might require improved models, such as transported Reynolds stress, to capture the relevant flow features such as flow recirculation, the onset of flow separation and its reattachment, etc. These advanced models are known to cause difficulties such as numerical instability during simulations compared to the two-equation models. Furthermore, the two-equation models are commonly used because, (i) they are easy to implement and use, (ii) numerically stable and (iii) provide sufficiently accurate solution to guide the analysis of turbulent flames, provided the mean reaction rate and turbulence-flame interactions are modelled correctly. Two variants of two-equation model commonly used in turbulent premixed flame calculations are discussed briefly next.

k-ε Model
Despite the advancement in understanding and modelling of the turbulence, the standard k-ε model is still one of the most widely used models for engineering calculations. The attractiveness of this model is rooted in its simplicity and favourable numerical characteristics, more importantly, in its surprisingly good predictive capabilities over a fair range of flow conditions. It represents a reasonable compromise between accuracy and cost while dealing with various flows. Transport equations for the standard k-ε model are (Jones, 1994) ∂ρ k ∂t where μ t is the turbulent eddy viscosity calculated as μ t = ρC μ k 2 / ε and Sc m is the turbulent Schmidt number for quantity m. The standard values of these model parameters are C μ = 0.09, C ε1 =1.44 and C ε2 =1.92. The pressure-dilation term can become important in turbulent premixed flames and it is modelled as (Zhang & Rutland, 1995) where τ is the heat release parameter defined as the ratio of temperature rise across the flame front normalised by the unburnt mixture temperature. The average value of the Favre (Jones, 1994).

Shear stress transport k-ω model
The Reynolds stress is closed using Eq. (10) in the approach also, but the eddy viscosity is obtained in a different manner as described below. The shear stress transport k-ω model proposed by Menter (1994) aims to combine the advantages in predictive capability of both the standard k-ε model in the free shear flow and the k-ω model, originally proposed by Wilcox (1988), in the near wall region. There are two important ingredients in this model. Firstly, a blending function, F 1 , is used to appropriately activate the k-ε model in free shear flow part and the k-ω model in near wall region of the flow. Secondly, the definition of eddy viscosity is modified to include the effects of the principal turbulent stress transport. In this method, the turbulent kinetic energy equation is very similar to Eq. (11) and ε equation is replaced by a transport equation for ω = ε/ k. This equation is obtained using Eqs. (11) and (12) and is written as whereα , β and σ ω 2 are model constants given by Menter (1994). The turbulent eddy viscosity is obtained using where F 2 is a function taking a value of one in boundary layer and zero in free shear region of wall bounded flows, and Ω = abs(∇× u). A number of other turbulence modelling approaches are discussed by Pope (2000) and interested readers can find the detail in there. The combustion sub-modelling is considered next.

Combustion sub-modelling
As noted earlier, the governing equations for thermo-chemical state of the reacting mixture can be reduced to one or two key scalars. The reaction progress variable is usually the key scalar for lean premixed flames. It is usually defined using either normalised fuel mass fraction or normalised temperature. Alternative definition using the sensible enthalpy have also been proposed in the literature (Bilger et al., 1991). Here, the progress variable is defined using temperature as c =( and T u , are the adiabatic and unburnt temperatures respectively. A transport equation for the instantaneous progress variable can be written as (Poinsot & Veynante, 2000) ∂ρ c ∂t By averaging the above equation, one gets the transport equation for c given as Eq. (9), which needs closure for the turbulent scalar flux, u ′′ j c ′′ , and the mean reaction rate,ω c . Many past studies (Echekki & Mastorakos, 2011;Libby & Williams, 1994;Swaminathan & Bray, 2011, see for example) have shown that the scalar flux can exhibit counter-gradient behaviour in turbulent premixed flames. The counter-gradient flux would yield a negative turbulent diffusivity and mainly arises when the local pressure forces accelerate the burnt and unburnt mixtures differentially due to their density difference. This phenomenon is predominant in low turbulence level, when u ′ /s 0 L is smaller than about 4, where thermo-chemical effects overwhelms the turbulence. Bray (2011) has reviewed the past studies on the scalar flux and has noted that the transition between gradient and non-gradient transport in complex flows deserves further investigation. However, it is quite common to model this scalar flux as a gradient transport in situations with large turbulence Reynolds number. The closure forω c in Eq. (9) is central in turbulent premixed flame modelling. Many closure models have been proposed in the past and they are briefly reviewed first before elaborating on the strained flamelet approach.

Eddy break-up model
This model proposed by Spalding (1976) is based on phenomenological analysis of scalar energy cascade in turbulent flames with Re ≫ 1a n dD a ≫ 1. The mean reaction rate is given byω = C EBU ρ ε c ′′ 2 / k,whereC EBU is the model parameter (Veynante & Vervisch, 2002). The large values of Re and Da implies that the combustion is in the flamelets regime where the flame front is thinner than the small scales of turbulence. In this regime, the variance can be written as c ′′ 2 = c(1 − c) and the reaction is assumed to be fast. This fast chemistry assumption does not hold in many practical situations and thus this model tends to over predict the mean reaction rate. Furthermore, this model does not consider the multi-step nature of combustion chemistry. A variant of this approach, known as eddy dissipation concept, is developed with provisions to include complex chemical kinetics (Ertesvag & Magnussen, 2000;Gran & Magnussen, 1996).

Bray-Moss-Libby model
This model (Bray, 1980;Bray & Libby, 1976;Bray & Moss, 1977) uses a reaction progress variable and its statistics for thermo-chemical closure. An elaborate discussion of this modelling can be found in many books (Bray, 2011;Libby & Williams, 1994, for example). The basic assumption in this approach is that the turbulent flame front is thin and the turbulence scales do not disturb its structure. This allows one to partition the marginal probability density function (PDF) of c into three distinct portions; fresh gases (0 ≤ c ≤ c 1 ), burnt mixtures ((1 − c 1 ) ≤ c ≤ 1) and reacting mixtures (c 1 ≤ c ≤ (1 − c 1 )) with probabilities α(x, t), β(x, t) and γ(x, t) respectively, obeying α + β + γ = 1. In the limit of c 1 → 0, the Reynolds PDF of c can be written as where the fresh gases and fully burnt mixtures are represented by the Dirac delta functions δ(c) and δ(1 − c) respectively. The progress variable is defined as c =( The interior portion, f (c; x, t), of the PDF represents the reacting mixture and it must satisfy 1 0 f (c; x, t) dc = 1. If the flame front is taken to be thin then γ ≪ 1whenDa≫ 1. Under this condition, it is straightforward to obtain α =(1 − c)/(1 + τ c) and β =(1 + τ) c/(1 + τ c) using ρ = ρ P dc with ρ/ρ =( 1 + τ c)/(1 + τ c). Now, the mean value of any thermo-chemical variable can be obtained simply as ϕ(x)= ϕ(x) P(c; x) dc =( 1 − c)ϕ u + cϕ b ,w h e r et h e Favre PDF is given by P = ρ P/ρ and the mean density is ρ = ρ u /(1 + τ c). Since the reaction rate is zero everywhere outside the reaction zones, its mean value, is proportional to γ.S i n c eγ has been neglected in the above analysis, alternative means are to be devised to estimate the mean reaction rate. One method is to treat the progress variable signal as a telegraphic signal (Bray et al., 1984). This analysis yields that the mean reaction rate is directly proportional to the frequency of undisturbed laminar flame front crossing a given location in the turbulent reacting flow. More detail of this analysis and its experimental verification is reviewed by Bray & Peters (1994). In another approach, the turbulent flame front is presumed to have the structure of unstrained planar laminar flame and this approach  gives a model for the mean reaction rate aṡ The symbol δ * L is a laminar flame thickness defined as δ * L = c(1 − c)/(1 + τc) dn,wheren is the distance along the flame normal. The small parameter ε,definedasε = 1 − c "2 /[ c(1 − c)], is related to γ. This implies that one must also solve a transport equation for the Favre variance, which is given in section 4 as Eq.(30). The Favre variance of c is equal to c(1 − c) when γ is neglected and the bimodal PDF is used. This simply means that ε = 0f o r Eq. (19) which also concurs our earlier observation on the mean reaction rate closure. The variance transport equation, Eq. (30), reduces to a second equation for c when γ is neglected and a reconciliation of these two transport equations led (Bray, 1979) to the conclusion that the mean reaction rate is directly proportional to the mean scalar dissipation rate, ǫ c = ρ D c (∇c" ·∇c")/ρ,andω where C m typically varies between 0.7 and 0.8 for hydrocarbon-air flames. The mean scalar dissipation rate is an unclosed quantity and if one uses a classical model, ǫ c = C d c "2 ( ε/ k),for it then one getsω which is very similar to the eddy break-up model discussed in section 3.1.

Flame surface density model
In this approach (Marble & Broadwell, 1977), the mean reaction rate is expressed as the product of reaction rate per unit flame surface area, ρ u s L , and the flame surface area per unit volume, Σ, which is known as the flame surface density (FSD). The straining, bending, wrinkling and contortion, collectively called as stretching, of the flame surface by turbulent eddies can influence the flame front propagation speed and thus it is quite useful and usual to write s L = s 0 L I 0 ,whereI 0 is known as the stretch factor, which is typically of order one. The FSD approach has been studied extensively and these studies are reviewed and summarised by Veynante & Vervisch (2002) and in the books edited by Libby & Williams (1994), Echekki & Mastorakos (2011) and Swaminathan & Bray (2011). Two approaches are normally used to model Σ; in one method an algebraic expression (Bray & Peters, 1994; is used and in another method a modelled differential equation is solved. A simple algebraic model proposed by  is given as with the three model parameters which are of order unity. When Da ≫ 1, suggesting that the generalised FSD scales with the laminar flame thickness rather than the turbulence integral length scale, Λ. Many earlier algebraic models discussed by Bray & Peters (1994) suggest a scaling with Λ. An algebraic FSD model has also been deduced using fractal theories by Gouldin et al. (1989). The unclosed transport equation for FSD was derived rigorously by Candel & Poinsot (1990) and Pope (1988) and this equation is written, when there is no flame-flame interaction, as where ··· s denotes the surface average. The three flux terms on the left hand side are due to the mean flow advection, turbulent diffusion and flame displacement at speed s d .T h e last three terms in the above equation require modelling and usually the propagation term is neglected in the modelling, which may not hold at all situations. The turbulent diffusion is usually modelled using gradient flux hypothesis. The term on the right hand side is the source or sink term for Σ due to the effects of turbulence on the flame surface. The quantity Φ is usually called as flame stretch which is a measure of the change in the flame surface area, A, and is given by (Candel & Poinsot, 1990) where k m is the mean curvature of the flame surface, e ij is the symmetric strain tensor and n i is the component of the flame normal in direction i. Turbulence, generally, has the tendency to increase the surface area implying that the average stretch rate is positive. However, the curvature term is negative because s d is negatively correlated to k m , while the tangential strain rate, a T , is positive Trouvé & Poinsot, 1994). Thus, to have Φ s > 0 the magnitude of the surface averaged tangential strain rate must be larger than or equal to the magnitude of 2s d k m s . The mean curvature of the flame surface was observed to be around zero in a number of studies (Baum et al., 1994;Chakraborty & Cant, 2004;Echekki & Chen, 1996;Gashi et al., 2005). However, fluctuation in the flame surface curvature contributes to the dissipation of the surface area. The modelling of a T is also typically done by splitting the contributions into mean and fluctuating fields and obtaining accurate models for the various contributions to the flame surface density has been the subject of many earlier studies Peters, 2000;Peters et al., 1998). The FSD method has also been developed for LES of turbulent premixed flames and these works are summarised by Vervisch et al. (2011) and Cant (2011).

G-equation, level set approach
A smooth function G,s u c ht h a tG < 0i nu n b u r n tm i x t u r e ,G > 0 in burnt mixture, and G = G o = 0 at the flame, is introduced for premixed combustion occurring in reaction-sheet, wrinkled flamelets, regime. A Huygens-type evolution equation can be written for the instantaneous flame element as (Kerstein et al., 1988;Markstein, 1964;Williams, 1985) ∂G ∂t where s G is the propagation speed of the flame element relative to the unburnt mixture. This propagation speed may be expressed in terms of s 0 L corrected for stretch effects using Markstein number. This number is a measure of the sensitivity of the laminar flame speed to the flame stretch (Markstein, 1964). Peters (1992;1999) has developed this approach for corrugated flamelets and thin reaction zones regimes of turbulent premixed combustion and also proposed (Peters, 2000) s G expressions suited to these regimes. Using the above instantaneous equation, Peters (1999) deduced transport equations for G and G "2 , which can be used in RANS simulations (Herrmann, 2006, see for example). The development and use of this method for LES has been reviewed by Pitsch (2006). A close relationship between the G field and the FSD is shown by Bray & Peters (1994) as where P (G; x) is the PDF of G and an approximate expression has been proposed for P(G o ; x) by Bray & Peters (1994).

Transported PDF approach
In this method, a transport equation for the joint PDF of scalar concentrations is solved along with equations for turbulence quantities. The transport equation for the joint PDF has been presented and discussed by Pope (1985). The attractive aspect of this approach is that the non-linear reaction rate is closed and does not require a model. However, the molecular flux in the sample space, known as micro-mixing, needs a closure model and many models are available in the literature. These models are discussed by Haworth & Pope (2011), Lindstedt (2011) and Dopazo (1994). The micro-mixing is directly related to the conditional dissipation. This dissipation rate, for example for the progress variable, is defined as N c |ζ = D∇c ·∇c|ζ ,whereζ is the sample space variable for c. The predictive ability of this method depends largely on the quality of the models used for the unclosed terms. The joint PDF equation is of (N + 4) dimensions for unsteady reacting flows in three spatial dimensions involving (N − 1) reactive scalars and temperature. This high dimensionality of the PDF transport equation poses difficulties for numerical solution. Also the molecular flux term will have a negative sign which precludes the use of standard numerical approaches such as the finite difference or the finite volume methods, and the Monte-Carlo methods (Pope, 1985) are generally used to solve the PDF transport equation. In this method, the computational memory requirement depends linearly on the dimensionality (number of particles used) of problem but one needs a sufficiently large number of particles to get a good accuracy.

Presumed PDF method
The marginal PDF of the key scalars are presumed to have a known shape, which is determined usually using computed values of the first two moments, mean and variance. A Beta function is normally used and it is given by where a and b are related to the first and second moments of the key scalar ϕ,w i t hs a m p l e space variable ζ,by The normalising factorβ is the Beta function (Davis, 1970), which is related to the Gamma function given byβ This presumed form provides an appropriate range of shapes: if a and b approach zero in the limit of large variance then the PDF resembles a bimodal shape of the BML PDF in section 3.2, which requires only the mean value, ϕ. In the limit of small variance, Eq. (27) develops a mono-modal form with an internal peak and it has been shown by Girimaji (1991) that this PDF behaves likes the Gaussian when the variance is very small. The variance equation is written as using the standard notations. The contributions from the scalar flux, dissipation rate and the chemical reactions denoted respectively by the last three terms in the above equation, need to be modelled. The modelling of scalar dissipation rate is addressed in many recent studies, which are discussed by Chakraborty et al. (2011). This quantity is typically modelled using turbulence time scale as has been noted in section 3.2 but, this model is known to be inadequate for premixed and partially premixed flames. A recent proposition by Kolla et al. (2009) includes the turbulence as well as laminar flame time scales and this model is given by where, β ′ is 6.7 and K * c /τ is 0.85 for methane-air combustion. The parameters C 3 and C 4 are given by The Karlovitz number, Ka, is defined as where t c is the chemical or flame time scale defined earlier, t η is the Kolmogorov time scale and ν is the kinematic viscosity. The Zeldovich thickness is related to the thermal thickness by δ 0 L /δ ≈ 2(1 + τ) 0.7 . In the presumed PDF approach, the reaction rate is closed aṡ for premixed flames and the functionω(c) is obtained using freely propagating laminar flame having the same thermo-chemical attributes as that of the turbulent flame. In this approach, it is implicit that the flame structure is undisturbed by the turbulent eddies. For partially premixed flames, one can easily extend the above model by including the dependence of the reaction rate on the mixture fraction asω(c, Z) and thus one must integrate over c and Z space to get the mean reaction rate after replacing the marginal PDF by the joint PDF, P(c, Z).I nt h is modelling practice, these two variables are usually taken to be statistically independent. More work is required to address the statistical dependence of Z and c, and its modelling. A closure model for the effects of chemical reaction on the variance transport, see Eq. (30), can be written asω As noted earlier, the closure in Eq. (34) assumes that the laminar flame structure is undisturbed by the turbulent eddies. The influence of fluid dynamic stretch can also be included in this approach as suggested by Bradley (1992) usinġ where P b is the burning rate factor, which can be expressed in terms of Markstein number (Bradley et al., 2005). Recently, Kolla & Swaminathan (2010a) proposed to use the scalar dissipation rate to characterise the stretch effects on flamelets for the following reasons, viz., (i) the chemical reactions produce the scalar gradient and thus the scalar dissipation rate in premixed and stratified flames and (ii) this quantity signifies the mixing rate between hot and cold mixtures, which are required to sustain combustion in premixed and partially premixed flames. This method is elaborated by Kolla & Swaminathan (2010a) and briefly reviewed in the next section.

Strained flamelet model
The closure for the mean reaction rate is given as (Kolla & Swaminathan, 2010a) where ψ is the sample space variable for the instantaneous scalar dissipation rate, N c = D(∇c · ∇c), of the progress variable. The marginal PDF, P(ζ) is obtained from Eq. (27) using the computed values of c and c "2 from Eqs. (9) and (30). The conditionally averaged reaction rate is given by The conditional PDF P(ψ|ζ) is presumed to be log-normal andω(ζ, ψ) is obtained from calculations of strained laminar flames established in opposed flows of cold reactant and hot products. This flamelet configuration seems more appropriate to represent the local scenario in turbulent premixed flames (Hawkes & Chen, 2006;Libby & Williams, 1982). The log-normal PDF is given by The mean, μ,andvariance,σ 2 ,ofln(ψ|ζ) are related to conditional mean N|ζ and variance of the scalar dissipation rate G 2 N via N|ζ = exp( μ + 0.5σ 2 ) and G 2 N = N|ζ 2 [exp(σ 2 ) − 1]. The conditional mean of scalar dissipation rate is related to the unconditional mean through where f (ζ) is the variation of N c normalised by its value at the location of peak heat release rate in unstrained planar laminar flame. A typical variation of f (ζ) is shown in Fig. 2. It has been shown by Kolla & Swaminathan (2010a) that f (ζ) is weakly sensitive to the stretch rate for ζ values representing intense chemical reactions and f (ζ) has got some sensitivity to the stretch rate in thermal region of the flamelet. Despite this, the variation shown in Fig. 2 is sufficiently accurate for turbulent premixed flame calculation and it must be also be noted that f (ζ) will strongly depend on the thermo-chemical conditions of the flamelet. The unconditional mean dissipation rate, ǫ c , is modelled using Eq. (31) and it is to be noted that the model parameters and their numerical values are introduced to represent the correct physical behaviour of ǫ c in various limits of turbulent combustion and thus they are not arbitrary.
The robustness of this model for ǫ c has been shown in earlier studies (Darbyshire et al., 2010;Kolla et al., 2009;Kolla & Swaminathan, 2010a;. The influence of Lewis number on this modelling is also addressed in a previous study (Chakraborty & Swaminathan, 2010).
The mean reaction rate can now be obtained using Eq. (37) for given values of c, c "2 and ǫ c and thus a three dimensional look up table can be constructed for use during turbulent flame calculations. However, care must be exercised to cover a range of fully burning flamelets to a nearly extinguished one. Such a turbulent flame calculation has been reported recently (Kolla & Swaminathan, 2010b) and this study aims to implement this approach in a complex, commercial type, CFD code and validate it by comparing the simulation results to the previously published results. Fig. 2. Typical variation of f (ζ) with ζ from unstrained planar laminar flame. The curve is shown for stoichiometric methane-air combustion at atmospheric conditions.

Partially premixed flame modelling
Partially premixed flame occurs if the fuel and oxidiser were mixed unevenly. As a result, the mixture fraction Z is required to describe the local mixture composition and it is defined following Bilger (1988) as where Z i is the mass fraction of element i with atomic mass W i .T h es u p e r s c r i p t sf and ox refer to reference states of the fuel and oxidiser respectively. The subscripts C,H and O refer to carbon, hydrogen and oxygen. The transport equations for the Favre mean mixture fraction, Z, given as Eq. (8), and its variance Z ′′2 ,givenby are usually solved in the presumed PDF approach to obtain the local mixing related information. The turbulent flux of the variance, u ′′ i Z ′′2 , is expressed using gradient hypothesis in the above equation. The scalar dissipation rate is modelled by assuming a constant ratio of turbulence to scalar time scales and this model is given as where C ξ is a model constant and it is typically unity. The strained flamelet modelling discussed in the previous section can be extended to turbulent partially premixed flame by considering this flame as an ensemble of strained premixed flamelets with mixture fraction ranging from the lean to rich flammability limits. Then, the mean reaction rate can be written asω = ω(ζ, ψ, ξ) P(ζ, ψ, ξ) dζ dψ dξ, where ξ is the sample space variable for Z. Modelling the joint PDF P(ζ, ψ, ξ) is a challenging task and as a first approximation, it is common to assume that ζ and ξ are statistically independent leading to P(ζ, ψ, ξ)=P(ζ, ψ)P(ξ). One can then follow the approach suggested by Eq. (37) to model P (ζ, ψ).T h e m a r g i n a l P D F ,P(ξ), can be obtained using a Beta PDF. Now, the flamelet library will have five controlling parameters instead of three for the purely premixed case. The assumption of statistical independence of Z and c may not be valid and can be easily removed by including the covariance, Z ′′ c ′′ , by solving its transport equation. However, no reliable modelling is available yet to close this equation. One can also extend the unstrained flamelet model, given in Eq. (34), to partially premixed flames by including the mixture fraction in this equation. This is same as Eq. (44) after removing ψ and the associated integral from it.

Assessment using DNS data
A priori assessment of the unstrained and strained flamelet modelling for partially premixed flames is discussed in this subsection. The DNS data for a hydrogen turbulent jet lifted flame (Mizobuchi et al., 2002) is used for this analysis. The time averaged reaction rateω c of the progress variable c, which is defined using the equilibrium value of H 2 Om a s sf r a c t i o na s c = Y H 2 O /Y eq H 2 O has been extracted from the DNS data. Figure 3 presents the radial variation ofω c at two axial positions. The radial distance, r, is normalised using the fuel jet nozzle diameter d.Th ev a luesofω c computed using the unstrained and strained flamelet modellings are also shown in this figure. The means and variances required to construct the PDFs are obtained from the DNS data for this analysis. In figure 3(a), it is clear that unstrained flamelet model agrees well with the DNS results in the region close to the centreline, but it generally over predicts the mean reaction rate for r/d > 1. The strained flamelet model under predictṡ ω c for r/d < 2.5 while giving a good agreement for r/d > 2.5. At a downstream location, figure 3(b) shows a similar trend where unstrained flamelet model over predictsω c while strained flamelet model gives a reasonable agreement. It is likely that the strain effects are important and need to be included to give correct mean reaction rate depending on axial and radial positions and unstrained flamelet model is insufficient. A note of caution is that the strained flamelet model used in this assessment only includes the strain effects for rich mixture and it is constructed with only 12 rich flamelets up to the fuel rich flammability limit. Further work needs to be done to include the strain effects for lean flamelets and to examine the effects of using more than 12 fuel rich flamelets. Exploring a way to combine the unstrained and strained flamelets for partially premixed flames in a unified modelling framework need to be addressed. Whether this approach would be sufficient or a completely different approach would be required, is an open question. Also, the cross dissipation, ǫ cZ = ρD(∇Z" ·∇c")/ρ, can play important role in these kind of closure modelling for partially premixed flames (Bray et al., 2005). It is clear that more works need to be done for partially premixed flames.

Model implementation
The implementation of the strained flamelet model into a commercial CFD software (for example, FLUENT), which can handle complex geometries that are common in industrial scenarios are discussed in this section. The flow and turbulence models available in the software are utilised to provide the required information for combustion calculation. Additional transport equations for Z, Z" 2 , c and c" 2 are included as user defined scalars (UDS). A transport equation for total enthalpy, h, is also included to obtain spatial temperature distribution using Eq. (7). Various sources and sinks appearing in these transport equations, Eqs. (8), (9), (30) and (42), are calculated using user defined functions (UDFs) and the turbulent transports are modelled using gradient hypothesis. The mean density is obtained using the equation of state. Some discrepancies in the modelling of the Reynolds stress in the Fluent for the mean momentum and the k-ε equations were noted and corrected for the results reported in this study. Instructions to incorporate the UDS transport equations and UDFs are provided in the theory and user guides of FLUENT. It is to be noted that the choice of the progress variable is guided broadly by the flame configuration. Progress variable definitions based on temperature or species mass fraction are popular choices for most premixed combustion calculations as noted earlier in this chapter. These choices are equally applicable for open as well as enclosed flames. However, if there are heat losses to the boundary then the fuel or product mass fractions can be used to define the progress variable. A prudent decision on the choice of the progress variable can go a long way in obtaining accurate CFD predictions of flame related quantities.
Once an appropriate progress variable is chosen, the flamelets reaction rate,ω c (ζ, ψ),i s calculated using unstrained and strained laminar flames. An arbitrarily complex chemical kinetics mechanism can be used for these calculations and GRI-3.0 is used for the flames computed and discussed in this chapter. The PREMIX and OPPDIF suites of Chemkin software is used for the flamelet computations. As noted earlier, reactant-to-product configuration is used for the OPPDIF calculations. These flamelets reaction rates are then used in Eq. (37) to obtain the mean reaction rate,ω, as explained in section 4. This mean reaction rate,ω ′′ c ′′ (see Eq. 35), c p , Δh 0 f for the mixture, and Y i are tabulated for 0 ≤ c ≤ 1, These tabulated values are read during a CFD calculation for the computed values of c, c" 2 and ǫ c in each computational cell. The converged fields of these three quantities can then be used to obtain the species mass fractions, Y i , from the tables as a post-processing step. For a purely premixed flames, there is no need to solve for Z and Z" 2 and for partially premixed flames one must solve for these two quantities.

Sample results
Pilot stabilised Bunsen flames  of stoichiometric methane-air mixture with three Reynolds number, based on bulk mean jet velocity and nozzle diameter, of 52000, 40000, and 24000 designated respectively as F1, F2 and F3 have been computed using the strained flamelet model. These flames are in the thin reaction zones regime of turbulent combustion. As noted earlier, the implementation of the strained flamelet model in a commercial CFD software is validated here by comparing the results with those published in an earlier study (Kolla & Swaminathan, 2010b). It is to be noted that the earlier study used a research type CFD code with HLPA (hybrid linear parabolic approximation) discretisation schemes (Zhu, 1991) and TDMA solver. A pressure correction based technique was used in that study while the current study uses a density based method with Roe scheme (Roe, 1981) and a second order accurate upwind discretisation scheme available in Fluent. The model constant C ε1 in the k-ε model is changed from its standard value of 1.44 to 1.52 to correct for the round jet anomaly (Pope, 1978). The turbulent Schmidt numbers for the scalar transport equations are taken to be unity and the turbulent Prandtl number for the enthalpy equation is 0.7 for this study. Mean axial velocity profiles at the nozzle exit measured and reported by Chen et al. (1996) are used as the boundary condition for the inlet velocity. The profiles of RMS of turbulent velocity fluctuations along with longitudinal length scale reported in the experimental study are used to specify the boundary conditions for k and ε. The numerical values for the various model parameters for turbulence and combustion models, turbulent Schmidt and Prandtl numbers, and the boundary conditions used in this study are consistent with those used by Kolla & Swaminathan (2010b).  Chen et al. (1996) and previously published results (--) of Kolla & Swaminathan (2010b).
The computational results for the F1 and F3 flames are compared with previously published results along with experimental measurements here. Figure 4 shows the normalised mean axial velocity and turbulent kinetic energy with radial distance for three axial locations for the flame F1. The distances are normalised by the nozzle diameter D.
The mean axial velocity profiles computed in this study compare well with the experimental measurements and previously published values. The centreline values of k/k o from the Fluent simulation are slightly over predicted compared to the experimental values, but they are almost the same as the values computed by Kolla & Swaminathan (2010b). However, the Fluent calculation over predicts the shear generation of turbulent kinetic energy, which is evident in the higher peak values seen in Fig. 4. The results shown here are grid independent and the differences between the Fluent and previous results are acceptable, given the difference in the numerical schemes, solutions methods and solvers used. Radial variation of the normalised Reynolds mean temperature, c =( T − T u )/(T b − T u ) and the fuel mass fraction are plotted for three axial locations in Fig. 5. The Fluent results are compared to the calculations of Kolla & Swaminathan (2010b) and the experimental measurements of Chen et al. (1996). The mean methane mass fraction shows good comparison with experimental data for flame F1 and the centreline values computed using Fluent agree with previously published values indicating that the flame length is predicted accurately. The computed values of peak mean temperatures at x/D = 2.5 is consistently higher than the experimental measurements, However, the peak mean temperature agrees well for x/D = 8.5 and is under predicted for x/D = 10.5. These trends are consistent with those reported by Kolla & Swaminathan (2010b). Note that the Fluent solution predicts a higher rate of turbulent diffusion of mean progress variable c, which results in lower mean temperature for r/D > 1.0. This could be explained by the higher peak values of turbulence quantities computed by Fluent. Higher values of turbulence would result in increase in the turbulent diffusivity, thus increasing the rate of turbulent diffusion of passive scalars. Also, the Fluent code seem to over predict the rate at which the ambient air is entrained into the reacting jet compared to the solution of Kolla & Swaminathan (2010b).  , --published results of (Kolla & Swaminathan, 2010b). Figure 6 shows the typical values of c and Y CH 4 computed using Fluent for flame F3 at x/D = 8.5 as a function of normalised radial distance, r/D. The experimental measurements and the results computed earlier are also shown for comparison. This flame has lower Reynolds and Karlovitz numbers compared to F1 and hence thermo-chemical effects are dominant compared to turbulence effects. The experimental data clearly shows that the peak value of mean temperature in F3 is larger compared to F1 (cf. Figs. 5 and 6). The relative role of the turbulence and thermo-chemistry is supposed to be naturally captured by the scalar dissipation rate based modelling of turbulent premixed flames, which is reflected well in the results shown in Fig. 6. There is some under prediction of the mean temperature in the calculations using Fluent compared to the previous results, which is due to, as noted earlier, over prediction of the entrainment rate. However, the agreement is good for r/D ≤ 1.0 and the trend is captured correctly for r/D > 1 for the flame F3.

Summary and future scope
In this chapter, a brief overview of various combustion modelling approaches to simulate lean premixed and partially premixed flames is given. The focus is limited to RANS framework because of its high usage in industry currently. The strained flamelet formulation developed recently is discussed in some detail and important details in implementing this model into a commercial CFD code are discussed. The results obtained for pilot stabilised turbulent Bunsen flames using Fluent with strained flamelet model are compared to experimental measurements and earlier CFD results. These published CFD results (Kolla & Swaminathan, 2010b) are obtained using another CFD code employing different numerical schemes and solver methodologies. A good comparison among the Fluent and previous CFD results and the experimental measurements is observed. These comparisons, gives good confidence on the implementation of the strained flamelets model and the associated source and sink terms in the commercial CFD code, Fluent. This initial work serves as a foundation for further studies of lean premixed, partially premixed combustion in industry relevant combustor geometries and, turbulence and thermo-chemical conditions using this modelling framework. Also, this implementation provides opportunities to study self induced combustion oscillations, interaction of flame and sound, interaction of flame generated sound waves with combustor geometries, etc., since a compressible formulation is used in the implementation. The influence of non-unity Lewis number on this combustion modelling framework is yet to be addressed.