State-of-the-art in premixed combustion modeling using flamelet generated manifolds

Flamelet based chemical reduction techniques are very promising methods for efficient and accurate modeling of premixed flames. Over the years the Flamelet Generated Manifold (FGM) technique has been developed by the Combustion Technology Group of Eindhoven University of Technology. Current state-of-the-art of FGM for the modeling of premixed and partially-premixed flames is reviewed. The fundamental basis of FGM consists of a generalized description of the flame front in a (possibly moving) flame-adapted coordinate system. The basic nature of the generalized flamelet model is that effects of strong stretch in turbulent flames are taken into account by resolving the detailed structure of flame stretch and curvature inside the flame front. The generalized flamelet model, which forms the basis on which FGM is built, is derived in Part I. To be able to validate numerical results of flames obtained with full chemistry and obtained from FGM, it is important that the generalized flamelet model is analyzed further. This is done by investigating the impact of strong stretch, curvature and preferential diffusion effects on the flame dynamics as described by the local mass burning rate. This so-called strong stretch theory is derived and analyzed in Part I, as well as multiple simplifications of it, to compare the strong stretch theory with existing stretch theories. The results compare well with numerical results for flames with thin reaction layers, but described by multiple-species transport and chemistry. This opens the way to use the generalized flamelet model as a firm basis for applying FGM in strongly stretched laminar and turbulent flames in Part II. The complete FGM model is derived first and the use of FGM in practice is reviewed. The FGM model is then validated by studying effects of flame stretch, heat loss, and changes in elements, as well as NO formation. The application to direct numerical simulations of turbulent flames is subsequently studied and validated using the strong stretch theory. It is shown that the generalized flamelet model still holds even in case of strong stretch and curvature effects, at least as long as the reaction layer is dominated by reaction and diffusion phenomena and not perturbed too much by stretch related perturbations. The FGM model then still performs very well with a low number of control variables. Turbulent flames with strong preferential diffusion effects can also be modeled efficiently with an FGM model using a single additional control variable for the changes in element mass fractions and enthalpy. Finally FGM is applied to the modeling of turbulent flames using LES and RANS flow solvers. For these cases, the flame front structure is not resolved anymore and unresolved terms need to be modeled. A common approach to include unresolved turbulent fluctuations is the presumed probability density function (PDF) approach. The validity of this FGM-PDF approach is discussed for a few test cases with increasing level of complexity.


Introduction
Detailed numerical modeling of reacting flows has gained a continuous growth of interest in the last few decades. This is related to the fact that the development of new and the improvement of existing combustion equipment is getting more and more important. It has become strikingly clear that we cannot continue with the emission of undesired polluting combustion-byproducts in the way we have been doing in the last century, as this might lead to the destruction of our complete ecosystem. Engines, gas turbines and industrial furnaces play a major role in this emission and it is therefore of the utmost importance to improve the combustion properties of these systems significantly in the near future to avoid a further pollution of the atmosphere.
However, the improvement of combustion processes is a very difficult task: the improvement of one aspect often leads to a deterioration of other combustion aspects due to the complicated nature of combustion systems. For example, the shift of combustion technology during the last few decades toward much leaner low-NO x premixed combustion processes in many applications has generally evolved in a reduction of the flame stability (leading e.g., to flame oscillations, noise, CO-emissions and even flame quenching). To be able to avoid all these problems, adaptations of combustion systems were needed, and to assist these studies, combustion CFD has become a very important tool. Another example is related to the demand for more sustainable combustion which has introduced a broad range of new fuels, some of which show really different combustion phenomena, for example hydrogen enrichment of fossil fuels from bio-sources introduces preferential diffusion phenomena which may have a huge impact on the structure and dynamics of underlying combustion processes. It is clear that combustion CFD is also very important to guide these developments. Therefore, further improvements by ad-hoc measures are insufficient and we have to rely on an in-depth understanding of the processes to be able to improve these systems significantly. Detailed experiments have to be carried out and detailed numerical modeling is needed to meet this goal. However, the numerical modeling of combustion systems is also very challenging from a scientific point of view. The interaction of the fluid flow, turbulence, chemical reactions and thermodynamics in reacting flows is of exceptional complexity. At the moment it has become within reach to model the most important physical aspects in detail, but this is still limited to small academic combustion problems. The modeling of the full details of practical combustion equipment will remain prohibitive in the next few decades, because of current and future limitations in computing power.
This problem asks for special treatments with respect to the modeling of flames. An important way to tackle this problem is by making use of the fact that the chemical time and length scales in most flames are very small. This idea can be exploited in different smart ways to reduce the number of equations to be solved, leading to an enormous reduction in computing effort. In the last decades two main routes have been followed using this idea in combustion science to model the detailed dynamics and structure of chemically reacting flows: chemical reduction techniques and laminar flamelet models.
Chemical reduction techniques -such as conventional reduction [1], Intrinsic Low-Dimensional Manifolds (ILDM) [2] and Computational Singular Perturbation (CSP) [3] -are based on the idea that most of the chemical time scales in the system are very small. If all transport processes are neglected, a time-scale analysis can be performed and the fastest time scales are assumed to be in steady-state. Mathematically, this means that all variables can be stored in a database as a function of a few controlling variables and during run-time only the equations for the controlling variables are solved. Large savings in computing time are reported by most methods with a small loss in accuracy.
Laminar flamelet methods [4] are based on the idea that flame structures are much thinner than most scales of the distortions in the flow, also implying that the chemical reactions are very fast compared to all other time scales. For that reason, the internal structure of the flame front is almost frozen while it moves around in the flow field. The dynamics of the thin flame front is then predicted by using a kinematic equation for the propagation of the flame front, the mixture fraction equation for the mixing and a computational fluid dynamics (CFD) solver for the flow [4].
In 1999, van Oijen and de Goey proposed a combination of these approaches, referred to as the Flamelet Generated Manifold (FGM) technique [5,6]. The FGM approach is more accurate than chemical reduction techniques like ILDM and CSP in colder flame parts, because it also takes transport effects into account in the reduction of chemical kinetics. The correct flame dynamics and propagation speed of the flame are only captured by taking into account the balance between convection, diffusion and reaction, which is essentially not taken into account in the chemical reduction techniques like ILDM/CSP, because they are based on chemical kinetics only. In hot flame parts, chemistry dominates and ILDM/ CSP methods are accurate, but transport processes start playing an important role in the colder flame parts where chemical time scales increase, taken into account only with flamelet based methods like FGM. The FGM method is quite promising, both because of the improvement in reduction efficiency as well as in the accuracy of the model. This method is based on a detailed analysis of (partially) premixed flames in the so-called laminar flamelet and thin reaction zones combustion regime. The basis of the FGM approach is the set of strongly stretched flamelet equations and the related strong stretch  Superscripts 0 stretchless reference condition [-] theory, developed in the last decade. This theoretical model is based on the idea that the most important aspects of the dynamics inside the internal structure of the flame front should be taken into account.
In this paper, we will review the essential ingredients of the strong stretch theory and the FGM method, show results, and indicate the relation between the models and with other existing methods. Various other methods have been developed to couple diffusion and chemistry in reduction methods. Bongers et al. [7] introduced the Phase-Space ILDM (PS-ILDM) method, which analyzes the time scales of a linearized form of the flamelet equations. Reaction Diffusion Manifolds (REDIM) were introduced by Bykov and Maas [8]. Like PS-ILDM, REDIM is an extension of the ILDM framework by including diffusion terms. Furthermore, the FGM method is very similar to the Flamelet Prolongation of ILDM (FPI) technique, which was introduced by Gicquel et al. [9] in 2000. Parallel to the development of FGM, the FPI approach has been applied in many numerical combustion studies (e.g., Refs. [10][11][12]). While FGM and FPI were at first developed for premixed flames and only later applied to non-premixed flames [13,14], the Flamelet/Progress Variable (FPV) approach was introduced by Pierce and Moin [15] for application in non-premixed flames. The FPV method is based on the non-premixed flamelet model of Peters [4,16], but uses a reaction progress variable to describe non-equilibrium chemistry. Pitsch, Ihme and coworkers further developed the FPV method focusing on turbulent non-premixed flames (e.g., Refs. [17][18][19][20]). Both FPI and FPV describe the chemical kinetics with a single reaction progress variable apart from variables to account for mixing and non-adiabatic effects, while FGM has been developed to be used with multiple reaction control variables [6,13,21].
The main goal of this paper is to form one consistent picture of the FGM approach by reviewing older papers, filling the gaps with new derivations and especially to discuss the main physical consequences, which were less obvious in our previous papers. This paper contains two parts. In Part I the set of strongly stretched flamelet equations is derived on which the FGM model is based. This set of equations is further used to derive and validate the strong stretch theory. This theory is the basic formalism for describing the effects of turbulence generated stretch and curvature on the flame propagation rate and is therefore the essential link to the FGM method described in Part II. Inclusion of strong stretch is the main focus because the FGM approach is intended to predict turbulent (partially) premixed combustion in the corrugated flamelet and thin reaction zones regime of the combustion regime diagram [16].

Introduction to strong stretch theory
Partially premixed flames of hydrocarbon-air mixtures generally propagate as thin layers with a thickness δf mainly determined by the balance of diffusion and reaction in the flame front. Such a thin layer contains an even thinner reaction layer with a thickness δr (with δ δ r f ) mainly determined by chemical kinetics. Turbulent flames can then be regarded as thin reaction-diffusion layers wrinkled by the turbulence, while the internal structure and the propagation speed are perturbed by flame stretch and flame curvature due to the turbulent structures. In most applications, flames are wrinkled by all turbulent scales up to the smallest Kolmogorov eddies with scale η, which can be of the order of the flame thickness δf. We therefore consider partially-premixed flames of which the smallest flow scales η satisfy η δ r while η can be of the same order of magnitude as the flame thickness, typically δ η δ f f ≤ ≤10 . This corresponds to flames in the so-called corrugated flamelet regime or the thin reaction zones regime. For this reason we can assume effectively that the reaction layer is infinitely thin (δ r → 0), while the diffusive layer is resolved.
It is important to realize that most existing theories to predict the effects of flame stretch and curvature on the flame propagation speed and flame structure, which have been developed over the years [22][23][24], assume that length and time scales of the perturbations are large compared to the flame scales. These so-called weak stretch theories effectively assume that the perturbations are homogeneous on the scale of the flame front and are therefore suited for the wrinkled flamelet regime ( η δ δ f r ), but making predictions about the effects of turbulence in the regime of main practical interest, where η δ ∼ f , probably unreliable. Still, these theories have been 'extrapolated' to the more practical regimes of turbulent premixed flames by most combustion researchers, simply because a theory for the case with η δ f was lacking. For that reason, the authors developed an extended theory, referred to as 'strong stretch theory', indicating that the Karlovitz number Ka = τf/τη is no longer small but can be of the order 1, with τf = δf/sL the flame time and τη = η/u ′ the flow time scale, where sL is the laminar burning speed and u ′ is the velocity fluctuation of the surrounding (turbulent) flow.
In the regime with η ≈ δf and δr→0, it is obvious that the flow perturbations cannot be assumed homogeneous inside the diffusive layer. The flame stretch rate and the flame curvature should therefore be resolved inside the flame front. This is the main difference between the existing weak stretch theories and our strong stretch theory. We will not assume that δf→0 thus we cannot treat the flame as an interface as in other models. We will treat the reaction layer as an interface though. There are more differences between the weak and strong stretch theories: weak stretch theories are derived for single-step kinetics between a few major species having a large activation energy while the strong stretch theory is valid for (realistic) multiple-species transport and chemistry as we will see, as long as the reaction chemistry is fast enough so that δ δ r f . This enables application to (numerical) studies of flames using realistic transport and chemistry models, like FGM. To resolve the structure of the diffusive layer, we consider a progress variable and we will evaluate the kinematic evolution of a series of isocontours of the progress variable, thereby taking the detailed interaction of fluid dynamics and transport processes inside the flame front into account. These aspects are the main ingredients of the strong stretch theory developed over the last decades [25], which forms the basis of the FGM method presented in Part II and further. In Part I, we will review and derive the theory step-by-step.
We will first present the equations governing (partially) premixed flames in Section 2.1.2. We then introduce a flame-adapted coordinate system and transform the equations to these curvilinear coordinates in Section 2.1.3. After introducing the mass-based flame stretch rate in Section 2.1.4, the transport equations are split in a kinematic equation describing flame propagation and a set of flamelet equations for the internal structure and propagation speed of the resolved structure in Section 2.1.5. Using integral analysis we then analyze the influence of strong stretch on the flame propagation rate in Section 2.2, while the theory is validated in Section 2.3.

Governing equations
In this section we summarize a mathematical model for laminar, premixed flames. Typically, we consider the combustion of hydrocarbons in air. The species in the flame are numbered 1 through Ns. Species Ns is nitrogen (N2) and is present in abundance. The governing equations for such flames are the conservation equations of mass, momentum and energy of the gas mixture, the balance equations of mass for the species involved and the thermal and caloric equations of state [26,27]. These can be written in the following form: ∂ ∂ ( )+ ∇⋅( )= −∇ + ∇⋅ + The independent variables in Eq. (1) are the density ρ, the (massweighted average) flow velocity v, the pressure p, the specific enthalpy h, the species mass fractions Yi and the temperature T. The variables ρ, Yi, T and h are the combustion variables, whereas v and p are the flow variables. Other variables/constants in Eq. (1) are the stress tensor T , the gravitational acceleration g, the thermal conductivity λ, the specific heat at constant pressure cp, the preferential diffusion flux of enthalpy Jh, the diffusion velocities Vi, the reaction rates ωi, the ambient pressure p amb , the universal gas constant R, the species molar masses M i s, , the average molar mass M and the specific enthalpies of formation h i,ref at a reference temperature T ref . In the sequel of this paper we will neglect gravity. Note that only the mass balance equations (1d) for the first Ns − 1 species are required, since the mass fraction Y Ns follows from the constraint ∑ = = i N i Y 1 1 s . In the remainder we assume that the fluid is an ideal Newtonian fluid. Furthermore, using the Hirschfelder-Curtiss approximation, we obtain for the diffusion velocities Vi the relations: are the mixture averaged diffusion coefficients and Lei the Lewis numbers of the species concerned. Here we replaced the gradient of mole fractions by a gradient of mass fractions, assuming that the gradient of mean molar mass is small [28]. For premixed flames burning in air, this is a valid assumption because M is mainly determined by the abundant N2. Furthermore, thermal diffusion is neglected in Eq. (2). Soret and Dufour effects are second-order diffusion effects, which can be neglected in most premixed hydrocarbon-air flames. Since thermal diffusion mainly affects the transport of light species, it only plays an important role for fuel mixtures with very large quantities of H2 [29]. The effect of this assumption is investigated in Section 3.4.1.2.
For the heat flux vector q we only take into account conduction and enthalpy transport by diffusion, i.e., where hi are the species specific enthalpies. Combining the relations in Eqs. (2) and (3) we obtain the following expression for the preferential diffusion flux of enthalpy s , . Finally, the low-Mach number approximation holds, i.e., p is set to the constant ambient pressure in the equation of state (1e) and the term Dp/Dt is neglected in the enthalpy equation (1c); see Ref. [30].
In the following of this paper we frequently need the element mass fractions Zj. Assuming that the elements are numbered 1 through Ne, the element mass fractions are defined by , , , where J Z j is the mass diffusion flux of element j given by for more details see Refs. [25,31]. These conservation equations are, of course, not independent of the set (1d). Note that the balance equations (7) are only required for the first Ne − 1 elements, since ∑ = = j N j Z 1 1 e .

