CFD Online Logo CFD Online URL
Home > Wiki > Combustion


From CFD-Wiki

(Difference between revisions)
Jump to: navigation, search
m (Description of Real Mixtures)
m (Governing Equations for Chemically Reacting Flows)
Line 994: Line 994:
their transport properties relatively to each other. They are also coupled to the Navier-Stokes equations (mass, momentum, and energy) thanks to density, velocity, pressure and temperature taking a role in their diffusion and convection fluxes and their source terms.
their transport properties relatively to each other. They are also coupled to the Navier-Stokes equations (mass, momentum, and energy) thanks to density, velocity, pressure and temperature taking a role in their diffusion and convection fluxes and their source terms.
The transport equation for the mass fraction <math> Y_k </math> of <i>k-th</i> species is
Due to these complicated behaviours, it is rather chosen to speak in terms of
diffusion velocities <math> V^j_k </math> such that the transport equation for the mass fraction <math> Y_k </math> of the <i>k-th</i> specie in a complex mixture
takes the form:
\frac{\partial}{\partial t} \left( \rho Y_k \right) +
\frac{\partial}{\partial t} \left( \rho Y_k \right) +
\frac{\partial}{\partial x_j} \left( \rho u_j Y_k\right) =  
\frac{\partial}{\partial x_j} \left( \rho u^j Y_k\right) = -
\frac{\partial}{\partial x_j} \left( \rho D_k \frac{\partial Y_k}{\partial x_j}\right)+ \dot \omega_k
\frac{\partial}{\partial x_j} \left( \rho V_k^j Y_k\right)+ \dot \omega_k
<math> \dot \omega_k </math> denotes the specie reaction rate.
where Ficks' law is assumed for scalar diffusion with <math> D_k </math> denoting the species difussion coefficient, and <math> \dot \omega_k </math> denoting the species reaction rate.
A non-reactive (passive) scalar (like the mixture fraction <math> Z </math>) is goverened by the following transport equation
A non-reactive (passive) scalar (like the mixture fraction <math> Z </math>) is goverened by the following transport equation
Line 1,103: Line 1,103:
== Reaction mechanisms ==
== Reaction mechanisms ==

Revision as of 03:08, 14 May 2009

The power of Fire, or Flame, for instance, which we designate by some trivial chemical name, thereby hiding from ourselves the essential character of wonder that dwells in it as in all things, is with these old Northmen, Loke, a most swift subtle Demon of the brood of the Jötuns... From us too no Chemistry, if it had not Stupidity to help it, would hide that Flame is a wonder. What is Flame?

Carlyle on Heroes Odin and Scandinavian Mythology., [1].


What is combustion -- physics versus modelling

Combustion phenomena consist of many physical and chemical processes which exhibit a broad range of time and length scales. A mathematical description of combustion is not always trivial, although some analytical solutions exist for simple situations of laminar flame. Such analytical models are usually restricted to problems in zero or one-dimensional space.


In a general manner, we shall try to stick to the official standards, the nomenclature for heat transfer, that may be found in an issue of Journal of Heat Transfer [2]

u' Turbulent integral RMS velocity

Fundamental Aspects

Main Specificities of Combustion Chemistry

Combustion can be split into two processes interacting with each other: thermal, and chemical.

The chemistry is highly exothermal (this is the reason of its use) but also highly temperature dependent, thus highly self-accelerating. In a simplified form, combustion can be represented by a single irreversible reaction involving 'a' fuel and 'an' oxidizer:

 \frac{\nu_F}{\bar M_F}Y_F + \frac{\nu_O}{\bar M_O} Y_O \rightarrow Product + Heat

Although very simplified compared to real chemistry involving hundreds of species (and their individual transport properties) and elemental reactions, this rudimentary chemistry has been the cornerstone of combustion analysis and modelling.

The most widely used form for the rate of the above reaction is the Arrhénius law:

 \dot\omega = \rho^2 A \left ( \frac{Y_F}{\bar M_F} \right )^{n_F} \left (\frac{Y_O}{\bar M_O}\right )^{n_O} \exp^{-T_a/T}

 T_a is the activation temperature, high in combustion, consistently with the temperature dependence. This is where the high non-linearity in temperature is modelled. A is the pre-exponential constant. One of the interpretation of the Arrhénius law comes from gas kinetic theory: the number of molecules whose kinetic energy is larger than the minimum value allowing a collision to be energetic enough to trig a reaction is proportional to the exponential term introduced above divided by the square root of the temperature. This interpretation allows one to think that the temperature dependence of A is very weak compared to the exponential term. A is eventually considered as constant. The reaction rate is also naturally proportional to the molecular density of each of the reactant. Nonetheless, the orders of reaction  n_i are different from the stoichiometric coefficients as the single-step reaction is global, not governed by collision for it represents hundreds of elementary reactions. If one goes into the details, combustion chemistry is based on chain reactions, decomposed into three main steps: (i) generation (where radicals are created from the fresh mixture), (ii) branching (where products and new radicals appear from interaction of radicals with reactants), and (iii) termination (where radicals collide and turn into products). The branching step tends to accelerate the production of active radicals (autocatalytic). The impact is nevertheless small compared to the high non-linearity in temperature. This explains why single-step chemistry has been sufficient for most of the combustion modelling work up to now.

The fact that a flame is a very thin reaction zone separating, and making the transition between, a frozen mixture and an equilibrium is explained by the high temperature dependence of the reaction term, modelled by a large activation temperature, and a large heat release (the ratio of the burned and fresh gas temperatures is about 7 for typical hydrocarbon flames) leading to a sharp self-acceleration in a very narrow area. To evaluate the order of magnitude of the quantities, the terms in the exponential argument are normalized:

 \beta=\alpha\frac{T_a}{T_s} \qquad \alpha=\frac{T_s-T_f}{T_s}

\beta is named the Zeldovitch number and \alpha the heat release factor. Here,  T_s has been used instead of  T_b, the conventional notation for burned gas temperature (at final equilibrium).  T_s is actually  T_b for a mixture at stoichiometry and when the system is adiabatic, i.e. this is the reference highest temperature that can be obtained in the system.  T_f is the ambient temperature of the fresh gases. That said, typical value for \beta and \alpha are 10 and 0.9, giving a good taste of the level of non-linearity of the combustion process with respect to temperature. Actually, the reaction rate is rewritten as:

 \dot\omega = \rho^2 A \left ( \frac{Y_F}{\bar M_F} \right )^{n_F} \left (\frac{Y_O}{\bar M_O}\right )^{n_O} 

where the non-dimensionalized temperature is:


The non-linearity of the reaction rate is seen from the exponential term:

  •  {\mathcal O}(\exp^{-\beta}) for \theta far from unity (in the fresh gas)
  •  {\mathcal O}(1) for \theta close to unity (in the reaction zone close to the burned gas whose temperature must be close to the adiabatic one  T_s ), more exactly  1-\theta \sim {\mathcal O}(\beta^{-1})
Temperature Non-Linearity of the Source Term: the Temperature-Dependent Factor of the Reaction Term for Some Values of the Zeldovitch and Heat Release Parameters
Note that for an infinitely high activation energy, the reaction rate is piloted by a \delta(\theta) function. The figure, beside, illustrates how common values of \beta around 10 tend to make the reaction rate singular around \theta of unity. Two set of values are presented: 
\beta = 10 and \beta = 8. The first magnitude is the representative value while the second one is a smoother one usually used to ease numerical simulations. In the same way, two values for the heat release \alpha, 0.9 and 0.75, are explored. The heat release is seen to have a minor impact on the temperature non-linearity.

Transport Equations

Additionally to the Navier-Stokes equations, at least with variable density, the transport equations for a reacting flow are the energy and species transport equations with their appropriate boundary conditions that may be of Dirichlet or Neumann type. In usual notations, the specie i transport equation is written as:

\frac{D \rho Y_i}{D t} = \nabla\cdot \rho D_i\vec\nabla Y_i - \nu_i\bar M_i\dot\omega

and the temperature transport equation:

\frac{D\rho C_p T}{Dt} = \nabla\cdot \lambda\vec\nabla T + Q\nu_F \bar M_F \dot\omega

The diffusion is modelled thanks to Fick's law that is a (usually good) approximation to the rigorous diffusion velocity calculation. Regarding the temperature transport equation, it is derived from the energy transport equation under the assumption of a low-Mach number flow (compressibility and viscous heating neglected). The low-Mach number approximation is suitable for the deflagration regime (as it will be demonstrated below), which is the main focus of combustion modelling. Hence, the transport equation for temperature, as a simplified version of the energy transport equation, is usually retained for the study of combustion and its modelling.

Note: The species and temperature equations are not closed as the fields of velocity and density also need to be computed. Through intense heat release in a very small area (the jump in temperature in typical hydrocarbon flames is about seven and so is the drop in density in this isobaric process, and the thickness of a flame is of the order of the millimetre), combustion influences the flow field. Nevertheless, the vast majority of combustion modelling has been developed based on the species and temperature equations, assuming simple flow fields, sometimes including hydrodynamics perturbations.

Low-Mach Number Equations

In compressible flows, when the motion of the fluid is not negligible compared to the speed of sound (which is the speed at which the molecules can reorganize themselves), the heap of molecules results in a local increase of pressure and temperature moving as an acoustic wave. It means that, in such a system, a proper reference velocity is the speed of sound and a proper pressure reference is the kinetic pressure. A contrario, in low-Mach number flows, the reference speed is the natural representative speed of the flow and the reference pressure is the thermodynamic pressure. Hence, the set of reference quantities to characterize a low-Mach number flow is given in the table below:

Density \rho_o A reference density (upstream, average, etc.)

Velocity U_o A reference velocity (inlet average, etc.)

Temperature T_o A reference temperature (upstream, average, etc.)

Pressure (static) P_o=\rho_o \bar r T_o From Boyle-Mariotte

Length L_o A reference length (representative of the domain)

Time L_o/U_o

Energy C_p T_o Internal energy at constant reference pressure.  C_p must also be chosen in a reference thermodynamical state.

The equations for fluid mechanics properly adimensionalized can be written:

Mass conservation:

 \frac{D\rho}{Dt} =0


\frac{D\rho\vec U}{Dt}=-\frac{1}{\gamma M}\vec\nabla P+\nabla\cdot\frac{1}{Re}\bar\bar\Sigma

Total energy:

\frac{D\rho e_T}{Dt}-\frac{1}{RePr}\nabla\cdot\lambda\vec\nabla T=-\frac{\gamma-1}{\gamma}\nabla\cdot P\vec U
+\frac{1}{Re}M^2(\gamma-1)\nabla\cdot \vec U\bar\bar\Sigma + \frac{\rho}{C_pT_o(U_o/L_o)}Q\nu_F \bar M_F\dot\omega


\frac{D\rho Y}{Dt}-\frac{1}{ScRe}\nabla\cdot\rho D\vec\nabla Y=-\frac{\rho}{U_o/L_o}\nu\bar M\dot\omega

State law:

 P=\rho T

The low-Mach number equations are obtained considering that  M^2 is small. 0.1 is usually taken as the limit, which recovers the value of a Mach number of 0.3 to characterize the incompressible regime.

Considering the energy equation, in addition to the terms with  M^2 in factor in the equation, the total energy reduces to internal energy as:  e_T = T/\gamma+M^2(\gamma-1)\vec U^2/2 . Moreover, the work of pressure is considered as negligible because the gradient of pressure is negligible (low-Mach number approximation is indeed also named isobaric approximation) and the flow is assumed close to a divergence-free state. For the same reason, volumic energy and enthalpy variations are assumed equal as they only differ through the addition of pressure. Hence, redimensionalized, the low-Mach number energy equation leads to the temperature equation as used in combustion analysis:

\frac{D\rho C_p T}{Dt} = \nabla\cdot \lambda\vec\nabla T + Q\nu_F \bar M_F \dot\omega

The Damköhler Number

A flame is a reaction zone. From this simple point of view, two aspects have to be considered: (i) the rate at which it is fed by reactants, let call  \tau_d the characteristic time, and the strength of the chemistry to consume them, let call the characteristic chemical time  \tau_c . In combustion, the Damköhler number, Da, compares these both time scales and, for that reason, it is one of the most integral non-dimensional groups:


If Da is large, it means that the chemistry has always the time to fully consume the fresh mixture and turn it into equilibrium. Real flames are usually close to this state. The characteristic reaction time,  (Ae^{-T_a/T_s})^{-1} , is estimated of the order of the tenth of a ms. When Da is low, the fresh mixture cannot be converted by a too weak chemistry. The flow remains frozen. This situation happens with ignition or misfire, for instance.

The picture of a deflagration lends itself to a description based on the Damköhler number. A reacting wave progresses towards the fresh mixture through preheating of the upstream closest layer. The elevation of the temperature strengthens the chemistry and reduces its characteristic time such that the mixture changes from a low-Da region (far upstream, frozen) to a high-Da region in the flame (intense reaction to equilibrium).

Conservation Laws

The processus of combustion transforms the chemical enthalpy into sensible enthalpy (i.e. rises the temperature of the gases thanks to the heat released). Simple relations can be drawn between species and temperature by studying the source terms appearing in the above equations:

 \frac{Y_F}{\nu_F\bar M_F} - \frac{Y_O}{\nu_O\bar M_O} = \frac{Y_{F,u}}{\nu_F\bar M_F} - \frac{Y_{O,u}}{\nu_O\bar M_O}
 Y_F + \frac{Cp T}{Q} = Y_{F,u} + \frac{CpT_u}{Q}

Those coupling functions are named Schwalb-Zeldovich variables.

Hence  T_b = T_u + \frac{Q Y_{F,u}}{Cp} ,  Y_{O,b} = Y_{O,u} - \frac{\nu_O \bar M_O}{\nu_F \bar M_F} Y_{F,u} = Y_{O,u} -sY_{F,u} and  Y_{F,b} = 0 . Here, the example has been taken for a lean case.

As mentioned in Sec. Main Specificities, the stoichiometric state is used to non-dimensionalize the conservation equations:

 Y_i^* = Y_{i,u}^* - \theta\frac{Y_{F,s}}{Y_{i,s}}s^{\delta_{i,O}} \qquad ; \qquad Y_i^*=Y_i/Y_{i,s}.

A comprehensive form of the reaction rate can be reconstituted to understand the difficulty of numerically resolving the reaction zone:

 \dot\omega = B \prod_{i=O,F}(Y_{i,u}^*-\theta\frac{Y_{F,s}}{Y_{i,s}}s^{\delta_{i,O}} )^{n_i} \exp{\left ( -\beta\frac{1-\theta}{1-\alpha(1-\theta)}\right)}

where  B stands for all the constant terms present in this reaction rate, plus density.

Source Term versus Temperature

For the stoichiometric case and a global order of two, the reaction rate is graphed versus the reduced temperature for different values of the heat release and Zeldovitch parameter. A high value of \beta makes the reaction rate very sharp, versus temperature. It means that reaction is significant beyond a temperature level (sometimes called ignition temperature) that is close to one (the exponential term above is non-negligible for 1-\theta \sim \beta^{-1}). The heat release has qualitatively the same impact but not so strong. Transposed to the case of a flame sheet, it effectively shows that the reaction exists only in a fraction of the thermal thickness of the flame (the region close to the flame that the latter preheats, hence, where the reduced temperature rises from 0 to 1 here) where the temperature deviates few from the maximal one (density can be assumed as constant and equal to its burned-gas value). Numerically capturing such a sharp reaction zone can be costly and the lower values of \beta and \alpha as presented here are usually preferred whenever possible.

Most problems in combustion involve turbulent flows, gas and liquid fuels, and pollution transport issues (products of combustion as well as for example noise pollution). These problems require not only extensive experimental work, but also numerical modelling. All combustion models must be validated against the experiments as each one has its own drawbacks and limits. In this article, we will address the modeling fundamentals only.

The combustion models are often classified on their capability to deal with the different combustion regimes.

Three Combustion Regimes

Depending on how fuel and oxidizer are brought into contact in the combustion system, different combustion modes or regimes are identified. Traditionally, two regimes have been recognized: the premixed regime and the non-premixed regime. Over the last two decades, a third regime, sometime considered as a hybrid of the two former ones to a certain extend, has risen. It has been named partially-premixed regime.

The Non-Premixed Regime

Sketch of a diffusion flame

This regime is certainly the easiest to understand. Everybody has already seen a lighter, candle or gas-powered stove. Basically, the fuel issues from a nozzle or a simple duct into the atmosphere. The combustion reaction is the oxidization of the fuel. Because fuel and oxidizer are in contact only in a limited region but are separated elsewhere (especially in the feeding system) this configuration is the safest. Diffusion of oxidant and fuel has to occur simultaneously to reaction to sustain combustion, the flame being a surface of separation of fuel and oxidant streams. The non-premixed flame has some other advantages. By controlling the flows of both reactants, it is (theoretically) possible to locate the stoichiometric interface, and thus, the location of the flame sheet. Moreover, the strength of the flame can also be controlled through the same process. Depending on the width of the transition region from the oxidizer to the fuel side, the species (fuel and oxidizer) feed the flame at different rates. This is because the diffusion of the species is directly dependent on the unbalance (gradient) of their distribution. A sharp transition from fuel to oxidizer creates intense diffusion of those species towards the flame, increasing its burning rate. This burning rate control through the diffusion process is certainly one of the reasons of the alternate name of such a flame and combustion mode: diffusion flame and diffusion regime.

