OpenFOAM Combustion Modeling — Premixed, Non-Premixed, and Spray

#+CATEGORY: openfoam # :PROPERTIES: # :ID: uuid-openfoam-combustion # :END:

OpenFOAM Combustion Modeling :: Premixed, Non-Premixed, and Spray Combustion

This note documents OpenFOAM's combustion modelling capabilities. Combustion in OpenFOAM is one of the most challenging applications: it couples turbulence with chemical kinetics, heat release with buoyancy, species transport with radiation, and finite-rate chemistry with turbulent mixing. The computational cost of detailed-chemistry combustion can be 100–1000× larger than simple incompressible flow because the chemistry source term =\dot{\omega}_k/ is computed at every cell at every timestep from a stiff system of ODEs.

Combustion modelling is a classification problem: is the reaction rate limited by mixing (non-premixed), by chemistry (premixed), or by both? The answer determines which model family is appropriate. This maps to Conservation Laws (species and energy conservation) and Turbulence Models (turbulence-chemistry interaction).

**Table of Contents**

• Combustion Regimes :: Premixed vs non-premixed vs partially-premixed
• Non-Premixed Combustion :: Flamelet, FPV, EDC models
• Premixed Combustion :: XiFoam, flame surface density
• Partially-Premixed Combustion :: EDC, PaSR hybrid models
• Chemistry Solvers :: VODE, SUNDIALS, CHEMKIN integration
• Spray Combustion :: Lagrangian droplets, evaporation, ignition
• Model Selection :: When to use which combustion model
• See Also
• References

Combustion Regimes :: Premixed vs Non-Premixed vs Partially-Premixed

Combustion regimes are classified by how fuel and oxidiser are mixed before reaction:

The Damköhler number =Da = \tau_t / \tau_c/ (turbulence time scale / chemical time scale) is the dimensionless parameter that determines whether a combustion regime is premixed or non-premixed:

• Da >> 1: Chemistry is fast compared to turbulence — reaction is mixing-limited (non-premixed)
• Da << 1: Chemistry is slow compared to turbulence — reaction is chemistry-limited (premixed)
• Da ~ 1: Both mixing and chemistry are important (partially-premixed)

The Karlovitz number =Ka = (\tau_c / \tau_\eta)^{-1/2}/ (chemical time / Kolmogorov time) determines the flame structure within the non-premixed regime:

• Ka << 1: Flamelets are not disrupted — thin reaction zone (flamelet model appropriate)
• Ka >> 1: Flamelets are thickened or extinguished — distributed reaction zone (PDF methods needed)

Non-Premixed Combustion :: Flamelet, FPV, and EDC Models

Non-premixed combustion occurs when fuel and oxidiser are injected separately and mix during the reaction process (e.g., a Bunsen burner, a gas turbine combustor, a ship engine exhaust). The flame is a diffusion flame — the reaction occurs at the stoichiometric surface where the local fuel/air ratio equals the stoichiometric ratio.

**Flamelet Model**

The flamelet model assumes the turbulent flame is a collection of laminar flamelets — thin reaction zones embedded in the turbulent flow. The flame structure is solved as a 1D laminar flame (the "flamelet") parameterised by the mixture fraction =Z/ and the scalar dissipation rate =\chi = -2D |\nabla Z|^2/:

// Transport equations for mixture fraction and its variance
// (required by the flamelet model)
MULES
{
    type MULES;
    flux phi;
    alpha phiAlpha;
    deltaCoeff 1;
}

// Scalar dissipation rate
chi = 2 * alpha * |grad(Z)|^2

The mixture fraction =Z/ is a conserved scalar (not affected by reaction) defined as:

= Z = \frac{s Y_F - Y_O + Y_{O,in}}{s Y_{F,in} + Y_{O,in}} /

where =s/ is the stoichiometric mass ratio, =Y_F/ is the fuel mass fraction, and =Y_O/ is the oxidiser mass fraction. At stoichiometric mixture (=Z_{st}/), the flame temperature is at its maximum (adiabatic flame temperature).