Flame adapted coordinates and flame kinematics
In this subsection we present our flame adapted coordinate system; see also Ref. [25]. The rationale behind our definition is the observation that premixed flames have a local one-dimensional structure, i.e., the gradients of many combustion variables are quite large along certain curves in the flame front, whereas in surfaces perpendicular to these curves the variables are virtually constant. Therefore, we consider a premixed flame as the x, t for some suitable combustion variable Y , taking the values Y u and Y b in the unburnt and burnt gases, respectively. The variable Y is called the progress variable and can be, say, one of the species mass fractions or the temperature. We assume that ∇ ≠ Y 0 everywhere in the flame. Consequently, we can identify iso- , which we refer to as flame surfaces. We define the curvilinear coordinate system ξ ξ ξ ξ = ( ) i.e., coordinate surfaces ξ 1 = const coincide with flame surfaces and the ξ 1 -coordinate lines are orthogonal to these. We refer to the ξ 1coordinate lines as flamelet paths. In each flame surface, (ξ 2 ,ξ 3 ) is a curvilinear coordinate system, not necessarily orthogonal, and is not specified unless stated otherwise. Note that the flame coordinates are not defined in the burnt and unburnt gas mixtures. For time-dependent flames the iso-surfaces will move in the spatial domain, resulting in a time dependent coordinate system, i.e., ξ = ξ(x,t). Moreover, we introduce the variable τ as the time corresponding to the curvilinear coordinate system. Obviously, τ = t. In Section 2.1.5, we will consider the transformation from the laboratory system (x,t) to the curvilinear flame adapted coordinate system (ξ,τ) and we reformulate the conservation laws in terms of these coordinates. It is obvious that the moving coordinate system and the kinematic behavior of the flame are closely related. Note that due to unsteady flow, the flame surfaces move in the spatial domain with velocity vf given by v v n n f d : , where sd is the displacement speed, and n the (local) unit normal vector to each flame surface, which we assume to be directed toward the burned gas mixture, i.e., −sdn is the velocity of the flame surfaces relative to the flow. In our previous papers (e.g., Refs. [25,32,33]), we referred to this velocity as the local laminar burning velocity using symbol sL. In order to align with literature, we use here the symbol sd and refer to it as the displacement speed. The laminar flame speed sL is a global property of the flame front and is defined as the displacement speed at the unburnt side of a onedimensional stretchless flame, s s L d u = , 0 . As we intend to resolve the internal structure of the flame front, it is regarded as a collection of flame surfaces and the motion of each surface is described by the (local) kinematic condition Alternatively, differentiating the first relation in Eq. (9) with respect to τ, we find where x x = ∂ ∂τ is the velocity of the moving system. Combining the conditions in Eqs. (11) and (12) Y 0, i.e., the flame velocity vf and the velocity x of the moving coordinate system have the same component normal to the flame surfaces. On the other hand, the tangential component of the velocity x of the moving coordinate system is not determined by Eq. (9), and for this we simply choose a v x , , a α being the contravariant unit vectors in the curvilinear coordinate system (see The consequence of this choice is that the curvilinear coordinate system moves with the flame velocity vf with respect to the laboratory system. This motion is described by: which follows readily if we substitute Eq. (10) in Eq. (11). If the conservation laws are reformulated in terms of this moving coordinate system, the flame front in terms of Y ξ 1 ( ) does not change any more, although the internal structure (in terms of other combustion variables) may change due to unsteady effects (described by partial derivatives with respect to τ), flame curvature and flame stretch (see Fig. 1).

Mass-based flame stretch
Traditionally, the flame stretch rate is defined as the relative rate of change of an area element d d d S = σ ξ ξ , (14) due to flow straining and flame motion. It is clear that this is a very useful definition in case of an infinitely thin flame sheet. This concept needs extension when the diffusive layer is resolved. For a flame of finite thickness, enclosed between the iso-surfaces Y Y 1 ) are defined as curves crossing the entire flame from the unburnt side to the burnt side in the direction of the normal n, i.e., they satisfy an initial value problem of the form where C is some function depending on the flamelet coordinate ξ 1 only, which we assume to increase when going from unburnt to burnt boundary of the flame.
It is important to investigate how much the mass flow rate changes along a flamelet path. For this reason, de Goey et al. [34] introduced the mass-based stretch rate K, defined as the fractional rate of change of the mass M(t) contained in an infinitesimal small control volume V(t) moving in the flame with velocity vf, i.e., . ρ (16) Applying Reynolds' transport theorem we find We now have the following expression for K Combining this relation with the continuity equation (1a) and substituting vf = v − sdn we get where m = ρsd is the mass flux relative to the flame surfaces. This relation gives rise to the following physical interpretation. Consider a patch S 1 t ( ) on the flame surface Y Y x, t ( ) = 1 and the patch S 2 t ( ) stating that stretch is responsible for the change of mass along flame paths.
Using the notation and nomenclature introduced in Appendix A, we can rewrite the expression (18) for K in the curvilinear coordinate system. If we take q = ρ and u v x = = f in relation (98e), see Appendix A, we can derive from Eq. (18) the following expression for the stretch rate where g is the Jacobian of the transformation x ξ, given by : . x (22) an area element on the flame surfaces and ds = h1dξ 1 is the arc length element along flame paths. Substituting Eq. (22) in Eq. (21) we find is the fractional rate of change of the density of a fluid parcel moving with the flame. The first two terms in Eq. (23) represent the fractional rate of change of the volume V(t). As indicated in Ref. [35], the change in volume can be regarded as the sum of changes in flame surface area and flame thickness. The flame thickness changes due to unsteady phenomena as well as local flow straining and curvature effects. This contribution is obviously not present in the traditional definition (14) of Kσ for an infinitesimally thin flame interface. Alternatively, starting from Eq. (19) and relation (98b) we find = ⋅ = a n . Substituting the expressions for m 1 and g in Eq. (24) and applying the relation g g 11 11 1 = , this equation can be rewritten as once more confirming that K determines the change of mass flux along the flamelet paths.

Flamelet equations
In this subsection we will reformulate the combustion equations in terms of the flame adapted coordinates. Let us start with the equation for the progress variable Y . Combining this equation with the expression for ρK in Eq. (18) and the kinematic condition (11) and writing all differential operators in ξ-coordinates, we obtain the quasi-one-dimensional equation In the derivation of Eq. (27) we additionally used the relations m α = 0, g 1 0 α = and ∂ ∂ = Y ξ α 0 for α = 2,3. Substituting again the expressions for m 1 , g and ds = h1dξ 1 as arc-length in Eq. (27) and applying the relation g g 11 11 1 = , we obtain the flamelet equation for Y , see Ref. [25], Equation (28) describes the convection-diffusion-reaction balance of Y along flamelet paths, however, with an additional source term −σρKY describing the loss/gain of mass due to flow along the flame surfaces.
The same procedure can be applied to the scalar conservation equations (1d), (1c) and (7). The result is the following set of quasione-dimensional differential equations, referred to as the Strongly Stretched Flamelet Equations (SSFE), Note that the flamelet equations for mass and for the progress variable are steady. All other equations do have a time derivative indicating the unsteady behavior of the corresponding variables with respect to Y . If a scalar is steady in the flame adapted coordinate system, it means that this variable has the same transient behavior as Y . Note that the time derivative in the flame adapted coordinate system can be written as where we introduced v f i as the velocity of the iso-surfaces of Yi, which follows from the kinematic condition for Yi, similar to Eq. (11). Equation (30) expresses that the time derivative arises because the local iso-surfaces of Yi move with respect to the iso-surfaces of Y . Furthermore, the terms Qh, Q Yi and Q Z j , which describe transport in the flame surfaces, arise because the local iso-surfaces of h, Yi and Zj generally are not parallel to the iso-surfaces of Y . These terms read where ∇ ⋅ t and ∇t denote the internal divergence and gradient operator, respectively, restricted to the flame surfaces; see Appendix A.
The Q-terms in Eq. (29) describe conduction/diffusion along the flame surfaces and are presumably small. In order to investigate this issue more precisely we have to make the flamelet equations (29) and the relations (31) dimensionless. As an example consider equation (29d) together with relation (31b). We consider combustion in the corrugated flamelet regime and assume that the surrounding flow, which distorts the flame layer, is characterized by a velocity where q is an arbitrary variable/parameter in Eq. (39) or Eq. (43). The resulting dimensionless equation reads (omitting all asterisks): where In the corrugated flamelet regime [16], Ka and ε are both much smaller than 1. In this regime, the time derivative, the stretch term and the Q-term can be neglected. In the thin reaction zones regime Ka~1, the Q-term is still much smaller than the stretch term and is negligible. However, the time derivative and stretch term are of order 1. In practice the order of magnitude of the time derivative is smaller because it is related to a difference in velocities (Eq. 30).
In flames that are not fully premixed or adiabatic, gradients of Zj and h may occur in the flame surfaces. Such gradients will lead to diffusion along the flame surfaces and may result in a nonnegligible Q-term. When the typical length scale of the gradient is comparable to or smaller than the flame thickness (e.g., Z Z j j ,u f ∇ δ ), then the Q-term cannot be neglected. In Part II we will show that the magnitude of the Q-term is proportional to the scalar dissipation rate in the flame surface and that this may result in non-negligible Q-terms in partially-premixed flames and flames that stabilize on a burner by heat loss. A quantitative analysis of all terms in flamelet equation (29d) is given in Ref. [36] for the turbulent expanding flame kernel shown in Fig. 1. For this premixed adiabatic case, the Q-terms were found to be negligible compared to the other terms in the flamelet equation. In the MFGM method [37] these terms are included in the manifold by computing multidimensional flamelets.
The analysis shows that the transformation to the flame adapted coordinate system leads to a quasi-one-dimensional system of equations (29), describing the internal structure of the flame as function of the arc-length s along flamelet paths, which are not straight lines. The mass burning rate m is the eigenvalue of this one-dimensional system. If this so-called 'flamelet' system would be solved together with the kinematic equation (13) for the set of flame isocontours Y in the range Y Y u b , ( ) this would be equivalent to solving the original set of transport equations. For infinitely thin flames, a single kinematic equation for a reference iso-plane would suffice. This is known as the G-equation model [16], where a kinematic equation is solved for a single iso-surface G = G0, where G is some arbitrary scalar field. For flames with a resolved diffusive layer an extension is necessary. However, as the numerical handling of this set of (loosely coupled) kinematic equations is very tedious [38], we will solve the transport equation (26) for the control variable Y instead (together with the flamelet system for the internal structure) to track the resolved flame volume. If the kinematic equation (13) is multiplied by ρ and added to Eq. (28) we recover the original transport equation (26). This is thus again equivalent to solving a set of kinematic equations. This approach is the starting point of the FGM method, presented in Part II.

Integral analysis of flame dynamics
The SSFE set of equations forms the basis of the FGM formalism for the modeling of flames in the flamelet and thin reaction zones regime and is described further in Part II. This set governs the impact of small and large scale fluctuations in the flow on the flame dynamics, viz. the influence of the K(s) and σ(s)-fields on the mass burning rate. To analyze the effects of stretch, curvature and preferential diffusion further, a strong stretch theory has been derived. This theory makes use of the fact that reaction layers are mostly thinner than all other scales and in the case of an infinitely thin reaction zone, the SSFE set can be used to derive analytical expressions for the influence of perturbations on the mass burning rate. These expressions in turn can be used to validate FGM which solves the SSFE set numerically without assuming that the reaction layer is infinitely thin. This strong stretch theory is derived and analyzed in this section and approximations and comparisons with other stretch theories and with numerical computations using detailed chemistry are presented in Section 2.3. Thus, we will analyze the solution of the SSFE system (29). In Section 2.2.1 we first neglect the Q-terms and the stretch terms (and insert Q = K = 0). The solution for the mass burning rate will be referred to as the mass burning rate of 'stretchless flames'. Note that this is the mass burning rate of a curved flame, since the flamelet path is generally not straight and the flame contours are curved (σ ≠ 1). This does not necessarily represent a physical flame because important stretch terms have been neglected. This procedure will just lead to a mathematical expression for the mass burning rate of stretchless flames and is a first step to derive the physically relevant stretched mass burning rate. Subsequently, in Section 2.2.2, we will analyze how the mass burning rate of stretchless flames changes into the full stretched mass burning rate if the right-hand side stretch terms are included.
We will solve the system using the so-called integral analysis method, first introduced by Law and coworkers [39]. This means that the SSFE system is integrated along flamelet paths connecting the unburnt and burnt mixtures. The final expression for the mass burning rate is the same as the expression obtained using matched asymptotic expansions [40,41], in case of an infinitely fast reaction taking place at a flame sheet.

Mass burning rate of stretchless flames
Let us first derive an expression for the mass burning rate of stretchless flames, i.e., for K = 0. Integration of the conservation equations for mass, enthalpy and element mass fractions of the set (29) from the unburnt to the burnt gas mixture gives: where it has been used that the diffusive fluxes all vanish in the unburnt and burnt gases. Equations (34) simply indicate that the mass, enthalpy and element composition are conserved over the flame area, i.e., h h b u 0 = and Z Z j j , , b u 0 = , as expected. This means that the stoichiometry of the burnt mixture near the reaction zone is equal to the stoichiometry of the unburnt mixture while the temperature equals the adiabatic temperature of the flame (as for a planar unstretched flame without heat loss).
Let us now turn to the mass burning rate m b 0 of the stretchless flame, which can be computed from the quasi-one-dimensional equation for Y : where we introduced the terms A , D and S: Inserting the approximation (38) The validity of the approximation in Eq. (38) is related to the fact that the reaction layer is much thinner than the preheating zone, so that the major contribution to the integral over A D where we used Eq. (34a). If the expressions for the integrals are substituted into Eq. (39), we obtain the following approximation for the mass burning rate: , σ σ λ ω (41) This approach to compute m b 0 gives the same result as the large activation energy asymptotics (LAEA) approach [27,42]. In that case, a scaling procedure is carried out for the two (preheating and reaction) zones in the flame and after expanding all terms in series of a small parameter ε inversely proportional to the activation energy, i.e., ε ∝ 1/Ea, solutions in the different zones are matched for the different orders of the parameter ε. A similar procedure can be carried out with Eq. (35) as follows. The solution in the preheating zone follows from Eq. (35) if we set S s ( ) = 0 , leading to expression (39) for the diffusive flux D b of Y at the edge of the diffusive layer and the reaction layer; cf. Eq. (38). Convective transport is small and neglected in the reaction zone. So after setting A s ( ) = 0 and multiplying Eq. (35) with D s ( ) and integrating, we find the expression (40a) for the diffusive flux D b at the edge between the two zones. Matching them gives Eq. (41). Physically, this can be interpreted as follows. The square root of the right-hand side integral (39) is proportional to the total amount of Y per unit flame surface area and per unit time produced or consumed in the reaction layer. Equation (39) thus expresses that the total amount of Y produced in the reaction layer leads to a diffusive flux D b of Y at the edge of the reaction layer. This flux of Y into the preheating zone has to be distributed somehow between convection and diffusion, which is expressed by Eq. (40a). Matching them fixes the mass burning rate needed for that.
It should be realized that the right hand side of Eq. (41) is in general not equal to the adiabatic mass burning rate m b,1 0 of a flat adiabatic stretchless flame, due to the factor σ σ b ( ) 2 . For an infinitely thin reaction layer, though, this factor becomes unity and this expression does reduce to the expression for the adiabatic mass burning rate of a planar flame. That curvature does not influence the mass burning rate of a flame with infinitely thin reaction layer is also clear from a physical point of view: the production or consumption of Y in the reaction layer leads to a flux D b of Y at the edge of the reaction layer directed toward the preheating zone. If this reaction layer is infinitely thin, this flux cannot be influenced by curvature. Convective and diffusive processes are responsible for the distribution of this 'flux' in the preheating zone. Without 'losses' (like flame stretch) in the preheating zone, this flux is equal to Eq. (40a) and independent of the structure of the preheating zone. This has already been shown by de Goey [32] who proved that the laminar burning velocity is not influenced by transport processes inside the preheating zone, as long as no losses are present. Equation (40a) expresses 'conservation' of Y in the preheating zone. It should be realized that Eq. (41) also shows that the way the distribution inside the preheating zone takes place is not relevant, since the expression depends on the diffusivity near the reaction layer but does not depend on the diffusivity inside the preheating zone where the source term is virtually 0. We can conclude from this that the mass burning rate of an unstretched flame is not dependent on curvature, i.e., flat, spherical or cylindrical flames all have the same m b,1 0 at the reaction layer in case of an infinitely thin reaction layer. As the mass burning rate of a curved flame changes through the preheating zone with m s m s 0 0 σ σ , this does not hold at other positions in the flame front, like in the unburnt gases. This coincides with observations from numerical calculations using detailed chemistry performed by Groot [38], who concluded that the mass burning rate is independent of curvature only at the 'inner layer' coinciding with the point of highest heat release rate.
Note that the expression Eq. (41) for m b 0 only depends on the unburnt state Y u due to the 'friction' factor 1 0 Y Y b u − , which has to be conquered by the flame, in the sense that Y has to change from Y u to Y b 0 by convection, diffusion and reaction. Furthermore, we also explicitly emphasized that the mass burning rate m b 0 also depends on the enthalpy h b 0 and element composition Z j,b 0 in the equilibrium state. This follows from the observation that the factor 1 Le  (41) then indicates that the integral over Y is independent of the initial composition and depends only on the equilibrium state. The local equilibrium state is completely described by the pressure (assumed to be constant), total enthalpy and element composition in the reaction layer, i.e., by h b 0 and Z j,b 0 for j N = 1 2 , , , … e . Though the FGM method is based on the assumption of the existence of a lower-dimensional manifold (see Part II), the 1D manifold assumption used here is only needed for the present integral analysis.