Because a diffusion flame is fully determined by the inter-penetration of the fuel and oxidizer streams, it is convenient to introduce a tracer of the state of the mixture. This is the role of the mixture fraction, usually called Z or f. Z is usually taken as unity in the fuel stream and is null in the oxidizer stream. It varies linearly between this two bounds such that at any point of a frozen flow the fuel mass fraction is given by  Y_F=ZY_{F,o} and the oxidizer mass fraction by  Y_O = (1-Z)Y_{O,o} .  Y_{F,o} and  Y_{O,o} are the fuel and oxidizer mass fractions in the fuel and oxidizer streams, respectively. The mixture fraction posses a transport equation that is expected to not have any source term as a tracer of a mixture must be a conserved scalar. First, the fuel and oxidizer mass fraction transport equations are written in usual notations:

\frac{D \rho Y_F}{D t} = \nabla\cdot \rho D\vec\nabla Y_F - \nu_F\bar M_F\dot\omega
\frac{D \rho Y_O}{D t} = \nabla\cdot \rho D\vec\nabla Y_O - \nu_O\bar M_O\dot\omega

The two above equations are linearly combined in a single one in a manner that the source term disappears:

\frac{D \rho (\nu_O\bar M_O Y_F-\nu_F\bar M_F Y_O)}{D t} = \nabla\cdot \rho D\vec\nabla (\nu_O\bar M_O Y_F-\nu_F\bar M_F Y_O)

The quantity (\nu_O\bar M_O Y_F-\nu_F\bar M_F Y_O) is thus a conserved scalar, one recognizes one of the Schwalb-Zeldovich quantities introduced in Sec. Conservation Laws. One also remarks that this transport equations combination is made in the equidiffusional approximation, i.e. all the scalars, including temperature, diffuse at the same rate. The last step is to normalize it such that it equals unity in the pure fuel stream ( Y_F=Y_{F,o} and  Y_O=0 ) and is null in the pure oxidizer stream ( Y_F=0 and  Y_O=Y_{O,o} ). The resulting normalized passive scalar is the mixture fraction:

Z=\frac{\Phi(Y_F/Y_{F,o})-(Y_O/Y_{O,o})+1}{\Phi+1}\qquad \Phi=\frac{\nu_O\bar M_O Y_{F,o}}{\nu_F\bar M_F Y_{O,o}}

governed by the transport equation:

\frac{D \rho Z}{D t} = \nabla\cdot \rho D\vec\nabla Z

The stoichiometric interface location (and thus the approximate location of the flame if the flow is reacting) is where \Phi(Y_F/Y_{F,o})-(Y_O/Y_{O,o}) vanishes (or Y_F and  Y_O are both null in the reacting case). This leads to a stoichiometry definition:


As the mixture fraction qualifies the degree of inter-penetration of fuel and oxidizer, the elements originally present in these molecules are conserved and can be directly traced back to the mixture fraction. This has led to an alternate defintion of the mixture fraction, based on element conservation. First, the elemental mass fraction  X_{j} of element j is linked to the species mass fraction  Y_i :

 X_j = \sum_{i=1}^{n} \frac{a_{i,j} \bar M_j}{\bar M_i} Y_i

where  a_{i,j} is a matrix counting the number of element j atoms in specie molecule named i and n is the number of species in the mixture. The group pictured by the summation above is a linear combination of  Y_i . Because the transport equations of species mass fraction, few lines earlier, are also linear, a transport equation for the elemental mass fraction can be written:

\frac{D \rho X_j}{D t} = \nabla\cdot \sum_{i=1}^n\frac{a_{i,j}\bar M_j}{\bar M_i}\rho D_i\vec\nabla Y_i

For mass is conserved, the linear combination of the source terms vanishes. Furthermore, by taking the same diffusion coefficient  D_i for all the species, the elemental mass fraction transport equation has exactly the same form as the specie transport equation (except the source term). Notice that the assumption of equal diffusion coefficient was also made in the previous definition of the mixture fraction and is justified in turbulent combustion modelling by the turbulence diffusivity flattening the diffusion process in high Reynolds number flows. Hence, the elemental mass fraction transport equation has the same structure as the mixture fraction transport equation seen above. Properly renormalized to reach unity in the fuel stream and zero in the oxidizer stream, the elemental mass fraction is a convenient way of determining the mixture fraction field in a flow. Indeed, it is widely used in practice for this purpose.

Simplified Diffusion Flame Solution

Because Schwalb-Zeldovitch variables rid off any source, they allow a first approximation of the flame description independent of the kinetics. From the definition of Z, the Schwalb-Zeldovich variables have the following expression:

