These notes synthesize material from F.A. Williams’s Combustion Theory, J.D. Buckmaster and G.S.S. Ludford’s Theory of Laminar Flames, A. Liñán and F.A. Williams’s Fundamental Aspects of Combustion, J.W. Dold and A.K. Kapila’s lectures on asymptotic combustion theory, and supplementary material from Cambridge DAMTP Combustion notes, Caltech ACM/Ae combustion lecture materials, and Princeton mechanical and aerospace engineering course notes.
Chapter 1: Fundamentals of Reactive Flows
1.1 Conservation Laws for Reacting Mixtures
Combustion is the study of exothermic chemical reactions coupled to fluid motion and heat transfer. The mathematical framework begins with the conservation laws of continuum mechanics, augmented by species transport equations and chemical source terms. Unlike inert fluid mechanics, the governing equations for reactive flows include production and destruction rates for each chemical species, and the energy equation contains a volumetric heat release term that drives the entire dynamics.
We consider a mixture of \(N\) chemical species with mass fractions \(Y_i\), \(i = 1, \ldots, N\), where \(Y_i = \rho_i / \rho\) is the ratio of the partial density of species \(i\) to the mixture density \(\rho\). The mass fractions satisfy the constraint \(\sum_{i=1}^N Y_i = 1\).
Definition 1.1 (Conservation Equations for Reactive Flow). The governing equations for a compressible reactive mixture are:
\[
\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{u}) = 0.
\]\[
\frac{\partial (\rho \mathbf{u})}{\partial t} + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u}) = -\nabla p + \nabla \cdot \boldsymbol{\tau} + \rho \mathbf{g}.
\]\[
\frac{\partial (\rho Y_i)}{\partial t} + \nabla \cdot (\rho Y_i \mathbf{u}) = -\nabla \cdot \mathbf{J}_i + \dot{\omega}_i, \quad i = 1, \ldots, N.
\]\[
\frac{\partial (\rho e)}{\partial t} + \nabla \cdot (\rho e \mathbf{u}) = -\nabla \cdot \mathbf{q} - p \nabla \cdot \mathbf{u} + \boldsymbol{\tau} : \nabla \mathbf{u} + \dot{Q},
\]
where \(\boldsymbol{\tau}\) is the viscous stress tensor, \(\mathbf{J}_i\) is the diffusive mass flux of species \(i\), \(\dot{\omega}_i\) is the mass rate of production of species \(i\) by chemical reaction, \(\mathbf{q}\) is the heat flux, and \(\dot{Q}\) is the volumetric heat release rate.
The species diffusive fluxes are commonly modeled by Fick’s law, \(\mathbf{J}_i = -\rho D_i \nabla Y_i\), where \(D_i\) is the mass diffusivity of species \(i\) into the mixture. The heat flux follows Fourier’s law, \(\mathbf{q} = -\lambda \nabla T\), where \(\lambda\) is the thermal conductivity. These constitutive models, while approximate for multicomponent mixtures, capture the essential physics of combustion at a level suitable for asymptotic and analytical work. For detailed transport, the full Stefan-Maxwell equations must be used, but these complications are deferred to computational treatments.
1.2 Arrhenius Reaction Rate and the Damköhler Number
The chemical source term \(\dot{\omega}_i\) is what distinguishes combustion from inert fluid mechanics. For a single irreversible reaction of the form \(\text{Fuel} + \text{Oxidizer} \to \text{Products} + \text{Heat}\), the overall reaction rate is typically modeled using the Arrhenius law, one of the most important empirical laws in chemistry. Svante Arrhenius proposed in 1889 that the rate constant depends exponentially on the inverse of temperature, reflecting the fact that only molecules with sufficient kinetic energy can overcome the activation energy barrier and react.
Definition 1.2 (Arrhenius Reaction Rate). For a single-step irreversible reaction with fuel mass fraction \(Y_F\) and oxidizer mass fraction \(Y_O\), the Arrhenius reaction rate is
\[
\dot{\omega} = B \rho^{n} Y_F^{a} Y_O^{b} \exp\left(-\frac{E_a}{R_u T}\right),
\]
where \(B\) is the pre-exponential factor (frequency factor), \(n\) is the overall reaction order (\(n = a + b\)), \(E_a\) is the activation energy, \(R_u\) is the universal gas constant, and \(T\) is the absolute temperature.
The exponential factor \(\exp(-E_a / R_u T)\) is extraordinarily sensitive to temperature. For typical hydrocarbon flames, \(E_a / R_u \approx 15{,}000\text{--}25{,}000\) K, so that a modest increase in temperature (say from 1500 K to 2000 K) increases the reaction rate by several orders of magnitude. This extreme temperature sensitivity is the source of both the power and the mathematical difficulty of combustion theory: it creates thin reaction zones (flame sheets), sharp fronts, and stiff differential equations.
Definition 1.3 (Damköhler Number). The Damköhler number is the ratio of the flow time scale \(\tau_{\text{flow}}\) to the chemical time scale \(\tau_{\text{chem}}\):
\[
\text{Da} = \frac{\tau_{\text{flow}}}{\tau_{\text{chem}}}.
\]
When \(\text{Da} \gg 1\), chemistry is fast compared to transport and the reaction zone is thin. When \(\text{Da} \ll 1\), transport dominates and the mixture is nearly homogeneous (well-stirred reactor limit).
The Damköhler number is the central dimensionless group in combustion. It controls the regime of burning: large Da gives thin flame sheets amenable to asymptotic analysis, while small Da gives distributed reaction zones requiring different mathematical tools. Much of the asymptotic theory in Chapters 2–5 exploits the large-Da limit, which corresponds physically to the large activation energy limit.
Example 1.4 (Estimating the Damköhler Number). Consider a premixed methane-air flame stabilized on a Bunsen burner with nozzle diameter \(d = 1\) cm and flow velocity \(U = 2\) m/s. The flow time scale is \(\tau_{\text{flow}} = d / U = 5 \times 10^{-3}\) s. The chemical time scale for stoichiometric methane-air at atmospheric pressure is approximately \(\tau_{\text{chem}} \sim (\kappa / S_L^2) \approx 0.2 / (40)^2 \approx 1.25 \times 10^{-4}\) s (using \(\kappa \approx 0.2\) cm\(^2\)/s and \(S_L \approx 40\) cm/s). The Damköhler number is therefore \(\text{Da} = \tau_{\text{flow}} / \tau_{\text{chem}} \approx 40\), which is large enough that the flame is a thin sheet compared to the burner diameter. If one reduces the burner diameter to \(d = 0.5\) mm (a micro-combustor), then \(\text{Da} \approx 2\), and the flame thickness becomes comparable to the device size — the thin flame sheet approximation begins to break down, and finite-rate chemistry effects become important.
1.3 Adiabatic Flame Temperature
When a premixed fuel-air mixture burns in an adiabatic, constant-pressure process, the final temperature is determined entirely by thermochemistry. This adiabatic flame temperature provides a fundamental upper bound on the temperature attainable in a combustion system and serves as a reference scale for nondimensionalization.
Definition 1.5 (Adiabatic Flame Temperature). For a premixed mixture initially at temperature \(T_u\), the adiabatic flame temperature \(T_{ad}\) is defined by the energy balance
\[
T_{ad} = T_u + \frac{Q Y_{F,u}}{c_p},
\]
where \(Q\) is the heat of combustion per unit mass of fuel, \(Y_{F,u}\) is the initial fuel mass fraction, and \(c_p\) is the (constant) specific heat capacity of the mixture.
For a stoichiometric methane-air mixture at standard conditions, \(T_{ad} \approx 2230\) K. For hydrogen-air, \(T_{ad} \approx 2380\) K. These temperatures are high enough to strongly ionize alkali metals in the combustion products and to emit visible radiation, which is why flames glow.
Example 1.6 (Computing \(T_{ad}\) for Lean Methane-Air). Consider a lean methane-air mixture with equivalence ratio \(\phi = 0.7\). At standard conditions, stoichiometric methane-air has \(Y_{F,u}^{\text{st}} \approx 0.055\) and \(Q \approx 50\) MJ/kg of fuel, with \(c_p \approx 1.2\) kJ/(kg\(\cdot\)K) and \(T_u = 300\) K. For a lean mixture, the fuel mass fraction scales as \(Y_{F,u} = \phi \, Y_{F,u}^{\text{st}} \approx 0.0385\). Then
\[
T_{ad} = 300 + \frac{50{,}000 \times 0.0385}{1.2} \approx 300 + 1604 = 1904 \text{ K}.
\]
This is approximately 330 K below the stoichiometric value, illustrating the strong dependence of flame temperature on equivalence ratio. The reduction has practical consequences: the thermal NO\(_x\) production rate (which depends exponentially on temperature through the Zeldovich mechanism) drops by roughly two orders of magnitude, which is why lean-burn strategies are used for emissions control in gas turbine engines.
A remarkable simplification of the reactive flow equations arises when the Lewis number \(\text{Le} = \lambda / (\rho c_p D)\) equals unity, meaning that thermal diffusivity equals mass diffusivity. In this case, certain linear combinations of temperature and species mass fractions satisfy source-free (non-reactive) transport equations. This observation, due independently to Shvab (1948) and Zeldovich (1949), reduces the number of independent equations that must be solved and is the starting point for much of classical flame theory.
Definition 1.8 (Shvab-Zeldovich Variables). Define the coupling functions (or Shvab-Zeldovich variables) by
\[
\beta_1 = \frac{T - T_u}{T_{ad} - T_u} + \frac{Y_F}{Y_{F,u}}, \qquad \beta_2 = \frac{Y_O}{\nu Y_{F,u}} - \frac{Y_F}{Y_{F,u}},
\]
where \(\nu\) is the stoichiometric mass ratio of oxidizer to fuel. When \(\text{Le} = 1\), both \(\beta_1\) and \(\beta_2\) satisfy the sourceless convection-diffusion equation
\[
\rho \frac{\partial \beta}{\partial t} + \rho \mathbf{u} \cdot \nabla \beta = \nabla \cdot (\rho D \nabla \beta).
\]
The physical content of the Shvab-Zeldovich formulation is that the chemical source term, which makes the species and energy equations nonlinear and coupled, can be eliminated from all but one equation. The surviving equation retains the Arrhenius nonlinearity, but the remaining coupling functions are governed by linear (or at least source-free) equations. This decomposition is the mathematical foundation of the Burke-Schumann flame sheet model for diffusion flames and of the Zeldovich-Frank-Kamenetskii analysis for premixed flames.
Example 1.9 (Dimensionless One-Step Model). Consider a one-dimensional, steady, premixed flame with constant density, unit Lewis number, and a single irreversible reaction. Nondimensionalize temperature as \(\theta = (T - T_u)/(T_{ad} - T_u)\) and fuel mass fraction as \(y = Y_F / Y_{F,u}\). The steady equations in the flame-attached frame become
\[
\frac{d\theta}{dx} = \frac{d^2\theta}{dx^2} + \text{Da} \cdot y \exp\left(\frac{-\beta(1 - \theta)}{1 - \alpha(1 - \theta)}\right),
\]
\[
\frac{dy}{dx} = \frac{d^2 y}{dx^2} - \text{Da} \cdot y \exp\left(\frac{-\beta(1 - \theta)}{1 - \alpha(1 - \theta)}\right),
\]
where \(\beta = E_a(T_{ad} - T_u)/(R_u T_{ad}^2)\) is the Zeldovich number and \(\alpha = (T_{ad} - T_u)/T_{ad}\) is the heat release parameter. Adding the two equations yields \((\theta + y)'' - (\theta + y)' = 0\), confirming that the coupling function \(\theta + y\) satisfies a source-free equation. With boundary conditions \(\theta + y \to 1\) as \(x \to -\infty\), one obtains \(\theta + y = 1\) everywhere, reducing the system to a single equation for \(\theta\).
This reduction to a single nonlinear ODE is the starting point for the asymptotic analysis of premixed flame structure in Chapter 2.
Chapter 2: Premixed Laminar Flames
2.1 The Flame as a Traveling Wave
A premixed laminar flame is a self-sustaining combustion wave that propagates through a homogeneous mixture of fuel and oxidizer. Mathematically, it is a traveling wave solution of a reaction-diffusion system. The flame converts cold, unburned reactants ahead of it into hot, burned products behind it, and propagates at a definite speed — the laminar flame speed — that is an intrinsic property of the mixture.
The existence and speed of such traveling waves is one of the classical problems of applied mathematics. It connects combustion theory to the broader theory of traveling waves in reaction-diffusion equations, including the KPP-Fisher equation from mathematical biology. The key difference is that combustion waves have an ignition-temperature nonlinearity (the Arrhenius term is negligible at low temperatures), which selects a unique wave speed rather than a continuous family.
Definition 2.1 (Premixed Flame as Traveling Wave). A premixed laminar flame is a traveling wave solution \(\theta(x - S_L t)\) of the reaction-diffusion equation
\[
\frac{\partial \theta}{\partial t} = \kappa \frac{\partial^2 \theta}{\partial x^2} + f(\theta),
\]
where \(\theta\) is the nondimensional temperature, \(\kappa\) is the thermal diffusivity, and \(f(\theta)\) is the reaction rate. In the traveling wave frame \(\xi = x - S_L t\), the equation becomes
\[
\kappa \theta'' + S_L \theta' + f(\theta) = 0,
\]
with boundary conditions \(\theta(-\infty) = 0\) (cold unburned gas) and \(\theta(+\infty) = 1\) (hot burned gas), where primes denote \(d/d\xi\).
The laminar flame speed \(S_L\) appears as an eigenvalue of this boundary value problem. For the Arrhenius nonlinearity \(f(\theta) = \text{Da} \cdot (1 - \theta) \exp(-\beta(1 - \theta)/(1 - \alpha(1-\theta)))\), the traveling wave exists and is unique (up to translation). The eigenvalue nature of \(S_L\) is a reflection of the fact that the flame speed is not imposed externally but is determined by the internal balance between diffusion and reaction.
2.2 The KPP-Fisher Equation and Combustion
The connection between combustion and mathematical biology deserves a brief digression. Kolmogorov, Petrovsky, and Piskunov (1937) and, independently, Fisher (1937) studied the equation \(\theta_t = \kappa \theta_{xx} + f(\theta)\) with \(f(0) = f(1) = 0\), \(f(\theta) > 0\) for \(\theta \in (0,1)\), as a model for gene propagation. They proved that traveling wave solutions exist for all speeds \(S_L \geq 2\sqrt{\kappa f'(0)}\) and that the minimum speed is selected by initial conditions with compact support.
Theorem 2.2 (KPP Minimum Speed Selection). For the reaction-diffusion equation \(\theta_t = \kappa \theta_{xx} + f(\theta)\) with \(f\) of KPP type (i.e., \(f(0) = 0\), \(f(1) = 0\), \(f(\theta) > 0\) for \(0 < \theta < 1\), and \(f(\theta) \leq f'(0)\theta\)), traveling wave solutions \(\theta(\xi) = \theta(x - S_L t)\) connecting \(\theta = 0\) to \(\theta = 1\) exist if and only if
\[
S_L \geq S_{\min} = 2\sqrt{\kappa f'(0)}.
\]
Moreover, compactly supported initial data evolve toward the traveling wave with speed \(S_{\min}\).
Proof. (Sketch.) In the traveling wave ODE \(\kappa \theta'' + S_L \theta' + f(\theta) = 0\), linearize near the cold boundary \(\theta = 0\) by setting \(f(\theta) \approx f'(0)\theta\). The characteristic equation is \(\kappa r^2 + S_L r + f'(0) = 0\), with roots
\[
r = \frac{-S_L \pm \sqrt{S_L^2 - 4\kappa f'(0)}}{2\kappa}.
\]
For the solution to remain non-negative, the roots must be real, requiring \(S_L^2 \geq 4\kappa f'(0)\), i.e., \(S_L \geq 2\sqrt{\kappa f'(0)}\). When \(S_L = S_{\min}\), the roots coincide and the approach to \(\theta = 0\) is algebraically slower (critical damping), which is selected by steep initial data. \(\square\)
In combustion, however, the Arrhenius nonlinearity has a crucial difference: \(f'(0) = 0\) (or effectively zero for large activation energy) because the reaction rate is negligibly small at the cold boundary temperature. This means the KPP minimum speed formula gives \(S_{\min} = 0\), which is physically meaningless. The Arrhenius nonlinearity is of “ignition type,” and the wave speed is uniquely selected rather than being the minimum of a continuous family. Determining this unique speed is the central problem of premixed flame theory.
2.3 The Zeldovich-Frank-Kamenetskii (ZFK) Analysis
The first successful determination of the laminar flame speed was achieved independently by Zeldovich and Frank-Kamenetskii in 1938, using what we would now call matched asymptotic expansions in the large activation energy limit. Their analysis exploits the fact that for large Zeldovich number \(\beta \gg 1\), the reaction rate is confined to a thin layer near the maximum temperature, and the flame structure naturally separates into a broad preheat zone (where diffusion balances convection) and a narrow reaction zone (where diffusion balances reaction).
Definition 2.4 (Zeldovich Number). The Zeldovich number is the dimensionless activation energy
\[
\beta = \frac{E_a (T_{ad} - T_u)}{R_u T_{ad}^2}.
\]
It measures the sensitivity of the reaction rate to temperature changes near the flame temperature. For typical hydrocarbon-air flames, \(\beta \approx 8\text{--}15\).
The Zeldovich number is the large parameter in the asymptotic expansion. When \(\beta \gg 1\), the Arrhenius exponential acts like a delta function concentrated at the adiabatic flame temperature, and the flame has a well-defined two-zone structure.
Theorem 2.5 (ZFK Flame Speed Formula). In the limit \(\beta \to \infty\) with all other parameters fixed, the laminar flame speed for the one-step irreversible Arrhenius model satisfies
\[
S_L^2 = \frac{2 \kappa B}{\beta^2} \exp\left(-\frac{E_a}{R_u T_{ad}}\right) \cdot I(\alpha, n),
\]
where \(I(\alpha, n)\) is an order-one integral that depends on the heat release parameter \(\alpha\) and the reaction order \(n\). For a first-order reaction with \(\alpha = O(1)\),
\[
S_L^2 \sim \frac{2 \kappa B}{\beta^2} \exp\left(-\frac{E_a}{R_u T_{ad}}\right).
\]
Proof. We construct the solution by matched asymptotic expansions.
Outer region (preheat zone): For \(\xi < 0\) (ahead of the flame), the reaction rate is exponentially small and the temperature satisfies the convection-diffusion equation \(\kappa \theta'' + S_L \theta' = 0\), whose solution is \(\theta = \exp(S_L \xi / \kappa)\).
Outer region (burned gas): For \(\xi > 0\) (behind the flame), all fuel is consumed and \(\theta = 1\).
\[
f(\theta) = \text{Da} \cdot \frac{\phi}{\beta} \exp\left(-\frac{\beta(1/\beta)\phi}{1 - \alpha\phi/\beta}\right) \approx \text{Da} \cdot \frac{\phi}{\beta} \exp(-\phi),
\]\[
\frac{d^2\phi}{d\eta^2} = -\Lambda \phi \, e^{-\phi},
\]\[
\frac{1}{2}\left(\frac{d\phi}{d\eta}\right)^2\bigg|_{\eta \to +\infty} = \Lambda \int_0^\infty \phi \, e^{-\phi} \, d\phi = \Lambda.
\]\[
\frac{S_L}{\kappa} = \frac{\beta}{\delta} \sqrt{2\Lambda},
\]
and determining \(\delta\) and \(\Lambda\) in terms of original parameters gives the stated formula. \(\square\)
The ZFK formula reveals the essential physics: the flame speed is exponentially sensitive to the adiabatic flame temperature (through the Arrhenius factor) and algebraically sensitive to the activation energy (through the \(\beta^{-2}\) prefactor, which reflects the thinness of the reaction zone). For stoichiometric methane-air, this gives \(S_L \approx 40\) cm/s, in good agreement with experiment.
Lemma 2.6 (Monotonicity of the Flame Profile). In the traveling wave frame, the temperature profile \(\theta(\xi)\) of the ZFK flame is strictly monotonically increasing from \(\theta = 0\) at \(\xi = -\infty\) to \(\theta = 1\) at \(\xi = +\infty\). The fuel mass fraction \(y(\xi) = 1 - \theta(\xi)\) is correspondingly monotonically decreasing. There are no overshoots, undershoots, or oscillations.
Proof. From the phase-plane formulation, the traveling wave ODE \(\kappa \theta'' + S_L \theta' + f(\theta) = 0\) can be written as the system \(\theta' = p\), \(\kappa p' = -S_L p - f(\theta)\). In the phase plane \((\theta, p)\), the trajectory connects the saddle point \((0, 0)\) to the stable node \((1, 0)\). Since \(f(\theta) \geq 0\) and \(S_L > 0\), the equation \(\kappa p' = -S_L p - f(\theta)\) ensures that \(p\) cannot change sign along the trajectory (if \(p = 0\) at an interior point, then \(p' = -f(\theta)/\kappa < 0\), contradicting the requirement that the trajectory reach \((1, 0)\) from above). Hence \(p = \theta' > 0\) throughout. \(\square\)
This monotonicity is a consequence of the maximum principle for parabolic equations and is specific to the one-dimensional, single-step reaction model. Multi-step chemistry and multi-dimensional effects can produce non-monotone temperature profiles (e.g., superadiabatic temperatures in rich hydrogen flames due to differential diffusion).
Example 2.8 (Flame Thickness Estimate). The preheat zone thickness is \(\delta_f \sim \kappa / S_L\). For a methane-air flame with \(\kappa \approx 0.2\) cm\(^2\)/s and \(S_L \approx 40\) cm/s, we obtain \(\delta_f \approx 0.05\) mm. The reaction zone is thinner by a factor of \(\beta \approx 10\), giving a reaction zone thickness of \(\sim 5 \, \mu\)m. This extreme thinness of the reaction zone is what makes the flame sheet approximation so effective: the reaction zone is indeed negligibly thin compared to any macroscopic scale in the flow.
Example 2.9 (Effect of Pressure on Flame Speed). The ZFK formula predicts the pressure dependence of the laminar flame speed. For a second-order reaction (\(n = 2\)), the pre-exponential factor scales as \(B \rho^2 \propto p^2\), while the thermal diffusivity scales as \(\kappa \propto 1/p\) (since \(\kappa = \lambda / (\rho c_p)\) and \(\lambda\) is nearly pressure-independent). Therefore
\[
S_L^2 \propto \kappa \cdot B \rho^n \propto p^{-1} \cdot p^2 = p,
\]
giving \(S_L \propto p^{1/2}\) for a second-order reaction. More generally, \(S_L \propto p^{(n-2)/2}\). For methane-air, the experimental pressure exponent is approximately \(-0.5\), suggesting an effective overall reaction order of about \(n \approx 1\). This negative pressure exponent means that high-altitude flames are faster (per unit mass) but weaker (per unit volume) than sea-level flames — a result of practical importance in aircraft engine design.
2.4 Flame Structure via Matched Asymptotics
The formal construction of the ZFK solution illustrates the method of matched asymptotic expansions in a physically compelling context. The outer solutions describe the preheat zone and the equilibrium zone, while the inner solution resolves the thin reaction zone. The matching conditions enforce continuity and determine the eigenvalue (flame speed).
Example 2.11 (Comparison of ZFK Prediction with Experiment). For a stoichiometric methane-air flame at 1 atm, the ZFK formula (Theorem 2.5) with parameters \(\kappa = 0.22\) cm\(^2\)/s, \(B = 2.1 \times 10^{14}\) s\(^{-1}\), \(E_a/R_u = 24{,}400\) K, and \(T_{ad} = 2230\) K gives
\[
S_L = \sqrt{\frac{2 \times 0.22 \times 2.1 \times 10^{14}}{10^2} \exp(-24{,}400/2230)} \approx \sqrt{9.24 \times 10^{9} \times 1.85 \times 10^{-5}} \approx 41 \text{ cm/s}.
\]
The measured value is \(S_L \approx 36\)--\(40\) cm/s (depending on the experimental technique and the precise definition of flame speed). The 5--15% discrepancy is typical of the ZFK formula applied with fitted global kinetic parameters and illustrates both the power of the asymptotic approach (the correct order of magnitude and functional dependences are captured) and its limitations (quantitative accuracy requires careful calibration or multi-step chemistry).
Chapter 3: Activation Energy Asymptotics
3.1 The Large Activation Energy Limit
Activation energy asymptotics (AEA) is the systematic exploitation of the limit \(\beta \to \infty\) to simplify the nonlinear reaction-diffusion equations of combustion. The Arrhenius exponential, when \(\beta\) is large, acts as a “switch” that is essentially zero below a critical temperature and extremely large near the flame temperature. The mathematical consequence is that the reaction is confined to a thin zone whose thickness scales as \(1/\beta\) relative to the diffusion scale, and the rest of the domain is chemically inert.
Definition 3.1 (Activation Energy Asymptotics). The method of
activation energy asymptotics consists of:
- Expanding the solution in powers of \(1/\beta\) (or \(\varepsilon = 1/\beta\)) in the outer (chemically frozen) regions.
- Introducing a stretched inner variable \(\eta = (\xi - \xi_f)/\varepsilon\) near the reaction sheet location \(\xi_f\).
- Solving the inner and outer problems order by order and matching them in the overlap region.
The leading-order outer solution gives the flame sheet approximation, and the inner solution resolves the reaction zone structure.
The AEA approach is due primarily to the work of Yakov Zeldovich and David Frank-Kamenetskii in the Soviet combustion school, and was elevated to a systematic asymptotic methodology by Amable Liñán in his seminal 1974 Acta Astronautica paper on diffusion flames, and by John Buckmaster and Gregory Ludford in their rigorous mathematical treatments.
3.2 Distinguished Limits
When multiple small or large parameters are present, the asymptotic structure depends on how they are related. In combustion, the Damköhler number, the Zeldovich number, and the Lewis number departure from unity all compete, and different scalings lead to different flame behaviors. These different scalings are called distinguished limits.
Definition 3.3 (Distinguished Limit). A
distinguished limit is a scaling relationship between two or more parameters such that competing effects remain in balance at leading order in the asymptotic expansion. In combustion, the principal distinguished limits relate the Damköhler number \(\text{Da}\) to the Zeldovich number \(\beta\):
- Premixed flame regime: \(\text{Da} = O(\beta^2 \exp(\beta))\), ensuring the flame speed is \(O(1)\).
- Near-extinction regime: \(\text{Da} = O(\beta^2 \exp(\beta)) \cdot (1 - \delta/\beta)\) with \(\delta = O(1)\).
- Ignition regime: \(\text{Da} = O(1)\), so that the Damköhler number is not exponentially large.
Each distinguished limit leads to a different canonical inner problem. The premixed flame regime recovers the ZFK analysis. The near-extinction regime introduces a turning point in the inner problem, leading to extinction. The ignition regime gives the Semenov or Frank-Kamenetskii thermal explosion problems. This classification, due to Liñán, provides a complete map of the parameter space for one-step reaction-diffusion systems.
Example 3.4 (Identifying the Distinguished Limit). Consider a counterflow diffusion flame with \(\beta = 10\), strain rate \(a\), and a pre-exponential factor \(B\). The raw Damköhler number is \(\text{Da} = B/a\). In the premixed flame distinguished limit, we write \(\text{Da} = C \beta^2 e^{\beta}\) with \(C = O(1)\). For \(\beta = 10\), this gives \(\text{Da} = 100 \, e^{10} \cdot C \approx 2.2 \times 10^6 \, C\). If the physical Damköhler number is \(\text{Da} = 10^7\), then \(C \approx 4.5\), confirming that this flame is in the premixed regime. Near extinction, \(C\) drops toward a critical value of order unity, and below that the flame cannot exist. The practical value of this classification is immediate: it tells the analyst which inner equation to solve and what physical behavior to expect, without performing a full numerical computation.
3.3 Inner and Outer Expansions
We now carry out the formal asymptotic expansion for a general one-dimensional steady flame, illustrating the systematic procedure that applies to all activation energy asymptotic problems.
Theorem 3.5 (Outer Solution Structure). Consider the steady reaction-diffusion equation
\[
\theta'' + c \, \theta' + \text{Da} \cdot (1 - \theta) \exp\left(\frac{-\beta(1-\theta)}{1 - \alpha(1-\theta)}\right) = 0,
\]
on \(-\infty < x < \infty\) with \(\theta(-\infty) = 0\), \(\theta(+\infty) = 1\). In the limit \(\beta \to \infty\), the outer solution (away from the reaction zone at \(x = 0\)) is, to all algebraic orders in \(\beta\),
\[
\theta^{(-)}_{\text{outer}}(x) = e^{cx}, \quad x < 0, \qquad \theta^{(+)}_{\text{outer}}(x) = 1, \quad x > 0.
\]
The reaction rate is transcendentally small (\(O(e^{-\beta})\)) in the outer regions and is concentrated entirely in the inner region.
Proof. In the outer region \(x < 0\), we have \(\theta < 1\) and \(1 - \theta = O(1)\), so the Arrhenius factor satisfies \(\exp(-\beta(1-\theta)/\ldots) \leq \exp(-c'\beta)\) for some \(c' > 0\). This is beyond all algebraic orders in \(\beta\), so the reaction term vanishes from the outer equation. The remaining equation \(\theta'' + c\theta' = 0\) with \(\theta(-\infty) = 0\) gives \(\theta = Ae^{cx}\). Matching to \(\theta \to 1\) as \(x \to 0^-\) gives \(A = 1\). Similarly, in \(x > 0\), the fuel is fully consumed and \(\theta = 1\) identically. \(\square\)
The outer solution has a discontinuous first derivative at \(x = 0\): the slope jumps from \(c\) (on the left) to \(0\) (on the right). This jump is produced by the delta-function-like reaction rate in the inner zone and carries information about the total heat release. The inner solution must smooth out this jump over the thin reaction zone.
Example 3.6 (Matching at the Reaction Sheet). To illustrate the matching procedure concretely, consider the outer preheat-zone solution \(\theta^{(-)}(x) = e^{cx}\) for \(x < 0\). As \(x \to 0^-\), we have \(\theta^{(-)} \to 1\) and \(d\theta^{(-)}/dx \to c\). In the inner zone, set \(\eta = \beta x\) and \(\theta = 1 - \phi/\beta\). The inner temperature deficit \(\phi\) must satisfy \(\phi \to 0\) as \(\eta \to -\infty\) (matching onto the burned side) and
\[
\frac{d\phi}{d\eta} \to -\frac{c}{\beta} \quad \text{as} \quad \eta \to +\infty,
\]
which, after rescaling, becomes the condition that the inner gradient matches the outer preheat-zone slope. This matching condition is what determines the eigenvalue \(c\) (and hence the flame speed \(S_L\)) in terms of the Damköhler number and activation energy. The inner equation itself,
\[
\frac{d^2\phi}{d\eta^2} = -\Lambda \phi \, e^{-\phi},
\]
is a second-order ODE that can be reduced to first order by the standard energy integral method. Multiplying both sides by \(d\phi/d\eta\) and integrating from \(\eta = -\infty\) to \(+\infty\) yields the eigenvalue relation \(\Lambda = c^2/(2\beta^2)\), completing the determination of the flame speed.
3.4 The Flame Sheet Approximation and Delta-Function Model
At leading order in the large-\(\beta\) limit, the flame reduces to a surface (a point in one dimension) across which temperature is continuous but its gradient jumps. The reaction rate, integrated across the thin reaction zone, appears as a delta function source in the outer equations. This is the flame sheet approximation, which is the workhorse of analytical combustion theory.
Definition 3.8 (Flame Sheet Approximation). In the flame sheet approximation, the governing equations outside the reaction zone are the source-free transport equations, and the effect of the reaction is captured by jump conditions across the flame sheet \(\Sigma\):
\[
[\theta]_{\Sigma} = 0, \qquad \left[\frac{\partial \theta}{\partial n}\right]_{\Sigma} = -\frac{S_L}{\kappa},
\]
where \([f]_{\Sigma} = f^+ - f^-\) denotes the jump across \(\Sigma\) and \(n\) is the normal to the flame sheet. The fuel mass fraction satisfies \(Y_F = 0\) on the burned side of the sheet.
3.5 Liñán’s Diffusion Flame Analysis
Amable Liñán’s 1974 analysis of the counterflow diffusion flame is one of the masterpieces of asymptotic combustion theory. He identified four distinct regimes — the premixed flame regime, the diffusion flame regime, the partial burning regime, and the ignition regime — each corresponding to a different distinguished limit of the Damköhler number. We shall return to the diffusion flame aspects in Chapter 4, but the mathematical structure of the inner problem is best introduced here.
Theorem 3.10 (Liñán's Canonical Inner Problem). In the diffusion flame regime with \(\beta \gg 1\), the inner reaction-zone equation, after appropriate stretching and nondimensionalization, takes the canonical form
\[
\frac{d^2\phi}{d\eta^2} = \delta_L(\phi - \gamma|\eta|) \exp(-\phi),
\]
where \(\phi\) is the rescaled temperature perturbation, \(\eta\) is the stretched inner coordinate, \(\gamma\) is the stoichiometric parameter measuring the asymmetry of the flame, and \(\delta_L\) is the reduced Damköhler number. The boundary conditions are \(\phi \sim \gamma|\eta|\) as \(\eta \to \pm\infty\).
The parameter \(\delta_L\) (Liñán’s reduced Damköhler number) serves as the bifurcation parameter. For \(\delta_L\) above a critical value \(\delta_E(\gamma)\), the problem has a solution with \(\phi\) everywhere bounded — corresponding to a burning diffusion flame. Below \(\delta_E(\gamma)\), no such solution exists — the flame has been extinguished. The transition is a saddle-node bifurcation, and the turning point corresponds to the extinction Damköhler number.
Example 3.11 (Symmetric Diffusion Flame, \(\gamma = 0\)). When \(\gamma = 0\) (equal diffusivities and stoichiometric symmetry), Liñán's equation simplifies to
\[
\frac{d^2\phi}{d\eta^2} = \delta_L \, \phi \, \exp(-\phi).
\]
This is equivalent (by the substitution \(\psi = d\phi/d\eta\)) to the phase-plane system \(\phi' = \psi\), \(\psi' = \delta_L \phi e^{-\phi}\). The extinction condition is obtained by requiring that the integral
\[
\int_0^{\phi_{\max}} \phi \, e^{-\phi} \, d\phi = \frac{1}{2\delta_L}\left(\frac{d\phi}{d\eta}\bigg|_{\eta \to \infty}\right)^2
\]
be satisfied simultaneously with the matching conditions. For the symmetric case, extinction occurs at \(\delta_E \approx e = 2.718\ldots\), a result noted by Liñán as "a happy numerical coincidence."
Chapter 4: Diffusion Flames
4.1 The Burke-Schumann Flame Sheet Model
Diffusion flames occur when fuel and oxidizer are initially separated and must mix by molecular diffusion before they can react. The canonical example is a candle: the fuel (wax vapor) issues from the wick, and the oxidizer (oxygen in air) surrounds it. The flame sits at the surface where fuel and oxidizer meet in stoichiometric proportions. In the infinite-rate chemistry limit (\(\text{Da} \to \infty\)), fuel and oxidizer cannot coexist and the flame is an infinitely thin sheet — the Burke-Schumann flame sheet.
S.P. Burke and T.E.W. Schumann introduced this model in 1928 in their pioneering paper on diffusion flames. Despite its simplicity, the model captures the essential geometry and structure of non-premixed combustion and remains the starting point for all theoretical work on diffusion flames.
Definition 4.1 (Mixture Fraction). The mixture fraction \(Z\) is a conserved scalar defined so that \(Z = 1\) in pure fuel and \(Z = 0\) in pure oxidizer. For a one-step reaction with \(\text{Le} = 1\), it satisfies the source-free equation
\[
\rho \frac{\partial Z}{\partial t} + \rho \mathbf{u} \cdot \nabla Z = \nabla \cdot (\rho D \nabla Z).
\]
The mixture fraction is related to the Shvab-Zeldovich variable \(\beta_2\) by a linear transformation.
Definition 4.2 (Stoichiometric Surface). The stoichiometric surface is the level set
\[
Z(\mathbf{x}, t) = Z_{st},
\]
where \(Z_{st} = 1/(1 + \phi)\) and \(\phi = \nu Y_{F,\infty} / Y_{O,\infty}\) is the equivalence ratio. In the Burke-Schumann limit, the diffusion flame sits on this surface.
The mixture fraction formulation is the foundation of modern computational approaches to non-premixed combustion. Because \(Z\) satisfies a linear equation, its statistics can be computed (or modeled) independently of the chemistry. The flame structure is then obtained by mapping from \(Z\)-space to physical space. This is the basis of the flamelet model, discussed in Chapters 7 and 8.
Example 4.3 (Computing the Stoichiometric Mixture Fraction). For methane (\(\text{CH}_4\)) burning in air, the stoichiometric reaction is \(\text{CH}_4 + 2\text{O}_2 \to \text{CO}_2 + 2\text{H}_2\text{O}\). The stoichiometric mass ratio is \(\nu = 2 \times 32 / 16 = 4\) (4 kg of \(\text{O}_2\) per kg of \(\text{CH}_4\)). With air containing 23.2% oxygen by mass, the stoichiometric air-to-fuel ratio is \(4 / 0.232 \approx 17.2\). The stoichiometric mixture fraction is
\[
Z_{st} = \frac{1}{1 + \nu Y_{F,\infty}/Y_{O,\infty}} = \frac{1}{1 + 4 \times 1 / 0.232} \approx 0.055.
\]
This small value means that the flame sits close to the oxidizer side of the mixing layer — a universal feature of hydrocarbon-air diffusion flames. For hydrogen (\(\nu = 8\)), \(Z_{st} \approx 0.028\), even closer to the air side, which is why hydrogen diffusion flames are extremely thin and susceptible to local extinction by turbulent straining.
Theorem 4.4 (Burke-Schumann Solution). In the limit \(\text{Da} \to \infty\) with \(\text{Le} = 1\), the temperature and species fields are piecewise linear functions of the mixture fraction:
\[
\theta(Z) = \begin{cases} \dfrac{Z}{Z_{st}}, & 0 \leq Z \leq Z_{st}, \\[6pt] \dfrac{1 - Z}{1 - Z_{st}}, & Z_{st} \leq Z \leq 1, \end{cases}
\]
\[
Y_F(Z) = \begin{cases} 0, & 0 \leq Z \leq Z_{st}, \\[6pt] Y_{F,\infty} \dfrac{Z - Z_{st}}{1 - Z_{st}}, & Z_{st} \leq Z \leq 1, \end{cases}
\]
\[
Y_O(Z) = \begin{cases} Y_{O,\infty} \dfrac{Z_{st} - Z}{Z_{st}}, & 0 \leq Z \leq Z_{st}, \\[6pt] 0, & Z_{st} \leq Z \leq 1. \end{cases}
\]
Fuel and oxidizer do not coexist: \(Y_F Y_O = 0\) everywhere.
Proof. Since \(\text{Da} \to \infty\), the reaction rate is infinitely fast and the reaction is limited by the rate at which fuel and oxidizer are brought together by diffusion. The coupling function \(\beta_2 = Y_O / (\nu Y_{F,\infty}) - Y_F / Y_{F,\infty}\) satisfies the same source-free equation as \(Z\) and is therefore a function of \(Z\) alone. With \(\beta_2 = (Z_{st} - Z)/Z_{st}\) and the constraint \(Y_F Y_O = 0\) (since any coexistence at infinite Da would trigger infinite reaction), the fuel and oxidizer profiles follow. The temperature profile is then obtained from the coupling function \(\theta + Y_F / Y_{F,u} = Z/Z_{st}\) on the oxidizer side and \(\theta + Y_O/(\nu Y_{F,u}) = (1-Z)/(1-Z_{st})\) on the fuel side. \(\square\)
The piecewise-linear structure in \(Z\)-space corresponds to a smooth field in physical space (since \(Z\) itself is smooth), but with a kink in the profiles at the flame sheet \(Z = Z_{st}\). The temperature reaches its maximum at the stoichiometric surface, a fundamental feature of diffusion flames.
Example 4.5 (Burke-Schumann Flame Height in a Co-Flow Jet). Consider a round fuel jet of radius \(R\) issuing into a co-flowing air stream. In the Burke-Schumann limit, the flame sits on the surface \(Z(r, x) = Z_{st}\), where \(x\) is the axial distance and \(r\) is the radial coordinate. For a uniform velocity \(U\) and constant diffusivity \(D\), the mixture fraction satisfies
\[
U \frac{\partial Z}{\partial x} = D \left(\frac{\partial^2 Z}{\partial r^2} + \frac{1}{r}\frac{\partial Z}{\partial r}\right),
\]
whose solution is a series of Bessel functions. The flame height — the axial distance at which \(Z = Z_{st}\) on the centerline \(r = 0\) — scales as
\[
L_f \sim \frac{U R^2}{D} \cdot \frac{1}{\ln(1/Z_{st})}.
\]
For methane in air, \(Z_{st} \approx 0.055\), so \(\ln(1/Z_{st}) \approx 2.9\). With \(U = 1\) m/s, \(R = 5\) mm, and \(D = 2 \times 10^{-5}\) m\(^2\)/s, one obtains \(L_f \approx 0.43\) m. The proportionality \(L_f \propto U R^2 / D\) explains why diffusion flame height increases with fuel velocity (unlike premixed flames, whose cone height is set by the flame speed). Burke and Schumann's 1928 prediction of this scaling was confirmed experimentally and remains a cornerstone of non-premixed combustion theory.
4.2 Counterflow Diffusion Flame
The counterflow configuration is the most important canonical geometry for studying diffusion flame structure. Two opposing streams — one of fuel, one of oxidizer — impinge on each other, creating a stagnation-point flow. The flame forms in the mixing layer between the streams. The great virtue of this geometry is that the problem reduces to a system of ODEs in the single variable \(Z\) (or equivalently in the coordinate normal to the flame).
Definition 4.7 (Counterflow Configuration). In the counterflow diffusion flame, fuel and oxidizer approach from opposite sides of a stagnation plane with strain rate \(a\) (units s\(^{-1}\)). The flow field is \(\mathbf{u} = (-ax, ay)\) in two dimensions (or the axisymmetric analogue). The mixture fraction satisfies
\[
-ax \frac{dZ}{dx} = D \frac{d^2Z}{dx^2},
\]
whose solution is \(Z(x) = \frac{1}{2}\text{erfc}(x / \sqrt{2D/a})\). The flame thickness scales as \(\delta_f \sim \sqrt{D/a}\).
The strain rate \(a\) controls the Damköhler number: \(\text{Da} \sim (B/a) \exp(-E_a / R_u T_{ad})\). As the strain rate increases (stronger flow), the Damköhler number decreases, the flame weakens, and eventually extinction occurs. The counterflow experiment is the standard method for measuring laminar flame properties and extinction strain rates.
Example 4.9 (Extinction Strain Rate). For a methane-air counterflow flame, typical parameters give an extinction strain rate of \(a_{\text{ext}} \approx 400\) s\(^{-1}\), corresponding to \(\text{Da}_{\text{ext}} = O(1)\) (in the appropriately rescaled sense). Using Liñán's analysis from Chapter 3 with \(\beta \approx 10\), the extinction condition is \(\delta_L = \delta_E(\gamma)\), which translates to
\[
a_{\text{ext}} \sim \frac{B}{\beta^3} \exp\left(-\frac{E_a}{R_u T_{ad}}\right) \cdot g(\gamma, Z_{st}),
\]
where \(g\) is a known function of the stoichiometric parameters. The \(\beta^3\) factor (rather than \(\beta^2\) as in the premixed case) arises because the diffusion flame reaction zone is thinner by an additional factor of \(\beta^{-1/2}\) compared to the premixed flame reaction zone.
4.3 The S-Curve: Flame Response Diagram
The relationship between the maximum flame temperature and the Damköhler number traces out a characteristic S-shaped curve when plotted, revealing the multiplicity structure of diffusion flame solutions. This S-curve is one of the most important conceptual diagrams in combustion theory.
Theorem 4.10 (S-Curve Structure). For a diffusion flame with large but finite activation energy, the maximum temperature \(T_{\max}\) as a function of the Damköhler number \(\text{Da}\) exhibits an S-shaped curve with three branches:
- Upper branch (burning): A strongly burning flame with \(T_{\max}\) close to \(T_{ad}\). This branch exists for \(\text{Da} > \text{Da}_E\) (the extinction Damköhler number) and is stable.
- Middle branch (unstable): An intermediate solution that is unstable to perturbations. This branch connects the turning points.
- Lower branch (frozen): A weakly reacting or non-reacting state with \(T_{\max}\) close to the mixing temperature. This branch exists for \(\text{Da} < \text{Da}_I\) (the ignition Damköhler number) and is stable.
The turning points \(\text{Da}_E\) (extinction) and \(\text{Da}_I\) (ignition) are saddle-node bifurcations.
The S-curve explains the hysteresis observed in combustion: a flame that is already burning can survive at Damköhler numbers below the ignition threshold, and an unburned mixture requires a Damköhler number above the ignition threshold to be lit. Between \(\text{Da}_E\) and \(\text{Da}_I\), both the burning and frozen states are stable, and the system is bistable. This bistability is the mathematical basis of ignition and extinction phenomena.
Chapter 5: Ignition and Extinction
5.1 Thermal Explosion Theory
The theory of thermal explosions predates modern combustion theory and addresses the most basic question: under what conditions will an exothermic reaction run away to thermal explosion? The competition is between heat generation (which increases exponentially with temperature through the Arrhenius law) and heat loss (which increases at most linearly with temperature through Newtonian cooling or conduction). When heat generation exceeds heat loss at all temperatures, no steady state exists and a thermal explosion occurs.
Definition 5.1 (Semenov Model). The Semenov model of thermal explosion considers a spatially uniform reactive system with Newtonian heat loss:
\[
\rho c_v \frac{dT}{dt} = Q B \exp\left(-\frac{E_a}{R_u T}\right) - \frac{hS}{V}(T - T_w),
\]
where \(h\) is the heat transfer coefficient, \(S/V\) is the surface-to-volume ratio, and \(T_w\) is the wall temperature. Steady states satisfy
\[
Q B \exp\left(-\frac{E_a}{R_u T}\right) = \frac{hS}{V}(T - T_w).
\]
Nikolai Semenov’s (1928) analysis of this algebraic equation is elegant in its simplicity. The heat generation curve (left side) is a convex exponential, while the heat loss line (right side) is a straight line through \(T_w\). When the line is tangent to the exponential, the system is at the critical condition for explosion. For subcritical heat loss, two steady states exist (a stable low-temperature state and an unstable intermediate state); for supercritical heat loss, no steady state exists and the system explodes.
Theorem 5.2 (Semenov Criticality Condition). In the large activation energy limit, the critical condition for thermal explosion in the Semenov model is
\[
\delta_{\text{cr}} = \frac{1}{e},
\]
where the Semenov number is
\[
\delta = \frac{Q B V E_a}{h S R_u T_w^2} \exp\left(-\frac{E_a}{R_u T_w}\right).
\]
For \(\delta < 1/e\), a stable steady state exists; for \(\delta > 1/e\), thermal explosion occurs.
Proof. Nondimensionalize the temperature as \(\theta = E_a(T - T_w)/(R_u T_w^2)\), so that the Arrhenius factor becomes (in the Frank-Kamenetskii approximation) \(\exp(\theta/(1 + \varepsilon\theta)) \approx e^{\theta}\) for \(\varepsilon = R_u T_w / E_a \ll 1\). The steady state equation becomes \(\delta \, e^{\theta} = \theta\). The critical condition is tangency: \(\delta \, e^{\theta} = \theta\) and \(\delta \, e^{\theta} = 1\) simultaneously, giving \(\theta_{\text{cr}} = 1\) and \(\delta_{\text{cr}} = e^{-1}\). \(\square\)
The Semenov analysis assumes spatial uniformity — the entire reactive mass is at one temperature. David Frank-Kamenetskii (1939) extended this to the more realistic case of a spatially distributed system where heat conduction within the material competes with heat generation.
Example 5.3 (Semenov Criticality for a Batch Reactor). Consider a well-stirred batch reactor of volume \(V = 1\) L containing a reactive gas at \(T_w = 600\) K, with heat transfer coefficient \(h = 100\) W/(m\(^2\cdot\)K) and surface area \(S = 0.06\) m\(^2\) (approximately a 10 cm sphere). Taking \(Q B = 10^{12}\) W/m\(^3\), \(E_a / R_u = 15{,}000\) K, the Semenov number is
\[
\delta = \frac{Q B V E_a}{h S R_u T_w^2} \exp\left(-\frac{E_a}{R_u T_w}\right) = \frac{10^{12} \times 10^{-3} \times 15{,}000}{100 \times 0.06 \times 600^2} \exp(-25).
\]
Computing: the prefactor is \(15 \times 10^{12} / (2.16 \times 10^6) \approx 6.94 \times 10^6\), and \(e^{-25} \approx 1.39 \times 10^{-11}\), giving \(\delta \approx 9.7 \times 10^{-5}\). Since \(\delta \ll 1/e \approx 0.368\), this reactor is far from the explosion limit and will reach a stable steady state with a temperature rise of approximately \(\Delta T \approx \delta \times R_u T_w^2 / E_a = 9.7 \times 10^{-5} \times 24 \approx 0.002\) K above the wall temperature. Reducing the cooling (say, by insulating half the surface) doubles \(\delta\), but the system remains safe. Only when the gas temperature or pre-exponential factor increases enough to bring \(\delta\) to \(1/e\) does criticality occur.
5.2 Frank-Kamenetskii Theory
Definition 5.4 (Frank-Kamenetskii Problem). The Frank-Kamenetskii problem is the steady-state heat equation in a bounded domain \(\Omega\) with Dirichlet boundary conditions and Arrhenius heat release:
\[
\lambda \nabla^2 T + Q B \exp\left(-\frac{E_a}{R_u T}\right) = 0 \quad \text{in } \Omega, \qquad T = T_w \quad \text{on } \partial\Omega.
\]
In the Frank-Kamenetskii approximation (\(\varepsilon = R_u T_w / E_a \ll 1\)), with \(\theta = E_a(T - T_w)/(R_u T_w^2)\), this becomes
\[
\nabla^2 \theta + \delta_{FK} \, e^{\theta} = 0 \quad \text{in } \Omega, \qquad \theta = 0 \quad \text{on } \partial\Omega,
\]
where \(\delta_{FK} = (Q B E_a / (\lambda R_u T_w^2)) L^2 \exp(-E_a / R_u T_w)\) is the Frank-Kamenetskii parameter and \(L\) is the characteristic dimension of \(\Omega\).
This is a semilinear elliptic PDE with exponential nonlinearity — the Bratu-Gelfand problem. It has a rich mathematical structure: for \(\delta_{FK}\) below a critical value \(\delta_{\text{cr}}\), at least one (and sometimes two) solutions exist; above \(\delta_{\text{cr}}\), no solution exists and thermal explosion occurs.
Theorem 5.5 (Critical Frank-Kamenetskii Parameters). The critical values \(\delta_{\text{cr}}\) for the Frank-Kamenetskii problem depend on the geometry of \(\Omega\):
\[
\delta_{\text{cr}} = \begin{cases} 0.878 & \text{infinite slab (1D)}, \\ 2.000 & \text{infinite cylinder (2D)}, \\ 3.322 & \text{sphere (3D)}. \end{cases}
\]
For the slab, the exact solution can be written in terms of hyperbolic functions. For the sphere, \(\delta_{\text{cr}}\) must be determined numerically.
Proof. (Slab case.) In one dimension, \(\theta'' + \delta e^{\theta} = 0\) on \([-1, 1]\) with \(\theta(\pm 1) = 0\). By symmetry, \(\theta'(0) = 0\) and \(\theta(0) = \theta_m\) (the maximum). Multiplying by \(\theta'\) and integrating from \(0\) to \(x\),
\[
\frac{1}{2}(\theta')^2 = \delta(e^{\theta_m} - e^{\theta}).
\]
With the substitution \(u = e^{(\theta_m - \theta)/2}\), one obtains \(\theta(x) = \theta_m - 2\ln\cosh(\sqrt{\delta e^{\theta_m}/2}\, x)\). The boundary condition \(\theta(1) = 0\) gives the implicit relation \(\theta_m = 2\ln\cosh(\sqrt{\delta e^{\theta_m}/2})\). The critical condition is found by requiring that this equation has a degenerate (tangent) solution, yielding \(\delta_{\text{cr}} = 0.878\ldots\) \(\square\)
Example 5.7 (Critical Size of a Reactive Stockpile). Consider a spherical stockpile of organic material (e.g., wood chips) at ambient temperature \(T_w = 300\) K, with thermal conductivity \(\lambda = 0.1\) W/(m\(\cdot\)K), heat of reaction \(Q = 5 \times 10^5\) J/kg, pre-exponential factor \(B = 10^8\) s\(^{-1}\), activation energy \(E_a/R_u = 10{,}000\) K, and density \(\rho = 500\) kg/m\(^3\). The Frank-Kamenetskii parameter is
\[
\delta_{FK} = \frac{\rho Q B E_a}{\lambda R_u T_w^2} R^2 \exp\left(-\frac{E_a}{R_u T_w}\right) = \frac{500 \times 5 \times 10^5 \times 10^8 \times 10{,}000}{0.1 \times 300^2} \, R^2 \, e^{-33.3}.
\]
Setting \(\delta_{FK} = 3.322\) (the spherical critical value) and solving for \(R\) gives the critical radius. After evaluating the exponential (\(e^{-33.3} \approx 3.7 \times 10^{-15}\)), one obtains \(R_{\text{cr}} \approx 4\) m. Stockpiles with radius exceeding this value will undergo spontaneous thermal runaway. This type of calculation, refined with measured kinetic parameters, is routinely used in warehouse safety assessments and in determining safe storage dimensions for reactive chemicals.
Example 5.8 (Explosion Limits of Hydrogen). The classical explosion limits of hydrogen-oxygen mixtures can be partially understood through thermal explosion theory, though the second limit (around 1 atm, 400--500°C) is actually a chain-branching phenomenon rather than a purely thermal one. The first and third limits are thermal: the first limit occurs at low pressure where mean free paths are long and heat loss to the walls is efficient (\(\delta_{FK} < \delta_{\text{cr}}\)); the third limit occurs at high pressure where the thermal runaway overwhelms conduction. The second limit requires the detailed chemistry of the \(\text{H}_2\text{-}\text{O}_2\) system (chain branching by \(\text{H} + \text{O}_2 \to \text{OH} + \text{O}\) vs. chain termination by \(\text{H} + \text{O}_2 + \text{M} \to \text{HO}_2 + \text{M}\)) and cannot be explained by single-step thermal theory.
5.3 Ignition Delay and the Critical Damköhler Number
Ignition is a transient phenomenon: a reactive mixture at low temperature evolves under the influence of chemical heat release until either a steady weakly-reacting state is reached or a thermal runaway occurs. The time to runaway is the ignition delay.
Definition 5.9 (Ignition Delay). The ignition delay \(t_{\text{ign}}\) is the time required for a homogeneous reactive mixture at initial temperature \(T_0\) to reach thermal runaway. In the Semenov model,
\[
t_{\text{ign}} = \frac{\rho c_v R_u T_0^2}{Q B E_a} \exp\left(\frac{E_a}{R_u T_0}\right) \cdot g(\delta),
\]
where \(g(\delta)\) is an order-one function of the Semenov number that diverges logarithmically as \(\delta \to \delta_{\text{cr}}\).
The Arrhenius dependence of the ignition delay on initial temperature — \(t_{\text{ign}} \propto \exp(E_a / R_u T_0)\) — is one of the most robust predictions of combustion theory and is confirmed by shock-tube experiments across a wide range of fuels. The shock tube, invented by Vieille in 1899 and refined for combustion kinetics studies by Glick, Squire, and Hertzberg in the 1950s, remains the primary experimental apparatus for measuring ignition delays. By rupturing a diaphragm between high-pressure driver gas and low-pressure test gas, a shock wave is generated that instantaneously heats the test mixture to a precisely controlled temperature and pressure. The reflected shock from the end wall creates a region of stagnant, heated gas — essentially a batch reactor at known initial conditions — and the ignition delay is measured optically (by CH\(^*\) or OH\(^*\) chemiluminescence) or by the pressure rise associated with thermal runaway. An Arrhenius plot of \(\ln t_{\text{ign}}\) versus \(1/T_0\) yields a straight line with slope \(E_a / R_u\), providing one of the most common methods for determining global activation energies.
Example 5.10 (Ignition Delay for n-Heptane). Shock-tube measurements of n-heptane ignition delay at pressures of 2--50 atm and temperatures of 700--1300 K reveal two distinct regimes. At high temperatures (\(T > 1000\) K), the ignition delay follows the simple Arrhenius law \(t_{\text{ign}} \propto \exp(E_a / R_u T)\) with an apparent activation energy \(E_a / R_u \approx 15{,}000\) K. At intermediate temperatures (700--900 K), a region of negative temperature coefficient (NTC) behavior is observed: the ignition delay actually increases as temperature increases. This NTC behavior arises from the competing effects of low-temperature chain-branching pathways (alkylperoxy radical chemistry) and cannot be captured by a single-step Arrhenius model. The NTC region is critically important for engine knock in internal combustion engines and for homogeneous charge compression ignition (HCCI) combustion strategies.
Proposition 5.12. In the Frank-Kamenetskii approximation, the maximum temperature perturbation \(\theta_m\) in a reactive body at the critical condition satisfies \(\theta_m = 1\) regardless of geometry, corresponding to a dimensional temperature rise of \(\Delta T_{\text{cr}} = R_u T_w^2 / E_a\). For a wall temperature of \(T_w = 500\) K and activation energy \(E_a/R_u = 20{,}000\) K, this gives \(\Delta T_{\text{cr}} = 12.5\) K — a remarkably small temperature rise triggers the explosion. This universality of \(\theta_m = 1\) is a consequence of the linearization of the Arrhenius exponential (the Frank-Kamenetskii approximation) and holds only asymptotically for large \(\beta\).
5.4 The Connection Between Ignition and Explosion Theory
The ignition delay analysis connects the transient Semenov model to the broader framework of chemical reactor theory. In the language of dynamical systems, the ignition event corresponds to the trajectory of the ODE system crossing from the basin of attraction of the stable low-temperature steady state to the unbounded domain leading to thermal runaway. The separatrix between these regions is the unstable intermediate steady state on the S-curve.
5.5 Extinction by Stretch and Multiplicity
Extinction of a diffusion flame by aerodynamic stretching is intimately connected to the S-curve of Chapter 4. As the stretch rate increases, the Damköhler number decreases along the upper branch until the extinction turning point is reached, at which the flame jumps abruptly to the lower (frozen) branch.
The connection between the S-curve turning point and the dynamic flame response provides the bridge between steady-state multiplicity analysis and unsteady flame behavior.
Theorem 5.14 (Extinction Criterion). For a counterflow diffusion flame in the large activation energy limit, extinction occurs when the reduced Damköhler number
\[
\delta_L = \frac{4\text{Da}}{\beta^3} \cdot F(Z_{st})
\]
falls below the critical value \(\delta_E(\gamma)\) determined by Liñán's canonical inner problem (Theorem 3.10). The function \(F(Z_{st})\) depends on the stoichiometric mixture fraction and captures the effect of flame location within the mixing layer.
Chapter 6: Detonation Waves
6.1 Rankine-Hugoniot Relations for Reactive Flow
A detonation is a supersonic combustion wave in which a shock wave triggers chemical reactions. Unlike the subsonic deflagrations (laminar flames) of Chapters 2–5, detonations propagate at velocities of 1500–3000 m/s and involve enormous pressure jumps. The mathematical framework for analyzing detonation waves begins with the Rankine-Hugoniot relations — the conservation laws applied across a discontinuity — augmented by an energy release term from the chemical reaction.
Definition 6.1 (Reactive Rankine-Hugoniot Relations). For a steady, one-dimensional reactive wave propagating at velocity \(D\) into a quiescent gas at state \((\rho_0, p_0, T_0)\), the conservation laws across the wave give:
Mass: \(\rho_0 D = \rho_1 (D - u_1) \equiv m\),
Momentum: \(p_0 + \rho_0 D^2 = p_1 + \rho_1(D - u_1)^2\), or equivalently \(p_1 - p_0 = m^2(v_0 - v_1)\),
Energy: \(h_0 + \frac{1}{2}D^2 = h_1 + \frac{1}{2}(D - u_1)^2 + q\),
where \(v = 1/\rho\) is the specific volume, \(h\) is the specific enthalpy, \(u_1\) is the product gas velocity in the lab frame, and \(q\) is the heat released per unit mass by the chemical reaction.
The reactive Rankine-Hugoniot curve differs from the inert Hugoniot by being shifted upward in the \((v, p)\)-plane by an amount proportional to the heat release \(q\). This shift is what makes detonation possible: it allows the final state to lie at a higher pressure and lower specific volume than the inert shocked state.
6.2 Chapman-Jouguet Theory
David Chapman (1899) and Émile Jouguet (1905) independently proposed a criterion for selecting the detonation velocity: the detonation propagates at the minimum velocity consistent with the Rankine-Hugoniot relations and the entropy condition. This minimum velocity corresponds to a tangency condition on the Hugoniot curve.
Theorem 6.2 (Chapman-Jouguet Condition). The Chapman-Jouguet (CJ) detonation corresponds to the point on the reactive Hugoniot curve where the Rayleigh line (the straight line connecting initial and final states in the \((v,p)\)-plane) is tangent to the Hugoniot. At this point:
- The flow velocity relative to the wave in the product gas equals the local sound speed: \(|D - u_1| = c_1\) (sonic condition).
- The detonation velocity is minimized: \(D = D_{CJ}\).
- The entropy of the product gas is maximized among all points on the Hugoniot reachable by a Rayleigh line through the initial state.
Proof. The Rayleigh line passing through \((v_0, p_0)\) with slope \(-m^2\) is \(p - p_0 = -m^2(v - v_0)\). The Hugoniot curve \(H(v, p; q) = 0\) relates final states consistent with conservation. At a tangency point, the Rayleigh line touches the Hugoniot, so
\[
\frac{dp}{dv}\bigg|_H = \frac{dp}{dv}\bigg|_R = -m^2.
\]
From thermodynamic identities applied to the Hugoniot, one can show that the slope of the Hugoniot at any point equals \(-c^2/v^2\) where \(c\) is the local sound speed, provided the thermodynamic derivative is taken at constant composition. Thus at the tangency point, \(m^2 = c_1^2/v_1^2 = \rho_1^2 c_1^2\). Since \(m = \rho_1(D - u_1)\), this gives \((D - u_1)^2 = c_1^2\), the sonic condition. That this tangency point minimizes \(D\) follows from the geometry: Rayleigh lines with smaller slope (smaller \(m\), hence smaller \(D\)) do not intersect the Hugoniot. \(\square\)
For a stoichiometric hydrogen-air mixture, the CJ detonation velocity is approximately 1970 m/s, with a CJ pressure of about 15.6 atm and a CJ temperature of about 2950 K. For stoichiometric methane-air, the CJ velocity is approximately 1800 m/s.
Example 6.4 (CJ Detonation Velocity for an Ideal Gas). For an ideal gas with constant \(\gamma\), the CJ detonation velocity can be obtained in closed form. Let \(\tilde{q} = 2(\gamma^2 - 1)q / (c_0^2)\) where \(c_0 = \sqrt{\gamma p_0 / \rho_0}\) is the initial sound speed. Then
\[
\frac{D_{CJ}}{c_0} = \sqrt{1 + \frac{\tilde{q}}{2}} + \sqrt{\frac{\tilde{q}}{2}}.
\]
For strong detonations (\(\tilde{q} \gg 1\)), this simplifies to \(D_{CJ} \approx c_0 \sqrt{\tilde{q}} \propto \sqrt{q}\), showing that the detonation velocity scales as the square root of the heat release. For a typical hydrocarbon with \(q \approx 3\) MJ/kg and \(\gamma = 1.25\), this gives \(D_{CJ} \approx 1800\) m/s.
6.3 ZND Detonation Structure
The Chapman-Jouguet theory treats the detonation as a discontinuity and says nothing about its internal structure. The ZND model, developed independently by Zeldovich (1940), von Neumann (1942), and Döring (1943), resolves this structure. In the ZND model, the detonation consists of a leading inert shock wave that compresses and heats the gas, followed by a reaction zone where the chemical energy is released. The gas enters the reaction zone at the von Neumann spike conditions (the inert Hugoniot state) and evolves along the Rayleigh line toward the CJ point as the reaction proceeds to completion.
Definition 6.5 (ZND Detonation Model). The
ZND model describes the internal structure of a detonation as:
- A leading shock that compresses the gas from the initial state \((p_0, v_0, T_0)\) to the von Neumann spike state \((p_N, v_N, T_N)\) on the inert Hugoniot.
- An induction zone where the temperature is high enough for the reaction to begin but the heat release is still negligible.
- A reaction zone where the reaction progress variable \(\lambda\) increases from 0 to 1 and the state slides along the Rayleigh line from the von Neumann point to the CJ point.
The governing equations in the shock-attached frame are
\[
m \frac{dv}{dx} = \frac{dp}{dx}, \qquad m \frac{du}{dx} = -\frac{dp}{dx}, \qquad m \frac{d\lambda}{dx} = \dot{\omega}(\rho, T, \lambda),
\]
supplemented by the energy equation and an equation of state.
Theorem 6.6 (Thermicity and Sonic Plane). Within the ZND reaction zone, the flow Mach number relative to the wave satisfies
\[
\frac{1 - M^2}{M^2} \frac{dM^2}{dx} = (\gamma + 1) \sigma \dot{\omega},
\]
where \(\sigma\) is the thermicity coefficient, defined as the fractional rate of change of the sound speed due to chemical heat release. The flow enters the reaction zone with \(M < 1\) (subsonic relative to the shock) for a CJ detonation and accelerates to \(M = 1\) at the end of the reaction zone. For an overdriven detonation (\(D > D_{CJ}\)), the flow remains subsonic throughout.
The sonic condition at the end of the CJ reaction zone is a generalization of the Chapman-Jouguet tangency condition. It implies that acoustic disturbances from behind the detonation cannot catch up with the wave front, making the CJ detonation an information-separating structure. This has profound implications for detonation stability and for the interaction of detonations with rarefaction waves.
6.4 Deflagration-to-Detonation Transition
One of the most dramatic phenomena in combustion is the deflagration-to-detonation transition (DDT), where a slow laminar flame (propagating at tens of cm/s) accelerates and abruptly transitions to a detonation (propagating at km/s). This transition involves a complex interplay of flame instabilities, turbulent flame acceleration, shock-flame interaction, and local explosion centers.
Definition 6.8 (Deflagration-to-Detonation Transition). DDT is the transition from a subsonic deflagration to a supersonic detonation, typically occurring in confined geometries (tubes). The process involves several stages:
- Flame acceleration: The flame wrinkles and accelerates due to hydrodynamic and thermal-diffusive instabilities and interaction with obstacles.
- Shock formation: The accelerating flame pushes compression waves ahead that coalesce into a shock.
- Hot spot formation: Shock reflections or flame-shock interactions create localized regions of high temperature and pressure.
- Detonation onset: A local explosion in a hot spot generates a detonation wave that overtakes the leading shock.
Example 6.10 (Run-Up Distance in a Tube). The distance from the ignition point to the location where a detonation first forms is called the run-up distance \(L_{\text{DDT}}\). For stoichiometric hydrogen-air in a smooth tube of diameter \(d = 5\) cm, experimental measurements give \(L_{\text{DDT}} \approx 1\)--\(2\) m. Adding obstacles (such as Shchelkin spirals with blockage ratio \(\sim 0.4\)) reduces the run-up distance to approximately 0.3--0.5 m by promoting turbulent flame acceleration. The run-up distance correlates roughly with the detonation cell size \(\lambda_c\) as \(L_{\text{DDT}} / d \propto d / \lambda_c\): mixtures with smaller cells (more reactive, more unstable detonation structure) have shorter run-up distances. For hydrogen-air, \(\lambda_c \approx 1.5\) cm, so \(d/\lambda_c \approx 3.3\), and the correlation gives \(L_{\text{DDT}} / d \approx 20\)--\(40\), consistent with observation.
6.5 Stability of Detonation Waves
Real detonation waves are not the steady, one-dimensional structures of the ZND model. They exhibit complex multi-dimensional instabilities that manifest as cellular patterns on the detonation front, observable as diamond-shaped patterns on soot foils placed in the path of the wave.
Theorem 6.11 (Linear Stability of ZND Detonation). The one-dimensional linear stability of a ZND detonation to planar perturbations is governed by the normal mode analysis. Perturbations of the form \(\tilde{f}(x) e^{st + iky}\) (where \(s\) is the growth rate and \(k\) is the transverse wavenumber) satisfy a boundary value problem in the reaction zone with a radiation condition at \(x \to -\infty\). The eigenvalue \(s(k)\) determines stability:
- If \(\text{Re}(s) > 0\) for some \(k\), the ZND detonation is linearly unstable.
- For one-dimensional perturbations (\(k = 0\)), the detonation is unstable when the activation energy is sufficiently large: specifically, when the stability parameter \(\chi = E_a \Delta T / (R_u T_N^2)\) (where \(\Delta T = T_N - T_0\)) exceeds a critical value of order unity.
The stability analysis of detonation waves was carried out by Erpenbeck (1964) using a Laplace transform approach and later by Lee and Stewart (1990) using a normal mode formulation. The key finding is that detonation stability is controlled primarily by the ratio of the activation energy to the post-shock temperature: high activation energies produce unstable detonations with complex cellular patterns, while low activation energies produce stable or weakly unstable detonations.
Example 6.13 (Detonation Cell Sizes for Common Fuels). Measured detonation cell sizes in stoichiometric fuel-air mixtures at standard conditions provide a practical database for safety engineering:
\[
\begin{aligned}
\text{Hydrogen-air:} \quad &\lambda_c \approx 1.5 \text{ cm}, \\
\text{Acetylene-air:} \quad &\lambda_c \approx 0.6 \text{ cm}, \\
\text{Ethylene-air:} \quad &\lambda_c \approx 2.5 \text{ cm}, \\
\text{Propane-air:} \quad &\lambda_c \approx 5.5 \text{ cm}, \\
\text{Methane-air:} \quad &\lambda_c \approx 30 \text{ cm}.
\end{aligned}
\]
The very large cell size for methane-air explains why methane detonations are difficult to initiate and sustain in laboratory-scale facilities. Acetylene, with its small cell size and high reactivity (due to the triple bond), is the most detonation-prone common hydrocarbon — a fact of considerable importance in the chemical process industry, where acetylene decomposition detonations are a recognized hazard.
Chapter 7: Flame Instabilities and Turbulent Combustion
7.1 Darrieus-Landau Hydrodynamic Instability
A planar premixed flame is unconditionally unstable to long-wavelength perturbations due to the density jump across the flame front. This remarkable result, obtained independently by Darrieus (1938) and Landau (1944), shows that perfectly planar flames cannot exist in nature — they are always wrinkled by the hydrodynamic instability, even in the absence of any other perturbing mechanism.
Theorem 7.1 (Darrieus-Landau Instability). Consider a planar premixed flame separating unburned gas (density \(\rho_u\)) from burned gas (density \(\rho_b < \rho_u\)). Subject the flame to a small sinusoidal perturbation with wavenumber \(k\). Treating the flame as a discontinuity propagating at the constant laminar burning velocity \(S_L\) relative to the unburned gas, the growth rate of the perturbation is
\[
\sigma = k S_L \left(\frac{-r + \sqrt{r^3 + r^2 - r}}{r + 1}\right),
\]
where \(r = \rho_u / \rho_b > 1\) is the density ratio (typically \(r \approx 5\text{--}8\) for hydrocarbon-air flames).
Proof. The linearized analysis proceeds by solving the incompressible Euler equations (potential flow) on either side of the perturbed flame sheet \(x = \xi(y,t)\) with \(\xi = \hat{\xi} e^{iky + \sigma t}\). On the unburned side (\(x < \xi\)), the velocity potential satisfies \(\nabla^2 \phi_u = 0\) with \(\phi_u \to -S_L x\) as \(x \to -\infty\). On the burned side (\(x > \xi\)), similarly with \(\phi_b \to -(S_L/r)x\) as \(x \to +\infty\). The jump conditions at the flame are: (i) the flame propagates at speed \(S_L\) relative to the unburned gas, and (ii) mass and momentum are conserved. Substituting the perturbation ansatz and linearizing gives a quadratic equation for \(\sigma/k\), whose positive root is the stated formula.
Since \(r > 1\), the discriminant \(r^3 + r^2 - r = r(r^2 + r - 1) > 0\) for all \(r > 1\), and \(\sqrt{r^3 + r^2 - r} > r\) (which can be verified by squaring), so \(\sigma > 0\) for all wavenumbers \(k > 0\). The flame is unconditionally unstable. \(\square\)
The Darrieus-Landau instability has a simple physical mechanism: when the flame bulges toward the unburned gas, the streamlines diverge ahead of the bulge (reducing the approach velocity and causing the flame to advance further) and converge behind it (increasing the approach velocity and causing the flame to recede further). The density jump amplifies this effect because the gas accelerates through the flame.
Example 7.2 (Growth Rate of the Darrieus-Landau Instability). For a stoichiometric methane-air flame with \(r = \rho_u/\rho_b \approx 6\) and \(S_L = 40\) cm/s, the Darrieus-Landau growth rate for a perturbation with wavelength \(\lambda = 1\) cm (wavenumber \(k = 2\pi/\lambda \approx 628\) m\(^{-1}\)) is
\[
\sigma = k S_L \left(\frac{-6 + \sqrt{6^3 + 6^2 - 6}}{6 + 1}\right) = 628 \times 0.4 \times \left(\frac{-6 + \sqrt{246}}{7}\right).
\]
Since \(\sqrt{246} \approx 15.7\), we obtain \(\sigma \approx 251 \times (9.7/7) \approx 348\) s\(^{-1}\). The corresponding e-folding time is \(\tau = 1/\sigma \approx 2.9\) ms, which is comparable to the flow time scale \(\lambda / S_L = 25\) ms. For perturbations at 10 times shorter wavelength (\(\lambda = 1\) mm), the growth rate is 10 times larger, giving an e-folding time of 0.29 ms. This rapid growth at short wavelengths is precisely what the Markstein correction suppresses: the curvature-dependent burning velocity introduces a stabilizing term that grows as \(k^2\), eventually overwhelming the destabilizing \(k\)-linear growth for \(k > k_c \sim 1/\mathcal{L}\).
7.2 Thermal-Diffusive Instability
A second fundamental instability of premixed flames arises from differential diffusion of heat and mass. When the Lewis number \(\text{Le} = \lambda/(\rho c_p D)\) differs from unity, the flame can be unstable even without the density jump. This thermal-diffusive (or Turing-type) instability was first analyzed by Barenblatt, Zeldovich, and Istratov (1962) and leads to the formation of cellular flame patterns.
Theorem 7.4 (Thermal-Diffusive Instability). In the limit of large activation energy \(\beta \gg 1\) and constant density (eliminating the Darrieus-Landau mechanism), a planar premixed flame is unstable to perturbations with wavenumber \(k\) when
\[
\text{Le} < \text{Le}_{\text{cr}} = 1 - \frac{2}{\beta} + O(\beta^{-2}).
\]
The instability is oscillatory (pulsating) when it first appears if the effective Lewis number is above a second threshold, and steady (cellular) otherwise. The most unstable wavenumber scales as \(k_{\max} \sim \beta / \delta_f\).
Proof. (Sketch.) Perturb the flame position and the inner reaction-zone solution of the ZFK analysis. The perturbation \(\hat{\xi} e^{iky + \sigma t}\) modifies the curvature of the flame sheet, which changes the local burning velocity through the Markstein effect. In the inner zone, the temperature perturbation satisfies a modified version of Liñán's canonical equation with a term proportional to \(k^2 (1 - \text{Le})\). The solvability condition for the inner problem gives the dispersion relation
\[
\sigma = S_L^2 k^2 \delta_f \left(\frac{\beta(\text{Le}_{\text{cr}} - \text{Le})}{2} - k^2 \delta_f^2\right) + O(\beta^{-1}).
\]
The flame is unstable (\(\sigma > 0\)) when \(\text{Le} < \text{Le}_{\text{cr}}\) for wavenumbers in the range \(0 < k < k_c\), where \(k_c \propto \sqrt{\beta(\text{Le}_{\text{cr}} - \text{Le})} / \delta_f\). \(\square\)
Example 7.5 (Critical Lewis Number for Hydrogen-Air). For a lean hydrogen-air flame with \(\beta \approx 8\), the critical Lewis number from Theorem 7.4 is
\[
\text{Le}_{\text{cr}} = 1 - \frac{2}{8} = 0.75.
\]
The effective Lewis number of lean hydrogen-air is approximately \(\text{Le} \approx 0.3\), which is far below the critical value. The instability growth rate from the dispersion relation scales as \(\sigma \propto \beta(\text{Le}_{\text{cr}} - \text{Le}) \approx 8 \times 0.45 = 3.6\), which is large, explaining the violent cellular and pulsating instabilities observed in lean hydrogen flames. By contrast, a lean propane-air flame has \(\text{Le} \approx 1.8 > \text{Le}_{\text{cr}}\), so it is stable to thermal-diffusive perturbations and exhibits only the hydrodynamic Darrieus-Landau instability with its characteristically smoother, larger-scale wrinkling.
The physical mechanism of the thermal-diffusive instability is intuitive: when \(\text{Le} < 1\), mass diffuses faster than heat. At a convex perturbation of the flame (bulging toward the unburned gas), the enhanced mass flux of fuel increases the local burning rate, but the defocused heat flux cools the flame less effectively. The net effect is that the perturbation grows. For \(\text{Le} > 1\), the opposite occurs and the flame is stable to this mechanism (though still unstable to Darrieus-Landau). Hydrogen, with \(\text{Le} \approx 0.3\), is strongly susceptible to thermal-diffusive instability, which is why hydrogen flames exhibit spectacular cellular and pulsating patterns.
7.3 Markstein Length and the Sivashinsky Equation
The Markstein length \(\mathcal{L}\) characterizes the linear response of the laminar burning velocity to flame stretch and curvature. It was introduced by George Markstein (1951) as a phenomenological parameter and later derived from first principles using activation energy asymptotics.
Definition 7.6 (Markstein Length). The local burning velocity \(S_n\) of a weakly curved and strained flame is related to the unstretched burning velocity \(S_L\) by
\[
S_n = S_L - \mathcal{L} \kappa - \mathcal{L}_s \frac{1}{A}\frac{dA}{dt},
\]
where \(\kappa = \nabla \cdot \mathbf{n}\) is the flame curvature (with \(\mathbf{n}\) the normal pointing toward the unburned gas), \(dA/dt\) is the rate of change of flame surface area (stretch), and \(\mathcal{L}\), \(\mathcal{L}_s\) are the curvature and strain Markstein lengths, respectively. In the large activation energy limit,
\[
\mathcal{L} = \delta_f \left(\frac{\beta(\text{Le} - 1)}{2} + C(r)\right),
\]
where \(C(r)\) is a contribution from the Darrieus-Landau instability that depends on the density ratio.
The combined effects of hydrodynamic instability, thermal-diffusive instability, and curvature stabilization are captured by the Sivashinsky equation, a nonlinear evolution equation for the flame front position.
Theorem 7.7 (Sivashinsky Equation). In the limit of weak thermal expansion (\(r - 1 = O(\varepsilon)\), \(\varepsilon \ll 1\)), the flame front position \(\xi(y, t)\) evolves according to the Michelson-Sivashinsky equation:
\[
\frac{\partial \xi}{\partial t} + \frac{1}{2}|\nabla \xi|^2 = \nu \nabla^2 \xi + \mathcal{I}[\xi],
\]
where \(\nu\) is an effective Markstein diffusivity and \(\mathcal{I}[\xi]\) is the Landau-Darrieus operator, a nonlocal pseudo-differential operator whose Fourier symbol is \(|k|\):
\[
\widehat{\mathcal{I}[\xi]}(k) = |k| \hat{\xi}(k).
\]
Gregory Sivashinsky derived this equation in 1977 by combining the Darrieus-Landau instability with a weakly nonlinear expansion. The equation captures the essential dynamics of wrinkled flames: the \(|\nabla\xi|^2\) term is the Huygens (geometric) nonlinearity that generates cusps, the Laplacian provides short-wavelength stabilization, and the nonlocal term drives the long-wavelength instability. The Sivashinsky equation has been the subject of extensive mathematical analysis and is known to produce solutions with fractal-like wrinkling of the flame front.
Example 7.8 (Cellular Flame Patterns). Consider a planar flame in a two-dimensional channel of width \(L\). The Sivashinsky equation with periodic boundary conditions admits steady cellular solutions. The number of cells \(n\) is selected by the competition between the Darrieus-Landau instability (which favors many cells) and the curvature stabilization (which penalizes short wavelengths). For a methane-air flame with \(\delta_f \approx 0.5\) mm, \(r \approx 6\), and \(\mathcal{L} \approx 1\) mm, a channel of width \(L = 10\) cm supports approximately \(n \approx L/(2\pi\mathcal{L}) \approx 16\) cells. This is consistent with experimental observations of cellular flames in Hele-Shaw cells.
7.4 Turbulent Flame Speed and Flamelet Models
When the flow is turbulent, the flame front is wrinkled by turbulent eddies across a range of scales, greatly increasing the flame surface area and hence the overall burning rate. The turbulent flame speed \(S_T\) — the effective propagation velocity of the turbulent flame brush — is one of the most important quantities in practical combustion engineering.
Definition 7.9 (Turbulent Flame Speed). The turbulent flame speed \(S_T\) is defined as the volume of unburned gas consumed per unit time per unit area of the mean flame surface:
\[
S_T = \frac{1}{A_{\text{mean}}} \int_{\text{flame}} S_L \, dA,
\]
where the integral is over the wrinkled flame surface and \(A_{\text{mean}}\) is the projected area.
Theorem 7.10 (Damköhler's Scaling Laws). Gerhard Damköhler (1940) proposed two limiting scaling laws for the turbulent flame speed:
- Large-scale turbulence (turbulent eddies much larger than the flame thickness): the flame is wrinkled but locally retains its laminar structure. The turbulent flame speed is
\[
\frac{S_T}{S_L} \approx 1 + \frac{u'}{S_L},
\]
where \(u'\) is the turbulent velocity fluctuation (rms). The wrinkling increases the flame area by a factor \(\sim u'/S_L\).
- Small-scale turbulence (turbulent eddies much smaller than the flame thickness): the eddies enhance the effective diffusivity within the flame. The turbulent flame speed is
\[
\frac{S_T}{S_L} \approx \sqrt{1 + \frac{D_T}{D}},
\]
where \(D_T \sim u' l_T\) is the turbulent diffusivity and \(l_T\) is the integral length scale.
These scaling laws, while approximate, frame the problem correctly. The first regime (now called the wrinkled flamelet regime or corrugated flamelet regime) is the basis of the flamelet model, in which the turbulent flame is viewed as an ensemble of laminar flamelets convected and distorted by the turbulent flow. The second regime (the distributed reaction zone regime) is less common in practice because the Kolmogorov scale is usually larger than the flame thickness for most practical flames.
Definition 7.11 (Flamelet Model). The flamelet model (Peters, 1984) assumes that the turbulent flame consists of an ensemble of thin, locally one-dimensional, laminar flame structures embedded in the turbulent flow. In the flamelet regime, the species and temperature fields are parameterized by the mixture fraction \(Z\) and the scalar dissipation rate \(\chi = 2D|\nabla Z|^2\):
\[
\rho \frac{\partial Y_i}{\partial t} = \frac{\rho \chi}{2} \frac{\partial^2 Y_i}{\partial Z^2} + \dot{\omega}_i.
\]
This is the unsteady flamelet equation, which replaces the full three-dimensional reactive Navier-Stokes equations by a family of one-dimensional problems parameterized by \(\chi\).
Chapter 8: Computational Combustion
8.1 Stiff ODE Systems from Chemical Kinetics
The numerical simulation of combustion flows is challenging primarily because of the extreme stiffness of the chemical kinetics. A detailed mechanism for methane combustion (GRI-Mech 3.0) involves 53 species and 325 reactions, with time scales ranging from nanoseconds (radical reactions) to milliseconds (slow oxidation steps) — a stiffness ratio of \(10^6\) or more. The mathematical character of this stiffness and the numerical methods required to handle it are the starting point of computational combustion.
Definition 8.1 (Stiffness in Chemical Kinetics). A system of ODEs \(d\mathbf{c}/dt = \mathbf{R}(\mathbf{c})\), where \(\mathbf{c} = (c_1, \ldots, c_N)^T\) is the vector of species concentrations and \(\mathbf{R}\) is the vector of net production rates, is stiff if the Jacobian \(J = \partial \mathbf{R}/\partial \mathbf{c}\) has eigenvalues \(\lambda_i\) with \(\text{Re}(\lambda_i) < 0\) and
\[
S = \frac{\max_i |\text{Re}(\lambda_i)|}{\min_i |\text{Re}(\lambda_i)|} \gg 1.
\]
The stiffness ratio \(S\) measures the separation of the fastest and slowest chemical time scales.
Stiffness arises because combustion involves both fast radical equilibration reactions (e.g., \(\text{H} + \text{O}_2 \rightleftharpoons \text{OH} + \text{O}\), time scale \(\sim 10^{-8}\) s at flame temperatures) and slow overall reactions (e.g., CO oxidation, time scale \(\sim 10^{-3}\) s). An explicit time integrator would need to resolve the fastest scale, making the computation prohibitively expensive. Implicit methods are essential.
Theorem 8.2 (A-Stability Requirement). For the stiff ODE system \(\mathbf{c}' = \mathbf{R}(\mathbf{c})\) with Jacobian eigenvalues satisfying \(\text{Re}(\lambda_i) < 0\) for all \(i\), a numerical method is
A-stable if its stability region contains the entire left half of the complex plane. A necessary condition for an A-stable linear multistep method is (Dahlquist's second barrier):
- The method has order at most 2.
- Among all A-stable linear multistep methods of order 2, the trapezoidal rule has the smallest error constant.
Consequently, practical stiff solvers for combustion chemistry use either implicit Runge-Kutta methods (e.g., SDIRK, Radau IIA) or backward differentiation formulas (BDF) of order 1--5 with adaptive order and step size control.
Example 8.3 (VODE/CVODE for Combustion Chemistry). The VODE solver (Variable-coefficient ODE solver, Brown, Byrne, and Hindmarsh, 1989) and its successor CVODE (from the SUNDIALS package) are the workhorses of combustion chemistry integration. They use BDF methods of order 1--5 with Newton iteration for the implicit solve at each step, and adaptive step size and order selection based on local error estimates. For a hydrogen-air autoignition problem at 1000 K and 1 atm (9 species, 21 reactions), CVODE typically takes \(\sim 200\) time steps to integrate from initial conditions to steady state, compared to \(\sim 10^6\) steps that would be required by an explicit method at the same accuracy. The key cost at each step is the LU factorization of the \(N \times N\) Jacobian matrix, which costs \(O(N^3)\) operations. For large mechanisms (\(N > 100\)), sparse linear algebra is essential.
8.2 Operator Splitting Methods
In multi-dimensional reactive flow simulations, the full system couples the Navier-Stokes equations (fluid mechanics) with the species transport equations (reaction-diffusion). Solving these as a monolithic system is possible but expensive. Operator splitting methods decompose the problem into simpler sub-problems — typically a flow step, a diffusion step, and a reaction step — that can be solved with specialized methods.
Definition 8.4 (Strang Splitting). Given the system \(\partial \mathbf{u}/\partial t = \mathcal{L}_1(\mathbf{u}) + \mathcal{L}_2(\mathbf{u})\), where \(\mathcal{L}_1\) represents transport and \(\mathcal{L}_2\) represents chemistry, the
Strang splitting advances the solution from \(t^n\) to \(t^{n+1} = t^n + \Delta t\) by:
- Advance \(\mathcal{L}_2\) by \(\Delta t/2\) (half-step of chemistry).
- Advance \(\mathcal{L}_1\) by \(\Delta t\) (full step of transport).
- Advance \(\mathcal{L}_2\) by \(\Delta t/2\) (half-step of chemistry).
This achieves second-order accuracy in \(\Delta t\) (compared to first-order for Lie splitting, which omits the symmetric half-steps).
Example 8.6 (Splitting Error in a Methane-Air Flame). To quantify the splitting error, consider a one-dimensional premixed methane-air flame computed with the GRI-Mech 3.0 mechanism. With Strang splitting and time step \(\Delta t = 1 \, \mu\)s, the flame speed computed by the split scheme differs from the monolithic (unsplit) solution by approximately 0.5%. Increasing \(\Delta t\) to \(10 \, \mu\)s increases the splitting error to about 5%, primarily because the transport step displaces the flame front by a distance comparable to the reaction zone thickness during one step, causing the chemistry sub-step to operate on an inconsistent temperature profile. A useful rule of thumb is that the splitting time step should satisfy \(\Delta t \lesssim \delta_r / S_L\), where \(\delta_r \sim \delta_f / \beta\) is the reaction zone thickness. For methane-air, \(\delta_r \approx 50 \, \mu\)m and \(S_L \approx 40\) cm/s, giving \(\Delta t \lesssim 1.25 \, \mu\)s, consistent with the observed accuracy threshold.
8.3 Adaptive Mesh Refinement for Flame Fronts
The thin reaction zone of a flame (thickness \(\sim 0.01\)–\(0.1\) mm) embedded in a flow domain of macroscopic size (centimeters to meters) creates a scale separation that makes uniform grids impractical. Adaptive mesh refinement (AMR) dynamically refines the computational grid in regions of steep gradients (the flame front) while keeping a coarse grid elsewhere.
Definition 8.7 (Block-Structured AMR). In block-structured AMR (Berger and Oliger, 1984; Berger and Colella, 1989), the computational domain is covered by a hierarchy of nested, logically rectangular grid patches at different refinement levels \(\ell = 0, 1, \ldots, L\). Each level \(\ell\) has mesh spacing \(\Delta x_\ell = \Delta x_0 / r^\ell\), where \(r\) is the refinement ratio (typically 2 or 4). The grid hierarchy is updated at each time step based on error indicators (typically gradients of temperature, density, or species mass fractions).
Example 8.8 (AMR for Detonation Simulation). Consider a two-dimensional detonation propagating in a channel of width 10 cm. The ZND reaction zone length is approximately 0.1 mm, requiring at least 10 grid points across the zone for adequate resolution. With a uniform grid, this would require \(\Delta x = 0.01\) mm and a total of \(10^8\) grid cells — prohibitively expensive. With three levels of AMR (refinement ratio 4), the base grid has \(\Delta x_0 = 2.56\) mm (approximately 400 cells across the channel), and the finest level has \(\Delta x_3 = 0.01\) mm, covering only the detonation front and its vicinity. The total number of cells is reduced by a factor of \(\sim 100\) compared to the uniform fine grid, making the simulation feasible. Codes such as AMReX (formerly BoxLib) and AMROC implement this strategy for reactive flow simulations.
8.4 Reduced Chemistry Mechanisms
For turbulent combustion simulations where the chemistry must be integrated millions of times per time step, detailed mechanisms with hundreds of species are too expensive. Reduced mechanisms aim to represent the essential combustion chemistry with a small number of species and reactions, preserving key quantities such as the flame speed, extinction strain rate, and ignition delay.
Definition 8.9 (Quasi-Steady-State Approximation). A species \(c_i\) is in quasi-steady state (QSS) if its net production rate is much smaller than its individual production and consumption rates:
\[
|\dot{\omega}_i| \ll \dot{\omega}_i^+ \quad \text{and} \quad |\dot{\omega}_i| \ll |\dot{\omega}_i^-|,
\]
where \(\dot{\omega}_i^+\) and \(\dot{\omega}_i^-\) are the production and consumption rates. The QSS approximation sets \(\dot{\omega}_i = 0\) and solves for \(c_i\) algebraically in terms of the remaining species, eliminating \(c_i\) from the ODE system.
Definition 8.10 (Partial Equilibrium Approximation). A reaction \(j\) is in partial equilibrium if its forward and reverse rates are nearly equal: \(r_j^f \approx r_j^r\). The partial equilibrium condition \(r_j^f = r_j^r\) provides an algebraic relation among the species concentrations, reducing the number of independent variables.
The systematic construction of reduced mechanisms from detailed ones is a major research area pioneered by Peters, Warnatz, Smooke, and others. The intrinsic low-dimensional manifold (ILDM) method of Maas and Pope and the computational singular perturbation (CSP) method of Lam and Goussis provide systematic, automatable frameworks for identifying fast and slow subspaces and constructing reduced descriptions.
Example 8.12 (Four-Step Reduced Mechanism for Methane). Peters's four-step reduced mechanism for methane-air combustion replaces the 53-species GRI-Mech with four global steps:
\[
\begin{aligned}
&\text{I:} \quad \text{CH}_4 + 2\text{H} + \text{H}_2\text{O} \to \text{CO} + 4\text{H}_2, \\
&\text{II:} \quad \text{CO} + \text{H}_2\text{O} \to \text{CO}_2 + \text{H}_2, \\
&\text{III:} \quad \text{H} + \text{H} + \text{M} \to \text{H}_2 + \text{M}, \\
&\text{IV:} \quad \text{O}_2 + 3\text{H}_2 \to 2\text{H}_2\text{O} + 2\text{H}.
\end{aligned}
\]
This mechanism involves only 7 species (\(\text{CH}_4\), \(\text{O}_2\), \(\text{CO}\), \(\text{CO}_2\), \(\text{H}_2\text{O}\), \(\text{H}_2\), \(\text{H}\)) plus the inert \(\text{N}_2\). The rate expressions are derived from the detailed mechanism using QSS and partial equilibrium approximations. The reduced mechanism reproduces the laminar flame speed of the detailed mechanism to within 5% for equivalence ratios 0.6--1.4 and pressures 1--40 atm, a remarkable compression from 325 reactions and 53 species to 4 reactions and 7 species.
8.5 Direct Numerical Simulation
Direct numerical simulation (DNS) resolves all scales of the turbulent reactive flow — from the largest energy-containing eddies down to the Kolmogorov scale and the flame thickness — without any turbulence or combustion modeling. DNS is the gold standard for understanding turbulent combustion physics, but its computational cost is staggering.
Theorem 8.13 (Cost Scaling of Reactive DNS). For a turbulent premixed flame in a cubic domain of side \(L\) with integral-scale Reynolds number \(\text{Re} = u' L / \nu\), the number of grid points required scales as
\[
N \sim \left(\frac{L}{\eta}\right)^3 = \text{Re}^{9/4},
\]
where \(\eta = L \, \text{Re}^{-3/4}\) is the Kolmogorov length scale, provided that the flame thickness \(\delta_f > \eta\) (the flamelet regime). If \(\delta_f < \eta\), additional resolution of the flame structure is needed and the cost increases further to \(N \sim \text{Re}^{9/4} (\eta / \delta_f)^3\). With \(N_s\) chemical species, the storage requirement is \(\sim (N_s + 5) N\) floating-point numbers, and the time-stepping cost is \(\sim N_s^3 N \cdot T / \Delta t\) operations (dominated by the chemistry solve if \(N_s\) is large).
Example 8.14 (State-of-the-Art Reactive DNS). As of 2024, the largest reactive DNS computations use grids of approximately \(10^{10}\)--\(10^{11}\) points with simplified chemistry (hydrogen or syngas, 9--12 species) and achieve Reynolds numbers of \(\text{Re} \approx 500\)--\(1000\). These simulations run on tens of thousands of GPU cores for weeks. A landmark computation by Chen et al. (Sandia) in the 2020s simulated a turbulent hydrogen-air jet flame at \(\text{Re} \approx 10{,}000\) with 11 species and a grid of \(1.3 \times 10^{10}\) points — the largest reactive DNS ever performed. Such simulations generate petabytes of data and provide unprecedented insight into flame-turbulence interactions, local extinction and re-ignition events, and the validity of flamelet and other modeling assumptions.
8.6 The Steady Laminar Flamelet Model in Practice
We close by connecting the theoretical flamelet concept from Chapter 7 to its computational implementation.
Definition 8.15 (Flamelet Table). A flamelet table (or flamelet library) is a pre-computed database of flame structure solutions indexed by the mixture fraction \(Z\) and scalar dissipation rate \(\chi_{st}\) (the scalar dissipation evaluated at the stoichiometric surface). For each value of \(\chi_{st}\), the steady flamelet equations
\[
\frac{\rho \chi}{2} \frac{d^2 Y_i}{dZ^2} + \dot{\omega}_i = 0, \qquad \frac{\rho \chi}{2} \frac{d^2 T}{dZ^2} - \frac{1}{c_p}\sum_i h_i \dot{\omega}_i = 0,
\]
are solved with detailed or reduced chemistry to obtain \(Y_i(Z; \chi_{st})\) and \(T(Z; \chi_{st})\). The flamelet table \(\{Y_i(Z, \chi_{st})\}\) is then used as a lookup in the turbulent flow simulation.
Example 8.17 (Flamelet Table Construction). For methane-air combustion using GRI-Mech 3.0 (53 species, 325 reactions), a typical flamelet table is constructed as follows. The counterflow diffusion flame equations are solved for 50 values of \(\chi_{st}\) ranging from near-zero (equilibrium limit) to \(\chi_{\text{ext}}\) (the extinction value). Each solution is discretized on 200 points in \(Z\)-space. The resulting table has dimensions \(53 \times 200 \times 50\) and occupies approximately 4 MB of storage. The table construction takes minutes on a single workstation using a flamelet code such as FlameMaster (Peters) or Cantera. During the CFD simulation, table lookups (with bilinear interpolation in \(Z\) and \(\chi_{st}\)) cost negligible time compared to the flow solver, making the flamelet approach both accurate and efficient for non-premixed turbulent combustion.
