化学種輸送方程式
Theory and Physics
Overview
Professor, the species transport equation is the foundation of combustion CFD, right?
Exactly. The species transport equation describes how the mass fraction of each chemical component (CH4, O2, CO2, H2O, CO, NO, etc.) in a fluid changes. It is the foundation for all combustion CFD models (EDC, Flamelet, Species Transport + Finite Rate Chemistry, etc.).
Governing Equations
Please explain the species transport equation.
The transport equation for the mass fraction $Y_i$ of species $i$ can be written as follows.
The meaning of each term is as follows.
- Left-hand side, 1st term: Temporal change (unsteady term)
- Left-hand side, 2nd term: Convective transport
- $\mathbf{J}_i$: Diffusion flux (Fick's law: $\mathbf{J}_i = -\rho D_{i,m} \nabla Y_i$)
- $R_i$: Source/sink term due to chemical reaction
- $S_i$: Other sources (e.g., spray evaporation)
Reaction Source Term
How is the reaction source term $R_i$ expressed?
$R_i$ is the sum of contributions from all reactions and is given by the following equation.
Here, $k_{f,r}$ is the forward reaction rate constant of Arrhenius type, $\nu'$, $\nu''$ are the stoichiometric coefficients of reactants and products, and $[C_j]$ is the molar concentration of species $j$.
Do you consider the reverse reaction?
Reverse reactions are also important under conditions close to equilibrium. The reverse reaction rate constant $k_{b,r}$ is obtained from the equilibrium constant $K_{eq,r}$ as $k_{b,r} = k_{f,r}/K_{eq,r}$. In high-temperature combustion, the dissociation of CO2 ($\text{CO}_2 \rightleftharpoons \text{CO} + \frac{1}{2}\text{O}_2$) becomes important as a reverse reaction.
Diffusion Model
How do you model the diffusion flux?
There are three levels.
| Model | Computational Cost | Accuracy | Application |
|---|---|---|---|
| Fick's Law ($D_{i,m}$) | Low | Moderate | RANS Standard |
| Modified Fick's Method (Mass Conservation Correction) | Low | Good | Fluent Standard |
| Stefan-Maxwell Equations | High | Highest | When accuracy for light species (H2, He) is critical |
In turbulent flows, turbulent diffusion dominates over molecular diffusion, right?
Correct. The effective diffusion coefficient in a turbulent flow is $D_{\text{eff}} = D_{i,m} + D_t$, where $D_t = \mu_t/(\rho\,Sc_t)$ is the turbulent diffusion. $Sc_t \approx 0.7$ is a common value in combustion analysis. In high-Re turbulence, $D_t >> D_{i,m}$, so the choice of molecular diffusion model doesn't have much impact in RANS. However, in LES or laminar flames, the accuracy of molecular diffusion becomes important.
Mass Fraction Constraint
When there are $N_s$ species, solving $N_s - 1$ transport equations is sufficient. The last species is determined algebraically from the constraint $\sum Y_i = 1$. Typically, the component with the largest mass fraction (e.g., N2) is determined algebraically.
So the species transport equation is composed of three elements—convection, diffusion, and reaction—and all combustion CFD models are built upon this equation.
Exactly. Models like flamelet or PDF models are merely efficiency improvements that use table lookups instead of directly solving this equation. The essence lies in this transport equation.
Fick's Law of Diffusion—How a 1855 Law Became the Foundation of Combustion CFD
The diffusion law ($J = -D \nabla c$) published by Adolf Fick in 1855 was originally derived from experiments on salt dissolving in water. It wasn't until the early 20th century that it was found applicable to chemical species transport in gas mixtures. Extensions to multicomponent mixtures (Stefan-Maxwell model) and simplifications (Hirschfelder approximation) accumulated, leading to the modern species transport equation. To think that a 170-year-old saltwater experiment formula now operates as the fundamental equation for combustion simulations in gas turbines and furnaces today makes one appreciate anew the universality of physical laws.
Physical Meaning of Each Term
- Temporal term $\partial(\rho\phi)/\partial t$: Imagine the moment you turn on a faucet. At first, water comes out in an unstable, spluttering manner, but after a while, the flow becomes steady, right? This "period of change" is described by the temporal term. The pulsation of blood flow with each heartbeat, or the flow fluctuation each time an engine valve opens and closes—all are unsteady phenomena. So what is steady-state analysis? It's looking only at "after sufficient time has passed and the flow has settled down"—in other words, setting this term to zero. Since computational cost drops significantly, solving first in steady-state 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 the flow becomes faster, this term rapidly strengthens, making control difficult. This is the root cause of turbulence. A common misconception: "Convection and conduction are similar things" → 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, after a while it naturally mixes, right? 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. When viscosity is large, the diffusion term becomes strong, and the fluid moves in a "sluggish" manner. In low Reynolds number flows (slow, viscous), diffusion dominates. Conversely, in high-Re 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. "Where there is a pressure difference, flow is generated"—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. If results go wrong immediately after switching to compressible analysis, confusion between absolute/gauge pressure might be the cause.
- 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, if you forget to include buoyancy, the fluid doesn't move at all—a physically impossible result where warm air doesn't rise in a room with the heater on in winter.
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 flows (shock capturing required), free surface flows (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 coefficient $\nu = \mu/\rho$ [m²/s] |
| Reynolds number $Re$ | Dimensionless | $Re = \rho u L / \mu$. Criterion 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
Please tell me what to watch out for in the numerical solution of the species transport equation.
There are three main challenges in the numerical solution of species transport: (1) Suppression of numerical diffusion, (2) Guarantee of mass conservation, (3) Handling of stiff reaction source terms.
Spatial Discretization Schemes
What spatial discretization scheme do you recommend for species?
To resolve sharp fronts of species (like flame fronts), schemes with low numerical diffusion are necessary.
| Scheme | Accuracy | Numerical Diffusion | Stability | Recommendation |
|---|---|---|---|---|
| First Order Upwind | 1st order | Large | High | Initial convergence only |
| Second Order Upwind | 2nd order | Moderate | Good | RANS standard |
| QUICK | 3rd order | Small | Somewhat unstable | Use with caution |
| Central Difference | 2nd order | None | Unstable (oscillations) | For LES momentum equations |
| Bounded Central Difference | 2nd order | Small | Good | Recommended for LES |
So Bounded Central Difference is recommended for LES.
Correct. Central Difference has zero numerical diffusion but produces non-physical oscillations (Gibbs phenomenon). The bounded version suppresses oscillations with limiters while maintaining low diffusion. In Fluent, Bounded Central Differencing is available as an option.
Mass Conservation
What is the mass conservation problem?
When solving $N_s - 1$ transport equations and determining the last species by $Y_{N_s} = 1 - \sum_{i=1}^{N_s-1} Y_i$, numerical errors from each equation can accumulate, sometimes causing $Y_{N_s}$ to become negative.
Countermeasures include:
- Species Bounding: Clip each $Y_i$ to [0, 1] (Fluent standard)
- Flux-Corrected Transport (FCT): High-order, conservative scheme
- Calculate N2 last: Determine the component with the largest mass fraction algebraically (errors are less noticeable)
Operator Splitting
Please explain Operator Splitting for the reaction source term.
It's a technique that splits the species transport equation into "transport part" and "reaction part" and solves them alternately.
Strang splitting has 2nd-order accuracy, but splitting error becomes problematic when $\Delta t$ is large. Typical settings:
What is the mass conservation problem?
When solving $N_s - 1$ transport equations and determining the last species by $Y_{N_s} = 1 - \sum_{i=1}^{N_s-1} Y_i$, numerical errors from each equation can accumulate, sometimes causing $Y_{N_s}$ to become negative.
Countermeasures include:
Please explain Operator Splitting for the reaction source term.
It's a technique that splits the species transport equation into "transport part" and "reaction part" and solves them alternately.
Strang splitting has 2nd-order accuracy, but splitting error becomes problematic when $\Delta t$ is large. Typical settings:
| Parameter | RANS | LES |
|---|---|---|
| CFD Timestep | -- (Steady) | $10^{-5}$ - $10^{-6}$ s |
| ODE Substep | Automatic (CVODE internal) | Automatic |
| Splitting Method | Sequential (Fluent) | Strang (OpenFOAM) |