\left \{
Y_O-sY_F & = & Y_{O,o}-Z(sY_{F,o}+Y_{O,o}) \\
\frac{Q Y_F}{C_p} + T & = & T_{O,o} + Z(T_{F,o}-T_{O,o}+\frac{Q Y_{F,o}}{C_p})

Two simplified limit cases for the diffusion flame equation can be built from those relations:

  • In the frozen flow, the chemical source term is null everywhere. Hence,  Y_F and  Z follow

strictly the same linear transport equation. Then, Y_F = Y_{F,o} Z as anticipated above. Follow:  Y_O=Y_{O,o}(1-Z) ,  T=T_{O,o}+Z(T_{F,o}-T_{O,o}) .

  • The equilibrium solution, also with zero source term, the 'historical' diffusion flame description laid out by Burke and Schumann in 1928:
    • Rich side:  Y_O = 0 ,  Y_F = -Y_{O,o}/s + Z(Y_{F,o}+Y_{O,o}/s),  T=T_{F,o}+(1-Z)(T_{O,o}-T_{F,o}+Y_{O,o}\frac{Q}{sC_p})
    • Lean side:  Y_F = 0 ,  Y_O = Y_{O,o} - Z(sY_{F,o}+Y_{O,o}),  T=T_{O,o}+Z(T_{F,o}-T_{O,o}+Y_{F,o}\frac{Q}{C_p})

From the value of Z at stoichiometry above, the adiabatic flame temperature is obtained:

T_s=T_{O,o}+Z_s(T_{F,o}-T_{O,o}+Y_{F,o} \frac{Q}{C_p})

Remark 1: The Burke-Schumann solution for a diffusion flame, despite its age, is still one of the most used model, under some refreshments such as Mixed-Is-Burned, or Equilibrium combustion model where the burned value of the quantities is obtained by Gibbs function extrema instead of the full adiabatic consumption. The reason is that the concentration of the species and temperature can be calculated from the mixture fraction field, which is in turn described in terms of statistics in the turbulent flow.

Remark 2: In the Burke-Schuman solution, as the chemistry is infinitely fast, the flame is reduced to an interface separating the rich and lean sides. Therefore, quantities must be continuous at the flame and their gradients observe the jump conditions, from Sec. Conservation Laws:

\left \{
s\rho D_F||\vec\nabla Y_F|| & = & \rho D_O||\vec\nabla Y_O|| \\
\frac{Q}{C_p}\rho D_F||\vec\nabla Y_F|| & = & \lambda ||\vec\nabla T||
\right .

Remark 3 From Remark 2, it is seen that the reactants fluxes must reach the flame in stoichiometric proportion (the fluxes and not the values themselves). Therefore, the relations developed in Remark 2 are true even without the equidiffusional approximation, whereas the relation based on  Z , no. On the other hand, the heat released by chemistry is conducted away from the flame thanks to the temperature gradients on both sides.

Dissipation Rate and Non-Equilibrium Effect

A very important quantity, derived from the mixture fraction concept, is the scalar dissipation rate, usually noted:  \chi . In the above introduction to non-premixed combustion, it has been said that a diffusion flame is fully controlled through: (i) the position of the stoichiometric line, dictating where the flame sheet lies; (ii) the gradients of fuel on one side and oxidizer on the other side, dictating the feeding rate of the reaction zone through diffusion and thus the strength of combustion. According the the mixture fraction definition, the location of the stoichiometric line is naturally tracked through the Z_s iso-line and it is seen here how the mixture fraction is a convenient tracer to locate the flame. In the same manner, the mixture fraction field should also be able to give information on the strength of the chemistry as the gradients of reactants are directly linked to the mixture fraction distribution. The feeding rate of the reaction zone is characterized with the help of the inverse of a time. Because it is done through diffusion, it must be obtained through a combination of mixture fraction gradient and diffusion coefficient (dimensional analysis):

 \chi_s = \frac{(\rho D)_s}{\rho_s} ||\vec \nabla Z||^2_s \qquad (\rho D)_s = \left ( \frac{\lambda}{C_p} \right )_s

where the subscript s refers to quantities taken effectively where the reacting sheet is supposed to be, close to stoichiometry. This deduction of the scalar dissipation rate, scaling the feeding rate of the flame, is obtained here through physical arguments and is to be derived from equations below in a more mathematical manner. Note that the transport coefficient for the mixture fraction is identified to the one for temperature. This notation is usually used in the literature to emphasize that the rate of temperature diffusion (that is commensurable to the rate of species diffusion) is the retained parameter (as introduced in the following approach highlighting the role of the scalar dissipation rate). Depending if the diffusion process is free (mixing layer, boundary layer) or controlled (stagnation), the dissipation rate and then the non-premixed combustion may be of unsteady or steady nature.

Ignition / Burning / Extinction Curve

Let consider a mixing layer between fuel and oxidizer whose strain (and then intensity of reactant inter-diffusion) is carefully controlled. Let start from a flame already existing. Let also use the Damköhler number as defined above with the inverse of the scalar dissipation rate being the characteristic mechanical time of flow of reactants feeding the reaction zone. For high Da (low dissipation rate), the reactants diffuse slowly. The reaction is not very intense but the chemistry is fully achieved such that the maximum temperature is reached. Now, let decrease slowly (to avoid any unsteady effect on the chemistry) the Da. On one hand, the rate of feeding of flame through diffusion increases, and so does the reaction rate (that said, the reaction rate normalized by the dissipation rate decreases as the Da). On the other hand, the chemistry may not have the time to `eat' every reactant molecules that begin to leak through the flame: the fully sensible value is not realized and a lower temperature is achieved. With still increasing the dissipation rate, this mechanism leads to lower temperature in the flame zone down to a level that cannot be sustained by combustion (that strongly depends on temperature). The diffusion flame leaves the diffusion-controlled burning regime and extinguishes suddenly. It is said to be quenched. This is experienced in real life when blowing off a small wood fire. Slightly blowing increases the transfer between reactants and strengthens the reaction but blowing too much extinguishes the fire.

S-Curve Diffusion Flame.

If now the start state is a frozen mixing layer between fuel and oxidant, at low Da (high dissipation rate), the flow remains frozen. When increasing Da, there will be a point where the chemistry will self-accelerate and the flame will light up through a sudden increase in temperature. Because the starting temperature is low, the ignition Da number is higher that the extinction Da above, exhibiting an hysteresis phenomenon.

When looking at the trend of the maximum temperature only versus Da, a curve with a shape in `S' appears, named as S-Curve for diffusion flames (see beside figure).

Flamelet Equation

In turbulent combustion modelling, the flamelet regime is identified as a regime that saves the integrity of the flame structure (Sec. The Wrinkled Regime). Turbulent eddies do not enter the structure and only contord the flame at a large scale. Therefore, it is usual to call the flamelet equation, the equation that describes the diffusion flame structure when not perturbed by turbulence. The model problem is usually retained as the counterflow diffusion flame. The strain imposed to the flame mimicks the one due to inhomogeneities in a real turbulent flow.

Sketch of the Stretched Mixing Layer Created by the Counter-Flow Configuration.

A counter-flow diffusion flame is basically made of two ducts, front-to-front, one issuing fuel at a mass fraction  Y_{F,o} (may be diluted by an inert) and temperature  T_{F,o}, and the other issuing oxidizer at a mass fraction  Y_{O,o} (the inert is usually di-nitrogen in air such that  Y_{O,o} =.23310 is a value commonly encountered) and temperature  T_{O,o}. A stretched mixging layer develops at the location the flows oppose. This stretch can be modulated by the flow rate from the feeding ducts. The diffusion flame sits around the iso-surface corresponding to stoichiometry. The figure beside sketches this widely used configuration, with the coordinate frame origin taken at the stagnation location in the middle of the mixing layer.

The Howarth-Dorodnitzyn transform and the Chapman approximation are applied to the above mixture fraction transport equation. In the Chapman approximation, the thermal dependence of (\rho D) is approximated as \rho^{-1}. The Howarth-Dorodnitzyn transform introduces \rho in the space coordinate system:  \vec\nabla \rightarrow \rho\vec\nabla \quad ; \quad \nabla\cdot\rightarrow \rho\nabla\cdot. The effect of these both mathematical operations is to `digest' the thermal variation of quantities such as density or transport coefficient. Hence, the mixture fraction equation comes in a simpler mathematical shape:

\rho\vec U \cdot \vec\nabla Z = \rho_s ( \rho D )_s \Delta Z

Here the references are taken in the stoichiometric zone (s subscript).

For the sake of simplicity, we shall use the potential flow result for this counter-flow problem with two identical jets although the presence of the flame (the associated density change) jeopardizes this theory, strictly speaking. The velocity solution is thus:

\left \{
u & = & -a\int\frac{dx}{\rho} \\
v & = & \frac{a}{2}\int\frac{dr}{\rho}

a is named the stretch of the mixing layer and is the inverse of a time.

In such a flow, the transport equation for the mixture fraction may be written as a purely 1-D problem (assuming constant properties):

-ax\frac{dZ}{dx} = \rho_s (\rho D)_s \frac{d^2Z}{dx^2}

from which a `length' scale emerges for the mixing layer:  \sqrt{\rho_s (\rho D)_s /a} . Non-dimensionalized, with 
\eta = x \sqrt{a/\rho_s (\rho D)_s} :

\eta \frac{dZ}{d\eta} + \frac{d^2 Z}{d \eta^2} = 0

The boundary conditions may be approximated as  \eta \rightarrow +\infty\quad ; \quad Z\rightarrow 1 and  \eta \rightarrow -\infty\quad ; \quad Z\rightarrow 0 . Note that this approximation requires a thin layer compared to the counter-flow distance, and thus a relatively high strain rate. Then, the solution for  Z :


The scalar dissipation rate is deduced from this equation as:

\chi = \frac{1}{2\pi}\exp\left( - [erf^{-1} (2Z)]^2 \right) \frac{a\rho^2}{\rho_s^2}

Because combustion is highly temperature-dependent, T is certainly the scalar to which attention must be paid for. The temperature equation in the Low-Mach Number regime (Sec. Low-Mach Number Equations) is written below in steady-state:

\rho\vec U \cdot \vec\nabla T = \nabla\cdot \frac{\lambda}{C_p}\vec\nabla T + \frac{Q}{C_p}\nu_F \bar M_F \dot\omega

After the Chapman approximation and the Howarth-Dorodnitzyn transform:

\frac{\rho}{\rho_s}\vec U \cdot \vec\nabla T = \left ( \frac{\lambda}{C_p}\right )_s \nabla\cdot  \vec\nabla T + \frac{Q}{\rho\rho_s C_p}\nu_F \bar M_F \dot\omega

Here the references are taken in the flame, i.e. close to the stoichiometric line (s subscript).

In a non-premixed system, strictly speaking,  \vec U , is not really relevant as the flame is fully controlled by the diffusion process. Notwithstanding, in practice, non-premixed flames must be stabilized by creating a strain in the direction of diffusion. This is the reason why the velocity is left in the equation. Because a diffusion flame is fully described by the mixture fraction field, a change in coordinate can be applied:

 \left (
r \\
\end{array}\right ) \longrightarrow 
\left (
r \\
\end{array}\right )

The Jacobian of the transform is given as:

 \left [
\frac{\partial r}{\partial r} & \frac{\partial Z}{\partial r} \\
\frac{\partial r}{\partial x} & \frac{\partial Z}{\partial x} 
\end{array}\right ] = 
\left [
1 & 0 \\
0 & (\rho l_d)^{-1} 
\end{array}\right ]

Note that the diffusive layer of thickness l_d is defined as the region of transition between fuel and oxidizer and is thus given by the gradient of Z along the x direction. This transform is applied to the vectorial operators:

 \nabla\cdot = \nabla_r\cdot+\nabla_x\cdot=\nabla_r\cdot+\nabla_x\cdot Z\nabla_Z\cdot
 \vec\nabla = \vec\nabla_r+\vec\nabla_x=\vec\nabla_r+\vec\nabla_x Z\nabla_Z\cdot

With this transform, the above temperature equation looks like:

\frac{\rho}{\rho_s}\vec U \cdot (\vec\nabla_r T + \vec\nabla Z \nabla_Z\cdot T) = \left (\frac{\lambda}{C_p}\right )_s \left ( \nabla_r\cdot \vec\nabla_r T +  \nabla_x\cdot Z \nabla_Z\cdot \vec\nabla_x Z \nabla_Z\cdot   T + \nabla_x\cdot Z \nabla_Z\cdot  \vec\nabla_r  T + \nabla_r\cdot \vec\nabla_x Z \nabla_Z\cdot  T \right ) + \frac{Q}{\rho\rho_s C_p}\nu_F \bar M_F \dot\omega

As mentioned above, the velocity and the variation along the tangential direction to the main flame structure r are not supposed to play a major role. To be convinced, a variable scaling may be done, considering that the reaction zone extends over a small fraction  \varepsilon of the diffusion thickness  l_d around stoichiometry. Then the convective term of the above equation is  {\mathcal O} ( \varepsilon^{-1} ) and the diffusive term is  {\mathcal O} ( \varepsilon^{-2} ) . By emphasizing the role of the gradient of Z along x as a key parameter defining the configuration the following equation is obtained:

0 = \left (\frac{\lambda}{C_p}\right )_s ||\vec\nabla_x Z||^2 \Delta_Z  T + \frac{Q}{\rho\rho_s C_p}\nu_F \bar M_F \dot\omega

This equation (sometimes named the flamelet equation) serves as the basic framework to study the structure of diffusion flames. It highlights the role of the dissipation rate with respect to the strength of the source term and shows that the dissipation rate calibrates the combustion intensity.


The flamelet equation may be used to explore the behaviour of the mixing-layer for not too high strain rate. As mentioned earlier, the Damk&oumlaut;hler may be high enough to allow a self-ignition. In this regime, the following assumptions may be made:

  • the exponential argument in the reaction term,  T_a/T, is linearized around the frozen temperature value (indeed, as the activation temperature is large, a small change in the temperature from its frozen value is magnified in the reaction rate magnitude):  T_a/T = T_a/T_o-T_a/T_o^2(T-T_o) = T_a/T_o - \phi . It can be seen that a small change of temperature leading to  T_a/T_o^2(T-T_o)\quad {\mathcal O} (1) impacts the exponential term of the reaction rate.
  • the reactants may be approximated by their frozen flow values:  Y_F=Y_{F,o}Z-\frac{Cp T_o^2}{QT_a}\phi\approx Y_{F,o}Z \quad ; \quad Y_O=Y_{O,o}(1-Z)-s\frac{Cp T_o^2}{QT_a} \phi \approx Y_{O,o}(1-Z) , under the hypothesis of large activation energy.

The flamelet equation is recast as: 
\Delta_Z\phi = - Z^{n_F}(1-Z)^{n_O} e^{\phi}\overbrace{\frac{Q\nu_F}{\rho_s\lambda_s}\frac{A}{\nabla^2_x Z \bar M_F^{n_F-1}\bar M_O^{n_O}} \frac{T_a}{T_o^2} e^{-\frac{T_a}{T_o}}Y_{F,o}^{n_F}Y_{O,o}^{n_O}}^{\mathcal D}
In this expression \mathcal D is a Damk\"ohler number.

This expression must be solved numerically but, with the help of simplifications, an analytical solution may be found. It must be noticed that this analytical solution has no quantitative prediction capabilities. However, it is sufficient to illustrate the critical Damk\"ohler number for ignition for pedagogical purpose.

First, the equation is integrated from Z_s to  Z 'somewhere in the rich side'.

\int_{Z_s}^{Z} dZ \Delta_Z\phi = \left [\frac{d\phi}{dZ} \right]_{Z_s}^{Z} = - \int_{Z_s}^{Z} dZ Z^{n_F}(1-Z)^{n_O} e^{\phi}{\mathcal D}

 \phi may be approximated as, on the asymptote  \hat{\phi}(Z-1)/(Z_s-1) with a maximum  \hat{\phi} on the stoichiometric line, where ignition is supposed the strongest. Then, on the stoichiometric line, the derivative of  \phi vanishes and, sufficiently far from the stoichiometric line it tends to  \hat{\phi}/(Z_s-1) . To keep the solution simple, we choose the integration limit to be on the border of the region where  \phi remains within 10% of its maximum such that it may be approximated as constant on the right hand side. The value of  Z that fits this trade-off of being sufficiently far from  Z_s to allow the derivative of \phi to reach its asymptote and sufficiently close to prevent \phi from decreasing is estimated to 40% of Z_s-Z. This estimation, quite disputable, is based on temperature profile observed in diffusion flames at low Da when the diffusive effect strongly smoothes it.

Ignition Da curve in a cold diffusion layer.

The end equation has the form:

\hat{\phi}e^{-\hat{\phi}}=-{\mathcal D} \overbrace{(Z_s-1)\left (\frac{Z^2-Z_s^2}{2}-\frac{Z^3-Z_s^3}{3}\right )}^{a}

This makes \mathcal D a function of \hat{\phi}, maximum temperature in the ignition zone at stoichiometry, that has a maximum for \hat{\phi}=1. This gives the bottom tipping point of the S-curve presented above. For the academical case  Z_s = .5 and  n_F=n_O=1 , this gives  \mathcal D = 16. which is in a satisfactory agreement with the full numerical solution. It is interesting to note that the numerical solution follows the physics process of ignition of the diffusion layer by slowly increasing the Damk\"ohler number. Therefore, it cannot catch properly the turning point and the doubled-value relation. Instead, once the maximum Damk\"ohler number corresponding to ignition is reached, the reduced temperature soars sharply. Further increasing D in the forbidden range, no solution can be found. Physically, for larger values of the Damköhler number, the configuration jumps from this nearly frozen situation to an equilibrium situation with a well developed diffusion flame, see the S-Curve above. The simplified equation for ignition is thus unadapted and a specific solution for a developed diffusion flame is now to be presented.

Diffusion Burning Process

From Sec. Simplified Solution of Diffusion Flame, the temperature jump at stoichiometry is identified as  Q Y_{F,o} Z_s/Cp .

Although chemistry in a diffusion flame is usually very fast, one expects a very residual amount of reactants in the reaction zone laying along the stoichiometric line. It thus means that equilibrium is not reached in practice and that the temperature is somewhat lower than the equilibrium temperature. Thermodynamics is moderated by kinetics effects or, in other words, the curves T(Z) and Y(Z) depend on the dissipation rate and are not simply made of straight lines meeting at stoichiometry as in the simplified Burke-Schumann solution for diffusion flames.

To a certain extend, the flamelet equation Sec. Flamelet Equation may be used to present a solution for the non-premixed flame when the flame is already ignited, and the chemistry strength controlled by diffusion (we have stressed the author that an interesting aspect of diffusion flames is the ability to control the reaction rate through the level of mixing in the neighbourhood of the flame). As already emphasized, due to a strong non-linearity in temperature, the structure of a burning flame can never be far from the Burke-Schumann simplified solution.

Non-dimensionalized, the temperature solution in the flamelet may be written as:

\theta=1+(Z-Z_s)\Delta T_o + \frac{Z_s - Z}{1-Z_s} -\varepsilon\Gamma

in the rich side of the mixing layer, and:

\theta=1+(Z-Z_s)\Delta T_o + \frac{Z-Z_s}{Z_s} -\varepsilon\Gamma

in the lean side. On recognizes the structure of the Burke-Schumann solution (Sec. Simplified Diffusion Flame Solution). \Delta T_o is the non-dimensionalized temperature difference between the reactants feed streams:

\Delta T_o = \frac{T_{F,o}-T_{O,o}}{Q Y_{F,o} Z_s/C_p}

Those solutions are amended by the term \varepsilon\Gamma with \Gamma\sim{\mathcal O}(1) in the reaction zone and vanishing on the sides. This amendment corresponds to the departure from the equilibrium solution due to kinetics effect. The quantity \varepsilon is taken small for the physical reasons explained above.

It is interesting to introduce \gamma=2(Z_s-1)(Z_s\Delta T_o+1)+1, combining the temperature difference between the reactant feeding streams and the stoichiometry of the combustion) as a synthetic parameter to inform about the geometry of the flame structure in the (Z,T) diagram (Burke-Schumann):

\theta=1-\frac{(Z-Z_s)}{2Z_s(1-Z_s)}(\gamma+1) -\varepsilon\Gamma

in the rich side of the mixing layer, and:

\theta=1-\frac{(Z-Z_s)}{2Z_s(1-Z_s)}(\gamma-1) -\varepsilon\Gamma

in the lean side. For \gamma=0 the slope in temperatue vs Z is the same on both sides of the flame. For  \gamma > 0 the rich side slope is the sharpest, and for  \gamma <0 this is the lean side. An interesting case is when |\gamma|>1: one for the feeding stream enters the flame zone with a temperature higher than the flame temperature. In that case, whatever the dissipation rate, the flame is stretched with a hot stream that sustains its temperature such that it is never extinguished.

Structure of a diffusion flame in temperature and species mass fractions vs mixture fraction burning under the diffusion-controlled regime.

The departure in temperature from equilibrium in the reaction zone is reflected into the species concentrations, as it corresponds to still unburned reactants. With the help of the above Conservation Laws one finds:

  • Lean side:  Y_F=\varepsilon\Gamma Y_{F,o}Z_s,  Y_O=Y_{O,o}(1-Z)-s(Y_{F,o}Z-\varepsilon\Gamma Y_{F,o}Z_s)
  • Rich side:  Y_F=Y_{F,o}Z+\frac{1}{s}(s Y_{F,o}Z_s\varepsilon\Gamma-Y_{O,o}(1-Z)) ,  Y_O=sY_{F,o}Z_s \varepsilon\Gamma

The picture beside illustrates the structure of a diffusion flame in its `normal' regime.

Recasting the Flamelet Equation above in terms of the reduced temperature, we find:

\chi_s\frac{\partial^2 \theta}{\partial Z^2} + \dot\omega\frac{\nu_F\bar M_F}{\rho\rho_s^2 Y_{F,o} Z_s} = 0

To describe the structure of the diffusion flame, the reduced mixture fraction is set:

 \xi = \frac{Z-Z_s}{2Z_s(1-Z_s)\varepsilon}

The utility of the reduced mixture fraction is to focus on the reaction zone. This reaction zone is supposed located on the stoichiometric line (this is why the reduced mixture fraction is centred on Z_s) and to be very thin (reason of the introduction of the magnifying factor \varepsilon). One may convince oneself by considering in the above equation linking \theta to  Z that one must stay close to the stoichiometric line if one wants the non-equilibrium effect to be sensible compared to the adiabatic drop of the temperature. The equation is also written in terms of \Gamma, the non-equilibrium departure rescaled by  \varepsilon and the source term as presented in Main Specificities of Combustion Chemistry is explicated:

\left \{
\frac{\chi_s}{4Z_s^2(1-Z_s)^2\varepsilon}\frac{\partial^2 \Gamma}{\partial \xi^2} & = & 
A \frac{\rho\nu_F}{\rho_s^2} \left ( \frac{Y_{F,o}}{\bar M_F} \right )^{n_F-1} \left ( \frac{Y_{O,o}}{\bar M_O} \right )^{n_O} 
e^{-\frac{\beta}{\alpha}} Z_s^{n_F+n_O-1} \varepsilon^{n_F+n_O} \Gamma^{n_F} \Phi^{n_O} (\Gamma-2\xi)^{n_O} 
e^{-\beta\varepsilon (\xi(\gamma-1)+\Gamma)} & Z<Z_s \\

\frac{\chi_s}{4Z_s^2(1-Z_s)^2\varepsilon}\frac{\partial^2 \Gamma}{\partial \xi^2} & = & 
A \frac{\rho\nu_F}{\rho_s^2} \left ( \frac{Y_{F,o}}{\bar M_F} \right )^{n_F-1} \left ( \frac{Y_{O,o}}{\bar M_O} \right )^{n_O} 
e^{-\frac{\beta}{\alpha}} Z_s^{n_F+n_O-1} \varepsilon^{n_F+n_O} (\Gamma+2\xi)^{n_F} \Phi^{n_O} \Gamma^{n_O} 
e^{-\beta\varepsilon (\xi(\gamma+1)+\Gamma)}  & Z>Z_s 

It is seen that one can form a Damköhler number as (see Sec. Damköhler Number for the general meaning of the Damköhler):

Da=\frac{A\rho\nu_F}{\chi_s\rho_s^2}\left(\frac{Y_{F,o}}{\bar M_F}\right)^{n_F-1}\left(\frac{Y_{O,o}}{\bar M_O}\right )^{n_O}e^{-\frac{\beta}{\alpha}}\Phi^{n_O} 4(1-Z_s)^2 Z_s^{n_F+n_O+1}

and that this number must scale with \varepsilon^{n_F+nu_O+1}. So, to have a thin reaction zone, the Damköhler number must be high. Physically, it means that when fuel and oxidizer meet close to the stoichiometric line they combine quickly such that the reaction can be completed in a very narrow region.

In the diffusion burning process, Da is sufficiently high to have \beta\varepsilon << 1 such that the departure from equilibirum for the temperature has no impact in the exponential term at the end of the RHS of the equations above. This exponential term approximates to unity (in the flame zone, so for small  \xi ). Physically, this is the Burke-Schumann solution when the chemistry is sufficiently fast such that the kinetics limitation (modelled by this exponential term in the reaction rate) is inexistent. This gives more details on the structure of the diffusion flame for the simplified solution presented above Sec. Simplified Diffusion Flame Solution. In particular, the relationship between Da and \varepsilon allows an estimation of the quantities in the flame zone.

As for the boundary conditions, in the diffusion regime, the chemistry is fast enough to deplete the reactants diffusing into the reaction zone. Then, no presence of fuel on the oxidizer side and vice versa is expected. In a concise manner, we must have:

\lim_{\xi \rightarrow \infty} \Gamma = \lim_{\xi \rightarrow \infty} \Gamma\pm 2\xi = 0

When, for a given Damköhler number, the Zeldovitch parameter is too high such that its product with  \varepsilon is no more negligible, the combustion regime enters an extinction process. Physically, the small departure from equilibrium is sufficient to impact the RHS exponential of the temperature equations above such that the non-linear temperature dependence is sensitive.  Da \varepsilon^{n_F+n_O+1} \sim {\mathcal O}(1) and  \beta \varepsilon \sim {\mathcal O}(1) lead to the construction of a reduced Damköhler number \delta = Da /\beta^{n_F+n_O+1}. This reduced Damköhler number is expected to be around unity in the extinction regime, meaning that the Da is still high enough to sustain combustion but that the non-linearity of combustion due to a high value of \beta makes it very dependent to the temperature such that a weakening Damköhler, further enhancing the departure (temperature drop) from equilibrium, will quench the flame.

The fact that the chemistry is not fast enough to `eat' the reactants entering the flame and turn them into full equilibrium is synonymous to say that, at extinction, we expect a leakage of one or both reactants across the flame such that a concentration of the order of \varepsilon may be found on the oxidizer side for the fuel and/or vice versa. This quantity is a function of \delta. When the leakage is too intense, the combustion is quenched. There is thus a quenching  \delta=\delta_e .

Depending on the value of \gamma, as said above, the temperature slope on either side of the flame may be different. When reaching extinction, it is expected that the leakage happens for the reactant coming from the side with the weaker temperature slope through the side with the stronger temperature slope. The maximum temperature zone is lower than the adiabatic flame temperature, as previously mentioned, and also shifted towards the side of the flame with the weaker temperature slope. Hence, depending of the value of  \gamma the flame at extinction is described by one of the two equations above. Fortunately, due to the symmetrical role of  \gamma , a synthetic correlation for the Da at extinction can be obtained from the numerical solution of the above set of equations (with  n_i = 1):

\delta_e=e(1-|\gamma| - (1-|\gamma|)^2 + .26(1-|\gamma|)^3+.055(1-|\gamma|)^4)

This correlation, together with its equation now in a generic form for each side of the flame:

\Gamma_{\xi\xi}=\delta\Gamma(\Gamma+2|\xi|)e^{-\Gamma}e^{-|\xi|(1-|\gamma|)} \qquad sign(\xi)=-sign(\gamma)

has an especial place in combustion theory and modelling. It has been first proposed in [3] in 1974 and is since used as a canonical equation for problems that may be reduced to. It is also at the heart of understanding of non-premixed combustion modelling such as the flamelet model that will be introduced later.

It is important to notice that due to the leakage of reactants at extinction, the boundary conditions must be weaker than for the diffusion controlled regime:

\lim_{\xi\rightarrow \infty} \Gamma_{\xi} = 0, 2\times (sign(\xi))

0 on the side where the flame lies at extinction (depending on \gamma as seen above and where the equation is written) and 2 for the other.

The idea is that, compared to the diffusion-controlled regime, the shape of the Burke-Schumann asymptotes is conserved but with a shift corresponding to the unburned leakage and the associated drop in sensible enthalpy. Then, compared to the diffusion-controlled regime the asymptotic boundary gradients are the same, but not the actual values.

In the case the flame is highly asymmetrical (|\gamma| \rightarrow 1 ) either because one of the reactant streams is preheated or because the stoichiometry is far from  Z_s = .5 , an analytical approach may be found for \delta_e. This approach allows us to illustrate qualitatively the mechanism of extinction described only with words so far. Let take the case \gamma < 0 , meaning that the slopes (temperature and species) are sharper on the oxidizer side. Then, the maximum temperature zone is shifted on the fuel side (it may be seen for instance even for a moderate departure from the symmetry as in the figure above) and a leakage of fuel is expected on the oxidizer side. In a first approximation, the leakage amount is approximated as  b/(1+\gamma), with  b\;{\mathcal O}(1). So, for highly unsymmetrical flames (\gamma\rightarrow -1), this leakage is considered important before extinction may occur. One may convince oneself by considering that the asymmetry is due to a preheating of the fuel stream. When \gamma\rightarrow -1, it means that the preheating is such that the fuel stream temperature is close to the adiabatic temperature of the system. One may imagine that the flame is fairly robust and may bear a high level of leakage before being quenched. Actually, when |\gamma|\approx 1-\varepsilon, by looking at the equation linking \theta, Z, and \Gamma, it is seen that an especial scaling happens between the non-equilibrium departure from the B-S solution and the drop of the B-S temperature away from stoichiometry. The non-premixed flame has thus the freedom to switch to a `premixed' regime on the side with the weakest temperature slope of the B-S solution and at a location corresponding to a strong chemistry. Hence, the fuel leakage on the oxidizer side is quantified as  \lim_{\xi\rightarrow -\infty} (\Gamma+2\xi) = b/(1+\gamma) (remember the maximum temperature location is expected on the fuel side and the corresponding equation is used for the diffusion flame as developed above). The new asymptote intersects the \xi axis at the location  b/2/(1+\gamma) , that is taken as the flame position at extinction (in other words, the dissipation has the effect of shifting with respect to the Burke-Schumann infinitely fast chemistry / equilibrium).

We use this shift to make a change of frame of reference, which is now centred around the flame position at extinction as just defined: \zeta = \xi - b/2/(1+\gamma) . The equation for the flame above becomes:

\Gamma_{\zeta\zeta} = \delta_e e^{-\frac{b}{2}} \left (\frac{b}{1+\gamma}\right)^{n_F} e^{-\Gamma-\zeta(1+\gamma)} \Gamma^{n_O}\left ( 1 + \frac{1+\gamma}{b} (\Gamma +2\zeta) \right )^{n_F}

With this shift, stronger boundary conditions may be devised:

  •  \zeta \rightarrow -\infty , we are on the oxidizer side of the flame in the new frame of reference and the fuel concentration must equals its leakage level, then \Gamma+2\zeta\rightarrow 0,
  •  \zeta \rightarrow +\infty , given the fact that we are in the limit case \gamma\rightarrow -1, it is expected that the oxidizer is consumed instead of leaking through the fuel side:  \Gamma\; ; \; \Gamma_{\zeta}\rightarrow 0.

This equation may be simplified to the first order in  1+\gamma:

\Gamma_{\zeta\zeta} = \delta_e e^{-\frac{b}{2}}\left ( \frac{b}{1+\gamma} \right)^{n_F} e^{-\Gamma} \Gamma^{n_O}

A first integral gives:

\left [\frac{\Gamma_{\zeta}}{2}\right ]_{-\infty}^{+\infty} = \delta_e e^{-\frac{b}{2}} \left [\frac{b}{1+\gamma}\right ]^{n_F} \int_{\Gamma(-\infty)}^{\Gamma(+\infty)}d\Gamma  e^{-\Gamma} \Gamma^{n_O}

With the boundary conditions:

  • \Gamma,\;\Gamma_{\zeta}\stackrel{\rightarrow}{+\infty}0
  • \Gamma \stackrel{\rightarrow}{-\infty} +\infty
  • \Gamma_{\zeta}\stackrel{\rightarrow}{-\infty} -2

one obtains:

2=\delta_e e^{-\frac{b}{2}}\left (\frac{b}{1+\gamma} \right )^{n_F} n_O!

The RHS is a measure of the integrated reaction rate across the flame, which is a function of the leakage term b. One retrieves some of the arguments developed previously. This RHS is positive and increases in value with b up to a maximum corresponding to  b=2n_F. Then it vanishes because the dissipation imposed to the flame is too strong for the chemistry. At this maximum, the minimum Damköhler that the flame can deal with is found, which is identified as the Damköhler number at extinction:

\delta_e = \frac{2 e^{n_F}}{n_O!\left (\frac{2 n_F}{1+\gamma}\right )^{n_F} }

For n_F=n_O=1 the first term of the correlation from [2] as provided above is recovered.

Effect of the Damköhler number on the temperature profile.

The figure beside illustrates the effect of the Damköhler number on the temperature profile across the flame (i.e. vs the mixture fraction). Those solutions have been obtained from the flamelet equation presented above with reactant orders and free stream mass fractions equal to unity. The stoichiometry is  Z_s = .5 and the reactant free stream temperatures are the same. Two values of \delta are chosen: 5, which is far from extinction and 2, which is relatively closer to extinction for this set of combustion parameters as seen below. It is retrieved that reducing the ratio chemistry strength / dissipation strength by reducing \delta leads to a larger departure from equilibrium in the flame zone, around  Z_s . The temperature profile is smoothed out in the flame.

Effect of the Damköhler number on the flamelet maximum temperature.

The figure on the right tracks the change in the maximal temperature in the flame zone for decreasing \delta . This figure illustrates well the non-linearity of combustion vs temperature: once the extinction value is passed, the temperature drop is as sharp as a cliff, meaning a sudden end of chemistry. Physically, this is what is experienced when blowing to much on a wood fire: as explained earlier, blowing reduces the Damköhler number. Below a threshold value, one experiences a sudden blow-off / extinction of the fire.

Reduced Damköhler number at extinction vs diffusion flame asymmetry factor.
The Damköhler number at extinction is given in function of \gamma in the figure beside. It is retrieved that when the asymmetry is sensitive, there is considerable improvement of the flame robustness as predicted.

Remark: To develop this extinction analysis, we have assumed that the Zeldovitch parameter was asymptotically high. This is what is commonly named Asymptotic Energy Analysis (AEA). It uses the fact that combustion, under normal conditions, is highly temperature dependent by assuming that the activation temperature of the combustion chemistry tends to infinity. For the numerical simulation of the same extinction solutions, we use finite values of the Zeldovitch parameter. When the combustion parameters are scaled to values representative of actual chemistry, it is seen that the discrepancies with the AEA may be substantial. In order to get close to the analytical development, the numerical simulations must use unreasonably high values of the Zeldovitch parameter. It must also be noticed that pushing up to the second order the perturbation analysis will bring a correction that is even worse since it will predict an infinitely robust flame at \gamma=0.

The Premixed Regime

Sketch of a premixed flame

In contrast to the non-premixed regime above, the reactants are here well mixed before entering the combustion chamber. Chemical reaction can occur everywhere and this flame can propagate upstream into the feeding system as a subsonic (deflagration regime) chemical wave. This presents lots of safety issues. Some situations prevent them: (i) the mixture is made too rich (lot of fuel compared to oxidizer) or too lean (too much oxidizer) such that the flame is close to its flammability limits (it cannot easily propagate); (ii) the feeding system and regions where the flame is not wanted are designed such that they impose strong heat loss to the flame in order to quench it. For a given thermodynamical state of the mixture (composition, temperature, pressure), the flame has its own dynamics (speed, heat release, etc) on which there is few control: the wave exchanges mass and energy through diffusion process in the fresh gases. On the other hand, those well defined quantities are convenient to describe the flame characteristics. The mechanism of spontaneous propagation towards fresh gas through the thermal transfer from the combustion zone to the immediate slice of fresh gas such that the ignition temperature is eventually reached for this latter was highlighted as early as by the end of the 19th century by Mallard and LeChatelier. The reason the chemical wave is contained in a narrow region of reaction propagating upstream is the consequence of the discussion on the non-linearity of the combustion with temperature in the Sec. Fundamental Aspects. It is of interest to compare the orders of magnitude of the temperature dependent term  \exp{(-\beta(1-\theta)/(1-\alpha(1-\theta)))} of the reaction source upstream in the fresh gas (\theta\rightarrow 0) and in the reaction zone close to equilibrium temperature (\theta\rightarrow 1 ) for the set of representative values:  \beta = 10 and \alpha=0.9. It is found that the reaction is about 10^{43} times slower in the fresh gas than close to the burned gas. It is known that the chemical time scale is about 0.1 ms in the reaction zone of a typical flame, then the typical reaction time in the fresh gas in normal conditions is about  10^{39} s . To be compared with the order of magnitude of the estimated Universe age:  1 0^{17} s . Non-negligible chemistry is only confined in a thin reaction zone stuck to the hot burned gas at equilibrium temperature. In this zone, the #Damköhler number is high, in contrast to in the fresh mixture. It is natural and convenient to consider that the reaction rate is strictly zero everywhere except in this small reaction zone (one recovers the Dirac-like shape of the reaction profile, provided that one can see the upstream flow as a region of increasing temperature towards the combustion zone and the downstream flow as in fully equilibrium).

As the premixed flame is a reaction wave propagating from burned to fresh gases, the basic parameter is known to be the progress variable. In the fresh gas, the progress variable is conventionally put to zero. In the burned gas, it equals unity. Across the flame, the intermediate values describe the progress of the reaction to turn into burned gas the fresh gas penetrating the flame sheet. A progress variable can be set with the help of any quantity like temperature, reactant mass fraction, provided it is bounded by a single value in the burned gas and another one in the fresh gas. The progress variable is usually named c, in usual notations:


It is seen that c is a normalization of a scalar quantity. As mentioned above, the scalar transport equations are assumed linear such that the transport equation for c can be obtained directly. Actually, the transport equation for T (Sec. Transport Equations) is linear if constant heat capacity is further assumed (combustion of hydrocarbon in air implies a large excess of nitrogen whose heat capacity is only slightly varying) and the progress variable equation is directly obtained (here for a default of fuel - lean combustion):

\frac{D\rho c}{Dt}=\nabla\cdot\rho D\vec\nabla c + \frac{\nu_F\bar M_F}{Y_{F,u}}\dot\omega \qquad \rho D = \frac{\lambda}{C_p}

The fact that the default or excess of fuel has been discussed above leads to the introduction of another quantity: the equivalence ratio. The equivalence ratio, usually noted \Phi, is the ratio of two ratios. The first one is the ratio of the mass of fuel with the mass of oxidizer in the mixture. The second one is the same ratio for a mixture at stoichiometry. Hence, when the equivalence ratio equals unity, the mixture is at stoichiometry. If it is greater than unity, the mixture is named rich as there is an excess of fuel. In contrast, when it is smaller than unity the mixture is named lean. The equivalence ratio presented here for premixed flames has little connection with the equivalence ratio introduced earlier regarding the non-premixed regime. Basically, the equivalence ratio as defined for non-premixed flames gives the equivalence ratio of a premixed mixture with the same mass of fuel and oxidizer. Moreover, the equivalence ratio as defined for a premixed mixture can be obtained based on the mixture fraction (it is thus the local equivalence ratio at a point in the non-homogeneous mixture described by the mixture fraction). From the definitions given above:


Premixed Flame Péclet Number

Earlier in this section, it has been said that a premixed flame posses its own dynamics, as a free propagating surface, and has thus characteristic quantities. For this reason, a Péclet number may be defined, based on these quantities. The Péclet number has the same structure as the Reynolds number but the dynamical viscosity is replaced by the ratio of the thermal conductivity and the heat capacity of the mixture. The thickness \delta_L of a premixed flame is essentially thermal. It means that it corresponds to the distance of the temperature rise between fresh and burned gases. This thickness is below the millimetre for conventional flames. The width of the reaction zone inside this flame is even smaller, by about one order of magnitude. This reaction zone is stuck to the hot side of the flame due to the high thermal dependency of the combustion reactions, as seen above. Hence, the flame region is essentially governed by a convection-diffusion process, the source term being negligible in most of it.

It is convenient to write the progress variable transport equation in a steady-state framework. The quantities at flame temperature (_f) are used to non-dimensionalize the equation:

 \overbrace{(\rho ||\vec S_L||)}^{||\vec M||}\frac{\vec S_L}{S_L}\vec\nabla c = (\rho D)_f \nabla\cdot (\rho D)^* \vec\nabla c

Note that the source term is neglected, consistently with what has been said above. This convection-diffusion equation makes appear a first approximation of a flame Péclet number:

 Pe_f = \frac{||\vec M|| \delta_L}{(\rho D)_f} \approx 1

From the Péclet number, it is possible to obtain an expression for the flame velocity (remembering that  \delta_L/S_{L,f} \approx \tau_c, vid. inf. Sec. Three Turbulent-Flame Interaction Regimes):

 S_{L,f}^2\approx \frac{(\rho D)_f}{\rho_f \tau_c}

For typical hydrocarbon flames, the speed is some tens of centimetres per second and the diffusivity is some  10^{-5} square metres per second. The chemical time in the reaction zone of about one tenth of a millisecond is recovered.

We anticipate here a result that will be obtained in the detailed analysis of the flame speed (vid. inf.), that is: the laminar flame speed is proportional to the square root of the diffusion coefficient and the reaction rate. This result is of interest to anticipate the trends of the premixed flame dynamics.

Details of the Premixed Unstrained Planar Flame

A plane combustion wave propagating in a homogeneous fresh mixture is the reference case to describe the premixed regime. At constant speed, it is convenient to see the flame at rest with a flowing upstream mixture. This is actually the way propagating flames are usually stabilized. In the frame of description, the physics is 1-D, steady with a uniform (the flame is said unstrained) mass flowing across the system. Two types of equations are thus sufficient to describe the problem, the temperature transport equation and the species transport equations, as in Sec. Transport Equations. The transport coefficients will be chosen as equal:  \rho D_i = \lambda / C_p (unity Lewis numbers). Suppose the 1-D domain is described thanks to a conventional (Ox) axis with a flame propagating towards negative x (this is the conventional usage), the boundary conditions are:

  • in the frozen mixture:
    • Y_i \rightarrow Y_{i,u} \qquad ; \qquad x\rightarrow-\infty
    • T \rightarrow T_u \qquad ; \qquad x\rightarrow-\infty
  • in the burned gas region supposed at equilibrium:
    • Y_i \rightarrow Y_{i,b} \qquad ; \qquad x\rightarrow+\infty
    • T \rightarrow T_b \qquad ; \qquad x\rightarrow+\infty
Y_{i,b} and T_b are obtained from Sec. Conservation Laws.

The quantities that have been mentioned just above (scalar and temperature profiles, mass flow rate through the system) are the solution to be sought.

According to the discussions above, the temperature transport equation in its full normalized form may be written as (lean /stoichiometric case):

 ||\vec M|| \frac{\partial \theta}{\partial x} = \frac{\partial \  }{\partial x}\frac{\lambda}{Cp} \frac{\partial \theta}{\partial x} + \frac{\nu_F \bar M_F B}{Y_{F,s}} \prod_{i=O,F}(Y_{i,u}^*-\theta\frac{Y_{F,s}}{Y_{i,s}}s^{\delta_{i,O}} )^{n_i} \exp{-\beta\frac{1-\theta}{1-\alpha(1-\theta)}}

This equation is further simplified by the variable change d\xi=||\vec M||/(\lambda/Cp)dx:

\frac{\partial \theta}{\partial \xi}=\frac{\partial^2 \theta}{\partial \xi^2} + \overbrace{\frac{\lambda/Cp}{||\vec M||^2}\frac{\nu_F \bar M_F B}{Y_{F,s}}}^{\Lambda} \prod_{i=O,F}(Y_{i,u}^*-\theta)^{n_i} \exp{-\beta\frac{1-\theta}{1-\alpha(1-\theta)}}

Although somewhat out of scope, the existence and unicity of the solution of this type of equation are usually demonstrated with the help of the Schauder Theorem and Maximum Principle. From the point of view of physicists and engineers, the solution that is found analytically is de facto considered as the unique solution of the equation.

Scenarii of Combustion Process in the Phase Portrait

In the frame moving with the flame, both phase variables are the reduced temperature and its gradient. To ease the reading with usual notations, it is written:  X_1 = \theta \quad ; \quad X_2 = \partial \theta / \partial \xi = \dot \theta . The system arises:

\dot X_1 & = & X_2 \\
\dot X_2 & = & X_2 - \varpi(X_1) 

with  \varpi being the full source term in the above equation.

In the frame moving with the flame, two singular nodes are found in the frozen flow  (X_1,X_2) = (0,0) and the equilibrium region  (\theta_b,0), i.e when  \varpi(X_1) vanishes.  x_1,x_2 are defined as small departures from the singular nodes such that the linearized system in their neighbourhood is:

\dot x_1 & = & l_{1,1} x_1 + l_{1,2} x_2 \\
\dot x_2 & = & l_{2,1} x_1 + l_{2,2} x_2  

provided:  l_{1,1} = 0 \quad l_{1,2} = 1 \quad l_{2,1} = -\varpi'_{X_1 = 0,\theta_b} \quad l_{2,2} = 1 . The characteristic polynom is, in usual notations:  s^2 - s + \varpi' such that the eigenvalues are:


A priori, those eigenvalues may be (i) real distinct, (ii) real identical, or (iii) conjugated complex. In the first case, the orbits in the phase diagram are organized, in the immediate neighbourhood of the singular node, with respect to the eigenvectors directions associated to the eigenvalues. The following task is to identify the nature of those eigenvalues and of the corresponding nodes. Because  0 < X_1 < \theta_b is bounded, complex eigenvalues are excluded as they would lead to a spiral node. This remark is important for the node on the cold side as it imposes a bound:

 \varpi'_{X_1=0} \le \frac{1}{4}

As the mass flow rate through the flame is included into \varpi, it imposes a minimum value on the flame speed to tackle with the cold boundary difficulty (rise of the chemical rate in the frozen flow). In this condition, it is an unstable node (improper in case of equality). On the other hand, because  \varpi'|_{X_1=\theta_b} is not positive, the node on the hot side is found as a saddle point. The overall scenario of combustion within the flame is thus an orbit leaving the cold node to join the hot node by branching on a trajectory compatible with the negative eigenvalue of the saddle.

Sketch of orbits for a combustion process across a premixed front. Dashed lines represent forbidden orbits from the physics. The red line describes the orbit expected in an idealized combustion process.

It must be noted that the associated eigenvectors are of the form:

\left (
t \\
t s^{\pm}  
\right ); \quad t\in \Re^*

that is, on the cold node, a positive departure on X_1 following any of the two eigendirections, leads to a consistent creation of positive temperature gradient, while on the hot node, only the stable direction will allow a consistent creation of a positive temperature gradient for any departure of the temperature towards region where it is inferior to \theta_b . Another remark is the structure of the eigendirections. The leaving directions on the cold node have a slope larger than the one of  \varpi while the stable direction of the hot node has a slope smaller than the one of  \varpi . It means that there is some point where the orbit must cross the profile of the chemical term versus temperature. For that temperature, the gradient equals the reaction rate through construction of the phase space. When looking back to the equation of the premixed flame, it happens in a region of inflexion for the temperature (the second order derivative must vanish). Furthermore, at this intersection, the orbit is horizontal (if the frame of reference for  (X_1, X_2) is Cartesian) due to the shape of the premixed flame equation above that can be recast into  X_{2,X_1}' = (X_2-\varpi)/X_2 . Close to the cold node, the orbits have a shape of parabola whose axis is the direction with the largest eigenvalue magnitude. Close to the hot node, the orbits have an hyperboloid shape with asymptots as the eigendirections. Now the ingredients are here to draw a sketchy scenario of the combustion in a premixed flame. It will be superimposed on the reaction rate graph studied in Sec. Fundamentals. Some typical orbits from the above analysis are drawn in the figure on the right. The basic geometrical arguments developed are reproduced. In particular, the dashed lines represent forbidden orbits by the physics (boundedness of X_1, irreversibility). Orbits must be travelled from left to right, corresponding to increasing free parameter  \xi . It is of integral importance to remember that, in combustion in conventional conditions, the source term is highly non-linear. Therefore, it is localized in a very thin sheet and, upstream,  \varpi and  \varpi' \rightarrow 0 . Only the most unstable eigendirection of the cold node is compatible as an orbit. The other trajectories, being paraboloidal, are tangent to the other direction that is flat at the limit. It physically means an elevation of temperature in bulk, that is contradictory to what is expected from a highly non-linear combustion term. Hence, the additional orbit in red is the one expected in idealized combustion conventional conditions.

Sketch of orbits for a combustion process across a premixed front, DNS simulations for the usual combinations of \beta and  \alpha . Note that the case with \alpha = 0 does not mean an athermal reaction but the mathematical simplification of the denominator of the exponential argument of the source term.

The following picture is the result of actual computations of the above 1-D flame equation (stoichiometry and  n_i = 1 ) with the help of a high-order (6) code. The already presented combinations of  \beta and  \alpha are used and the orbits are retrieved. For most of the cases, as predicted above, the system selects a solution leaving the cold node with the most unstable direction (identifiable with its slope close to unity for vanishing \varpi'_{X_1=0}). There is also an additional curve, for  \beta = 10 with the denominator of the exponential argument suppressed. This curve is remarkable as the solution selected by the system leaves the cold node in a manner fully controlled by the  \varpi'_{X_1=0} . This is a very singular solution, not expected in combustion in conventional conditions, as explained above. The purpose of this remark is to question the well-posedness of considering simplifying the exponential argument for \beta "sufficiently high", as it is usually proposed for this type of modelling. As observed, the dynamical system analysis demonstrates a switch in the nature of the solution selected. At the physical level, when the orbit follows the most unstable direction with a slope close to unity, it means that  X_2 "follows"  X_1 , which is a signature of a diffusion process. In other words, the preheating mechanism of a premixed flame propagation as proposed for more than one century is in work. On the other hand, when  X_2 is dependent on the evolution of  \varpi', it shows that "cold" chemistry drives the solution in the frozen flow and not the acknowledged mechanism of deflagration.

Flame Solution

As already mentioned, the flame system may be split into three zones. Upstream, the conventional mechanism of deflagration is supported by diffusion of heat. Downstream, the mixture is at equilibrium after combustion. In between, there exists the reaction layer. For large \beta the reaction layer is very thin such that it can be seen as a discontinuity between the fresh and burned gases. This is this difference in scales that introduces the use of the asymptotic method to resolve some of the flame characteristics, such as speed, time, heat region thickness, or reaction zone thickness.

The domain is partitioned, according to this zoning defined by the scales driving the physics with an outer domain, driven by large scales and an inner domain refining the description within the discontinuity. If the discontinuity (flame reaction zone) is at  \xi=0 , everywhere but 0, the equation is simplified as:

\frac{\partial^2 \theta}{\partial \xi^2} = \frac{\partial\theta}{\partial \xi}
  • For  \xi>0 , it is expected that the mixture has reached equilibrium chemistry, such that:  \forall \xi > 0, \quad \theta=\theta_b \quad ; \quad \partial \theta / \partial \xi =0 .
  • For  \xi < 0 , this is the preheat zone and the solution is  \theta = \theta_b e^{\xi} with  \theta reaching \theta_b at the disconstinuity and vanishing, together with its gradient, far upstream in the frozen mixture.

The solution for the species and the value of  \theta_b are obtained from Conservation Equations above. The 'big picture' is thus an exponential variation in the thermal thickness matched with a plateau in the downstream region, the line of matching being the discontinuity (flame) that has no thickness at this scale of description.

To refine the analysis in the discontinuity region, a magnifying factor  \varepsilon is used to stretch the coordinates:  \xi = \varepsilon \Xi . The inner solution is thus a slowly-varying function of  \Xi . Hence, in this inner region, the equation for the premixed flame becomes:

\frac{1}{\varepsilon}\frac{\partial \theta}{\partial \Xi}=\frac{1}{\varepsilon^2}\frac{\partial^2 \theta}{\partial \Xi^2} +\varpi

In order to stretch and 'look inside' a discontinuity, \varepsilon is very small. It yields two remarks:

  1. convection is negligible compared to diffusion. The heat losses from the reaction zone are essentially diffusion driven.
  2. The reaction zone is governed by a diffusion-reaction budget and the reaction term \varpi must be strong to balance the intense heat loss due to the sharp diffusion (the zone is very thin, hence the gradients are sharp).

The mechanism is thus different from the outer region that was convection-diffusion driven.

Each quantity is developed in a series of  \varepsilon . At the leading order, for \theta, in the lean case, the conservation relations (Sec. Conservation Laws) yield:

 \theta = 1 - \varepsilon \Gamma - (1 - Y^*_{F,u})

where \Gamma is the first-order development of the departure of \theta from the maximum value due to the incomplete combustion, and 1-Y_{F,u}^*=1-\theta_b is the reduction of temperature for non-stoichiometric cases. Injected into the above equation:

 \frac{1}{\varepsilon}\frac{\partial^2 \Gamma}{\partial \Xi^2} = \Lambda (\varepsilon \Gamma)^{n_F} (Y_{O,u}^*-Y_{F,u}^*\frac{Y_{F,s}}{Y_{i,s}}s^{\delta_{i,O}} +\varepsilon\Gamma\frac{Y_{F,s}}{Y_{i,s}}s^{\delta_{i,O}} )^{n_O}\exp{-\beta\frac{1-Y_{F,u}^*+\varepsilon\Gamma}{1-\alpha(1-Y_{F,u}^*+\varepsilon\Gamma)}}

Although the full develoment is not achieved, a number of scaling may be highlighted:

  1. because the temperature cannot be much below unity, Y_{F,u}^* must be close to 1  {\mathcal O}(\varepsilon) . For clarity, it is not expanded in an \varepsilon series.
  2. The denominator of the exponential argument simplifies to unity for small \varepsilon.
  3. To get a finite rate in the reaction zone, \varepsilon scales with \beta^{-1}.

The burning rate eigenvalue, \Lambda, is naturally expanded as:  \Lambda = \varepsilon^{-n_O-n_F-1}(\Lambda_0 + {\mathcal O}(\varepsilon)) . The low-order equation to be solved is:

\frac{d^2 \Gamma}{d \Xi^2} = \frac{1}{2}\frac{d(\Gamma_{\Xi})^{'2}}{d\Gamma} = \Lambda_0\exp{-\beta (1-Y_{F,u}^*)} \Gamma^{n_F} (\frac{Y_{O,u}^*-Y_{F,u}^*\frac{Y_{F,s}}{Y_{i,s}}s^{\delta_{i,O}} }{\varepsilon}+\Gamma\frac{Y_{F,s}}{Y_{i,s}}s^{\delta_{i,O}} )^{n_O}\exp{-\Gamma}
 \frac{d\Gamma}{d\Xi}(-\infty)=-\theta_b, \qquad \frac{d\Gamma}{d\Xi}(\infty)=0

The boundary conditions are obtained from the matching of the outer solutions on the right and left sides of the flame as written above (the outer solutions are reached at infinity for a very small magnifying factor \varepsilon). Once integrated with respect to those boundary conditions, the burning-rate eigenvalue (from which  \dot M is extracted) is obtained as:

 \Lambda_0 = \Bigg( 2\int_0^{\infty}d\Gamma\; \Gamma^{n_F} (\beta (Y_{O,u}^*-Y_{F,u}^*\frac{Y_{F,s}}{Y_{i,s}}s^{\delta_{i,O}} )+\Gamma\frac{Y_{F,s}}{Y_{i,s}}s^{\delta_{i,O}} )^{n_O}\exp{-\beta (1-Y_{F,u}^*)}\exp{-\Gamma} \Bigg )^{-1}

The RHS integral is not developed for clarity but presents no peculiar difficulties.

The development has been carried out at the first order in  \varepsilon . As soon as a second order development is attempted, some expressions are no more analytically tractable. On the other hand, a second order development allows introduce the temperature-dependent trends of some terms in  \Lambda . Physical results are retrieved such as a slight decrease of the speed for a positive sensitivity of transport parameters to temperature around equilibrium conditions.

Unstrained planar premixed flame speed with respect to fuel mass fraction (lean case) for a single global irreversible Arrh\'enius term. Symbols are obtained from a high-order DNS code. Continuous line is the theory exposed here. The usual combinations of \beta and  \alpha are used. Very high values of  \beta (and thus negligible effect of  \alpha are also presented to show the problem of slow convergence for finite value of the Zeldovitch parameter. Not that \alpha = 0 does not mean an athermal reaction but the mathematical simplification of the denominator of the exponential argument of the source term.

The image beside illustrates the response of the premixed flame speed with respect to fuel concentration (chosen as the limiting component here; equivalent findings are obtained for default of oxidizer) for different values of the chemical parameters  \beta and \alpha (as noted above, the pre-exponential constant  A impacts only the overall magnitude). The theoretical expression (continuous line) is tested against a 1-D high accuracy code with the given chemistry implemented. It is seen that:

  • the Zeldovitch parameter drives the drop for non-stoichiomeric mixtures,
  • the drop is relatively well modelled by the theoretical expression, and
  • the absolute magnitude converges slowly towards the theoretical one when increasing \beta in the code (effect of the finiteness of  \beta in the real case).
Typical profiles in 1-D premixed flame at stoichiometry. Representative value of global chemistry parameters. Specie profiles are simply the complementary to the temperature profile for simple chemistry.

The picture exhibiting profiles for different Zeldovitch and heat release parameters shows the factual impact: upstream, the exact exponential profile is recovered and corresponds to the pre-heating region (thermal thickness). In the reaction zone (just upstream of the extremum), the departure from the exact solution is due to the kinetic effect. This kinetic effect is more pronounced when  \beta is lower because the lower the Zeldovitch parameter, the lower the temperature reaction zone can be without leading to extinction. The flame takes this opportunity to maximize its transfer in heat and reactant with the cold zone. This is the physical understanding of an inverse dependence of the maximum flame speed with  \beta .

Modification of the Flame Speed with Curvature

Flame seen as an interface between fresh and burned gases. Its curved profile towards the burned side increases the transfers with the fresh side.

When the thickness of the flame (\lambda/Cp)_f/||\vec M|| is considered small compared to inhomogeneities existing in the flow, the flame can be reduced to an interface between fresh and burned gases. This interface may not be strictly plan in the general case. For instance, when the interface is curved towards the burned gas, it offers a larger opportunity for transfer of mass and heat with the fresh gas. As the ability of the chemistry to burn the coming matter is limited, the flame has thus to reduce its displacement speed. The figure beside provides a 2-D sketch of this situation that happens in contorded turbulent fields. Curvature effect is thus an ingredient appearing in combustion models.

To mathematically give the expression showing that the curvature influences the flame speed (for small curvature), the non-dimensionalized temperature equation across a planar premixed flame above is slightly recast:

||\vec M^c||||\vec\nabla_{\xi}\theta|| - \varpi= \nabla_{\xi}\cdot\vec\nabla_{\xi}\theta=-\nabla_{\xi}\cdot(||\vec\nabla_{\xi}\theta||\vec n) = - \vec\nabla_{\xi}||\vec\nabla_{\xi}\theta||\cdot\vec n - ||\vec\nabla_{\xi}\theta||\nabla_{\xi}\cdot\vec n

In this expression, \dot M^c is the flame mass flow rate when perturbed by a curvature normalized by the reference one developed above.  \vec n is the normal to the iso-temperature pointing towards the fresh gas, what is named the normal to the flame. From the expression developed above and for a slightly perturbed flame, the gradients, production and diffusion terms are very close to the unperturbed flame. Only the last term, the normal divergence, does not exist for a non-perturbed flame. Hence, the above equation may be simplified against the one for unperturbed flame presented earlier:

 ||\vec M^c|| = 1 - \nabla_{\xi}\cdot \vec n

This is the divergence of the flame normal that contains the information upon geometrical perturbation (curvature) that impacts the speed.

Local geometrical approximation of a flame surface

The key is thus to get an idea of the geometrical significance of the normal divergence. The normal divergence theorem says that this is the sum of the principal curvatures of the flame surface at the location considered. To give a good mental picture, the simplest configuration without loss of generality is to consider the flame surface approached by an osculatory revolution ellipsoid, as in the figure beside. The location of the approximation is the intersection of the Ox axis and the flame surface in red (where the ellipsoid is tangent to the flame). At this location, the normal divergence is the sum of the curvature of the basic ellipse (before its rotation around Oy) at its minor extremum, and the curvature of the circle corresponding to the rotation of this point around Oy when the 3-D shape is formed. Giving a lecture on conical coordinate systems is beyond the purpose but the interested reader may want to follow the corresponding steps to check this result: (i) write the divergence of a vector in the ellipsoid coordinate system, (ii) consider that the vector is the normal, i.e. it has only one constant component perpendicular to the local ellipsoidal iso-coordinate, to simplify the divergence expression, (iii) split the resulting terms into two parts by identifying the second derivative of the basic shape 2-D ellipse as one principal curvature, and the inverse of the radius of the circle corresponding to the rotation of the ellipse around Oy at the point where the flame surface is approached, in the figure this is simply the minor axis of the ellipse.

The expression is readily written as:

 ||\vec M^c|| = 1 - (R_1^{-1}+R_2^{-1})

where  R_1 and  R_2 are the two radius of curvature (non-dimensionalized by the reference flame thickness) local to the surface. For instance, they are the small axis of the basic ellipse and its radius of curvature at its minor maximum when the surface is approached by the osculatory ellipsoid as above.

Natural Instabilities of Premixed Flames

Another aspect accounted for in turbulent combustion modelling of premixed flame is the creation of flame surface (corrugation) by hydrodynamics instabilities. These instabilities have been mentioned by Darrieus as early as 1938. The motivation of discussing about this is also related to the framework of description used that is widely employed for combustion models development. This framework is named the hydrodynamics limit, where the flame is isolated as a zero-thickness interface in the flow, and has been first introduced just above. In this framework, any diffusive and energetic aspects disappear and the set of equations is limited to two incompressible Euler systems. One system in the fresh gases (with a constant density of cold mixture). One system in the burned gases (with a constant density of mixture at equilibrium).

To understand the basic properties of a premixed flame leading to the birth of instabilities, it is first important to realize that a premixed flame, in the hydrodynamic limit, behaves as a dioptre with a refractive index in the burned gases larger than in the fresh gases. For a flow with an angle of attack on the flame (i.e. a flame not strictly 1-D perpendicular to the flow), the tangential component of the flow speed relative to the flame surface is conserved while the normal component is accelerated by a factor corresponding to the ratio of the density in order to conserve mass flow across the interface. Hence, the streamlines are pushed towards the normal to the flame when crossing, that is similar to rays of light entering a refracting medium.

Basic mechanism explaining the unstable nature of a planar premixed flame.

If one considers a region of a premixed flame that bumps a little bit towards the fresh gas (see figure beside -- the same approach is symmetrically true for the bump towards the burned gas), the local stream tube slightly opens on the bumpy interface before being refracted and coming back to its original section at constant mass flow rate. Hence, just in front of the bump, the gas velocity in this stream tube decreases and does not oppose the flame motion, allowing the bump to increase in magnitude. This is the fundamental mechanism of such instabilities.

The equations sets are in usual notations:

\left \{
U_x + V_y & = & 0 \\
U_t + UU_x + VU_y & = & \frac{P_x}{\rho} \\
V_t +UV_x+VV_y & = & \frac{P_y}{\rho} 

 x is the coordinate used above along the flame path and y is in the tangential plan. u and v are the respective velocities. Those equations may be normalized by steady-state reference quantities such as flame speed, flame length, gauge pressure and density with respect to fresh gas properties and are written on both sides of the interface.

These sets of equations must match at the interface (flame surface) through conservation of mass flow rate and momentum:

\left \{
\rho_u (\vec U_u -\vec U_{f})\cdot \vec n & = &  \rho_b (\vec U_b -\vec U_{f})\cdot \vec n\\
(\rho_u \vec U_u (\vec U_u -\vec U_{f})+P_u)\cdot \vec n & = &  (\rho_b \vec U_b (\vec U_b -\vec U_{f})+P_b)\cdot \vec n

The subscripts  u and  b stand for unburned and burned sides, respectively. The subscript  f points the interface. In these equations, the flame motion  U_f is reintroduced because we want to track the instability movement over the mean position of the flame. This instability motion is described by the equation  x=F(y,t) such that the motion of the interface normal to itself is obtained as:

 \vec U_f\cdot\vec n = - \frac{F_t}{\sqrt{1+F^2_y}}

with a normal to the front:

 \vec n = \left (
-\frac{1}{\sqrt{1+F^2_y}} \\
\right )

As the instability motion is considered at its birth, i.e. when it is still small, the quantities are linearized with  \varepsilon as a small parameter:

 F=\varepsilon f
 \vec U = \vec {\mathcal U} + \varepsilon\vec u
 P = {\mathcal P} + \varepsilon p

The Eulerian system is reduced to (at the first order in \varepsilon, the order zero is simplified for homogeneous flow):

\left \{
u_x + v_y & = & 0 \\
u_t + {\mathcal U}u_x & = & \frac{p_x}{\rho} \\
v_t +{\mathcal U}v_x & = & \frac{p_y}{\rho} 

The one-sided equations for the jump conditions become (terms up to the first order in \varepsilon are kept):

\left \{
\rho (\vec U-\vec U_f)\cdot\vec n \times \sqrt{1+F^2_y} & = & -\rho {\mathcal U}-\rho\varepsilon u+\rho\varepsilon f_t \\
(\rho \vec U (\vec U -\vec U_{f})+P)\cdot \vec n \times \sqrt{1+F^2_y} & = &
\left \{
-\rho {\mathcal U}^2 -2\rho{\mathcal U}\varepsilon u+\rho{\mathcal U}\varepsilon f_t - {\mathcal P} -\varepsilon p \\
-\rho {\mathcal U}\varepsilon v + {\mathcal P}\varepsilon f_y

The following jump conditions are obtained (by recalling that unburned steady-state quantities are chosen as references):

  • From the first equation at the order 0, \rho_b U_b = \rho_u U_u \Leftrightarrow U_b=\rho_u/\rho_b
  • From the first equation at the order 1,  \rho_u f_t - \rho_u u_u = \rho_b f_t -\rho_b u_b
  • From the second equation at the order 0,  -\rho_u {\mathcal U}_u^2 - {\mathcal P}_u = -\rho_b {\mathcal U}_b^2 - {\mathcal P}_b \Leftrightarrow {\mathcal P}_b = 1 - \rho_u/\rho_b
  • From the second equation at the order 1,  -2\rho_u{\mathcal U}_u u_u +\rho_u {\mathcal U}_u f_t -p_u = 
-2\rho_b{\mathcal U}_b u_b +\rho_b {\mathcal U}_b f_t -p_b \Leftrightarrow p_b-p_u = 2 u_u-2 u_b
  • From the third equation at order 1 (no order 0),  -\rho_u {\mathcal U}_u v_u +{\mathcal P}_u f_y = -\rho_b {\mathcal U}_b v_b +{\mathcal P}_b f_y \Leftrightarrow v_b-v_u = f_y(1-\rho_u/\rho_b)

The solution for the linearized, autonomous Euler system is:

\left (
u \\ 
v \\ 
\right )
\left (
\bar u  \\ 
\bar v  \\ 
\bar p 
\right )
\exp{(\sigma x)}\exp{(\alpha t -iky)}

where one recognizes an account for the perturbation of the field in the  x direction (\sigma), the wave number of the instability following  y ,  k and the growth rate with time  \alpha . The eigenvalues of the system are used to determine the  x dependence of the solution  \sigma = - \alpha/U,\; k,\; -k, with the positive ones applying on the fresh side and the negative ones on the burned size to have a vanishing perturbation far from the flame.

The eigenmodes give the flow pertubations on either side of the flame:

\left ( 
u \\ v \\ p
\right )
a\left ( 
1 \\ -i \\ -1-\frac{\alpha}{k}
\right )
e^{kx+\alpha t - iky}

\left ( 
u \\ v \\ p
\right )
b\left ( 
1 \\ i \\ -1+\rho_b\frac{\alpha}{k}
\right )
e^{-kx+\alpha t - iky}
c\left ( 
1 \\ i\rho_b\frac{\alpha}{k} \\ 0
\right )
e^{-\rho_b \alpha x+\alpha t - iky}

The jump conditions above applied to this perturbation field yield the following system ( f has also been put into the harmonic form  f = \bar f e^{\alpha t-iky} ):

\left \{ 
\alpha \bar f -a & = & \rho_b \alpha \bar f -\rho_b (b+c) \\
a(1-\frac{\alpha}{k}) & = & b(1+\rho_b\frac{\alpha}{k}) +2c \\
a + b+ c\frac{\rho_b \alpha}{k} & = & k\bar f (\frac{1}{\rho_b}-1)
\right .

Additionaly, from the definition of the flame path  F , we have the kinematic relation  u_u = a = \alpha \bar f .

The above set of equations forms a system of four unknowns  a,\; b,\; c,\; \alpha/k for a given (but unknown) shape information on the flame bump  k\bar f. Solving for the growth rate gives:

\left (\frac{\alpha}{k} - \frac{1}{\rho_b}\right )\left(\frac{\alpha}{k} + 2 + \rho_b\frac{\alpha}{k} -\left (\frac{1}{\rho_b}-1\right )\left (\frac{\alpha}{k}\right )^{-1} \right )=0

This third-degree polynom has three solutions, namely the dispersion relation found by Darrieus, a stable mode, and the trivial solution (with no physical meaning), respectively: 
\frac{\alpha}{k} = \left ( \frac{1}{\rho_b + 1}\left (\sqrt{1+\frac{1-\rho_b^2}{\rho_b}}-1 \right ),\qquad \frac{1}{\rho_b + 1}\left (-\sqrt{1+\frac{1-\rho_b^2}{\rho_b}}-1 \right ),\qquad \frac{1}{\rho_b}\right )

Remark 1 By considering the shape of the eigenvectors giving the perturbation, one recovers that, upstream of the flame, the flow is potential (the rotational of the velocity vector is null and we have chosen a constant density flow, that is barotropic and divergence free), while downstream, additionally to another perturbation of the potential type, one finds a vorticity mode (the mode with the eigenvalue -\rho_b\alpha . The drift of vorticity at the crossing of the flame front is a known property of flames.

Remark 2 In more conventional literature, the strange trivial solution for the instability growth rate is swept under the rug. Given the pedagogical nature of this electronic documentation, we can dig a little bit as the appearance of this trivial, fool solution is a good example of some modelling issue. It is important to realize that a mathematical model and the physics it describes belong to different realities. Hence, the mathematical model will generate all the solutions that the mathematics can reach in its own space. Some of them are still connected to the physics. Some others, like this trivial solution, belongs only to the mathematical solution space, an indirect way of pointing out a model limitation. The model under scope here is the hydrodynamic limit. In this model, the domain is divided into two subdomains, one upstream of the flame interface, one downstream, that are put in relation with each other through a limited number of jump conditions. In reality, the physics connects these both domains much more tightly. The curious reader will have observed that the trivial solution makes the equation system for  a,\; b,\; c (i.e. for a fixed growth rate) undetermined, and so are the jump conditions. By generating this trivial solution, the mathematics decouples both domains. The information coming from upstream to downstream (thanks to the shape of the linearized Euler system), a solution is found only for the upstream domain, the downstream being not solved and remains undetermined in lack of information. This solution cannot happen in reality because the fresh and burned gases in a real system are connected by many aspects, and not only by the jump conditions.

Stretch / Compression of Premixed Flame

The last of the `big' three turbulent ingredients (including aforementioned curvature and instabilities) impacting the flame at the local level is the compression or stretch of a premixed flame due to inhomogeneities in the flow, likely to happen with turbulence.

Interpretation of flow inhomogeneities stresses on a 1-D premixed flame in terms of compression or stretch.

The physics is pictured in the beside figure. In the flame area, the mass flow rate along the main direction may decrease (or increase) with distance. In a frame of reference whose origin is the core of the flame, a stretch (compression) results due to the difference in mass flow rate entering and leaving the flame zone. Important Remark Here, we are interested in the inhomogeneity in the mass flow rate along one direction only. We do NOT write that the flame volume is a region of mass source/sink.

In the same manner as above for the curvature effect, we introduce a (small) inhomogeneity of the flow as  \vec\nabla\vec M such that the following equation can be substracted from the unperturbed flame equation.

||\vec M^{s/c}||||\vec\nabla_{\xi}\theta|| - \vec n\cdot \vec\nabla\vec M \cdot \vec n ||\vec\nabla_{\xi}\theta||   - \varpi= \nabla_{\xi}\cdot\vec\nabla_{\xi}\theta

Hence, for a small perturbation, the linearized dependence of the flame on stretch/compression is:

||\vec M^{s/c}|| = 1 + \vec n \cdot \vec\nabla\vec M \cdot \vec n

The Partially-Premixed Regime

Ideal sketch of a partially-premixed flame

This regime is somewhat less academics and has been recognized two decades ago. It is acknowledged as a hybrid of the premixed and the non-premixed regimes but the degree of interaction of these two modes of combustion to accurately describe a partially-premixed flame is still not well understood. It can be simply pictured by a lifted diffusion flame. Let us consider fuel issuing from a nozzle into the air. If the exit velocity is large enough, for some fuels, the flame lifts off the rim of the nozzle. It means that below the flame base, fuel and oxidizer have room to premix. Hence, the flame base propagates into a premixed mixture. However, it cannot be reduced to a premixed flame (although it is often simplified as this): (i) the mixing is not perfect and the different parts of the flame front constituting the flame base burn in mixtures of different thermodynamical states. This provides those parts with different deflagration capabilities such that the flame base has a complex shape. Indeed, it is convex, naturally leaded by the part burning at stoichiometry, unless `exotic' feeding temperatures are used. (ii) Because the mixture is not homogeneous, transfer of species and temperature driven by diffusion occurs in a direction perpendicular to the propagation of the flame base. Because the flame front is not flat, those transfers act as a connection vehicle across the different parts of the leading front. (iii) The unburned left downstream by the sections of the leading front not burning at stoichiometry diffuse towards each other to form a diffusion flame as described above. The connection of the leading front with the trailing diffusion flame has been evidenced as complicated and the siege of transfers of species and temperature. These two last items are the state-of-the-art difficulties in understanding those flames and do not appear in the models although it has been demonstrated they have a major impact and are certainly a fundamental characteristic of partially-premixed flames.

The partially-premixed flame is usually described using c and Z as introduced earlier. Because the framework is essentially non-premixed, the mixture fraction is primarily used to describe the flame. Regarding the head of the flame where partial-premixing has an impact, each part of the front is described with a local progress variable. The need of defining a local progress variable is that each section of the partially-premixed front has a different equivalence ratio leading to a different definition of c:


Three Turbulent-Flame Interaction Regimes

Combustion is a chemical reaction and turbulence is recognized as a powerful mean to enhance chemical reaction by accelerating fuel and oxidizer mixing. In the case of combustion, the impact is more subtle. Through the strong heat release in a narrow zone, a flame is able to accelerate and distord the streamlines. It may lead to flame-generated turbulence as the vortices created behind a premixed flame subject to Darrieus instabilities, vid. sup. Sec. Hydrodynamics Instabilities, or, a contrario, kill the existing eddy motion. Turbulence may also be diminished through a flame because gas viscosity increases with temperature. On the other hand, turbulence may curve a flame and disturb its internal structure, increasing the reaction rate and/or stretching the flame to extinction.

It may appear odd to try to describe here what is the reason of combustion modelling research: the interaction of the turbulence with chemistry. However, one of the first steps in building knowledge in turbulent combustion was the qualitative exploration of what might be the dynamics of a flame in a turbulent environment. This led to what is now known as combustion diagrams. As explained above, the premixed regime lends itself the easiest to such an approach as it exhibits natural intrinsic quantities which are not as objectively identifiable in the other combustion regimes. Note that these quantities may depend on the geometry of the flame: for instance turbulence can bend a flame sheet, leading to a change in its dynamics compared to the flat flame propagating in a medium at rest. In this section, the turbulence-flame interaction modes will be described for a premixed flame. Only remarks will be added regarding the non-premixed and partially-premixed regimes.

An integral quantity to assess the interaction between a premixed flame sheet and the turbulence is the Karlovitz number Ka. It compares the characteristic time of flame displacement with the characteristic time of the smallest structures (that are also the fastest) of the turbulence.

 Ka= \frac{\tau_c}{\tau_k}

\tau_c is the chemical time of the flame. To estimate it, it is necessary to come back to the above progress variable transport equation in a steady-state framework.

 \overbrace{(\rho S_L)}^{\dot M}\frac{\vec S_L}{S_L}\vec\nabla c = (\rho D)_f\nabla\cdot (\rho D)^* \vec\nabla c + \dot\omega_c

The premixed wave propagates at a speed S_L because it is fed by reactants diffusing inside the combustion zone and which are preheated because temperature diffuses in the reverse direction. The speed at which the flame progresses is thus related to the rate of species diffusion into the reaction zone which are then consumed. As the premixed flame is a free propagating wave whose speed of propagation is only limited by the chemical strength, the characteristic chemical time is based only on the diffusion and the mass flow rate experienced by the flame:

 \tau_c = \frac{\rho (\rho D)_f}{\dot M^2}

The smallest eddies are the ones being dissipated by the viscous forces. Their characteristic time is estimated thanks to a combination of the viscosity and the flux of turbulent energy to be dissipated (also called turbulent dissipation  \varepsilon=u'^3/l_t ):

\tau_k = \sqrt{\frac{l_t\nu}{u'^3}}

Thanks to those definitions of chemical and small structure times, it is possible to give another definition of the Karlovitz number:

 Ka=\left (\frac{\delta_L}{l_k} \right)^2

which is the square of the ratio between the premixed flame thickness and the small structure scale: Ka actually compares scales. To arrive to this latter result, the three following assumptions must be used: (i) the flame thickness is obtained thanks to the premixed flame Péclet number (vid. sup.); (ii) the turbulence small structure (Kolmogorov eddies) scale is given by:  l_k=(\nu^3/\varepsilon)^{1/4} following the same dimensional argument as for the estimation of its time; and (iii) scalar diffusion scales with viscosity.

Remark Regarding the Diffusion Flame

From what has been presented above, a diffusion flame does not have characteristic scales. Setting a turbulence combustion regime classification for non-premixed flames has still not been answered by research. Some laws of behaviour will only be drawn around the scalar dissipation rate which is the parameter of integral importance for a diffusion flame.

Indeed, the dynamics of a diffusion flame is determined by the strain rate imposed by the turbulence. As for the premixed flames, the shortest eddies (Kolmogorov) are the ones having the largest impact. The diffusive layer is thus given by the size of the Kolmogrov eddies:  l_d\approx l_k and the typical diffusion time scale (feeding rate of the reaction zone) is given by the characteristic time of the Kolmogorov eddies:  \tau_k^{-1}\approx \chi_s as the Reynolds number of the Kolmogorov structures is unity. Here, \chi_s is the sample-averaging of \chi based on (conditioned) stoichiometric conditions, where the flame is expected to be.

The Wrinkled Regime

Wrinkled flamelet regime

This regime is also called the flamelet regime. Basically, it assumes that the flame structure is not affected by turbulence. The flame sheet is convoluted and wrinkled by eddies but none of them is small enough to enter it. Locally magnifying, the laminar flame structure is maintained.

This regime exists for a Karlovitz number below unity (vid. sup.), i.e. chemical time smaller than the small structure time or flame thickness smaller than small structure scale. Notwithstanding, the laminar flame dynamics can be disrupted for  u'>S_L. In that case, although the flame structure is not altered by the small structures, it can be convected by large structures such that areas of different locations in the front interact. It shows that, even with a small Karlovitz number, the turbulence effect is not always weak.

The Corrugated Regime

Corrugated flamelet regime

The formal definition of a flame is the region of temperature rise. However, the volume where the reaction takes place is about one order of magnitude smaller, embedded inside the temperature rise region and close to its high temperature end. Hence, there exist some levels of turbulence creating eddies able to enter the flame zone but still large enough to not affect the internal reaction sheet. In other words, the flame thermal region is thickened by turbulence but the reaction zone is still in the wrinkled regime. This situation is called the Corrugated Regime.

Due to the structure of the Karlovitz number, once written in terms of length scales (vid. sup.), this situation arises for an increase of the Karlovitz number by two orders of magnitude compared to the value for the wrinkled regime. Hence, in the range  1 < Ka < 100 , the laminar structure of the reaction zone is still preserved but not the one of the preheat zone.

The Thickened Regime

Thickened regime

In this last case, turbulence is intense enough to generate eddies able to affect the structure of the reaction zone as well. In practice, it is expected that those eddies are in the tail of the energy spectrum such that their lifetime is very short. Their impact on the reaction zone is thus limited.

Obviously, Ka > 100. A topological description is of little relevance here and a well-stirred reactor model fits better.

Description of Real Mixtures

As explained above, combustion has been analyzed and modelled with the help of simple irreversible chemistry involving `a' fuel and `an' oxidizer. However, combustion is a chemical reaction highly complicated, except maybe for ozone and hydrogen in oxygen. A chemical mechanism involving more than thousand species and several hundreds of reactions between them is still unsuficient to describe the combustion of a simple hydrocarbon molecule as methane in air. This is due to the chemistry route from reaction to reaction that may depend on the conditions such as pressure and temperature. Also, species are unstable at high temperature, meaning that the combustion products are usually exposed to dissociations. Therefore, the full conversion of the chemical enthalpy into sensible enthalpy is not realized with substantial drop in combustion temperature. Furthermore, the calorimetry and transport properties of the species in the mixture are not unique: the heat capacities are temperature dependent and species diffuse at various speeds within each other, potentially leading to exotic phenomena for light molecules as hydrogen.

Governing Equations for Chemically Reacting Flows

Together with the usual Navier-Stokes equations for compresible flow (See Governing equations), additional equations are needed in reacting flows for species. In the earlier paragraphs dedicated to combustion fundamentals, we have seen that the description is reduced to simple transport equations. However, in real mixtures, numerous species with different transport properties may exist. Their transport equations are coupled together more or less directly by their rate of creation or disappearence (chemical reaction source/sink terms), and by their transport properties relatively to each other. They are also coupled to the Navier-Stokes equations (mass, momentum, and energy) thanks to density, velocity, pressure and temperature taking a role in their diffusion and convection fluxes and their source terms.

Due to these complicated behaviours, it is rather chosen to speak in terms of diffusion velocities  V^j_k such that the transport equation for the mass fraction  Y_k of the k-th specie in a complex mixture takes the form:

\frac{\partial}{\partial t} \left( \rho Y_k \right) +
\frac{\partial}{\partial x_j} \left( \rho u^j Y_k\right) = -
\frac{\partial}{\partial x_j} \left( \rho V_k^j Y_k\right)+ \dot \omega_k

 \dot \omega_k denotes the specie reaction rate.

A non-reactive (passive) scalar (like the mixture fraction  Z ) is goverened by the following transport equation

\frac{\partial}{\partial t} \left( \rho Z \right) +
\frac{\partial}{\partial x_j} \left( \rho u_j Z \right) = 
\frac{\partial}{\partial x_j} \left( \rho D \frac{\partial Z}{\partial x_j}\right)

where  D is the diffusion coefficient of the passive scalar.

RANS equations

In turbulent flows, Favre averaging is often used to reduce the scales (see Reynolds averaging) and the mass fraction transport equation is transformed to

\frac{\partial \overline{\rho} \widetilde{Y}_k }{\partial t} +
\frac{\partial \overline{\rho} \widetilde{u}_j \widetilde{Y}_k}{\partial x_j}=
\frac{\partial} {\partial x_j} \left( \overline{\rho D_k  \frac{\partial Y_k} {\partial x_j} } -
 \overline{\rho} \widetilde{u''_i Y''_k } \right)
+ \overline{\dot \omega_k}

where the turbulent fluxes  \widetilde{u''_i Y''_k} and reaction terms   \overline{\dot \omega_k} need to be closed.

The passive scalar turbulent transport equation is

\frac{\partial \overline{\rho} \widetilde{Z} }{\partial t} +
\frac{\partial \overline{\rho} \widetilde{u}_j \widetilde{Z} }{\partial x_j}=
\frac{\partial} {\partial x_j} \left( \overline{\rho D  \frac{\partial Z} {\partial x_j} } -
 \overline{\rho} \widetilde{u''_i Z'' } \right)

where  \widetilde{u''_i Z''} needs modelling. A common practice is to model the turbulent fluxes using the gradient diffusion hypothesis. For example, in the equation above the flux  \widetilde{u''_i Z''} is modelled as

\widetilde{u''_i Z''} = -D_t \frac{\partial \tilde Z}{\partial x_i}

where  D_t is the turbulent diffusivity. Since  D_t >> D , the first term inside the parentheses on the right hand of the mixture fraction transport equation is often neglected (Peters (2000)). This assumption is also used below.

In addition to the mean passive scalar equation, an equation for the Favre variance  \widetilde{Z''^2} is often employed

\frac{\partial \overline{\rho} \widetilde{Z''^2} }{\partial t} +
\frac{\partial \overline{\rho} \widetilde{u}_j \widetilde{Z''^2} }{\partial x_j}=
\frac{\partial}{\partial x_j}
\left(  \overline{\rho} \widetilde{u''_i Z''^2}  \right) -
 2 \overline{\rho} \widetilde{u''_i Z'' }
- \overline{\rho} \widetilde{\chi}

where  \widetilde{\chi} is the mean Scalar dissipation rate defined as  \widetilde{\chi} =
2 D \widetilde{\left| \frac{\partial Z''}{\partial x_j} \right|^2 } This term and the variance diffusion fluxes needs to be modelled.

LES equations

The Large eddy simulation (LES) approach for reactive flows introduces equations for the filtered species mass fractions within the compressible flow field. Similar to the #RANS equations, but using Favre filtering instead of Favre averaging, the filtered mass fraction transport equation is

\frac{\partial \overline{\rho} \widetilde{Y}_k }{\partial t} +
\frac{\partial \overline{\rho} \widetilde{u}_j \widetilde{Y}_k}{\partial x_j}=
\frac{\partial} {\partial x_j} \left( \overline{\rho D_k  \frac{\partial Y_k} {\partial x_j} } -
J_j \right)
+ \overline{\dot \omega_k}

where  J_j is the transport of subgrid fluctuations of mass fraction

J_j = \widetilde{u_jY_k} - \widetilde{u}_j \widetilde{Y}_k

and has to be modelled.

Fluctuations of diffusion coefficients are often ignored and their contributions are assumed to be much smaller than the apparent turbulent diffusion due to transport of subgrid fluctuations. The first term on the right hand side is then

\frac{\partial} {\partial x_j}  
\overline{ \rho D_k  \frac{\partial Y_k}  {\partial x_j} }
\frac{\partial} {\partial x_j}  
\overline{\rho} D_k  \frac{\partial \widetilde{Y}_k} {\partial x_j}  

Reaction mechanisms

Combustion is mainly a chemical process. Although we can, to some extent, describe a flame without any chemistry information, modelling of the flame propagation requires the knowledge of speeds of reactions, product concentrations, temperature, and other parameters. Therefore fundamental information about reaction kinetics is essential for any combustion model. A fuel-oxidizer mixture will generally combust if the reaction is fast enough to prevail until all of the mixture is burned into products. If the reaction is too slow, the flame will extinguish. If too fast, explosion or even detonation will occur. The reaction rate of a typical combustion reaction is influenced mainly by the concentration of the reactants, temperature, and pressure.

A stoichiometric equation for an arbitrary reaction can be written as:

\sum_{j=1}^{n}\nu' (M_j) = \sum_{j=1}^{n}\nu'' (M_j),

where  \nu denotes the stoichiometric coefficient, and M_j stands for an arbitrary species. A one prime holds for the reactants while a double prime holds for the products of the reaction.

The reaction rate, expressing the rate of disappearance of reactant i, is defined as:

RR_i = k \, \prod_{j=1}^{n}(M_j)^{\nu'},

in which k is the specific reaction rate constant. Arrhenius found that this constant is a function of temperature only and is defined as:

k= A T^{\beta} \, exp \left( \frac{-E}{RT}\right)

where A is pre-exponential factor, E is the activation energy, and \beta is a temperature exponent. The constants vary from one reaction to another and can be found in the literature.

Reaction mechanisms can be deduced from experiments (for every resolved reaction), they can also be constructed numerically by the automatic generation method (see [Griffiths (1994)] for a review on reaction mechanisms). For simple hydrocarbons, tens to hundreds of reactions are involved. By analysis and systematic reduction of reaction mechanisms, global reactions (from one to five step reactions) can be found (see [Westbrook (1984)]).

Infinitely fast chemistry

All combustion models can be divided into two main groups according to the assumptions on the reaction kinetics. We can either assume the reactions to be infinitely fast - compared to e.g. mixing of the species, or comparable to the time scale of the mixing process. The simple approach is to assume infinitely fast chemistry. Historically, mixing of the species is the older approach, and it is still in wide use today. It is therefore simpler to solve for #Finite rate chemistry models, at the overshoot of introducing errors to the solution.

Premixed combustion

Premixed flames occur when the fuel and oxidiser are homogeneously mixed prior to ignition. These flames are not limited only to gas fuels, but also to the pre-vaporised fuels. Typical examples of premixed laminar flames is bunsen burner, where the air enters the fuel stream and the mixture burns in the wake of the riser tube walls forming nice stable flame. Another example of a premixed system is the solid rocket motor where oxidizer and fuel and properly mixed in a gum-like matrix that is uniformly distributed on the periphery of the chamber. Premixed flames have many advantages in terms of control of temperature and products and pollution concentration. However, they introduce some dangers like the autoignition (in the supply system).


Turbulent flame speed model

Eddy Break-Up model

The Eddy Break-Up model is the typical example of mixed-is-burnt combustion model. It is based on the work of Magnussen, Hjertager, and Spalding and can be found in all commercial CFD packages. The model assumes that the reactions are completed at the moment of mixing, so that the reaction rate is completely controlled by turbulent mixing. Combustion is then described by a single step global chemical reaction

F + \nu_s O \rightarrow (1+\nu_s) P

in which F stands for fuel, O for oxidiser, and P for products of the reaction. Alternativelly we can have a multistep scheme, where each reaction has its own mean reaction rate. The mean reaction rate is given by

\bar{\dot\omega}_F=A_{EB} \frac{\epsilon}{k} 

where \bar{C} denotes the mean concentrations of fuel, oxidiser, and products respectively. A and B are model constants with typical values of 0.5 and 4.0 respectively. The values of these constants are fitted according to experimental results and they are suitable for most cases of general interest. It is important to note that these constants are only based on experimental fitting and they need not be suitable for all the situations. Care must be taken especially in highly strained regions, where the ratio of k to \epsilon is large (flame-holder wakes, walls ...). In these, regions a positive reaction rate occurs and an artificial flame can be observed. CFD codes usually have some remedies to overcome this problem.

This model largely over-predicts temperatures and concentrations of species like CO and other species. However, the Eddy Break-Up model enjoys a popularity for its simplicity, steady convergence, and implementation.

Bray-Moss-Libby model

Non-premixed combustion

Non premixed combustion is a special class of combustion where fuel and oxidizer enter separately into the combustion chamber. The diffusion and mixing of the two streams must bring the reactants together for the reaction to occur. Mixing becomes the key characteristic for diffusion flames. Diffusion burners are easier and safer to operate than premixed burners. However their efficiency is reduced compared to premixed burners. One of the major theoretical tools in non-premixed combustion is the passive scalar mixture fraction  Z which is the backbone on most of the numerical methods in non-premixed combustion.

Conserved scalar equilibrium models

The reactive problem is split into two parts. First, the problem of mixing , which consists of the location of the flame surface which is a non-reactive problem concerning the propagation of a passive scalar; And second, the flame structure problem, which deals with the distribution of the reactive species inside the flamelet.

To obtain the distribution inside the flame front we assume it is locally one-dimensional and depends only on time and the scalar coodinate.

We first make use of the following chain rules

\frac{\partial Y_k}{\partial t} = \frac{\partial Z}{\partial t}\frac{\partial Y_k}{\partial Z}

\frac{\partial Y_k}{\partial x_j} = \frac{\partial Z}{\partial x_j}\frac{\partial Y_k}{\partial Z}

and transformation

\frac{\partial }{\partial t}= \frac{\partial }{\partial t} + \frac{\partial Z}{\partial t} \frac{\partial }{\partial Z}

upon substitution into the species transport equation (see #Governing Equations for Reacting Flows), we obtain

\rho \frac{\partial Y_k}{\partial t} + Y_k \left[
\frac{\partial \rho}{\partial t} + \frac{\partial \rho u_j}{\partial x_j}
+ \frac{\partial Y_k}{\partial Z} \left[
 \rho \frac{\partial Z}{\partial t} + \rho u_j \frac{\partial Z}{\partial x_j} -
\frac{\partial}{\partial x_j}\left( \rho D \frac{\partial Z}{\partial x_j} \right)
\rho D \left( \frac{\partial Z}{\partial x_j} \frac{\partial Z}{\partial x_j} \right)
\frac{\partial^2 Y_k}{\partial Z^2} + \dot \omega_k

The second and third terms in the LHS cancel out due to continuity and mixture fraction transport. At the outset, the equation boils down to

\frac{\partial Y_k}{\partial t} = \frac{\chi}{2} \frac{\partial ^2 Y_k}{\partial Z^2} + \dot \omega_k

where  \chi = 2 D  \left( \frac{\partial Z}{\partial x_j} \right)^2 is called the scalar dissipation which controls the mixing, providing the interaction between the flow and the chemistry.

If the flame dependence on time is dropped, even though the field  Z  still depends on it.

 \dot \omega_k= -\frac{\chi}{2} \frac{\partial ^2 Y_k}{\partial Z^2}

If the reaction is assumed to be infinetly fast, the resultant flame distribution is in equilibrium. and  \dot \omega_k= 0. When the flame is in equilibrium, the flame configuration  Y_k(Z) is independent of strain.

Burke-Schumann flame structure

The Burke-Schuman solution is valid for irreversible infinitely fast chemistry. With a reaction in the form of

F + \nu_s O \rightarrow (1+\nu_s) P

If the flame is in equilibrium and therefore the reaction term is 0. Two possible solution exists, one with pure mixing (no reaction) and a linear dependence of the species mass fraction with  Z . Fuel mass fraction

Y_F=Y_F^0 Z

Oxidizer mass fraction


Where  Y_F^0 and  Y_O^0 are fuel and oxidizer mass fractions in the pure fuel and oxidizer streams respectively.

The other solution is given by a discontinuous slope at stoichiometric mixture fraction  Z_{st} and two linear profiles (in the rich and lean side) at either side of the stoichiometric mixture fraction. Both concentrations must be 0 at stoichiometric, the reactants become products infinitely fast.

Y_F=Y_F^0 \frac{Z-Z_{st}}{1-Z_{st}}

and oxidizer mass fraction

Y_O=Y_O^0 \frac{Z-Z_{st}}{Z_{st}}

Finite rate chemistry

Premixed combustion

Coherent flame model

Flamelets based on G-equation

Flame surface density model

Non-premixed combustion

Flamelets based on conserved scalar

Peters (2000) define Flamelets as "thin diffusion layers embedded in a turbulent non-reactive flow field". If the chemistry is fast enough, the chemistry is active within a thin region where the chemistry conditions are in (or close to) stoichiometric conditions, the "flame" surface. This thin region is assumed to be smaller than Kolmogorov length scale and therefore the region is locally laminar. The flame surface is defined as an iso-surface of a certain scalar  Z , mixture fraction in #Non premixed combustion.

The same equation used in #Conserved scalar models for equilibrium chemistry is used here but with chemical source term different from 0

 \frac{\partial Y_k}{\partial t} =
\frac{\chi}{2} \frac{\partial ^2 Y_k}{\partial Z^2} + \dot \omega_k.

This approach is called the Stationary Laminar Flamelet Model (SLFM) and has the advantage that flamelet profiles  Y_k=f(Z,\chi) can be pre-computed and stored in a dtaset or file which is called a "flamelet library" with all the required complex chemistry. For the generation of such libraries ready to use software is avalable such as Softpredict's Combustion Simulation Laboratory COSILAB [1] with its relevant solver RUN1DL, which can be used for a variety of relevant geometries; see various publications that are available for download. Other software tools are available such as CHEMKIN [2] and CANTERA [3].

Flamelets in turbulent combustion

In turbulent flames the interest is  \widetilde{Y}_k . In flamelets, the flame thickness is assumed to be much smaller than Kolmogorov scale and obviously is much smaller than the grid size. It is therefore needed a distribution of the passive scalar within the cell.  \widetilde{Y}_k cannot be obtained directly from the flamelets library  \widetilde{Y}_k \neq Y_F(Z,\chi) , where  Y_F(Z,\chi) corresponds to the value obtained from the flamelets libraries. A generic solution can be expressed as

\widetilde{Y}_k= \int  Y_F( \widetilde{Z},\widetilde{\chi}) P(Z,\chi) dZ d\chi

where  P(Z,\chi) is the joint Probability Density Function (PDF) of the mixture fraction and scalar dissipation which account for the scalar distribution inside the cell and "a priori" depends on time and space.

The most simple assumption is to use a constant distribution of the scalar dissipation within the cell and the above equation reduces to

\widetilde{Y}_k= \int Y_F(\widetilde{Z},\widetilde{\chi}) P(Z) dZ

 P(Z) is the PDF of the mixture fraction scalar and simple models (such as Gaussian or a beta PDF) can be build depending only on two moments of the scalar mean and variance, \widetilde{Z},Z''.

If the mixture fraction and scalar dissipation are consider independent variables,  P(Z,\chi) can be written as  P(Z) P(\chi). The PDF of the scalar dissipation is assumed to be log-normal with variance unity.

\widetilde{Y}_k= \int Y_F(\widetilde{Z},\widetilde{\chi}) P(Z) P(\chi) dZ d\chi

In Large eddy simulation (LES) context (see #LES equations for reacting flow), the probability density function is replaced by a subgrid PDF  \widetilde{P}. The same equation hold by replacing averaged values with filtered values.

\widetilde{Y}_k= \int Y_F(\widetilde{Z},\widetilde{\chi}) \widetilde{P}(Z) \widetilde{P}(\chi) dZ d\chi

The assumptions made regarding the shapes of the PDFs are still justified. In LES combustion the subgrid variance is smaller than RANS counterpart (part of the large-scale fluctuations are solved) and therefore the modelled PDFs are thinner.

Unsteady flamelets

Intrinsic Low Dimensional Manifolds (ILDM)

Detailed mechanisms describing ignition, flame propagation and pollutant formation typically involve several hundred species and elementary reactions, prohibiting their use in practical three-dimensional engine simulations. Conventionally reduced mechanisms often fail to predict minor radicals and pollutant precursors. The ILDM-method is an automatic reduction of a detailed mechanism, which assumes local equilibrium with respect to the fastest time scales identified by a local eigenvector analysis. In the reactive flow calculation, the species compositions are constrained to these manifolds. Conservation equations are solved for only a small number of reaction progress variables, thereby retaining the accuracy of detailed chemical mechanisms. This gives an effective way of coupling the details of complex chemistry with the time variations due to turbulence.

The intrinsic low-dimensional manifold (ILDM) method {Maas:1992,Maas:1993} is a method for in-situ reduction of a detailed chemical mechanism based on a local time scale analysis. This method is based on the fact that different trajectories in state space start from a high-dimensional point and quickly relax to lower-dimensional manifolds due to the fast reactions. The movement along these lower-dimensional manifolds, however, is governed by the slow reactions. It exploits the variety of time scales to systematically reduce the detailed mechanism. For a detailed chemical mechanism with N species, N different time scales govern the process. An assumption that all the time scales are relaxed results in assuming complete equilibrium, where the only variables required to describe the system are the mixture fraction, the temperature and the pressure. This results in a zero-dimensional manifold. An assumption that all but the slowest 'n' time scales are relaxed results in a 'n' dimensional manifold, which requires the additional specification of 'n' parameters (called progress variables). In the ILDM method, the fast chemical reactions do not need to be identified a priori. An eigenvalue analysis of the detailed chemical mechanism is carried out which identifies the fast processes in dynamic equilibrium with the slow processes. The computation of ILDM points can be expensive, and hence an in-situ tabulation procedure is used, which enables the calculation of only those points that are needed during the CFD calculation.

--Fredgauss 07:37, 25 August 2006 (MDT)

U. Maas, S.B. Pope. Simplifying chemical kinetics: Intrinsic low-dimensional manifolds in composition space. Comb. Flame 88, 239, 1992.

Ulrich Maas. Automatische Reduktion von Reaktionsmechanismen zur Simulation reaktiver Str¨omungen. Habilitationsschrift, Universit ¨at Stuttgart, 1993.

Conditional Moment Closure (CMC)

In Conditional Moment Closure (CMC) methods we assume that the species mass fractions are all correlated with the mixture fraction (in non premixed combustion).

From Probability density function we have

\overline{Y_k}= \int <Y_k|\eta> P(\eta) d\eta

where  \eta is the sample space for  Z .

CMC consists of providing a set of transport equations for the conditional moments which define the flame structure.

Experimentally, it has been observed that temperature and chemical radicals are strong non-linear functions of mixture fraction. For a given species mass fraction we can decomposed it into a mean and a fluctuation:

Y_k= \overline{Y_k} + Y'_k

The fluctuations  Y_k' are usually very strong in time and space which makes the closure of  \overline{\omega_k} very difficult. However, the alternative decomposition

Y_k=  <Y_k|\eta> + y'_k

where  y'_k is the fluctuation around the conditional mean or the "conditional fluctuation". Experimentally, it is observed that  y'_k<< Y'_k , which forms the basic assumption of the CMC method. Closures. Due to this property better closure methods can be used reducing the non-linearity of the mass fraction equations.

The Derivation of the CMC equations produces the following CMC transport equation where  Q \equiv <Y_k|\eta> for simplicity.

\frac{ \partial Q}{\partial t} + <u_j|\eta>  \frac{\partial Q}{\partial x_j} =
\frac{<\chi|\eta> }{2} \frac{\partial ^2 Q}{\partial \eta^2} + 
\frac{ < \dot \omega_k|\eta> }{ <\rho| \eta >}

In this equation, high order terms in Reynolds number have been neglected. (See Derivation of the CMC equations for the complete series of terms).

It is well known that closure of the unconditional source term   \overline {\dot \omega_k} as a function of the mean temperature and species ( \overline{Y}, \overline{T}) will give rise to large errors. However, in CMC the conditional averaged mass fractions contain more information and fluctuations around the mean are much smaller. The first order closure 
 < \dot \omega_k|\eta> \approx \dot \omega_k \left( Q, <T|\eta> \right)
is a good approximation in zones which are not close to extinction.

Second order closure

A second order closure can be obtained if conditional fluctuations are taken into account. For a chemical source term in the form  \dot \omega_k = k Y_A Y_B with the rate constant in Arrhenius form  k=A_0 T^\beta exp [-Ta/T] the second order closure is (Klimenko and Bilger 1999)

< \dot \omega_k|\eta> \approx  < \dot \omega_k|\eta >^{FO}

\left[1+ \frac{< Y''_A Y''_B |\eta>}{Q_A Q_B}+ \left( \beta + T_a/Q_T \right)
\frac{< Y''_A T'' |\eta>}{Q_AQ_T} + \frac{< Y''_B T'' |\eta>}{Q_BQ_T}
\right) + ...


where  < \dot \omega_k|\eta >^{FO} is the first order CMC closure and  Q_T \equiv <T|\eta> . When the temperature exponent  \beta or  T_a/Q_T are large the error of taking the first order approximation increases. Improvement of small pollutant predictions can be obtained using the above reaction rate for selected species like CO and NO.

Double conditioning

Close to extinction and reignition. The conditional fluctuations can be very large and the primary closure of CMC of "small" fluctuations is not longer valid. A second variable  h can be chosen to define a double conditioned mass fraction

Q(x,t;\eta,\psi) \equiv <Y_i(x,t) |Z=\eta,h=\psi  >

Due to the strong dependence on chemical reactions to temperature,  h is advised to be a temperature related variable (Kronenburg 2004). Scalar dissipation is not a good choice, due to its log-normal behaviour (smaller scales give highest dissipation). A must better choice is the sensible enthalpy or a progress variable. Double conditional variables have much smaller conditional fluctuations and allow the existence of points with the same chemical composition which can be fully burning (high temperature) or just mixing (low temperature). The range of applicability is greatly increased and allows non-premixed and premixed problems to be treated without ad-hoc distinctions. The main problem is the closure of the new terms involving cross scalar transport.

The double conditional CMC equation is obtained in a similar manner than the conventional CMC equations

LES modelling

In a LES context a conditional filtering operator can be defined and  Q therefore represents a conditionally filtered reactive scalar.

Linear Eddy Model

The Linear Eddy Model (LEM) was first developed by Kerstein(1988). It is an one-dimensional model for representing the flame structure in turbulent flows.

In every computational cell a molecular, diffusion and chemical model is defined as

\frac{\partial}{\partial t} \left( \rho Y_k \right)  =
\frac{\partial}{\partial \eta} \left( \rho D_k \frac{\partial Y_k}{\partial \eta }\right)+ \dot \omega_k

where  \eta is a spatial coordinate. The scalar distribution obtained can be seen as a one-dimensional reference field between Kolmogorov scale and grid scales.

In a second stage a series of re-arranging stochastic event take place. These events represent the effects of a certain turbulent structure of size  l , smaller than the grid size at a location  \eta_0  within the one-dimensional domain. This vortex distort the  \eta field obtain by the one-dimensional equation, creating new maxima and minima in the interval  (\eta_0, \eta + \eta_0) . The vortex size  l is chosen randomly based on the inertial scale range while  \eta_0  is obtained from a uniform distribution in  \eta  . The number of events is chosen to match the turbulent diffusivity of the flow.

PDF transport models

Probability Density Function (PDF) methods are not exclusive to combustion, although they are particularly attractive to them. They provided more information than moment closures and they are used to compute inhomegenous turbulent flows, see reviews in Dopazo (1993) and Pope (1994).

PDF methods are based on the transport equation of the joint-PDF of the scalars. Denoting  P \equiv P(\underline{\psi}; x, t) where  \underline{\psi} = ( \psi_1,\psi_2 ... \psi_N) is the phase space for the reactive scalars  \underline{Y} = ( Y_1,Y_2 ... Y_N) . The transport equation of the joint PDF is:

\frac{\partial <\rho | \underline{Y}=\underline{\psi}> P }{\partial t} + \frac{  
\partial <\rho u_j | \underline{Y}=\underline{\psi}> P }{\partial x_j} =
\sum^N_\alpha \frac{\partial}{\partial \psi_\alpha}\left[ \rho \dot{\omega}_\alpha P \right]
- \sum^N_\alpha \sum^N_\beta \frac{\partial^2}{\partial \psi_\alpha \psi_\beta}
\left[ <D \frac{\partial Y_\alpha}{\partial x_i} \frac{\partial Y_\beta}{\partial x_i} | \underline{Y}=\underline{\psi}> \right] P

where the chemical source term is closed. Another term appeared on the right hand side which accounts for the effects of the molecular mixing on the PDF, is the so called "micro-mixing " term. Equal diffusivities are used for simplicity  D_k = D

A more general approach is the velocity-composition joint-PDF with  P \equiv P(\underline{V},\underline{\psi}; x, t) , where  \underline{V} is the sample space of the velocity field  u,v,w . This approach has the advantage of avoiding gradient-diffusion modelling. A similar equation to the above is obtained combining the momentum and scalar transport equation.

The PDF transport equation can be solved in two ways: through a Lagrangian approach using stochastic methods or in a Eulerian ways using stochastic fields.


The main idea of Lagrangian methods is that the flow can be represented by an ensemble of fluid particles. Central to this approach is the stochastic differential equations and in particular the Langevin equation.


Instead of stochastic particles, smooth stochastic fields can be used to represent the probability density function (PDF) of a scalar (or joint PDF) involved in transport (convection), diffusion and chemical reaction (Valino 1998). A similar formulation was proposed by Sabelnikov and Soulard 2005, which removes part of the a-priori assumption of "smoothness" of the stochastic fields. This approach is purely Eulerian and offers implementations advantages compared to Lagrangian or semi-Eulerian methods. Transport equations for scalars are often easy to programme and normal CFD algorithms can be used (see Discretisation of convective term). Although discretization errors are introduced by solving transport equations, this is partially compesated by the error introduced in Lagrangian approaches due to the numerical evaluation of means.

A new set of  N_s scalar variables (the stochastic field  \xi ) is used to represent the PDF

P (\underline{\psi}; x,t) = \frac{1}{N} \sum^{N_s}_{j=1} \prod^{N}_{k=1} \delta \left[\psi_k -\xi_k^j(x,t) \right]

Other combustion models


The Multiple Mapping Conditioning (MMC) (Klimenko and Pope 2003) is an extension of the #Conditional Moment Closure (CMC) approach combined with probability density function methods. MMC looks for the minimum set of variables that describes the particular turbulent combustion system.


Derived from the #Eddy Dissipation Concept (EDC).


[1] Bone, W.A., and Townend, D.T.A. (1927), Flame and Combustion in Gases, Longmans, Green, and Co., Ltd.

[2] (1999), "", Journal of Heat Transfer, Vol. 121, No. 4, pp. 770-73.

[3] Liñan, A. (1974), "The Asymptotic Structure of Counter Flow Diffusion Flames for Large Activation Energies", Acta Astronautica, Vol. 1, No. 7/8, pp. 1007-1039.

  • Dopazo, C. (1993), "Recent development in PDF methods", Turbulent Reacting Flows, ed. P. A. Libby and F. A. Williams.
  • Fox, R.O. (2003), Computational Models for Turbulent Reacting Flows, ISBN 0-521-65049-6,Cambridge University Press.
  • Kerstein, A. R. (1988), "A linear eddy model of turbulent scalar transport and mixing", Comb. Science and Technology, Vol. 60,pp. 391.
  • Klimenko, A. Y., Bilger, R. W. (1999), "Conditional moment closure for turbulent combustion", Progress in Energy and Combustion Science, Vol. 25,pp. 595-687.
  • Klimenko, A. Y., Pope, S. B. (2003), "The modeling of turbulent reactive flows based on multiple mapping conditioning", Physics of Fluids, Vol. 15, Num. 7, pp. 1907-1925.
  • Kronenburg, A., (2004), "Double conditioning of reactive scalar transport equations in turbulent non-premixed flames", Physics of Fluids, Vol. 16, Num. 7, pp. 2640-2648.
  • Griffiths, J. F. (1994), "Reduced Kinetic Models and Their Application to Practical Combustion Systems", Prog. in Energy and Combustion Science,Vol. 21, pp. 25-107.
  • Peters, N. (2000), Turbulent Combustion, ISBN 0-521-66082-3,Cambridge University Press.
  • Poinsot, T.,Veynante, D. (2001), Theoretical and Numerical Combustion, ISBN 1-930217-05-6, R. T Edwards.
  • Pope, S. B. (1994), "Lagrangian PDF methods for turbulent flows", Annu. Rev. Fluid Mech, Vol. 26, pp. 23-63.
  • Sabel'nikov, V.,Soulard, O. (2005), "Rapidly decorrelating velocity-field model as a tool for solving one-point Fokker-Planck equations for probability density functions of turbulent reactive scalars", Physical Review E, Vol. 72,pp. 016301-1-22.
  • Valino, L., (1998), "A field montecarlo formulation for calculating the probability density function of a single scalar in a turbulent flow", Flow. Turb and Combustion, Vol. 60,pp. 157-172.
  • Westbrook, Ch. K., Dryer,F. L., (1984), "Chemical Kinetic Modeling of Hydrocarbon Combustion", Prog. in Energy and Combustion Science,Vol. 10, pp. 1-57.

External links and sources

My wiki