Mass burning rate of stretched flames
In this subsection we consider the change in the mass burning rate if stretch terms are included. As in the previous subsection, we first study the conservation equations for mass, enthalpy and element mass fractions.
In the same way, integration of the enthalpy equation and element mass fraction equations in Eq. (29) along the flamelet and using Eq.
Note that mass, enthalpy and elements are not conserved anymore, the 'change' being expressed by the integrals of mass change ρK, enthalpy change ρK(h − hu) and element mass fraction change ρK Z Z  , , with hT the thermal enthalpy. We also introduced normalized variables as follows The Karlovitz integrals express the relative change of temperature and species mass fractions over the flamelet path due to flame stretch. If a Karlovitz integral Kai is of order unity, this means that the 'loss' (positive integral) or 'gain' (negative integral) of Yi due to convection along the flame surfaces (flame stretch) is of the same order of magnitude as the change of Yi along the flamelet path.
Equations (47a) and (47b) describe the influence of preferential diffusion and flame stretch on the local enthalpy and element composition of the burnt mixture. These equations have been derived without assumptions about the magnitude of the Karlovitz integrals and therefore their validity is not limited to weak stretch. It should be noted from the definitions (45a) and (45b) of the Karlovitz integrals that the influence of flame stretch ρK in the preheating zone is effectively damped exponentially by factors like Y i and h T . These factors are equal to 1 in the burnt gases and tend to 0 exponentially for s → su. The enthalpy hb and the element mass fractions Z j,b in the burnt mixture determine the local stoichiometry and equilibrium composition in the burnt mixture. These quantities, following from Eqs. (47a) and (47b), have an important influence on the local mass burning rate mb of the stretched flame, because mb is determined to a large extend by the mixture composition and enthalpy in the reaction layer, close to the burnt mixture (see also Eq. (41)). The precise description of this influence will be studied hereafter.
Thus, let us now turn to the evaluation of the mass burning rate for stretched flames. The mass burning rate is again determined from Eq.
while D and S are still defined in Eq. (36). The derivation of the expression for the mass burning rate is analogous to the derivation presented in the previous subsection, with the expression for A s ( ) replaced by Eq. (48). For the integral over A in Eq. (40a) we so that we find for the mass burning rate )the expression for the stretchless mass burning rate as defined in Eq. (41), however, with h b 0 and Z j,b 0 replaced by their stretched counterparts hb and Z j,b . This term is again a function of all combustion variables Yi and T. Once again, the integral ) effectively runs over the thin reaction layer, where the source term is non-zero and where the composition approaches the equilibrium point, following the one-dimensional manifold described by Y i Y ( ) and T Y ( ). The second integral in the right-hand side of Eq. (50) indicates that the flux D b of Y into the preheating zone is not conserved inside the preheating zone. Flame stretch induces a convective 'loss' during the distribution of Y inside the diffusive layer. This transport 'loss' due to flame stretch is apparent from Eq. (49). Using Eq. (41) for m b 0 and the definition of the Karlovitz integral (45a) for Ka Y then gives: denotes the vector of element mass fractions.
So far we did not choose a specific variable for Y . The theory is independent of its choice but for the case of Y = h T , the expressions simplify. Therefore, we adopt this choice in the following, so that Eq. (51) reduces to If we finally use (52) in Eqs. (47a) and (47b) and introduce w i The set of equations (53) for Δhb, ΔZb and mb describes the enthalpy, element composition and mass burning rate of stretched and curved premixed flamelets and forms the basis of the subsequent analysis of the flamelet system. This set is also valid in case of strong stretch rates: Eqs. (47a) and (47b) follow rigorously from the flamelet equations and to derive Eq. (52) we only have assumed that the reaction layer is infinitely thin. However, the system (53) with hot containing higher-order terms in Ka i 0 , Ka T 0 . The super-script°indicates that stretchless expressions are taken, because we are considering lowest order stretch contributions. For the enthalpy/ element change we can now substitute Equations (54) and (55) are compared with theoretical expressions derived using LAEA in Section 2.3.1.

Physical interpretation
Although Eq. (53) describes the effects of curvature and flame stretch on the mass burning rate of premixed flames, it is hard to visualize the physical background of these expressions. We will therefore investigate some schematic examples by which it might be possible to understand these expressions in a better way. We will consider 1. a planar unstretched flame, 2. a spherical steady flame, 3. a planar stagnation flame with unity Lewis number, 4. a planar stagnation flame with non-unity Lewis number 5. a curved flame propagating in a homogeneous flow field of a mixture with non-unity Lewis number.
In all cases we will resolve the diffusive structure of the preheating layer, but the reaction layer will be treated as an interface as before. Our theory takes care of multiple-species transport effects through the use of a Lewis number for each species in the reaction mechanism. However, for the sake of simplicity, we consider in the following examples only the diffusion flux JF of one single species (fuel, mass fraction Y ) and the diffusion flux JT of thermal energy.

Planar unstretched premixed flame.
We consider a steady planar flame in a homogeneous flow field, like the one stabilized on a Bunsen burner through which a plug flow is issuing. As an example see Fig. 2a. The flow velocity can be split in a component perpendicular to the flame and a component parallel to the front. For the sake of simplicity we assume that the density is constant, so that the flow remains undisturbed by the flame. The flame is not stretched because the parallel flow component is constant so that ∂ ∂ = v ξ ξ 2 2 0 for a flame in a two-dimensional geometry. The consumption of fuel in the reaction layer leads to a flux from the reaction layer into the preheating layer (see Fig. 2a). This flux is transported further by diffusion and convection inside the diffusive layer. Consider a small volume element around the flamelet (red dashed area) in which this flux enters. (For interpretation of the references to color in this text, the reader is referred to the web version of this article.) As the convective flux of Y due to v ξ 2 into and out of the small volume element is the same, there is no loss or gain of fuel in the preheating zone and it is conserved. This means that we obtain the mass burning rate of an adiabatic planar stretchless flame, denoted by m b,1 0 . The kinematic equations of all flame isoplanes of Y s ( ) inside the diffusive layer have the same mass burning rate m b,1 0 , because ∂m/∂s = 0 (take σ(s) = 1 in Eq. 29a).

Steady spherical flame.
A spherical flame may be generated from a spark in a homogeneous mixture. This flame propagates outwards and grows in size. However, if a sink of mass would be positioned in the origin, the flame will become steady at a certain radius rf (see Fig. 2b). At this point, all iso-planes of Y are spherical and curved. However, for an infinitely thin reaction layer the amount of fuel consumed in the reaction sheet, and therefore also the flux Db into the diffusive layer is exactly the same as for a planar flame. If a small volume element is formed around the flamelet (red

Planar stagnation flame with unity Lewis number.
Planar premixed stagnation flames may be stabilized in a flow of two opposing nozzles, often applied to measure the adiabatic burning velocity of combustible mixtures (see e.g., Ref. [45]). A simple representation is given in Fig. 3a. For the case of unity Lewis number it is clear that the fuel mass fraction and the temperature diffuse with the same strength inside the preheating zone, indicated by the equal length of the diffusion flux vectors of fuel, JF, and thermal energy, JT. There is no reason for the enthalpy or mixture fraction to change and as a result, the flame attains the same adiabatic flame temperature and equilibrium composition as the planar flame of Fig. 2a. The consumption rate inside the reaction layer also remains unchanged and the flux Db of Y is the same as for the first example. There is no sign of preferential diffusion in this case. However, the distribution of fuel inside the preheating zone is different now. The flame is strained so that ∂ ∂ > v ξ ξ 2 2 0, which means that the parallel convective flux of Y in the small volume element around the flamelet is not conserved. As the flux Db into the volume remains the same, while mass is leaking out of the volume element, this requires a smaller convective mass transport mb along the flamelet path to distribute the fuel (see Eq. (49)). The mass burning rate thus decreases according to Eq. (52). In the remainder we will refer to this as the 'direct' effect of stretch in the flame. It exists even if the flame temperature and reactive mixture composition remain the same and is directly related to the loss of mass inside the diffusive layer, referred to as flame stretch.

Planar stagnation flame with non-unity Lewis number.
If we consider the same flame as in Fig. 3a, but now propagating in a mixture having a Lewis number not equal to 1, the situation changes as follows (see Fig. 3b). The amount of mass leaking away out of the volume element around the flamelet is now the same as in the previous case, so the 'direct' stretch effect is also present. However, due to the fact that diffusive fluxes are not the same, this leads to a change in element composition and enthalpy in the preheating zone. This can be understood by the simple example displayed in Fig. 3b, where two diffusive fluxes JF and JT have op-posite directions but a different magnitude. This leads to a change in element composition/enthalpy in the volume element. Diffusive fluxes are all aligned and not responsible for the change of total enthalpy/elements that the volume element and as such the flame front encounters. The local enthalpy increases inside the diffusive layer in the current example of Fig. 3b because there is a larger thermal energy flux JT into the diffusive layer than the energy carried by the flux JT out of the dashed area. Note that this unbalance of diffusive fluxes in stretchless flames leads to a local change in enthalpy/element composition inside the diffusive layer but these quantities remain conserved in a flamelet because there is no loss from the flamelet due to homogeneous flow. Since flame stretch is  related to an inhomogeneous flow perpendicular to the flamelet, the non-constant convective flux ρ ξ v 2 in the current stretched flame of Fig. 3b is responsible for the loss/gain of enthalpy/elements out of the volume element through the sides of the red dashed area. Enthalpy and elements are now not conserved along a flamelet. This means that the flame temperature changes and as a result also the amount of mass of fuel consumed in the front changes, which in turn changes the flux Db into the diffusive layer. The stretched flame in this example thus not only encounters a 'direct' effect of flame stretch but also an additional effect of 'preferential diffusion', and the mass burning rate behaves as described by Eqs. (54) and (55). 0 near the right flamelet, assuming the coordinate ξ 2 increases going from left to right along the flame surfaces. This again introduces a direct gain/ loss, previously indicated as the direct effect of flame stretch. As before, the diffusive fluxes have opposite directions but are of unequal strength, leading to a local change in enthalpy/elements in the diffusive layer. The enthalpy locally increases in the preheating zone in the current example as in the previous example. Again, the non-zero flow along the iso-planes of Y perturbs global conservation. Enthalpy and flame temperature increase in the left flamelet, while the opposite is true for the right flamelet. This flame thus encounters a 'direct' stretch effect and a 'preferential diffusion' effect, and the mass burning rate is again described by Eqs. (54) and (55).
The examples discussed in the previous subsections indicate how flame stretch influences phenomena taking place inside the diffusive layer of a flame front. It is clear that transport effects inside this layer dominate the physics. Transport effects inside the reaction layer are less important if this layer is much thinner than the preheating zone, simply because there is less time and space for transport phenomena to be of significance. We will validate these phenomena using numerical modeling in the following. However, it must be clear that this is not the first and only analysis of these phenomena, although this description is more detailed than most other theories, which are based on a flame sheet assumption and one-step chemistry (e.g., Refs. [46][47][48]). For some special cases it is possible to compare our extended strong stretch theory with existing theories derived by others. This will be one of the objectives of the first part of the next section.

Validation of strong stretch theory
In the previous section, we formulated general equations for computing the changes in enthalpy, element mass fractions and mass burning rate in strongly stretched premixed flames. The preheating zone had been resolved, while the reaction layer was reduced to a sheet. In this section, we will validate these expressions comparing them with existing theories and numerical modeling.

Comparison with large-activation energy asymptotics
In this subsection we will show that the strong stretch equations previously derived reduce to the equations available for weak stretch theories, based on large-activation energy asymptotics (LAEA). To be able to make a comparison possible, we will restrict the analysis to (a) weak stretch rates, (b) planar flames (σ(s) = 1), (c) a onestep irreversible reaction F P → , with one rate-determining lean species F and a single constant Lewis number Le Le Le F P = = , and (d) a constant coefficient λ/cp. For weak stretch, it can be shown that the stretch rate inside the flame zone is close to spatially homogeneous. For instance, the stretch rate of an expanding spherical flame K = 2vf/r, where the flame velocity vf(s) of the different isoplanes Y s ( ) is almost the same, while the radius r is large, much larger than the flame thickness for weak stretch. This means that we can assume that K(s) = Kb is constant.
Note that we then have a single element mass fraction Z Y Y Y = + = F P F ,u which is constant as F and P have identical Lewis numbers. This means that Y s Y s Δh c T T In the last equation we introduced the familiar Zeldovich number Ze Ea being the activation energy of the onestep reaction. The Zeldovich number in Eq. (57) indicates the sensitivity of the mass burning rate to changes in enthalpy (or temperature). It has been assumed that the mass burning rate m b 0 only changes as a function of the flame temperature T b 0 , given by m [49], as is often applied in LAEA theories.
Note that the expression (57) for m m b b 0 contains two contributions, viz. the term Ka T 0 related to the 'direct' stretch effect introduced in the previous section, and the term proportional to the Zeldovich number, which is related to the 'differential diffusion effect'. We will show that this last term vanishes if Le = 1, while the 'direct' stretch always remains unless K = 0.
We will consider three different theoretical cases in the following. The first case focuses on the (academic case of) a constant density ρ 0 (s) = ρb and the other two cases include variable density ρ 0 In the second case, we will look at the mass burning rate at the burnt side of the flame (s = sb). However, since there are also expressions available for the unburnt flame edge, this will be considered in the third case. This requires some extra attention because false interpretation and application of these expressions by some authors has been detected. This analysis also gives some directions on how to apply these expressions for interpreting experimental data of stretch induced phenomena.

Ka
Ka . ρ δ Then, Eq. (57) can be rewritten as: where M b is the so-called Markstein number, given by: In 1962 Barenblatt et al. [49] studied the dynamic flame response with constant density, leading to a Markstein number for the unburnt and burnt gases equal to Ze . To lowest-order in 1 0 Ze , Eqs. (60) and (61) are therefore equal.
2.3.1.2. Variable density case: burnt flame boundary. We will repeat the analysis above, but now for the case that the density ρ 0 (s) in the expression (45a) for the Karlovitz integrals varies as in a nonstretched deflagration wave with negligible pressure variation: where we introduced the thermal expansion coefficient . Now, substituting Eqs. (58) and (62) where the relation ρ ρ θ  Le Combining Eqs. (45b) and (64) In 1985 Clavin derived a similar expression using large activation energy asymptotics [40]. Clavin assumed in his analysis that Le − 1 is of order 1 0 Ze . To lowest order in the reciprocal Zeldovich number, Eq. (65) is therefore equal to the result of Clavin [40].

Variable density case: unburnt flame boundary.
For the mass burning rate mu in the unburnt mixture, an equation equivalent to Eq. (54) can be derived. This can be done by following the same procedure as for mb, which has been described in ten Thije Boonkkamp et al. [34,35]. This procedure then yields where Ka T * is given by , the normalized temperature with respect to the burnt mixture instead of the unburnt mixture, as used in Eq. (45a). We again take the stretch rate constant in the flame but use the unburnt flame edge as reference, i.e., K(s) = Ku. The preferential diffusion terms in Eq.
Using ρ 0 (s) from Eq. (62) where su denotes the position in the unburnt mixture where the unburnt flame boundary is chosen. Substituting Eqs. (70) and (64) in Eq. (68), taking the limit su → −∞ in the logarithmic term and retaining the term su/δf in Eq. (94), then gives the Markstein number M u u s ( ) in the unburnt gases: When this equation is evaluated at su = 0, it is equivalent to the equation of Clavin and Williams [41] up to order 1 0 Ze , if Le Ze dictates that ∂ ∂ =− m s K u u ρ . Deviations from this linear behavior arise in the preheating zone due to temperature changes which influence the density. This finally leads to a significant effect near the reaction layer. A similar effect is in principle present when using m m b b 0 in the burnt gases and extrapolating to s = 0. However, this effect is much smaller than in mu because the reaction layer is much thinner than the preheating zone. In case of an infinitely thin reaction layer, there is no error in m m b b 0 .

Comparison with numerical results
In this section we will validate the flame stretch theory by comparing the expressions for the mass burning rate with numerical simulations. First, methane-air flames are studied in a counterflow configuration. In order to study the behavior at strong stretch, the applied strain rate is increased up to the maximum level at which a steady flame exists. At higher strain rates the flame extinguishes. Second, we will investigate preferential diffusion effects in lean and rich propane-air flames. Different approximations in the theoretical strong stretch model will be assessed. Third, the mass burning rate of stretched lean methane-hydrogen-air flames is investigated to assess the theoretical model for fuels consisting of multiple components with different Lewis numbers.

Strongly stretched methane-air counterflow flames.
The objective of this section is to investigate the accuracy of the flame stretch theory for strongly stretched flames. This is done by comparing results obtained by the expression (52) with results from numerical simulations of back-to-back counterflow flames. Premixed lean methane-air (equivalence ratio ϕ = 0.7) counterflow flames are computed in the back-to-back configuration employing the GRI-Mech 3.0 reaction mechanism [50]. As a first step, we avoid preferential diffusion effects by setting the Lewis numbers to 1 for all species. For methane flames with an effective Lewis number close to 1, this is a valid assumption. Preferential diffusion effects are included in the analysis in the following sections. Steady flames are computed for a range of applied strain rates, a. Starting at a = − 100 1 s , the strain is increased until no steady solution can be obtained. The last steady flame is found at a = × − 1 65 10 3 1 . s . At higher strain rates the flame extinguishes.
The mass burning rate at the burnt side, mb, is plotted as a function of Karlovitz number Ka 0 in Fig. 6. For these steady flames, the  . As expected, this expression agrees with the numerical results for weak stretch ( Ka 0 0 1 < . ) but shows significant deviations at larger stretch rates. The present results demonstrate the validity of the theoretical expressions for mb derived by integral analysis of the set of SSFE. This basically confirms that the main underlying assumptions (i.e., a thin reaction layer and a one-dimensional manifold near chemical equilibrium) are valid. In Section 3.3.1 of Part II, we will investigate how accurate stretched flames are modeled by using FGM.

Lean and rich propane-air flames.
The aim of this section is to get more insight in the accuracy of the different steps, which are taken from the original linearized strong stretch expressions (54) and (55) to the final expressions like Eq. (57) focusing on preferential diffusion effects. This is done by comparing analytical results with numerical results for stretched lean (ϕ = 0.6) and rich (ϕ = 1.4) propane-air flames computed using the San Diego Mechanism [51]. The numerical mass burning rate mb is shown as a function of the Karlovitz integral Ka T 0 in Fig. 7 with black markers. To understand the impact of different assumptions on the Markstein number expressions, several assumptions in the derivation of Eq. (57) are relaxed. To start with, we consider the full expressions for the mass burning rate (54) and (55)  The enthalpy and element mass fraction changes can be written in terms of changes in temperature and mass fractions of the major species H2O, CO2 and O2 as shown by Groot et al. [38]. Neglecting the species' terms and taking only variations in Tb into account gives: For relation (57) taking ΔTb from the numerical solution of the fully stretched flame structure gives the green curves in the figure. This approximation is still accurate for both the lean and rich cases. Next, assuming one-step chemistry, the temperature change ΔTb can be expressed as the difference between the Karlovitz integrals of temperature and the deficient (lean) species: where 'lean' refers to the Karlovitz integral of the lean species (either fuel or oxidizer). This relation is plotted with red curves. It works well for the lean case but not for the rich, because the burnt gas contains large fractions of CO and H2. Finally, assuming analytical (exponential) profiles for ρ 0 (s), T 0 (s) and Y s i 0 ( ) gives Eq. (65). This expression is shown with blue lines and is a very poor representation of the true result. Note that this means that a much more accurate evaluation can be obtained substituting the one-dimensional numerical solution of the undistorted flame structure for ρ(s), T(s) and Yi(s) instead of approximate (exponential) solutions, like Eqs. (58) and (62).

Lean methane-hydrogen-air flames.
In the previous part of this section we have shown that the linearized strong stretch theory is able to reproduce numerically computed mass burning rates of propane-air flames. The full theory therefore seems to be accurate, but it would be desirable if more simplified versions of the theory could be used, which would make interpretations easier. It was found that rich propane-air flames are governed by multiple Lewis numbers, while lean flames are dominated by the Lewis number of the fuel. On the other hand, simple expressions found in the literature (e.g., Ref. [40]) can only predict the qualitative behavior of these flame parameters. Furthermore, we noticed that if we restrict preferential diffusion effects to temperature variations ΔTb in stretched propane-air flames, the theory is still reasonably accurate. In the following part of this section, we will study different parts of the theory in more detail and for other flames. It is investigated whether temperature changes are dominant also in other flames. It is shown that this is generally not the case, but that the combination of temperature changes and equivalence ratio changes yields much more accurate results. The contribution of different element and enthalpy changes to the change in mass burning rate will also be investigated. We will study methane-hydrogenair flames for this validation.
In previous sections we considered flame quantities at the burnt side of the flame (where fuel consumption has decreased to 5% of its maximum value). For a proper comparison with numerical results for methane-hydrogen-air flames computed with GRI-Mech 3.0 [50], it is more appropriate to use the inner layer inside the flame as point to evaluate the numerical and theoretical results. Equation (57) predicts the mass burning rate at the burnt side of a stretched flame. By integrating the continuity equation (29a) between these two locations, the mass burning rate mi at the inner layer (indicated by the subscript i) can then be determined. The inner layer is defined here as the point where the chemical source term of the progress variable has its maximum value. As shown in Ref. [52] with Δψ ψ ψ , , , … e ) at the inner layer and a modified Karlovitz integral, given by: where H is the Heaviside function. Predictions of mi from this theoretical expression are compared with results of a numerical simulation, in which the full transport equations are solved. The effects of further simplifications of the theory, which are used in most existing theoretical studies, are investigated as well. One of these simplifications that are mostly used is to assume that the unstretched mass burning rate does not depend on all four variables in ψi but only on flame temperature Ti through the Zeldovich number, as in the previous section. To study the effect of this, Δψi is projected on the one-dimensional subspace spanned by e T T = ∂ ∂ ψ i i . This vector is determined numerically from unstretched flame simulations for a range of Tu-values. Another simplified version of Eq. (54) is similar to the previous one, except it includes only variations in the equivalence ratio ϕ. In this case, the variations Δψi are projected on the one-dimensional subspace spanned by e φ ψ φ = ∂ ∂ i . This vector is determined numerically from unstretched flame simulations for a range of ϕ-values. The last simplified version is a combination of the previous two, i.e., it includes variations in temperature and equivalence ratio independently. In this case, the variations Δψi are projected on the two-dimensional subspace spanned by eT and eϕ.
In the results which follow, three different equivalence ratios ϕ = 0.6, 0.8, and 1.0 are considered. The equivalence ratio takes into account the overall stoichiometry of the methane-hydrogen mixture. For each stoichiometry, the mole fraction of hydrogen in the fuel is varied from pure methane, X H2 0 0 = . , to pure hydrogen, X H2 1 0 = . .  Table 1.
Methane-air flames. Before the influence of hydrogen addition is investigated, the pure methane case is considered. In Fig. 8, the mass burning rate mi is shown against the Karlovitz integral Ka i for different equivalence ratios. For all cases the mass burning rate decreases with positive stretch. The numerical results for mi (symbols) are compared with different theoretical expressions (lines). The theoretical expression (74) using all 4 components in Δψ predicts the mass burning rate within 0.5% for all cases (Fig. 8a). Simplified versions of this expression taking only changes in T or ϕ into account are clearly not able to predict the right behavior (Fig. 8b,c). For most cases large deviations are observed. On the other hand, the theoretical model taking both T and ϕ changes into account shows the same trend as the numerical results (Fig. 8d). Although the results are somewhat worse than that of the complete expression, the difference in mi with the numerical results is approximately 2.5%. Note that the influence of equivalence ratio fluctuations was smaller for propane-air flames, as shown in Section 2.3.2.2.
The explanation for this behavior can be found in Fig. 9. 0 of a lean flame is more sensitive to changes in the fuel concentration than in the amount of oxidizer. It can be seen in Fig. 9 that the changes in ZC and enthalpy have the largest contributions to Eq. (74), while ZO and ZH have negligible effects on mi for these methane-air flames. Furthermore, the effects increase with decreasing equivalence ratio. To a large extent, this is caused by the higher sensitivity at lower ϕ (cf. Table 1). The mass burning rate of lean flames is more sensitive to changes in ϕ than that of stoichiometric ones. Although the separate contributions are comparable to the direct effect of stretch (i.e., the second term in the r.h.s. of Eq. (74)), the sum of the contributions is much smaller because the changes in ZC and h have opposite effects. This is an effect noticed for the special case of methane-air flames, but this is not the case for ethane-air and propane-air flames (see also de Swart et al. [52]).
When the changes are projected on the space spanned by eT (see Fig. 10a), it is obvious that mainly enthalpy changes remain. Small non-zero values for ZH can be observed, since the numerical approximation of ∂ ∂ Z T H i u , 0 is small but not exactly 0. With this simplified model the predicted mi is too low. When Δψ is projected on the space spanned by eϕ (see Fig. 10b), mainly contributions of ZC and ZH are observed. Since the ratio of C and H is fixed in the fuel, they change proportionally. As a result, the effect of ZH is overestimated. The oxygen mass fraction ZO also changes, but because of the very low sensitivity it has a negligible contribution. The mass burning rate predicted by this version of the model is too high. When both effects are combined (see Fig. 10c), the most essential effects are covered. However, the contribution of ZC in the full version of the model is in this version replaced by the combined effect of ZC and ZH. In this simplified version, the ratio between ZC and ZH is fixed because changes in ϕ and T do not alter this ratio. Hydrogen addition. The influence of hydrogen addition is investigated next. It is well known that the burning velocity of unstretched flames increases when hydrogen is added at constant ϕ. Here, it is investigated how the mass burning rate of methanehydrogen flames changes due to flame stretch. Furthermore, this behavior is explained by analyzing the role of preferential diffusion.
In Fig. 11 the mass burning rate of stretched methane-hydrogen flames is shown for fuels of various equivalence ratio and a hydrogen content up to 60%. It can be seen that the scaled mass burning rate m m i i 0 of positively stretched flames increases when hydrogen is blended with methane. This effect is stronger for leaner flames. In Fig. 11, values of mi predicted by Eq. (74) are shown as well (lines).
For increasing hydrogen content, the difference with numerical values increases from 0.5% at X H2 0 0 = . to 2% at X H2 0 6 = . . The behavior of m m i i 0 observed in Fig. 11 is caused by changes in the enthalpy and element mass fractions Δψ. The separate contributions of these changes to mi are shown in Fig. 12 for various hydrogen contents. It can be seen that the contribution of ZC decreases with increasing X H2 . On the other hand, the contribution of ZH increases significantly when hydrogen is added to the fuel mixture. This effect increases for increasing hydrogen content up to X H2 0 6 = . . Markstein numbers. The Markstein numbers are determined from the previous results by fitting a third order polynomial to the  < . ). At higher X H2 the theoretical models follow the trend of the numerical results, but the quantitative agreement is less. The main reason is that for a hydrogen flame ( X H2 1 = ) at low temperatures the heat release ωT as a function of the scaled temperature T is non-zero in a broad region. This means that the assumption that the reaction layer is thin compared to the   preheating zone does not hold any more when reaching the pure hydrogen limit. In fact, the reaction layer shifts from the burnt side in the pure methane case to the leading edge of the flame for pure hydrogen. Neither the current theoretical model nor other existing models assuming a thin reaction layer take this behavior into account. The Markstein numbers presented in Fig. 13 correspond to the range of values given in literature. The behavior of M i as a function of X H2 was also observed in Ref. [53]. A quantitative comparison is, however, not very useful because different configurations and definitions are used in literature. The actual value of M is, e.g., sensitive to the position in the flame at which m is evaluated. However, the trends as function of ϕ and X H2 presented here do not change and the conclusions remain the same.
To conclude, we can say that changes in enthalpy and element mass fractions give rise to significant changes in mass burning rate. Predictions of mi using the full theory correspond well with numerical values. Simplified versions of the theory taking only changes in equivalence ratio or temperature into account are not sufficient to predict the numerical behavior of mi versus Ka i . Therefore, one should be careful applying theories or reduced numerical models that are based on one of these changes only. When variations in both temperature and equivalence ratio are included, the results are acceptable for the present case, though the large change in C to H ratio at the inner layer is not accounted for. This conclusion is not surprising since the burning velocity is to leading order dependent on temperature and equivalence ratio. This has been observed by others as well. However, unlike in previous studies, both separate effects are studied quantitatively here. Moreover, the theoretical model used in this study is able to predict these effects resulting in Markstein numbers that correspond both qualitatively and quantitatively with the numerical results for not too high amounts of hydrogen addition.
From the results presented here, one may also conclude that when highly diffusive hydrogen ( Le H2 0 3 = . ) is blended with methane ( Le CH4 0 97 = . ), the increase in ZH due to stretch becomes substantial. This effect can be so large that the Markstein number changes sign and becomes negative. This means that blending hydrogen to a methane-air flame reduces its diffusive-thermal stability possibly leading to sharp cusps, cellular structures and an increased flame surface area in turbulent flames as seen by many others (e.g., Refs. [54][55][56]). Thermal-diffusive instabilities might lead to increased flame stretch rates due to strong curvature effects. Although the present study is limited to weak stretch and considers a linearized version of the theoretical model, the original theory includes non-linear effects at high stretch rates.