The flamelet library (=P/\chi/ relationship) is pre-computed from detailed 1D laminar flame simulations using a chemistry solver (VODE or SUNDIALS). This library maps (=Z, \chi/) → (T, Y_k/) for each species k. During the CFD simulation, the flamelet library is interpolated at each cell based on the local =Z/ and =\chi/:

For RANS simulations with a k-epsilon or k-omega turbulence model, the mixture fraction variance is computed from turbulent kinetic energy:

= \tilde{Z''^2} = C_\mu^{1/2} k / \sigma_Z \times \tilde{\epsilon} / /

where =\sigma_Z/ is the turbulent Schmidt number for mixture fraction (typically ~0.7), and =C_\mu/ is the standard RANS constant (0.09).

**Flamelet Progress Variable (FPV) Model**

The FPV model extends the flamelet model by adding a progress variable =Y_P/ (e.g., mass fraction of CO2 or a lumped product species) in addition to the mixture fraction:

= (T, Y_k) = f(Z, Y_P, \chi) /

This adds one dimension to the flamelet library (now 3D instead of 2D), which captures finite-rate chemistry effects more accurately than the 2D flamelet model (which assumes equilibrium chemistry at each =\chi/ level). The progress variable equation:

= \frac{\partial (\rho \tilde{Y}_P)}{\partial t} + \nabla \cdot (\rho \vec{v} \tilde{Y}_P) = \nabla \cdot \left(\frac{\mu_t}{Sc_t} \nabla \tilde{Y}_P\right) + \tilde{\dot{\omega}}_P /

where =\tilde{\dot{\omega}}_P/ is the filtered reaction rate from the flamelet library.

**EDC (Eddy Dissipation Concept)**

The EDC model, originally from Magnussen & Hjertager (1976), assumes turbulence consists of fine structures (reaction zones) embedded in a bulk fluid. Reaction occurs only in the fine structures, which occupy a fraction =\gamma/ of the volume:

= \gamma = C_\gamma \frac{\nu}{\nu_t} \sqrt{\frac{\nu \epsilon}{\nu_t^2}} /

where =C_\gamma/ is a model constant (typically ~2.7). The fine-structure residence time is =\tau_c / 4/, and the reaction rate in the fine structures is computed from detailed chemistry:

= \tilde{\dot{\omega}}_k = \rho_f \dot{\omega}_{k,f} \frac{\gamma^2}{1+\gamma} \approx \rho_f \dot{\omega}_{k,f} \gamma^2 /

where =\rho_f/ is the fine-structure density and =\dot{\omega}_{k,f}/ is the fine-structure reaction rate (from detailed chemistry). EDC does NOT require pre-tabulated flamelets — it solves the chemical ODEs at each cell at each timestep. This is more expensive than flamelet models but more general (no assumption of thin-flame structure).

Premixed Combustion :: XiFoam and Flame Surface Density

Premixed combustion occurs when fuel and oxidiser are mixed before ignition — the flame front propagates through the mixture (e.g., spark-ignition engine, gas stove). The flame is modelled using Flame Surface Density (FSD):

// Flame Surface Density transport equation
\frac{\partial (\rho \Sigma)}{\partial t} + \nabla \cdot (\rho \vec{v} \Sigma) = \nabla \cdot (\rho D_\Sigma \nabla \Sigma) + S_T - S_D

where:
- \Sigma = \rho |\nabla c| is the flame surface density (m²/m³)
- D_\Sigma is the turbulent diffusivity of the flame surface
- S_T is the flame production term (turbulent stretching)
- S_D is the flame destruction term (curvature, extinction)

XiFoam solves the FSD equation coupled with the conservation equations. The turbulent burning velocity is modelled as:

= S_T = S_L f(\frac{u'}{S_L} / = S_L \cdot f(u'/S_L, l_t) /

where =S_L/ is the laminar burning velocity (function of temperature, pressure, and equivalence ratio), =u'/ is the turbulent intensity, and =l_t/ is the turbulent integral length scale. The function =f/ is model-dependent:

XiFoam is suitable for spark-ignition (SI) engines, lean premixed prevaporized (LPP) combustors, and other premixed combustion applications. For partially-premixed combustion (some mixing before ignition, some during), see the EDC model above. For engine simulations with moving geometry, use engineFoam (which couples XiFoam with dynamic mesh).

Partially-Premixed Combustion :: EDC and PaSR

Partially-premixed combustion is the most general regime — it includes cases where the mixture is neither perfectly premixed nor perfectly non-premixed. The Eddy Dissipation Concept (EDC) handles it without explicit regime classification. A more advanced model is the Partially-Stirred Reactor (PaSR) model:

**PaSR (Partially-Stirred Reactor)**

The PaSR model divides each cell into two regions: a reactive zone (fine structures) where chemistry occurs, and a non-reactive zone (bulk fluid) where no chemistry occurs. The fraction of the cell that is in the reactive zone is =\gamma/ (same as the EDC model):

= \gamma = C_\gamma \frac{\nu}{\nu_t} \sqrt{\frac{\nu \epsilon}{\nu_t^2}} /

in the reactive zone, detailed chemistry is solved. In the non-reactive zone, the composition is =\tilde{\phi}/ (the cell-averaged composition). The cell-averaged reaction rate is =\tilde{\dot{\omega}} = \gamma \dot{\omega}_{reactive} + (1-\gamma) \dot{\omega}_{non-reactive}/. Since =\dot{\omega}_{non-reactive} = 0/:

= \tilde{\dot{\omega}} = \gamma \dot{\omega}_{reactive} /

=PaSR/ is more general than EDC because it uses the cell-averaged conditions to select the reactive-zone composition, whereas EDC uses a fixed fine-structure composition. PaSR is available in newer OpenFOAM versions (ESI 2020+) through the combustion library.

Chemistry Solvers :: VODE, SUNDIALS, and CHEMKIN Integration

The chemical kinetics source term =\dot{\omega}_k/ is a system of stiff ODEs:

= \frac{d[Y_k]}{dt} = \frac{\dot{\omega}_k}{\rho} /

where the reaction rate =\dot{\omega}_k = \sum_j \nu_{kj} (k_{f,j} - k_{b,j}) \prod_i [X_i]^{\alpha_{ij}}/ is computed from the Arrhenius equation for each reaction j, with forward rate constant =k_{f,j} = A_j T^{\beta_j} \exp(-E_{a,j}/RT)/.

For detailed mechanisms (50–500 species, 1000–5000 reactions), this system can have timescales ranging from nanoseconds (radical formation) to milliseconds (bulk conversion). This stiffness requires implicit ODE solvers. OpenFOAM supports:

CHEMKIN-format mechanism files are the standard format for detailed chemical kinetics. A mechanism file contains:

// Mechanism file: mech.xml (CHEMKIN format)
ELEMENTS  12 C  H  O  N  Ar  N2  O2  OH  H  O  HO2  H2O  N
SPECIES   12 C H O2 N2 H2 H2O CO CO2 OH CH4  C2H6 AR
REACTION  24
...

The mechanism file specifies elements, species, and reactions. Each reaction contains: reactants, products, rate coefficients (Arrhenius parameters), and reaction type (third-body, pressure-dependent, etc.). The mechanism is parsed at solver startup, and the Jacobian matrix for the chemistry system is constructed. OpenFOAM's reactingFoam, XiFoam, and engineFoam all support CHEMKIN-format mechanisms.

For the reactingFoam solver, the thermophysical package is set via:

// thermophysicalProperties
thermoType
{
    type            heMulti mixture Therm;  // multi-mixture enthalpy
    mixture         polyMixture;            // multiple species
    package         basicThermo;
    equationOfState perfectGas;
    specie          specie;
    energy          enthalpy;
    equationOfState Hconst;                 // constant Cp (or: hConst)
    transport       const;                  // constant mu, k
    thermodynamics  hConst;                 // constant Cp per species
    diffusion       constant;               // constant diffusivity
}

Spray Combustion :: Lagrangian Droplets, Evaporation, Ignition

Spray combustion (e.g., diesel engines, gas turbine fuel injection) couples Lagrangian droplet tracking with turbulent non-premixed combustion. OpenFOAM's reactingTwoPhaseEulerFoam and DPMFoam with combustion support this:

The droplet evaporation model computes the droplet evaporation rate from the mass transfer equation:

= \frac{dm_p}{dt} = - \pi d_p Sh D_v \rho_g (Y_{v,s} - Y_{v,\infty}) /

and the heating rate:

= \frac{dT_p}{dt} = \frac{Nu k_g}{\rho_p c_p d_p} (T_g - T_p) + \frac{\dot{m}_p \Delta H_v}{\rho_p c_p \pi d_p^3 / 6} /

where =Nu/ is the Nusselt number, =Sh/ is the Sherwood number, and =\Delta H_v/ is the latent heat of vaporisation. Spray combustion is one of the most computationally expensive applications in CFD because it couples DPM (lagrangian), non-premixed combustion (flamelet or EDC), and turbulence (RANS or LES). Each timestep requires:

1. Solve the continuous phase Navier-Stokes equations (U, p, turbulence) 2. Solve the droplet equations (position, velocity, temperature, mass) — one ODE per parcel 3. Compute the interphase exchange terms (drag, heat transfer, mass transfer) 4. Feed the source terms back to the continuous phase 5. Evaluate the chemistry source term at each cell (stiff ODEs) 6. Repeat until converged

In practice, the chemistry step dominates the total computational cost for spray combustion — the stiff ODE solver for the chemistry source term requires many sub-timesteps per CFD timestep. For a detailed mechanism (~100 species), chemistry can account for 80–90% of the total computational cost. Subgrid-scale models (flamelet) are mandatory for production simulations — detailed chemistry is only feasible for code verification or DNS.

Turbulent Combustion Model Selection :: When to Use Which Model

The default recommendation for all non-premixed combustion in RANS simulations is the flamelet model with a finite-rate chemistry library. It offers the best balance between accuracy and computational cost for industrial applications. For LES combustion, the filtered PDF method or the Dynamic Flame Surface Density approach (in newer OpenFOAM versions) is preferred. For code verification, DNS requires fully-resolved detailed chemistry (which is limited to small domains, typically the size of a laminar flame).

See Also :: Related Notes

. Conservation Laws — species and energy conservation . Turbulence Models — turbulence-chemistry interaction . OpenFOAM Solver Selection — reactingFoam, XiFoam, engineFoam . OpenFOAM Numerical Schemes — schemes for species transport . OpenFOAM Case Setup — combustionProperties . OpenFOAM Multiphase Flows — spray combustion . OpenFOAM Dynamic Meshes — engineFoam dynamic mesh . ANSYS Combustion — non-premixed, premixed, EDC comparison

References

. OpenFOAM User Guide. ESI OpenCFD. (Combustion solver documentation.) . OpenFOAM Programming Guide. ESI OpenCFD. (Combustion library API.) . Pope, S.B. (2000). Turbulent Flows. Chapter 12 (Turbulent combustion). . Candel, S. (2002). "Combustion dynamics and control." Progress in Energy and Combustion Science, 28(6), 545-576. . Bray, K.N.C. (1980). "Prediction of autonomous flames with large chemical reaction rates." Combustion Theory, 63, 201-212. (Flamelet modelling.) . Magnussen, B.F., & Hjertager, B.H. (1976). "On mathematical modeling of turbulent combustion." 16th Symposium (International) on Combustion. (EDC model.) . Gicquel, L.Y.M., Pope, N., & Bray, K.N.C. (2000). "LES of turbulent flame structure using flamelet PDF." 19th Symposium (International) on Combustion. (FPV model.)