Chemical Reaction Kinetics Fundamentals
Theory and Physics
Overview
Teacher! Chemical kinetics is the fundamental theory underlying combustion simulations, right?
Exactly. Chemical kinetics is the discipline that mathematically describes how fast fuel and oxidizer react. In combustion CFD, if this reaction rate is not modeled correctly, results for flame temperature and exhaust gas composition will be completely off the mark.
Specifically, what kind of equations are central?
The core is the Arrhenius-type reaction rate equation. For each elementary reaction, the rate constant is expressed by three parameters: the pre-exponential factor $A$, the temperature exponent $n$, and the activation energy $E_a$.
Governing Equations
First, please explain the Arrhenius equation.
The reaction rate constant $k$ can be written as follows.
Here, $A$ is the pre-exponential factor (frequency factor), $n$ is the temperature exponent, $E_a$ is the activation energy [J/mol], and $R$ is the universal gas constant 8.314 J/(mol K). For example, in the basic hydrogen reaction H$_2$ + O$_2$ system, $E_a$ is on the order of tens of kJ/mol.
Does a larger activation energy mean the reaction is slower?
Yes. The larger $E_a$ is, the smaller the exponential term $\exp(-E_a/RT)$ becomes, and at low temperatures, the reaction hardly proceeds. Conversely, as temperature increases, the rate increases exponentially. This directly relates to combustion ignition delay and flame stability.
The time rate of change of the mass fraction $Y_i$ of chemical species $i$ can be written by summing the contributions of all elementary reactions as follows.
$\nu'$ and $\nu''$ are the stoichiometric coefficients for reactants and products, right?
Correct. $[X_j]$ is molar concentration, $M_{w,i}$ is molecular weight. Multiple elementary reactions intertwine to constitute the species production/consumption rate $\dot{\omega}_i$.
Stiff Chemical Reaction Systems
I've heard that time scales in combustion chemical reactions can be extremely different.
That's precisely the problem of Stiffness. For example, the detailed reaction mechanism for methane/air, GRI-Mech 3.0, contains 53 species and 325 elementary reactions, but the time constant for radical reactions is on the order of $10^{-9}$ seconds, while the main reactions are on the order of $10^{-3}$ seconds—a difference of over six orders of magnitude.
With such a large difference, ordinary explicit methods can't solve it, right?
Exactly. Explicit Euler or Runge-Kutta methods would require extremely small time steps, making the computation time impractical. Therefore, implicit methods (BDF methods, Rosenbrock methods, etc.) or tabulation techniques like ISAT (In-Situ Adaptive Tabulation), which will be covered in detail in the next article, are used.
Hierarchy of Reaction Mechanisms
Are there different levels of reaction mechanisms?
Reaction mechanisms have the following hierarchy.
| Category | Number of Species | Number of Reactions | Representative Examples | Applications |
|---|---|---|---|---|
| Global Single-Step | 2-5 | 1-2 | Westbrook-Dryer | Rough estimates / Initial studies |
| Reduced Mechanism | 10-30 | 20-100 | DRM-19, Lu-Law | 3D RANS/LES |
| Skeletal Mechanism | 30-100 | 100-500 | skeletal-iso-octane | Detailed RANS |
| Detailed Mechanism | 50-300+ | 300-3000+ | GRI-Mech 3.0, USC Mech II | 0D/1D, DNS |
Is it difficult to directly solve detailed mechanisms in 3D combustion LES?
Because the computational cost becomes enormous. Trying to solve a detailed mechanism with 100 species on a 3D mesh with 1 million cells would require solving a system of 100 ODEs in each cell every time step. Therefore, in practice, reduced mechanisms are used, or dimensionality is reduced using ISAT or FGM (Flamelet Generated Manifold).
I've gained a solid understanding of the basic theory of chemical kinetics. I see clearly that the stiffness problem is at its core.
Good. The details of numerical methods will be covered in the next article. If you make a mistake in setting the Arrhenius equation parameters, the ignition temperature can be off by hundreds of Kelvin, so always verify against experimental data.
Arrhenius's "Activation Energy" in 1889 — A Legacy for Modern CAE
The paper published by Svante Arrhenius in 1889 was considered a "peculiar theory among chemists" at the time. He expressed the empirical rule that "reaction rate increases sharply with temperature" with the simple formula $k = A \exp(-E_a/RT)$. Yet, over 130 years later, this formula remains at the heart of reaction rate calculations in combustion CFD. Each of the 325 reactions in GRI-Mech 3.0 has three parameters: A, Ea, and n. In other words, modern combustion simulation solves over 900 simultaneous equations based on the 1889 formula. Classical physical chemistry powers modern supercomputers—a story that makes you feel the continuity of science.
Physical Meaning of Each Term
- Time term $\partial(\rho\phi)/\partial t$: Imagine turning on a faucet. At first, water comes out in an unstable, spluttering manner, but after a while, it becomes a steady flow, right? This "period of change" is described by the time term. The pulsation of blood flow from a heartbeat, or the flow fluctuation each time an engine valve opens and closes—all are unsteady phenomena. So what is steady-state analysis? It looks only at "after sufficient time has passed and the flow has settled down"—meaning this term is set to zero. Since computational cost drops significantly, solving first with a steady-state approach is a basic CFD strategy.
- Convection term $\nabla \cdot (\rho \mathbf{u} \phi)$: What happens if you drop a leaf into a river? It gets carried downstream by the flow, right? This is "convection"—the effect where fluid motion transports things. Warm air from a heater reaching the far corner of a room is also because the "carrier," air, transports heat via convection. Here's the interesting part—this term contains "velocity × velocity," making it nonlinear. That is, as flow speed increases, this term rapidly strengthens, making control difficult. This is the root cause of turbulence. A common misconception: "Convection and conduction are similar" → They are completely different! Convection is carried by flow, conduction is transmitted by molecules. There's an order-of-magnitude difference in efficiency.
- Diffusion term $\nabla \cdot (\Gamma \nabla \phi)$: Have you ever put milk in coffee and left it? Even without stirring, it naturally mixes after a while. That's molecular diffusion. Now a question—honey and water, which flows more easily? Obviously water, right? Honey has high viscosity ($\mu$), so it flows poorly. Higher viscosity strengthens the diffusion term, and the fluid moves in a "thick" manner. In low Reynolds number flows (slow, viscous), diffusion dominates. Conversely, in high Re number flows, convection overwhelmingly dominates, and diffusion plays a supporting role.
- Pressure term $-\nabla p$: When you push the plunger of a syringe, liquid shoots out forcefully from the needle tip, right? Why? Because the plunger side is high pressure, the needle tip is low pressure—this pressure difference provides the force that pushes the fluid. Dam discharge works on the same principle. On a weather map, where isobars are tightly packed? That's right, strong winds blow. "Flow arises where there is a pressure difference"—this is the physical meaning of the pressure term in the Navier-Stokes equations. A point of confusion here: "Pressure" in CFD is often gauge pressure, not absolute pressure. When you switch to compressible analysis and suddenly get strange results, it might be due to mixing up absolute/gauge pressure.
- Source term $S_\phi$: Heated air rises—why? Because it becomes lighter (lower density) than its surroundings, so buoyancy pushes it upward. This buoyancy is added to the equation as a source term. Other examples: chemical reaction heat generated by a gas stove flame, Lorentz force acting on molten metal in a factory's electromagnetic pump... These are all actions that "inject energy or force into the fluid from the outside," expressed by source terms. What happens if you forget the source term? In a natural convection analysis, forgetting buoyancy means the fluid won't move at all—a physically impossible result, like turning on a heater in a winter room but the warm air doesn't rise.
Assumptions and Applicability Limits
- Continuum assumption: Valid for Knudsen number Kn < 0.01 (mean free path ≪ characteristic length)
- Newtonian fluid assumption: Shear stress and strain rate have a linear relationship (non-Newtonian fluids require viscosity models)
- Incompressibility assumption (for Ma < 0.3): Treat density as constant. For Mach numbers above 0.3, consider compressibility effects
- Boussinesq approximation (natural convection): Consider density changes only in the buoyancy term, using constant density in other terms
- Non-applicable cases: Rarefied gases (Kn > 0.1), supersonic/hypersonic flow (shock capturing required), free surface flow (VOF/Level Set, etc. required)
Dimensional Analysis and Unit Systems
| Variable | SI Unit | Notes / Conversion Memo |
|---|---|---|
| Velocity $u$ | m/s | When converting from volumetric flow rate for inlet conditions, pay attention to cross-sectional area units |
| Pressure $p$ | Pa | Distinguish between gauge and absolute pressure. Use absolute pressure for compressible analysis |
| Density $\rho$ | kg/m³ | Air: approx. 1.225 kg/m³@20°C, Water: approx. 998 kg/m³@20°C |
| Viscosity coefficient $\mu$ | Pa·s | Be careful not to confuse with kinematic viscosity $\nu = \mu/\rho$ [m²/s] |
| Reynolds number $Re$ | Dimensionless | $Re = \rho u L / \mu$. Indicator for laminar/turbulent transition |
| CFL number | Dimensionless | $CFL = u \Delta t / \Delta x$. Directly related to time step stability |
Numerical Methods and Implementation
Details of Numerical Methods
How are stiff ODE systems solved in CFD? Please explain the specific methods.
Numerical integration of chemical reactions in combustion CFD is broadly divided into three approaches: (1) Direct Integration (DI), (2) Tabulation methods, and (3) Operator splitting methods.
Direct Integration Method
First, please explain the direct integration method.
This method solves the ODE system for chemical species in each CFD cell every time step. To handle stiff systems, implicit multi-step methods are used.
Let's compare typical implicit solvers.
| Method | Order | Characteristics | Representative Implementation |
|---|---|---|---|
| BDF (Backward Differentiation Formula) | 1st-5th order | A-stable at high order | CVODE (SUNDIALS) |
| Rosenbrock Method | 2nd-4th order | Jacobian evaluated once | DASPK |
| SDIRK | 2nd-4th order | Diagonally Implicit RK | OpenFOAM standard |
| Exponential Integrators | Variable | Matrix exponential function | Research stage |
CVODE is famous. I've heard it's used in Fluent too.
Yes. In Ansys Fluent, a CVODE-based integrator is incorporated as the "Stiff Chemistry Solver." STAR-CCM+ also internally uses the SUNDIALS library. Evaluating the Jacobian matrix becomes a bottleneck, so automatic generation of analytical Jacobians is important.
ISAT (In-Situ Adaptive Tabulation)
How does ISAT work?
A method proposed by Pope in 1997, it records the input-output of once-calculated chemical reactions in a binary tree structure table. When a new composition point arrives, if the answer can be returned via linear approximation from existing records, direct integration is skipped.
What is the approximation accuracy?
The tolerance $\epsilon_{\text{tol}}$ is set by the user. Fluent's default is around $10^{-4}$. ISAT can often reduce direct integration calls by over 90%. However, table size can strain memory, and in large-scale parallel computing, each process holds an independent table, making memory efficiency a challenge.
Operator Splitting Method
What is the operator splitting method?
A method that solves fluid transport and chemical reactions in separate steps. Strang splitting is representative, alternately executing CFD transport steps and ODE integration for chemical reactions.
Related Topics
なった
詳しく
報告