Conclusions on Part I
A generalized formalism has been introduced to split the full set of conservation equations for describing the detailed behavior of flames into a kinematic equation for the motion and dynamics of the flame front and a set of SSFE equations for the internal structure of the flame structure. The full dynamics of the flame can be very well approximated if the SSFE model is further reduced into a model in which the flame structure is split into a preheating zone, where detailed transport phenomena including stretch, curvature and preferential diffusion take place, and a much thinner reaction layer in which only chemical reaction and molecular diffusion take place and balance each other. This is shown in terms of an advanced flame stretch theory for the mass burning rate of flames, which reduces to existing weak stretch theories under certain conditions. The theoretical formalism has shown to be an excellent prediction of the dynamics of numerically modeled flames using detailed transport and chemistry, even in the case of (strong) flame stretch, curvature, elemental/enthalpy effects. This opens the way to formulate flamelet-based reduction techniques in which the chemical source term of the flame is accurately predicted. In such techniques, it is important to include major effects related to molecular diffusion in the chemical source term evaluation, but the inclusion of stretch, and curvature is less important inside the reaction layer, while these effects have to be resolved by the model in the thicker preheating zone. The reproduction of the enthalpy and the elemental composition at the reaction layer is, however, important and has to be included. The FGM model, developed and described in Part II, is based on these ideas.

Introduction to Part II
In Part I of this paper, a mathematical formalism has been introduced to describe premixed combustion systems in terms of a kinematic equation for the flame propagation in combination with a flamelet system of equations (SSFE) describing the internal flame physics and chemistry. This method forms the basis of the Flamelet Generated Manifold (FGM) method for combustion processes taking place in thin active layers. The basics of the FGM approach and its application to laminar and turbulent premixed flames will be presented in this Part II. The principle idea of the method consists of concepts from both reduced chemistry as well as from laminar flamelet models. In a sense the FGM approach forms the bridge between these two methods, which were thus far regarded as two separate methods. In a number of previous publications, parts of the FGM approach have been introduced and applied. The FGM method was applied to premixed laminar flame first, see e.g., Refs. [6,21,57]. The applicability of FGM to partially premixed and nonpremixed laminar flames was also tested with success in Refs. [7,12,13,[58][59][60]. The FGM method has an important advantage compared to chemistry based reduction methods because the main influence of transport effects has been included in the solution of the reduced chemical mechanism.
A derivation of the FGM method with the relevant approximations and assumptions will be presented in this Part II. In Section 3.2, we will first present and derive the complete FGM model and compare it with other existing methods (Section 3.2.1) and describe how to use the model in practice by considering representative flamelets (Section 3.2.2). The numerical implementation of FGM and the coupling with a flow solver are described in Section 3.2.3. The method is further validated and applied step-by-step to various laminar combustion cases in Section 3.3. Later it is extended to DNS, LES and RANS models for turbulent flames in Section 2.4. The conclusions of Part II are given in Section 3.5. Special attention is given to the modeling of flame stretch and preferential diffusion effects.

Complete FGM model: from flamelet equations to manifold
The basic principles of the FGM approach are similar to those of reduced chemistry models that are based on dimension reduction techniques, such as ILDM [2] and CSP [3]. The idea is that the combustion process is described by a number of slow processes, while most of the fast processes are in a 'generalized' steady state. Suppose that the time scales are ordered, such that y y y = + +1 in general, but it should be realized that Zj and h are often not independent, so that one does not need to treat each conserved quantity as an independent degree of freedom. For instance, the element mass fractions are not independent since ∑ = Z j 1. Furthermore, for fully premixed adiabatic flames (with Le i = 1), h and Zj are all constant and do not have to be treated as degree of freedom. For nonpremixed flames in an adiabatic two-fluid stream, h and Zj are dependent on a single degree of freedom, referred to as the mixture fraction Z, which leaves us with only one single degree of freedom for the conserved variables. So the manifold with the lowest dimensionality for fully premixed flames is a one-dimensional manifold with y1 as degree of freedom. If the accuracy is not good enough also y2 or y2 and y3 can be used. For partially-premixed and nonpremixed flames, one has to treat at least y1 and Z1 as degrees of freedom (see Refs. [59,61]) and if a higher accuracy is desired one could add y2 (or y2 and y3).
There are different ways to define or approximate the solution of the manifold describing yk ( k r N N = + − 1, , … e ) as a function of yi (i = 1,. . .,r), h and Zj. Let us first look at the full most accurate manifold definition and consider approximations afterwards. To derive the exact manifold equations, we transform the combustion equations to a flame-adapted coordinate system as derived in Section 2.1.3. This thus leads to the set of strongly stretched flamelet equations SSFE (Eq. 29). This set of SSFE defines the inner 'flamelet' structure and the local 'manifold' is governed by the magnitude of the different convective, diffusive and reactive time scales, expressed in the flame-adapted coordinate system. So for any reduced chemistry or flamelet-based model, the SSFE set forms the basis. The exact solution of a combustion CFD problem would be obtained, if the local values of the flame stretch, curvature and Qi-fields along flamelet paths in the CFD problem would be derived and the SSFE set would be solved with those fields, to determine the inner structure in terms of y k ( k r N N = + − 1, , … e ). Although the focus has been on premixed flames in the derivation of the flamelet equations, the approach can also be applied to non-premixed flames. The flame coordinate system is then attached to iso-surfaces of the mixture fraction. This is explained in Appendix B.
There are different ways to combine the manifold with the transport equations. In the G-equation flamelet model (see Ref. [16]), the kinematic equation is solved for the motion of an infinitely thin flame front, while the local mass burning rate m of this front is derived from some approximation of the SSFE set. A more complete model would be to solve the full kinematic equation, Eq. (13), and locally derive the mass burning rate m from the full SSFE set. This would mean that the original set of transport equations is solved exactly. This method was used by Groot et al. for modeling the motion of spherically expanding premixed laminar flames, but they found it to be not very stable in a numerical sense [38].
Another method is to solve a set of transport equations for yk, k = 1, . . ., r for the r slowest processes together with those of h and Zj, like in so-called reduced chemistry models, e.g., ILDM [2], CSP [3], FGM [6], FPI [9], FPV [62], PS-ILDM [7] and REDIM [8]. In that case one does not use the full solution of the SSFE set, but depending on the method, some approximation of this SSFE set is proposed. How this SSFE set is approximated depends on the method involved.
There are several ways to derive reduced chemistry manifolds from the SSFE set Eq. (29). In case of ILDM or CSP, all diffusive, stretch and curvature effects in Eq. (29) are neglected and the reduced chemistry manifold is deduced from a time scale analysis of the remaining set of equations which only contains chemistry. The fastest N − Ne − r chemical time scales are assumed in chemical steady-state which defines the reduced chemistry manifold. In case of FGM, the transport and stretch/curvature terms are not neglected in the computation of the manifold. In practice, the FGM manifold is approximated from a series of independent 1D computations of a fully premixed or non-premixed flame (see next section).
To conclude, there are multiple ways to model flames based on the SSFE set. FGM is a promising technique for two reasons: First, it combines the accuracy of dimension reduction methods by solving a set of complete transport equations for y1,. . ., yr together with h and Zj, and second it does not only consider chemical steady-state of the SSFE set, but also takes into account most important transport effects. This is very important since flame structures are dominated by diffusion-reaction balances with a major influence of diffusion effects. This advantage of FGM was demonstrated by van Oijen et al. [6] by comparing a 1D FGM with a 1D ILDM for a premixed hydrogen-air flame. While at high temperatures close to chemical equilibrium the manifolds are similar, the FGM is more accurate in the colder parts of the reaction layer. In practice, it is too expensive to solve the full SSFE set (29) and we therefore approximate it by a representative set of solutions instead. How this is done is explained in the next subsection.

Practical FGM model: 'representative' flamelets
In the previous section we explained that the set of SSFE equations forms the basis of the manifold of the FGM method. A local strongly stretched flamelet in a CFD problem may be computed from this set on the basis of local stretch fields K(s) and curvature fields σ(s) as input. Direct computation of this set for each flamelet in the combustion problem is too involved and we therefore restrict to best approximations of this set in practice. Convective, diffusive and chemical source terms are always involved, but some transport, curvature and stretch terms are neglected, depending on the situation. The Qi-terms are always neglected, as explained in Part I.
The following procedure is used to derive the manifold from an approximation of the SSFE set. First, the conserved quantities h and Zj are considered. For the combustion system of interest, the number of independent 'mixing' parameters is evaluated. For instance, in case of a fully premixed system with unity Lewis numbers, Zj is constant, but h might be changing. In case of heat losses, h is changing and this is the only conserved scalar taken into account. To generate a 2D manifold that includes enthalpy variations, a series of 1D representative flamelets is solved with step-wise varying initial temperature to mimic the changes in h in the system. The 2D manifold is now parameterized by two scalars h and y1, while transport equations for h and y1 are solved during run-time. This is quite accurate for problems with constant Zj (See Section 3.3.2). However, in case of strong preferential diffusion effects, Zj and h are changing as well, even in case of a fully premixed combustion system. For that case one should also take into account variations in Zj and h as governed by the preferential diffusion phenomena. It is found that for some cases changes in Zj and h due to preferential diffusion are coupled, so that one may again define a 2D manifold with one parameter for these coupled variations in Zj and h. This is studied in Section 3.3.1. In case of a non-premixed system, at least one mixing parameter Z is needed to predict variations in Zj and h, but this could also increase if the elements and enthalpy are not coupled, e.g., in flames with more than two different inlet streams.
In any case, if the number of independent variations of Zj and h is known and parameterized, the principle method is to solve a set of 1D flames which mimic the combustion system best, by varying the inlet and/or boundary conditions according to the variations in conserved quantities. In case of premixed combustion systems, the flame curvature and stretch could be of interest. If that is the case, a series of 1D flamelets is solved with step-wise varying stretch and/ or curvature fields. In case of non-premixed or partially-premixed flames, a series of stagnation flames is solved, where the stagnation flow is governed by the stretch field. Again, boundary conditions are used which mimic the variations in Zj and h, but now on both sides of the flame. Contrary to the classical flamelet method, the solution is not parameterized in terms of mixture fraction Z and scalar dissipation rate χ = ∇ 2 2 D Z but by Z and y1. If auto-ignition and quenching phenomena are present in the combustion system, like in the case of engine combustion, the solution may have to contain unsteady flamelets in order to accurately predict these phenomena [63,64].
Apart from the procedure explained above, there is no principle way to generate the representative flamelets. The user has to analyze which flamelets are most representative of the system at stake. This can be regarded as a disadvantage of the FGM method compared to other dimension reduction techniques.

Implementation of FGM
The implementation of the FGM method consists of three main steps:

Computing representative flamelets 2. Building a look-up table 3. Coupling the table with a flow solver
These three steps will be explained in the following.

3.2.3.1.
Computing representative flamelets. The first step consists of computing flamelets which are representative for the combustion system to which the FGM will be applied. It is obvious that when combustion takes place under (non-)premixed conditions in the application, the most representative flamelet is a (non-)premixed flamelet. The flamelets are computed by solving the flamelet equations (29) under the appropriate conditions. For a premixed flame this means that the boundary conditions correspond to an inflow of the reactants at one side (s → −∞) and an outflow at the other side. When we furthermore assume that σ = 1, K = 0, Qi = 0, and that the time derivatives (∂/∂τ) are zero, the equations describe the most elementary 1D premixed flamelet: an adiabatic, flat, stretchless, steady flame. CHEM1D [65] is a code that can solve the flamelet equations numerically employing detailed chemistry and transport models. The solution consists of the profiles of mass burning rate m, enthalpy h, and species mass fractions Yi as functions of the spatial coordinate s. Note that for this special case m is constant and an eigenvalue of the system. So the 1D flamelet solution yields the thermochemical variables φ parameterized by one variable: the spatial coordinate, φ = φ(s). When more degrees of freedom are required, e.g., to account for changes in Zj, h and p, multiple flamelets have to be computed.
For instance, heat loss to the walls of the combustion chamber decreases the enthalpy h in the computational domain. In order to take this into account in the FGM tabulation process, laminar flamelets have to be solved for different values of enthalpy, introducing enthalpy as an additional degree of freedom. The enthalpy in an adiabatic premixed flame is constant and equal to the enthalpy of the burnt mixture, hb, apart from small local changes due to non-unity Lewis number effects. The enthalpy of the flamelets can be reduced by lowering the temperature of the reactants, dilute them with cold combustion products, computing burner-stabilized flames, by adding a radiation source term, or by rescaling the heat release source term. These different approaches are explained in Refs. [6,21,57,66,67]. With enthalpy as an additional parameter, the thermochemical state is parameterized by two variables: ϕ ϕ = ( ) s h , b . In general, an Nm-dimensional FGM can be generated from a series of one-dimensional flamelets with Nm − 1 parameters π:

Storage and retrieval.
During the CFD run, the code needs to retrieve thermochemical variables from the table for given values of the control variables. Therefore, a transformation needs to be performed: ϕ π π ϕ s y y N N , , , , with αi the weight factor of the mass fraction of species i. This choice should result in a monotonic function Y x ( ) for all hb in order to have a non-singular mapping. Therefore, linear combinations of reactants or products are commonly used. In the example shown in Fig. 14, we have chosen a combination of carbon dioxide and oxygen, with α CO2 22 7 = . and α O2 31 3 = − . . This choice results in a monotonous progress variable for the case studied here. Note that this choice is not unique; many other choices that result in a monotonous progress variable are possible. Although theoretically the exact choice is not important, it can be optimized to reduce numerical interpolation errors of the retrieval procedure (see e.g., Refs. [68,69]). Apart from monotonicity and numerical resolution, flame coverage is also relevant.
As an example, the heat release rate is shown in Fig. 14 as a function of the original parameters s and hb and as a function of y 1 = Y and y2 = h after transformation. The heat release rate decreases with decreasing enthalpy. Below a certain enthalpy level, the flame temperature is too low to sustain a steady flame. Since steady flamelets can not be computed for lower Tb, the tabulated data end at this level. In the application of the manifold, lower enthalpy levels might occur, which requires an extrapolation of the data in Fig. 14. Often linear extrapolation is employed, but a more accurate extrapolation is described in Ref. [6].
A multitude of methods has been used to store the data and to retrieve values from it, e.g., artificial neural networks [70], in-situ adaptive tabulation [71], k-d trees [72] and orthogonal polynomials [73], each with its own merits. The most common and straightforward method, however, is to store a discrete representation of the function ϕ = ( ) f y y N 1 , , … m on a structured mesh and using multi-dimensional linear interpolation for retrieval. In this case, the thermochemical variables are interpolated to a curvilinear mesh in y y N 1 , , …   with y i * the entry of control variable i. This may seem like a complicated and computationally expensive method. However, in most cases the control variables yi are smooth monotonic functions of ηi and only weakly dependent on the other ηj (j ≠ i). Therefore, the Jacobian matrix ∂yi/∂ηj is diagonally dominant and as a consequence the solution is typically obtained with only a few iterations when Newton's method (or similar method) is applied.
Often the solution method of Eq. (82) can be simplified even further. For instance, when unity Lewis numbers are applied, the enthalpy is constant in each flamelet. For a 2D manifold with y 1 = Y and y2 = h, this would mean that y2 is only a function of η2. Therefore, η2 is easily found and subsequently back substituted in Eq. (82a) to find η1. In this case, the Jacobian is a 2 × 2 upper triangular matrix and the solution can be obtained efficiently using back substitution. More generally in Nm dimensions, this holds when control variable yi is independent of ηj for all j < i. If in addition the grid points are distributed equally along the mesh lines, then the retrieval procedure becomes basically search free.
The main advantages of this method are its simplicity, fast retrieval, and the fact that one has full control over inter-and extrapolation errors. The large memory requirement is its mainand probably its only -drawback. A 3D table of 15 variables stored with double precision at 100 points in each dimension results in 15 100 8 12 10 3 7 × × = × bytes or approximately 120 MB. A 4D table would result in approximately 12 GB, which is more than the memory of a single CPU in many computer clusters. A higherorder interpolation procedure could reduce the required number of grid points, but it increases complexity and computational cost to retrieve a value from the table. A method to reduce the memory requirements by using a memory abstraction layer was introduced by Weise et al. [74].

Coupling with a flow solver.
After the manifold is stored in a lookup table, it can be linked to a standard CFD code. First, in the initialization phase, the database is loaded into memory. Then, the CFD code must solve transport equations for the control variables, together with the momentum and continuity equations. In manifold methods, the conservation equations for the control variables are derived by a projection of the full system onto the manifold (see e.g., Refs. [75,76]). The projection determines how processes that drive the composition of the manifold are forced back to the manifold by the fast chemical processes. In flamelet-based methods, this projection is usually not considered. It was shown by van Oijen [6,77] that the projection has only a small contribution in the FGM method and can be ignored. However, this omission causes the final result of an FGM calculation to depend on the choice of control variables.
Ignoring the projection, the transport equations for the control variables can be derived by taking the proper linear combinations of the original species equations (1d). For the progress variable, which can be rewritten as where ω Y is the progress variable source term given by The first term on the right hand side of Eq. (84) represents the fluxes due to preferential diffusion. In case of unity Lewis numbers, this term is zero. By applying the chain rule, this term can be rewritten in terms of the gradients of the control variables ∇yi. This can be simplified by assuming that the gradients of the control variables are not independent but are correlated as in the 1D flamelets: with ci a coefficient that is a function of the control variables. The transport equation then reads [6,77,78] ∂ ( ) The preferential diffusion coefficient D Y is stored in the FGM table. In a similar way, transport equations for the other control variables can be derived. For the enthalpy, it reads with D h the preferential diffusion coefficient [6]. The transport equations for the control variables are solved during run-time together with the momentum and continuity equations, while all other parameters (e.g., ρ, cp, λ, ω Y , D Y , D h , T) are retrieved from the FGM database. For low-Mach pressure based solvers, ρ and T can be taken directly from the manifold. For density based solvers, however, an energy equation should be solved that includes acoustic terms (contrary to Eq. 89). The temperature can then be computed from the energy equation and the pressure follows from the gas law. The implementation of FGM in fully compressible solvers is described in detail by de Swart et al. [78] and Vicquelin et al. [79].
Since the progress variable is a combination of species mass fractions, the boundary conditions for Y are straightforward to implement and follow from the definition in Eq. (79). The boundary conditions for h are somewhat more cumbersome because they are often not defined in terms of enthalpy itself but in terms of temperature. Therefore, an iterative procedure might be needed to obtain the enthalpy at a boundary. To find the enthalpy at a wall with constant temperature T wall , one has to solve T h T Y, ( )= wall (90) for h with Y given at the wall. Since T is a nearly linear function of h, this equation can be solved with only a few iterations. Alternatively, one can solve Eq. (90) in a pre-processing step for a given T wall and store the solution h wall as a function of Y . This has also been studied by Ketelheun et al. [80].
In the method described here, the species mass fractions are not required to solve the equations and they are only retrieved from the database during post-processing and visualization of the results. Several alternative methods to couple the manifold with a flow solver, in which all or major species are transported, were compared by Jha and Groth [81].

Concluding summary
Summarizing, we have explained how the strongly stretched flamelet equations can be used to derive reduced chemical models. FGM is a particular promising method, since it does not only consider chemical kinetics but also takes into account the main convective and diffusive effects in the SSFE. A complete FGM model is obtained when K and σ profiles along flamelet paths in the CFD problem are used to solve the set of SSFE. This would lead to an exact representation of the full set of equations apart from the Qi-terms, which were shown to be negligible. In practice, however, this is computationally too expensive and the set of SSFE is solved in a pre-processing step assuming representative conditions. These representative flamelets constitute a manifold which is parameterized by control variables and stored in a lookup table. To couple the FGM model with a flow solver, transport equations are solved for the control variables and the variables needed to solve these equations (source terms, diffusion coefficients, etc.) are retrieved from the lookup table. This coupling has been studied in many papers by TU/e, but also by many other groups, e.g., the EM2C group at Ecole Central Paris, CERFACS, CORIA and the University of Darmstadt. Results of some studies are summarized in the following parts. In the following section, the application of FGM in simulations of laminar flames is discussed.

Application of FGM in simulations of laminar flames
In the previous section the basics of FGM were explained. In this section, we will discuss the application of FGM in simulations of laminar flames. First, following the analysis of flame stretch effects in combination with preferential diffusion at the end of Part I, we will study how these effects can be modeled by using FGM. In Section 3.3.1, flame stretch effects will be systematically analyzed by considering different diffusion models. We will discuss how changes in enthalpy and element mass fraction caused by flame stretch can be accounted for by increasing the dimension of the manifold. In the subsequent two sections, we will study much larger changes in enthalpy and mixture fraction that are caused by the boundary conditions of the system. In Section 3.3.2, heat loss effects are investigated in simulations of burner-stabilized flames. It is explained how an FGM with a progress variable and enthalpy as control variables can be used to simulate these effects. Then, the application of FGM in partially premixed flames is discussed in Section 3.3.3. An FGM with Y and Z as control variables is used to simulate steady propagating triple flames in a stratified methane-air mixture and the results are compared with detailed chemistry simulations. Finally, in Section 3.3.4 we study the modeling of formation of nitrogen monoxide by using FGM. Since NO formation is a relatively slow chemical process, extra care is required for the use of FGM.

Modeling flame stretch effects with FGM
Flame stretch and curvature effects play an important role in the dynamics of flames because they can have a large influence on the mass burning rate as explained in Part I. Since the most elementary 1D FGM is constructed from a single stretchless flat flamelet, these effects are not present in the 1D manifold. However, in this section we will investigate how stretch effects -possibly in combination with preferential diffusion -are taken into account by using the FGM method. Moreover, we will discuss how an FGM can be extended with additional degrees of freedom to obtain a higher accuracy in stretched flame simulations.
In order to investigate the role of flame stretch, we have modeled stretched lean propane-air flames (ϕ = 0.8) with both detailed chemistry [51] and FGM. In Ref. [21] a similar study was performed for stretched methane-air back-to-back stagnation flames. Here propagating lean propane-air flames are studied, which are steady (in the flame reference system) and flat, but they are subjected to a constant stretch field K(s) = a, with a the strain rate, which is varied from 0 to 1000 1 s − . To analyze the effect of preferential diffusion, three different diffusion models will be considered: 1. Unity Lewis numbers are assumed for all species ( Le i = 1). In this case the element mass fractions and the enthalpy are constant: ΔZj = 0 and Δh = 0. Additional dimensions to account for these changes are therefore not needed. Note that the Le i = 1 assumption is used in many theoretical studies and turbulent combustion models. 2. All Lewis numbers equal 1.5 ( Le i = 1 5 . ). This is an artificial case, in which differential diffusion effects occur, since heat diffusion occurs at a different rate than mass diffusion. All species, however, have the same diffusion coefficient. Therefore, the element mass fractions are constant, but the enthalpy is not. To account for variations in h, the manifold is extended with an additional dimension. We will discuss how this can be done and how it affects the results. 3. All Lewis numbers are constant but not the same ( Le constant i = ). This is the most realistic model considered here. In this case, both the element mass fractions and the enthalpy will vary. In principle, they change independently from each other and an additional dimension is needed for each element mass fraction and the enthalpy.
In the following subsections we will discuss these cases subsequently.

Results for
Le i = 1. When unity Lewis numbers are applied, the element mass fractions and enthalpy are constant and do not have to be used as additional degrees of freedom. Therefore a 1D FGM is constructed with the reaction progress variable defined as in Eq. (80). In Fig. 16 the computed mass burning rate mb of these stretched flames is shown as a function of the dimensionless strain It can be seen that the mass burning rate of the flames computed with the FGM is almost the same as those computed with detailed chemistry. This is expected for small stretch rates Ka b < 0 1 . , but for higher Ka b the stretch terms in the flamelet equations are not negligible and larger deviations might have been expected.
The reason for the good agreement between FGM and detailed chemistry results can be found in Fig. 17, in which the chemical source term of the progress variable is plotted as a function of the progress variable for all stretch rates. The flamelet with a = 0 is shown with a red curve, all the other curves are blue. (For interpretation of the references to color in this text, the reader is referred to the web version of this article.) It can be observed that the chemical source term is hardly changed by the applied stretch field. Therefore, the source term in the FGM, which was generated with a = 0, can also be used for cases with moderate stretch. At very high stretch rates ( Ka b 1), which may occur in turbulent flames, the reaction layer might be affected by flame stretch and additional degrees of freedom are required for the FGM. The observation that the source term occurs in a thin layer and is hardly changed by stretch when there are no preferential diffusion effects was the very basis of the strong stretch theory and the derivation of Eqs. (41) and (51) in Part I. The small changes in source term can be accounted for by increasing the dimension of the manifold. This was shown in Ref. [21] where a 2D FGM was generated by changing the inlet composition of the flamelets while keeping h and Zj constant.

Results for
Le i = 1 5 . . In our second diffusion model, all Lewis numbers are assumed to be Le i = 1 5 . . In this case the enthalpy changes, but the element mass fractions are constant. This makes it a suitable case to investigate how variations in a conserved variable can be accounted for by the FGM method.
In Fig. 18  . the enthalpy has been decreased by 180 J/ g, which corresponds to a decrease in flame temperature Tb of approximately 100 K. This lower flame temperature has a direct effect on the reaction rates in the flame. This can be observed in Fig. 19, in which the source term of the progress variable is plotted as a function of the progress variable for flames with different stretch rates. Contrary to the unity Lewis number case shown in Fig. 17, the source term changes significantly. The lower temperature at high strain rates leads to a reduced source term. At the highest strain rate of a = − 1000 1 s , the maximum of the source term is reduced by 25% of the value at a = 0.
To account for this change in source term, the FGM is extended with enthalpy as an additional variable. To do so, a series of flamelets with different enthalpies has to be computed. Here we compare two ways to lower the enthalpy: The first manifold (FGM A) has been created by changing the temperature Tu of the reactants. Since it would require unrealistically low Tu to lower the enthalpy suffi-ciently, an alternative method is needed. This is done by computing 1D burner-stabilized flames [57]. By lowering the inlet gas velocity, the flame temperature and enthalpy are decreased. The second manifold (FGM B) was constructed from flamelets of which the reactants were diluted by cooled combustion products [6]. The reactants were diluted with carbon dioxide, water and nitrogen in such a way that the enthalpy is decreased but the element mass fractions are unchanged. A third manifold (FGM C) was created in the same way as FGM A, but the flamelets were computed using unity Lewis numbers. Therefore, each flamelet has a constant enthalpy and shows no preferential diffusion effects. However, when this manifold is used, the preferential diffusion terms in the transport equations (88) and (89) for the control variables are retained with Le i = 1 5 . . This third approach seems inconsistent but it is included in the comparison because the flamelet equations are often solved assuming unity Lewis numbers, since that simplifies the storage and retrieval procedure (see Section 3.2.3.2). Furthermore, it will show whether it is more important to include enthalpy as additional dimension or to include differential diffusion in the  = a m ρ δ 0 for the different chemistry models. The results are very similar, except for the 1D FGM where Δhb = 0. Due to the lower flame temperature and hence the lower chemical source term, the mass burning rate decreases more than for the unity Lewis number case. Since the 1D FGM does not include variations in enthalpy, it cannot predict this preferential diffusion effect. All 2D FGMs, however, include these effects and follow the trend of the detailed result much better. The small difference between FGM A and B indicates that the result is not very sensitive to the way the enthalpy is changed in generation of the manifold. It is also interesting to see that the manifold based on unity Lewis number flamelets predicts the preferential diffusion effects accurately. This demonstrates that it is more important to allow for changes in hb by retaining the preferential diffusion terms in the transport equations for the control variables, Eqs. (88) and (89).
FGM B and C overestimate the effect of preferential diffusion on the mass burning rate. This small deviation can be traced back to the larger decrease in enthalpy in the simulations with FGM B and C as can be seen in Fig. 18. While the enthalpy in the 1D FGM cannot change, the decrease in enthalpy is predicted by all 2D FGMs but slightly more accurate by FGM A.

Results for Le constant
The last and most realistic diffusion model considered here assumes constant but not equal Lewis numbers. As explained in Part I, the element mass fractions Zj and the enthalpy h will change due to preferential diffusion in combination with flame stretch. In principle, to account for all these variations we need Ne additional degrees of freedom for Ne − 1 independent element mass fractions and enthalpy. However, for weak stretch it was shown that these changes are linear functions of Ka (see Eq. 55). As a result, they are proportional to each other and can be parameterized by a single variable. In this case, one additional dimension for the manifold is sufficient to account for preferential diffusion effects. Therefore Δ with c the varying parameter and ΔZj and Δh obtained from Eq. (55) as explained in Ref. [21]. The enthalpy and element mass fractions of the unburnt mixture are changed by varying the mass fractions of the major species and temperature. The flamelet data are stored as a function of the progress variable and a second control variable, which is chosen to be ZC + ZH because it results in a unique mapping and these elements are changing the most. A second manifold (FGM B) was generated by solving the flamelet equations including a constant stretch term K = a for a series of stretch rates a. The applied flame stretch results in changes in enthalpy and element mass fractions. These changes are not linearly dependent but can be parameterized by a single control variable. The same control variables as for FGM A are used.
Both 2D FGMs are used to compute stretched 1D flames and the results are compared with simulations employing the full kinetic mechanism and a 1D FGM. In Fig. 21  Ka 0 1 = − . On the contrary, both 2D FGMs yield good agreement with the detailed calculations. With increasing Ka b the agreement becomes worse for FGM A, because the linear relation between the changes in elements and enthalpy is only valid for weak stretch. FGM B accurately predicts the mass burning rate because it includes the correct correlation between the changes in elements and enthalpy since it is based on flamelets with the same constant stretch field. However, in general, the stretch field is not constant and not known in the pre-processing step when the manifold is generated. Therefore, to get a higher accuracy at high stretch rates, probably more independent dimensions are required.
It is important to emphasize that in all three test cases with various transport models, good results are obtained by using FGM tables based on stretchless flamelets. This demonstrates that the direct effect of stretch on the chemical source terms is negligible, but that the indirect effect through changes in Z j,b and hb due to preferential diffusion is important.

Flame stabilization by heat loss
Changes in enthalpy not only happen due to preferential diffusion effects, but they can also occur due to heat transfer by radiation or convection. Premixed laminar flames are often stabilized on a burner by heat loss. The heat loss to the burner lowers the temperature of the flame and as a result its burning velocity is reduced. To account for these heat loss effects, the enthalpy is added as an  additional control variable to the FGM method. This is done by solving the flamelet equations for different values of the enthalpy as was discussed in Section 3.2.3. The FGM method was used in Ref. [6] to model a premixed laminar flame that stabilizes on a 2D slot burner. A methane-air flame with an equivalence ratio of ϕ = 0.9 was simulated by using the reaction mechanism of Smooke [82] and a 2D FGM while constant Lewis numbers were adopted. Here, the enthalpy of the flamelets in the manifold was changed by lowering the temperature of the reactants and by diluting the reactants with cold reaction products. The mass fraction of O2 was chosen as progress variable. Further details of the numerical implementation can be found in Ref. [6]. Fig. 22 shows isocontours of T and the mass fractions of O2, CO and H on a portion of the computational domain for both the detailed and reduced chemistry computations. The results obtained with the 2D manifold are in excellent agreement with the detailed chemistry computations. The location of the flame front is nearly the same in both results, which means that the FGM approach accurately describes the burning velocity and how it is affected by heat loss and flame stretch. Not only the position of the flame front is in good agreement, but the values of the mass fractions are as well. Also in the flame tip, where stretch and curvature are very important, the reduced computations appear to coincide with the detailed chemistry calculations. Note that preferential diffusion effects as discussed in the previous section are not very large for this methaneair flame with an effective Lewis number close to unity.
Near the burner wall, however, some differences between the FGM and detailed results can be observed. The CO mass fraction at the wall is much larger in the detailed results than in the FGM results, because in the detailed model CO can diffuse along flame surfaces toward the wall, while in FGM it can not. In Fig. 23 the mass frac-tion of CO is plotted as a function of the dimensionless tangential coordinate ξ 2 /δf in the flame surface Y = = Y O2 0 07 . , with ξ 2 = 0 at the wall. In the FGM approach, Y CO is a function of the control variables Y and h. In the manifold the CO concentration rapidly decays with decreasing enthalpy. When the flame temperature is below 1400 K, the CO mass fraction is less than 10% of its value at  adiabatic conditions. Since the temperature of the wall is 300 K, the CO mass fraction, which is retrieved from the manifold, is nearly zero for ξ 2 /δf < 1. In the detailed case, CO diffuses from the adiabatic zone along the ξ 2 -coordinate toward the burner wall. Diffusion in this direction is represented by the Q-terms (31) in the flamelet equations. At the wall, Q YCO kg m s = − − 1 1 3 1 .
, which is two orders of magnitude smaller than the chemical source term ω CO at adiabatic conditions. However, at the wall the convective terms and the chemical source term are negligible and the Q YCO -term is balanced by diffusion in flame normal direction, i.e., parallel to the burner surface. Therefore, the Q Y CO -term is actually of leading order in the flamelet equations and cannot be ignored in this region. The magnitude of the Q term can also be estimated from the FGM result. When λ/cp is assumed constant and Y Y h in which χ t,h is introduced as the scalar dissipation rate of enthalpy in the flame surface. With this expression it is possible to estimate the Q-term in an FGM simulation and to have an estimate of the validity of the 1D flamelet assumption. Another difference between the detailed and reduced chemistry computations is the computation time. The computation with detailed chemistry took approximately two weeks to converge, while the FGM results were obtained within a few hours. These figures are only a rough indication of the gain in computation time, because different numerical solvers, grids and initial fields have been used. A more systematic analysis of the computational gain of FGM is given in Ref. [77]. The reduced number of equations and the reduced stiffness of the set of equations are the main factors that cause the reduction in computing time.

Stratified or partially premixed flames
In many engineering applications of combustion, fuel and oxidizer are not perfectly mixed before entering the combustion chamber and combustion occurs under partially premixed conditions. Since fuel and oxidizer are not perfectly mixed in partially premixed flames, variations in the element mass fractions Zj and enthalpy h are inevitable. The element mass fraction of carbon in a methane/air mixture for instance varies between ZC = 0.75 in pure methane and ZC = 0 in pure air. In non-premixed flames, the local equivalence ratio, which is determined by the local element composition, is usually described by the mixture fraction Z, which is defined in such a way that it is a conserved scalar unchanged by chemical reaction. For hydrocarbon mixtures it can be expressed in terms of the element mass fractions ZC, ZH and ZO as follows [83]: where the subscripts fu and ox denote pure fuel and oxidizer quantities, respectively. The mixture fraction is scaled such that Z = 1 in the fuel stream and Z = 0 in the oxidizer stream. To account for variations in the mixture fraction due to mixing, Z can be added to the FGM as additional control variable. This is similar to the addition of h as control variable to account for non-adiabatic effects as described in the previous Section 3.3.2. Recently, the same burner geometry as investigated in the previous section was studied with a stratification of the fuel in the inlet [84]. A parabolic profile was adopted for the mixture fraction going from Z = 0.033 (ϕ = 0.6) at the wall to Z = 0.0713 (ϕ = 1.35) at the center of the inlet. By using a 3D FGM ( Y, , h Z) good agreement was found with detailed chemistry results.
Here we want to isolate the effects of fuel stratification from heat loss effects. Therefore, we discuss the application of FGM to partially premixed flames in the context of triple flames as studied in Ref. [59]. A triple flame is a flame structure generated by flame propagation in a partially premixed system. A schematic representation of such a triple (or tribrachial) flame structure is shown in Fig. 24. In a partially premixed field the mixture fraction determines the local equivalence ratio and thereby the value of the mass burning rate. Since premixed flame speeds reach the maximum for near stoichiometric conditions, a flame in a partially premixed field propagates preferentially along surfaces of stoichiometric mixture, i.e., near Z Z = st . On the fuel-lean side of such a surface there is a lean premixed flame front and on the fuel-rich side there is a rich premixed flame front, both propagating with a lower burning velocity than the leading edge of the flame, called the triple point. Behind the partially premixed flame front, two streams, one containing unburnt intermediates like CO and H2, and the other unburnt oxidant, come together and burn as a diffusion flame.
In Ref. [59] a manifold was constructed from 1D premixed flamelets. Since the mixture fraction in such flamelets is conserved, changes in the mixture fraction as they occur in partially premixed flames have to be taken into account by adding Z as control variable to the manifold. The procedure to add Z as extra control variable is similar to the one described before for adding h as variable to account for enthalpy changes. Because the triple flames considered here are adiabatic, the enthalpy is not used as additional control variable. In order to add the mixture fraction as additional control variable, the 1D flamelet equations are solved for different values of the initial mixture fraction Z Z x −∞ = =−∞ ( ). The value of Z −∞ is simply changed by varying the ratio between fuel and air in the initial mixture. Since unity Lewis numbers are used, there are no differential diffusion effects and Z Z = −∞ in each flamelet. The progress variable is chosen as a linear combination of CO 2 , H2 and H2O, which is monotonously increasing in each flamelet. The Z range of the manifold 0 028 0 082 . , .
[ ] is within the flammability limits and is wide enough to cover the variations in mixture fraction in the triple flame discussed here. In applications with much larger variations, the manifold needs to be extended outside the flammability ranges. A straightforward linear interpolation of the species mass fractions and enthalpy between the lean and rich flammability limit and the pure oxidizer and fuel composition, respectively, is commonly used and gives reasonable results (see, e.g., Refs. [85,86]). Linear extrapolation of the density leads to incorrect results because density is a non-linear function of the control variables. Therefore, it is better to calculate the density from the linear extrapolated mass fractions and enthalpy. Details of the numerical set-up of these triple flame simulations can be found in Ref. [59]. Fig. 25 shows iso-contours of temperature T and the species mass fractions Y CH4 , Y O2 , Y CO , YH and Y CH O 2 computed with detailed chemistry [82] and FGM. The lean and rich premixed flame branches of the triple flame (see Fig. 24) can clearly be identified in Fig. 25. The temperature rises from 300 K in the unburnt mixture to approximately 2000 K behind the premixed flame branches depending on the local stoichiometry. The highest temperatures of approximately 2150 K are found in regions of near-stoichiometric mixture. However, since heat is conducted away from the stoichiometric line, these temperatures are almost 100 K lower than the maximum equilibrium temperature of 2240 K.
The steady propagation speed of this flame structure computed by using detailed chemistry is 44 0 1 . cm s − . When FGM is used, the propagation speed is 44 7 1 . cm s − ; a difference less than two percent. For more details, the reader is referred to Ref. [59]. This small difference indicates that FGM yields accurate results for these triple flames. This is confirmed by the temperature and species profiles shown in Fig. 25. For most variables the results of the detailed and reduced chemistry computations agree very well. Since transport processes along the iso-surfaces of Y are not included in FGM, the CO formed on the fuel-rich side cannot diffuse along the premixed flame toward the fuel-lean side. Therefore, the CO concentration in the fuel-rich premixed flame branch is slightly overpredicted by FGM. This effect is very similar to the diffusion of CO toward the burner wall as discussed in the previous section. In this case, however, the variation along the flame surface is caused by a gradient in mixture fraction and, hence, the Q-term is proportional to the scalar dissipation rate χt of Z in the flame surface.
The mass fractions of H and other radicals are hard to predict accurately using conventional reduction methods based on chemical steady states [2,75]. However, when using FGM the results agree well with the results of detailed chemistry computations, though small differences can be observed. Again, as in the case of CO, diffusion along the ξ 2 -coordinate leads to a smaller peak value of H mass fraction in the detailed computations. By using a similar expression as Eq. (91) but now with the scalar dissipation of mixture fraction, the magnitude of Q YH is estimated to be 10% of the flame normal diffusion term, which explains the small difference between FGM and detailed chemistry. The effect of the scalar dissipation on the accuracy of FGM in the modeling of partially premixed counterflow flames has been investigated in Refs. [12,61].

Prediction of NO formation
Formation of nitrogen oxide and soot formation are relatively slow chemical process, which cannot be assumed in quasi-steady state and therefore need special treatment. For flamelet models the problem can be solved by using unsteady flamelets which are computed at runtime (see, e.g., Refs. [87,88]). This results in accurate NO predictions, but at an increased computational cost. When using precomputed tables as in FGM, two approaches to model the slow formation of NO can be distinguished: 1. Include NO in the definition of reaction progress variable Y 2. Solve a transport equation for NO with the source term from the look-up table With mixed results both approaches have been applied in simulations of turbulent flames [86,89,90]. It is not easy to draw conclusions from these simulations because the predictions are strongly influenced by the modeling of turbulence-chemistry interaction and the accuracy of the NO measurements is not always high enough. Therefore, both approaches were compared in simulations of unsteady one-dimensional laminar flames [91]. This made it possible to validate the reduced models against simulations with detailed chemistry. In the following section it is explained how the FGM method is used to model NO formation. After that the onedimensional test problem is described and the results and conclusions are presented.

Parametrization of NO chemistry in FGM.
Flamelet generated manifolds are constructed from solutions of steady laminar flamelet simulations with detailed chemistry. In accordance with the application, premixed methane-air flamelets are used here. Since the equivalence ratio is not fixed in the test problem, a series of premixed flamelets is computed for mixtures of different stoichiometry as explained in Section 3.3.3. This set of flamelets forms a twodimensional manifold, which can be parameterized by mixture fraction Z and a progress variable Y . In this case, the reference progress variable Y 0 is defined as Since the formation of NO continues in the burnt gas, the progress variable has almost reached its maximum value, while the mass fraction of NO is still increasing. This results in extremely large gradients in the database ( ∂ ∂ →∞ Y NO Y ) and, hence, in large interpolation errors during data retrieval. This effect can be observed in Fig. 26. In this figure the NO mass fraction is shown versus progress variable at stoichiometric conditions ( Z Z = st ). It can be seen that the mass fraction of NO increases sharply when the progress variable has nearly reached its maximum value Y Y = max 0 . Obviously, this choice of parametrization would result in huge interpolation errors. One way to solve this problem is to add the mass fraction of NO to the definition of the progress variable, viz.
The resulting parametrization is also shown in Fig. 26 for different values of α. It can be seen that the gradient ∂ ∂ Y NO Y becomes smaller for increasing α. Due to the very small values of Y NO compared to the major species, a rather large value of α is needed to get a significant contribution of Y NO in the progress variable.
A similar observation can be made for the chemical source term ω NO of NO, which is shown in Fig. 27. Although the effect is less dramatic than for the NO mass fraction, also interpolation errors in the source term may occur. Since this source term is used to solve the transport equation for NO, it also results in errors in the NO predictions.

Test problem and results.
The two approaches for NO modeling as described above have been applied to a simple test problem [91]. It consists of a one-dimensional, adiabatic, premixed methaneair flame, which is subjected to sudden changes in stoichiometry. The initial solution consists of a steady lean premixed flame ( Z Z = 0 5 . st ). At t = 0 the equivalence ratio of the inlet mixture is changed to rich conditions ( Z Z = 1 7 . st ). This is done by prescribing the composition of the unburnt mixture at the inlet of the computational domain (s = −1 cm). The flame is stabilized in the domain by adapting the inlet gas velocity in such a way that T s = ( ) = 0 1500 K . At t = 0.25 s the unburnt mixture composition is switched back to the initial lean condition.
In Fig. 28 the profiles of mixture fraction Z, temperature T, and NO mass fraction Y NO are shown as a function of space s and time t. These results have been computed using detailed chemistry, i.e., the GRI-Mech 3.0 reaction mechanism [50]. In the plot for mixture fraction, it can be seen that the inlet composition suddenly changes at t = 0 and t = 0.25 s. While this sudden change in mixture fraction travels with the gas velocity toward the flame front, it is smoothed by molecular diffusion. When this mixture fraction wave travels through the flame front, its wavelength is increased by the increase in gas velocity due to expansion of the reaction products. In the temperature plot, it can be seen that when the mixture fraction wave hits the flame front the temperature rises, because the stoichiometry changes from lean to stoichiometric. When the wave has passed the reaction front, the temperature decreases again since the conditions have changed from stoichiometric to rich. When the mixture fraction is changed back to lean conditions at t = 0.25 s, a similar but reversed process is observed.
The profiles of NO mass fraction closely follow the trend in the temperature profiles. When the mixture fraction is close to stoichiometry and the temperature is high, a lot of NO is produced. This results in a growing NO peak which is convected away from the flame front by the flow of burnt gases. This can also be observed in Fig. 29, in which snapshots of the NO profile are shown at three time instants. In the top figures the mixture fraction has just increased, in the middle figures a steady rich condition is obtained, and in the bottom figures the mixture fraction has just changed back to lean conditions.
In Fig. 29 the results of FGM simulations with the two different approaches and different values of α are shown as well. The results with α = 10 indicate that adding a small amount of NO to the progress variable does not prevent large interpolation errors. In fact, for this value of α the profile of NO is still not resolved by the grid spacing used in the look-up table. The interpolation errors result in an overestimation of the mass fraction of NO, but an underestimation of its chemical source term. For larger values of α, the results become more accurate. The steady situation at t = 0.25 s is then very well captured by both approaches. However, the unsteady parts of the simulation are more accurately predicted by the approach in which a transport equation for NO is solved. The direct lookup ap-proach gives only acceptable results for α = 10 3 . Note that in the limit α → ∞, both approaches are the same because then the progress variable Y is equivalent to Y NO . In order to include the dependency of the NO source term on concentration of NO, it is possible to split the source term in production and a consumption part. This was investigated by Ihme and Pitsch [18].

Application of FGM in simulations of turbulent flames
The previous section showed how FGM can be successfully applied to laminar flames. These were all flames that are moderately perturbed by steady flow and mildly unsteady effects, in one or two dimensions. In this section, we will discuss the application of FGM in simulations of turbulent flames. In the numerical modeling of turbulent flames, a distinction is made between direct numerical simulation (DNS), in which all turbulent scales of motion are resolved on the numerical grid, and Reynolds-Averaged Navier-Stokes (RANS) or Large-Eddy Simulation (LES) methods, in which the governing equations are filtered or averaged before they are solved. DNS is extremely costly, requiring tens of millions cpuhours, but it provides important insight in the details and fundamentals of turbulent combustion. Recent reviews on DNS of turbulent combustion are given by Chen [92] and Trisjono and Pitsch [93]. The computational burden for RANS and LES is considerably smaller than for DNS. However, the former require additional modeling of the unresolved terms in the governing equations. Reviews on models for LES of turbulent combustion have been given by, e.g., Janicka and Sadiki [94], Pitsch [95] and Gicquel et al. [96]. The application of FGM in simulations of turbulent premixed flames will Since in DNS all scales of turbulent motion are resolved, it is the ideal tool to investigate flameturbulence interaction. In various DNS studies the interaction of a flame front with turbulence has been investigated in canonical problems in order to clarify the effect of flame stretch and curvature on burning rates, to assess the accuracy of reduced chemical models, and to develop models for closing the unresolved terms in RANS and LES [93].
As discussed in Part I, the set of flamelet equations gives an accurate description of the flame structure when Ka < 1 and perturbations from 1D behavior are not too large. In turbulent flames this implies that the model is valid as long as the conditions are within the flamelet regime [16], in which the aerodynamic scales are larger than the thickness of the laminar flame front. This ensures that the combination DNS-FGM performs without potential loss of accuracy, even if the conditions are more severe than the laminar ones that were discussed in Section 3.3.
First implementations of DNS-FGM were presented in Refs. [33,97,98], in which it was analyzed whether the strategy correctly performed under three-dimensional time dependent conditions with a simple fuel (Le = 1) in a state-of-the-art com-pressible DNS code. Similar studies using a low-Mach code were reported in Refs. [99,100]. In these studies, premixed flame kernels in homogeneous isotropic turbulence were simulated. By varying the intensity of the velocity fluctuations, flame kernels in the thinreaction zones regime (T1) and in the distributed reaction zones regime (T2) were investigated [33]. The turbulent Reynolds numbers for cases T1 and T2 were Re . = ′ = u l s t L f 0 3 7 δ and 37, respectively. In Fig. 30, cross sections through the center of the flame kernels are shown for t = lt/u′. Three isotherms are shown to indicate the location of the unburnt side, the inner layer and the burnt side of the flame front. From the figures it is clear that in case T1 the turbulent flow structures are small and strong enough to disturb the preheat zone of the flame. The reaction zone, however, is not significantly affected by the flow structures. This is the typical behavior expected in the thin-reaction zones regime [16]. For case T2, it was found that not only the structure of the preheat zone but also that of the reaction zone was altered by the turbulent flow, which is characteristic for combustion in the distributed-reaction zones regime. The high turbulence level in case T2 was found to locally quench the flame, ultimately leading to global extinction.
The gray curves in Fig. 30 represent flamelet paths as defined in Eq. (15). These flame paths are 3D curves that are projected on the center plane and cross this plane at the location of the inner layer.   stretch rates. By using the actual profiles of density ρ, progress variable Y , flame surface σ and stretch rate K as shown in Fig. 31 in the Karlovitz integral (75), the mass burning rates of approximately 10 4 flame paths in the DNS were predicted by Eq. (74) with a mean deviation of 2% for case T1. When constant K and/or σ profiles were applied in the evaluation of the Karlovitz integrals, the difference between theory and numerical result increased severely, which indicates that these approximations in theoretical expressions should be used with great care [98]. For case T2, the deviation was higher but still a strong correlation between theory (Eq. 74) and DNS results was observed, which shows the strength of this theoretical model, since the dimensionless stretch rates are of order O 10 2

( )!
Due to the large stretch rates observed in these flames, one might wonder whether the 1D FGM model applied in the DNS is valid for these conditions. The fact that the flame stretch theory accurately predicts the mass burning rate may be interpreted as a positive sign because the FGM method is based on the same flamelet equations with similar assumptions. However, in order to assess the accuracy of the FGM model a comparison with simulations employing the full kinetic mechanism is desired. The most direct way would be by performing a DNS of the same case with detailed chemistry. Since such a DNS would require enormous computational power, which is not available to everyone, other approaches have been used. In Ref. [98], the FGM model was validated by comparing it with detailed chemistry solutions of the flamelet equations for the flame paths in the DNS results with the actual flame stretch rate K and surface σ as shown in Fig. 31. Neglecting the unsteady terms and the Q terms, the flamelet equations were solved for O 10 4 ( ) flame paths with a 1D flame code using both FGM and detailed chemistry.
The results of the detailed chemistry calculations, employing GRI-Mech 3.0 [50] without nitrogen chemistry, are shown in Fig. 32 as scatter plots of H2, O and CO in composition space. In the top row of this figure, it can be seen that the scatter is located in a narrow region around the 1D FGM but that deviations of approximately 10% occur. In the bottom row, conditional scatter plots are shown for Y = 0 6 . . In these plots, cross sections of a 2D FGM and a 2D ILDM are shown as well. The 2D FGM was created by solving the flamelet equations with a constant stretch rate K = a. A series of flamelets was computed for a m δ ρ  a) were then used to construct an FGM φ(y1, y2) with two control variables. It can be seen that, contrary to the 2D ILDM, the 2D FGM is in excellent agreement with the scatter of the detailed calculations. Apparently, for the unity Lewis number case studied here, the flamelets computed with a constant stretch rate occupy the same states in composition space as the flamelets with largely varying stretch rates.
Since the 2D FGM is a better approximation of the space accessed by the detailed calculations, it is expected that it will lead to a more accurate representation of the chemical source terms, ultimately leading to a more accurate modeling of the mass burning rate. This was tested in Ref. [98] by using the 2D FGM to solve the flamelets from the DNS. In Fig. 33, the results of the 1D and 2D FGM are both compared with results of the simulations with detailed chemistry. The mass burning rate at the inner layer is compared for all the flame paths in two scatter plots, in which the FGM result is plotted on the vertical axis and the detailed chemistry result on the horizontal axis. The average deviation between the 1D FGM results and the detailed results is 5%. The largest deviations occur for the lowest mass burning rate, because these flames are subjected to the highest stretch rates. The 2D FGM results show a better agreement with the detailed results. The average deviation is only 1%.

Preferential diffusion effects.
Although 3D DNS with detailed chemistry is extremely expensive and only possible with petascale computing power [92], 2D simulations with detailed chemistry are tractable on affordable computer clusters. Despite the fact that 2D turbulence behaves very differently from 3D turbulence, 2D DNS studies can yield very useful information on how vortical flow fields interact with flame fronts and how the internal structure of the flame front in terms of species and temperature profiles is affected by flame stretch and curvature. As an example, we discuss the validation of FGM in a 2D DNS study of lean premixed methanehydrogen-air flame kernels [101]. The fuel mixture consists of 60% methane and 40% hydrogen by volume and is premixed with air at an equivalence ratio of ϕ = 0.7. Apart from the fresh gas mixture and the dimensionality of the problem, the simulation setup is the same as applied in the studies discussed so far in this section [33,98]. However, as it was shown in Section 2.3.2.3, the presence of hydrogen with a Lewis number of 0.3 in the fuel leads to strong preferential diffusion effects, which affect the mass burning rate and hence the dynamics of the flame front. To capture the variations in enthalpy and element mass fractions caused by preferential diffusion effects, the manifold needs to be extended with additional dimensions as explained in Section 3.3.1. There, it was shown that one additional degree of freedom is enough to describe this effect for small Karlovitz numbers. Here we will investigate whether this also holds for the present fuel mixture of methane and hydrogen.
The Both manifolds are parameterized by a set of two control variables. The first control variable is the progress variable, which is here . with ϕi = Yi/Mi. The second control variable is chosen as a combination of element mass frac- H . This set of control variables results in a monotonic parametrization of both manifolds. The chemical source term of the progress variable in FGM A is shown in Fig. 34. The source term increases for higher values of Z , which corresponds to an increased burning rate at richer conditions. The two flamelet curves demonstrate that a positive stretch rate leads to an increase in Z which is in agreement with the results in Fig. 12.  Both manifolds are used to compute strained flames similar to the study in Section 3.3.1. Results obtained with the 2D manifolds are compared with results of calculations using detailed chemistry (GRI-Mech 3.0 [50]) and a 1D FGM. In Fig. 35 the mass burning rate is plotted versus Karlovitz integral. The results computed using detailed chemistry show the strong effect of preferential diffusion. These effects were analyzed in detail in Section 2.3.2.3. Fig. 35 also shows the results of calculations that employ a complex transport model, which solves the Stefan-Maxwell equations and includes thermal diffusion [29]. The results indicate that the simplified transport model with constant Lewis numbers (Eq. 2) is also valid for this fuel with high hydrogen content. The 1D FGM cannot predict the preferential diffusion effects because it does not include changes in hb and Z j,b . Both 2D FGMs, however, predict the correct trend in mass burning rate. This good agreement and the small difference between the two 2D FGMs indicate that one degree of freedom is sufficient to describe the changes in enthalpy and element mass fractions, and that it is not very important how they are imposed.
The 2D FGM A is used in DNS of an expanding flame kernel and compared with the 1D FGM and the detailed mechanism [50]. Snapshots of the mass fraction of the hydrogen radical, which is an important intermediate in hydrocarbon oxidation and an indicator for the fuel consumption rate, are shown in Fig. 36. The detailed result shows the increased mass fraction of H in parts of the flame front that are curved outward. These parts of the flame are positively stretched and become richer due to preferential diffusion effects (see Section 2.3.2.3). This local increase in stoichiometry leads to an increase in mass burning rate, which attributes to the cellular instability of lean flames with Le < 1 [102]. While the 1D FGM cannot predict the change in Zj, the 2D FGM reproduces the detailed chemistry results very well. The increased H mass fractions in the outward curved parts are very well predicted. At closer inspection, however, one can still find differences. Due to the variations of species mass fractions along isocontours of the progress variable, diffusion in this direction starts to play a role. Delhaye [36] performed a quantitative analysis of all terms in the flamelet equations (29) for a similar DNS of a lean methane-air mixture as shown in Fig. 1. He found that the Q-terms are mostly much smaller than the stretch and curvature terms, but in the part with the highest curvature, where the curvature radius is comparable to the thickness of the reaction layer δr, the Q-terms become comparable in magnitude. Since Delhaye investigated a methane flame with a near unity Lewis number, the variations in Zj along the flame front are not as large as in the present methane-hydrogen mixture. Therefore, the Q-terms are expected to be larger for the present case. For an accurate representation of these tangential diffusion fluxes, additional degrees of freedom have to be added to the manifold.
DNS modeling of lean methane-hydrogen-air mixtures with FGM is extensively described by de Swart et al. [103]. Three different methane-hydrogen ratios were studied in 2D simulations of flamevortex interaction and 3D simulations of statistically planar turbulent flames in the thin reaction zones regime. In this work, a similar flame path analysis was performed as described at the beginning of this section. Stretch and curvature profiles were computed and used to calculate the Karlovitz integrals. From the scatter plots of the scaled mass burning rate m/m 0 versus Karlovitz integrals, Markstein numbers M were obtained using linear least-squares fits. The Markstein numbers were found to decrease with increasing hydrogen content, in line with the results presented in Fig. 13. However, the change in M appeared to be much smaller than in the steady laminar flames. The changes in element mass fractions ΔZ j,b and enthalpy Δhb in flame paths in the DNS were not as large as in steady laminar flames with the same Karlovitz number. This difference can be attributed to the unsteadiness of K and σ. When these variables change in time, the element mass fractions will follow with a finite response time. This is demonstrated in Fig. 37, in which the evolution of ΔZ j,b is shown for a flat flame that is subjected to a sudden change in stretch rate. The flame is initially not stretched and ΔZ j,b = 0, but at t = 0 the stretch rate is increased to K = − 10 3 1 s and the element mass fractions start to change. However, it takes about two flame time scales τf = δf/sL before the steady values are reached. This finite response time causes the flame to act as a low band-pass filter, effectively averaging out the effect of the stretch rate fluctuations. Note that expressions (43b) and (43a) for ΔZ j,b and Δhb do not include this transient effect because the unsteady terms ∂/∂τ were ignored in their derivation. In the FGM approach, however, transport equations are solved for Zj and h if they are control variables, which implies that this unsteady effect can be accounted for.
The local variations in Zj and h lead to a change in mass burning rate m, which on its turn affects the propagation of the flame front. This effect is not clearly visible in Fig. 36 because only one eddyturnover time was simulated and in this time span the movement of the flame front is mainly determined by the velocity of the gas. The effect of preferential diffusion on the dynamics of the flame was investigated by Vreman et al. [104] by using FGM in DNS. Premixed turbulent flames of lean methane-hydrogen-air mixtures were simulated on a slot burner geometry that was based on the experiments of Filatyev et al. [105] and that was used in earlier FGM-DNS studies of methane combustion [106,107]. In this study the preferential diffusion effects were demonstrated by comparing a 1D FGM without changes in Zj and h, with a 2D FGM that does account for these changes. In Fig. 38, instantaneous contours of the progress variable source term are shown for both cases. The FGM model with preferential diffusion effects (Fig. 38b) shows an increase in the source term in regions of the flame front that are convex with respect to the unburnt mixture and vice versa. Therefore, the flame front becomes more wrinkled than in case of the 1D FGM. This increase in flame surface density leads to an increase in the averaged burning rate and a shorter averaged flame length. In Ref. [104], the turbulent burning velocity was found to increase by 30% when preferential diffusion was included. The change in local equivalence ratio does not only affect the mass burning rate but it also changes the NO production rate, because the flame temperature changes. The NO mass fraction was found to increase by approximately a factor two at the burnt side of the convex flame front in Fig. 38b. Variations in equivalence ratio at the flame front are not only caused by preferential diffusion, but they also arise often due to imperfect premixing of fuel and oxidizer. The effect of such a fuel stratification on the dynamics of the flame was investigated in a similar numerical setup with DNS-FGM by Ramaekers et al. [108].

RANS and LES
In the previous section, it was shown that the use of FGM in DNS yields results close to results of detailed chemistry calculations, but at a fraction of the computational cost. This makes it possible to perform DNS of high Reynolds number flows in real burners (see e.g., Moureau et al. [109]). In practice, however, a DNS is still too costly for most applications and RANS or LES approaches have to be used. Typically, two types of models are required in order to close the first moment equations of turbulent reacting flows with flamelet based tables. First, a fluid mechanical model must describe the unresolved stress and flux terms. Usually, eddy-viscosity and gradient transport assumptions are employed to close these terms [26]. Second, a closure method is needed for the mean values of the highly non-linear chemical terms, such as the averaged source term, density, temperature and species mass fractions.
First implementations of FGM in RANS of premixed turbulent flames are described in Refs. [12,110]. In these studies, 3D manifolds with Y , Z and h as control variables are used and the chemical closure problem is tackled by describing the control variables in a stochastic way. It is assumed that locally the probability of occurrence of a certain state is described by a presumed shape probability density function (PDF). Usually, the joint PDF of the control variables is written as a product of the marginal PDFs assuming statistical independence of the variables. Beta functions parameterized by the first two Favre averaged moments (mean and variance) are used to describe the marginal PDFs. The combination of a tabulated premixed flamelet with a beta PDF for RANS modeling of turbulent combustion was already introduced by Bradley et al. [111] in 1988. Cook and Riley [112] proposed to use the beta PDF for the modeling of unresolved terms in the context of LES. Since then it has become a standard closure technique in both LES and RANS, although various other PDF shapes have been proposed and used successfully as well [113,114].
The computational cost of LES is much larger than for RANS, because LES is intrinsically unsteady and requires a finer grid to resolve a sufficient part of the turbulent length scales in the flow. However, since a larger part of the turbulent fluctuations is resolved and the modeled terms are much smaller, LES has a superior accuracy compared to RANS. Together with the continuous increase in available computer power, this has led to considerable attention for LES modeling in recent years. FGM or other tabulated chemistry models have been combined with many different approaches to close the unresolved or sub-grid scale terms in LES. An overview of various approaches to couple tabulated chemistry with LES is given in a review by Fiorina et al. [115]. A non-exhaustive list of recently used methods for premixed flames features: The levelset or G-equation formalism [116,117], presumed PDF methods [107,118], Flame-Surface Density methods [119,120], Thickened Flame models [80,121], and Filtered Flamelet methods [107,[122][123][124][125]. Several of these LES approaches have been used recently by various groups to simulate the turbulent stratified flame experiments performed at TU Darmstadt [126,127]. A comparison of these LES results is described in a joint paper [128]. Good predictions were found for all the models that were based on premixed flamelet tables that included mixture fraction and enthalpy as additional degrees of freedom. Since the premixed flamelet models performed well, it is expected that the diffusion fluxes tangential to the flame front (Qterms in Eq. (31)) are relatively small for this mild stratified flame, although a quantitative analysis was not given in Ref. [128]. Including heat loss by adding enthalpy as a control variable was necessary to quench the flame near the cold burner tip leading to a lifted flame as it was observed in the experiments.
Another setup, in which the modeling of heat loss has been investigated, is the confined premixed jet experiment performed at the German Aerospace Center (DLR). This geometry was selected and analyzed experimentally in Refs. [129,130] because of the pro-nounced recirculation zone that arises due to the off-center position of the jet-nozzle exit. The strong recirculation of the hot product leads to a large heat loss to the walls of the combustion chamber, effectively lowering the flame temperature. RANS simulations of this setup were presented by Domenico et al. [131] and Donini et al. [132]. In Ref. [132], the FGM method was applied and both adiabatic and non-adiabatic manifolds were used to model the flame. As to be expected, the FGM with enthalpy as additional control variable gave much more accurate results. Excellent agreement was found between the non-adiabatic FGM results and measured temperature and velocity profiles in the region close to the nozzle. At farther distance (more than 10 nozzle diameters) the agreement was found to be less satisfactory. The differences were principally attributed to the RANS modeling, which was not able to predict the extent of the recirculation zone in an accurate way. More recently, Proch and Kempf [66] performed LES of the same setup combining FGM with an artificial thickened flame model. They also investigated the inclusion of heat loss and basically came to the same conclusion. Furthermore, they compared two different ways to lower the enthalpy of the flamelets. It was found that this choice has almost no effect on the result, which is in line with findings of laminar flame studies (see, e.g., Ref. [57], and Section 3.3.1). Although the LES predictions were more accurate than the RANS ones, the flame length was still not in perfect agreement with the measurements. Proch and Kempf attributed the underprediction of the flame length to flame stretch effects.
The last test case that will be discussed here is a turbulent swirling flame in a laboratory-scale gas turbine combustor developed at DLR. This modified version of a practical gas turbine combustor has been extensively investigated by Weigand, Meier and coworkers [133,134]. This burner setup provides a suitable test case for verification and validation of combustion models, given the challenging complexity of the flow and the availability of a comprehensive set of experimental measurements. Donini et al. [135] applied the FGM method in combination with a presumed PDF approach in an LES of this setup. Since the fuel (methane) is injected separately from the air, the mixture fraction has to be used as additional control variable. Additionally, enthalpy is used as control variable to account for heat losses due to gas radiation and convective cooling at the combustor walls. The complex flow behavior characterized by inner and outer recirculation zones and a precessing vortex core, as well as the stabilization of the flame are accurately predicted by the numerical model. A detailed discussion of the flow pattern and a comparison with measurements is given in Ref. [135]. In Fig. 39, the computed instantaneous and time-averaged distributions of the control variables are shown. High turbulence levels in the shear layers near the inlet lead to fast mixing of fuel and air. At the location of the flame front, they are almost completely premixed and burn in a premixed mode. This is confirmed by an analysis of the resolved scalar dissipation rate χ, which drops from approximately 10 3 near the inlet to values below 1 s −1 at the flame front. The scalar dissipation rate in the flame surface, for which χt ≤ χ holds, is therefore much smaller than the flame time scale sL/δf. The low scalar dissipation rate indicates that the related Q-terms in the flamelet equations can be ignored. Heat loss to the combustor walls leads to a significant enthalpy deficit in the outer recirculation, but the region where the flame stabilizes is nearly adiabatic. Strong enthalpy gradients can therefore be neglected in the FGM modeling of this case. However, when the flame front comes close to the wall, local quenching may occur and the enthalpy gradients will have the same length scale as the flame thickness. In that case, large Q-terms are expected as was demonstrated for the laminar flame in Section 3.3.2.
The FGM model was also used to predict NO concentrations in this combustor by solving a transport equation with the source term obtained from the look-up table as discussed in Section 3.3.4. The predicted NO level in the exhaust of the combustor was in excellent agreement with the measured value of 6 ppm. The influence of gas radiation on the NO emission was found to be small, but ignoring heat loss to the walls resulted in an overprediction of 50% due to an overestimation of flame temperature.

Concluding summary of Part II
In Part II, the basics of the FGM approach have been presented and its application in laminar and turbulent premixed flames has been discussed. The principle idea of the method consists of a combination of reduced chemistry and laminar flamelet models. It was shown that the FGM method can be derived from the strongly stretched flamelet equations taking the major convective and diffusive effects into account. A complete FGM model is obtained when the actual flame stretch and curvature profiles from the CFD calculation are included in the SSFE, which have to be solved at run time. In practice, this is computationally too expensive and the SSFE are solved in a pre-processing step assuming representative conditions. The representative flamelets constitute a manifold which is parameterized by control variables and stored in a database. Methods to store the manifold and to couple it to a flow solver are presented and discussed.
The FGM method was further validated and applied step-bystep to various laminar combustion cases. First, the modeling of stretch effects was investigated. It was shown that a manifold constructed from stretchless flamelets can accurately predict the mass burning rate of stretched laminar flames because the chemical source terms are hardly influenced by stretch. Only in combination with preferential diffusion, flame stretch leads to a change in chemical source terms, because the element mass fractions and enthalpy change. When these changes are included in the manifold by adding additional dimensions, the effect of preferential diffusion on the mass burning rate is predicted by FGM. While for low Karlovitz numbers one additional dimension is enough, more degrees of freedom are required at high stretch rates.
Second, the modeling of heat loss by using FGM was studied. It was shown that by using a 2D FGM with enthalpy as an additional control variable, the stabilization of a premixed flame on a cooled burner can be accurately modeled. Stretch and curvature effects in the tip of this methane-air flame were accurately predicted by the 2D FGM. Close to the burner wall, where the flame quenches, diffusion along iso-surfaces of the progress variable cannot be neglected leading to small inaccuracies.
Third, the FGM method was used to simulate triple flames in order to investigate its application in partially-premixed systems. To account for mixture fraction variations in the FGM method, Z was added to the manifold as an additional control variable. The agreement between the detailed chemistry and FGM results was very satisfactory. Since transport of species along iso-surfaces of Y is not included in the FGM computations, differences appear for species, whose concentration changes significantly along the premixed flame front. The differences between FGM and full chemistry results are small for the flames studied here because the length scale of the variations along the flame front is much larger than the premixed flame thickness. When the length scale of mixing ∇ − Z 1 is in the same order as δf, additional degrees of freedom will be required in the FGM. When ∇ − Z 1 δ f , it is more appropriate to use diffusion flamelets to generate the manifold (see Verhoeven et al. [60]). A combined approach that switches from a premixed to a non-premixed table depending on the local gradients was investigated by Ramaekers et al. [58]. The enhanced accuracy of this combined approach for partially premixed flames, however, did not outweigh the additional complexity of the implementation in their study.
Finally, the modeling of NO formation was discussed in the context of laminar flame simulations. The slow chemical time scales of NO formation (or other slow chemical processes) form a challenge for the FGM method because they result in very steep gradients in the manifold, which lead to severe interpolation errors when NO or its source term is retrieved from the database. Adding NO to the definition of the progress variable results in a smoother mapping of the NO mass fraction and its chemical source term and reduces the interpolation errors. When a transport equation is solved for NO and its source term is well resolved, the FGM results closely match the detailed results. In general, we can conclude that the FGM method can predict accurately and efficiently the structure and dynamics of laminar flames including emission formation.
In Section 3.4, the application and accuracy of the FGM method in simulations of turbulent flames were investigated. First, the application of FGM in DNS was discussed. Premixed expanding flame kernels were simulated and their interaction with the turbulent flow was analyzed. Turbulent eddies distorted the preheat zone of the flame and lead to large varying flame stretch rates. The accuracy of the FGM approach was assessed by reconstructing flame paths form the turbulent flame and solving the flamelet equations for these paths with detailed chemistry and FGM. It was shown that the mean error in the predicted mass burning rate was 5% for a 1D manifold and that by using a 2D manifold based on stretched flamelets, this error was reduced to 1%.
A direct comparison with detailed chemistry calculations was shown for a 2D DNS of a methane-hydrogen flame kernel. The large preferential diffusion effects caused by flame stretch and curva-ture in this flame require additional dimensions in the manifold to capture changes in enthalpy and element mass fractions. As in the laminar test cases, a 2D FGM reproduces the detailed results very well, while the 1D FGM cannot account for the changes in h and Z j. Small deviations in highly curved regions of the flame are attributed to large Q-terms due to large tangential diffusion fluxes. In order to improve the FGM modeling, these Q-terms might have to be included in the flamelet calculations. The relevance of these preferential diffusion effects was investigated in FGM-DNS of premixed methane-hydrogen-air combustion on a slot burner. The variations in stoichiometry lead to changes in local burning rate that in turn lead to more flame wrinkling and therefore a shorter averaged flame length with an increased turbulent flame speed.
For application in LES or RANS, the FGM method needs to be coupled with a turbulence model that accounts for unresolved terms in the governing equations. Many different approaches exist to close these terms and to couple FGM with LES or RANS. The presumed PDF approach is probably the most widely used method, but more advanced methods are being developed and are gaining interest. Several test cases were discussed, in which the role of fuel stratification and heat loss was investigated. It was found that LES with manifolds based on premixed flamelets with additional degrees of freedom to account for variations in mixture fraction and enthalpy, in general yields satisfactory agreement with measurements. The results obtained for the DLR gas-turbine combustor demonstrate the large potential of FGM in combination with LES for the modeling of combustion in real devices.

General conclusions and outlook
State-of-the-art of the Flamelet Generated Manifold (FGM) technique has been reviewed and consolidated. The FGM technique is one possible solution to the important problem to model the detailed structure of flames described by complex chemistry and transport processes in a much more efficient and fast way. This is not done by solving all transport equations at hand, but by reducing this set to a few equations describing the main progress in the flame, while detailed information on all other relevant parameters is stored in a database, prepared before the actual flame computation. The basic idea of the FGM technique is that the flame structure is a combination of a preheating zone where detailed transport phenomena including stretch, curvature and preferential diffusion effects take place and a much thinner reaction layer in which only chemical reaction and molecular diffusion take place and balance each other. The chemical source term, perturbed and balanced by molecular diffusion effects, is the key information stored in the FGM database. This chemical source is not very sensitive to stretch and curvature, but major effects due to changes in conserved variables like pressure, enthalpy and elemental composition at the reaction layer must be taken into account. This is invoked by creating a multiple-dimensional manifold, generated by choosing a set of appropriate flamelets in which the same or similar perturbations appear. A generalized flame stretch theory, based on these ideas, has been developed to show this and to analyze the flame dynamics. The FGM technique has been applied with success to many DNS studies of laminar and turbulent flames in which all phenomena in the preheating zone are resolved, while the source term has been prepared and stored on the basis of the above principles. Heat loss, fuel stratification, flame stretch and preferential diffusion effects have been accounted for. Diffusion transport along iso-surfaces of the progress variable is, however, not accounted for. This leads to inaccuracies in regions of intense cooling or mixing where large gradients of enthalpy and element mass fractions occur. It remains to be investigated how these effects can be included in future FGM approaches.
The application of the FGM model in 'unresolved' models using LES or RANS for the flow field is also studied. To cover the unresolved behavior, various approaches exist of which the presumed PDF method is probably the most applied one. This model has appeared to be very promising as well, but there are many open questions to be asked. Important questions are, for instance, what the shape and structure of the PDF should be and how to model accurately the propagation of a flame front which is much thinner than the mesh size. New ideas in the direction of this last problem are currently under investigation, e.g., by the preparation of flamelets, suitably filtered to recover the correct propagation and dynamics of the thin embedded fronts (see e.g., Refs. [107,109,[122][123][124]). Another question is how to include preferential diffusion effects in 'unresolved' models. In DNS the interaction between flame stretch and preferential diffusion was found to play an important role, but in LES this effect happens at the sub-grid scale and requires proper modeling.

A. Transformation rules
Let us consider the transformation rules belonging to the coordinate transformation from a Cartesian system to a moving curvilinear system x, a a a (97) For the case described by Eq. (9), the (symmetric) covariant and contravariant metric tensors are given by g g 11 11 Y denotes the derivative of Y with respect to ξ 1 . We assume that ′ > Y 0 . Using the Einstein summation convention, i.e., products containing repeated indices should be summed over the appropriate range of this index, we have for the differential operators in this paper the following expressions Another difference is that in the present case the transformation is locally orthogonal, while this is not the case for the method of Peters.