Use this Skill for combustion simulations with Cantera: mechanism loading, freely propagating flames, ignition delay, reactors, and sensitivity analysis.
Use this Skill for combustion simulations with Cantera: mechanism loading, freely propagating flames, ignition delay, reactors, and sensitivity analysis.
import cantera as ct
import numpy as np
import matplotlib.pyplot as plt
defadiabatic_flame_temperature(fuel, oxidizer="air", phi_range=None, mechanism="gri30.yaml",
T_initial=300.0, P=101325.0):
"""
Compute adiabatic flame temperature over a range of equivalence ratios.
Parameters
----------
fuel : str
Fuel species name (e.g., "CH4", "H2", "C2H5OH")
oxidizer : str
"air" or specific composition string
phi_range : array-like
Equivalence ratios to sweep
mechanism : str
Cantera mechanism file (.yaml or .cti)
Returns
-------
pd.DataFrame with phi, T_adiabatic, and selected species
"""if phi_range isNone:
phi_range = np.arange(0.5, 2.01, 0.05)
gas = ct.Solution(mechanism)
results = []
for phi in phi_range:
try:
gas.set_equivalence_ratio(phi, fuel, oxidizer)
gas.TP = T_initial, P
# Equilibrate at constant enthalpy and pressure (adiabatic combustion)
gas.equilibrate("HP")
results.append({
"phi": phi,
"T_adiabatic_K": gas.T,
"X_CO2": gas["CO2"].X[0],
"X_CO": gas["CO"].X[0],
"X_H2O": gas["H2O"].X[0],
"X_O2": gas["O2"].X[0],
})
except ct.CanteraError as e:
print(f" phi={phi:.2f}: {e}")
import pandas as pd
df = pd.DataFrame(results)
print(f"\nAFT for {fuel}/{oxidizer}:")
print(f" Max T_ad = {df['T_adiabatic_K'].max():.0f} K at φ = {df.loc[df['T_adiabatic_K'].idxmax(), 'phi']:.2f}")
return df
# Compare methane and hydrogen
fig, axes = plt.subplots(2, 1, figsize=(10, 8))
for fuel, color, label in [("CH4", "b", "Methane (GRI-3.0)"),
("H2", "r", "Hydrogen (GRI-3.0)")]:
aft_df = adiabatic_flame_temperature(fuel, phi_range=np.arange(0.5, 2.01, 0.05))
axes[0].plot(aft_df["phi"], aft_df["T_adiabatic_K"], f"{color}-o",
ms=4, label=label)
axes[1].plot(aft_df["phi"], aft_df["X_CO2"] * 100, f"{color}--", label=f"CO₂ ({label[:4]})")
axes[1].plot(aft_df["phi"], aft_df["X_CO"] * 100, f"{color}:", label=f"CO ({label[:4]})")
axes[0].axvline(1.0, color='k', linestyle='--', linewidth=0.8, label="φ=1 (stoich)")
axes[0].set_xlabel("Equivalence ratio φ"); axes[0].set_ylabel("T_adiabatic (K)")
axes[0].set_title("Adiabatic Flame Temperature"); axes[0].legend()
axes[1].set_xlabel("Equivalence ratio φ"); axes[1].set_ylabel("Mole fraction (%)")
axes[1].set_title("Product Species Mole Fractions"); axes[1].legend(fontsize=8)
plt.tight_layout()
plt.savefig("adiabatic_flame_temperature.png", dpi=150)
plt.show()
Step 2: Ignition Delay Time
import cantera as ct
import numpy as np
import matplotlib.pyplot as plt
defcompute_ignition_delay(fuel, phi=1.0, T_range=None, P=20*101325,
mechanism="gri30.yaml"):
"""
Compute ignition delay time via constant-volume homogeneous reactor.
Parameters
----------
T_range : array-like
Initial temperatures (K) to sweep (Arrhenius plot)
Returns
-------
dict: T_K → tau_ign_ms
"""if T_range isNone:
T_range = np.arange(1000, 1800, 50)
gas = ct.Solution(mechanism)
results = []
for T0 in T_range:
gas.set_equivalence_ratio(phi, fuel, "air")
gas.TP = T0, P
reactor = ct.IdealGasReactor(gas)
net = ct.ReactorNet([reactor])
# Integrate until OH peak (ignition criterion)
t_end = 0.01# 10 ms max
dt = 1e-6
t, OH_prev = 0.0, gas["OH"].X[0]
t_ign = Nonewhile t < t_end:
t = net.advance(t + dt)
OH_cur = reactor.thermo["OH"].X[0]
if OH_cur > 10 * OH_prev and OH_cur > 1e-4:
t_ign = t
break
OH_prev = max(OH_cur, OH_prev)
tau = (t_ign * 1000) if t_ign else np.nan
results.append({"T_K": T0, "tau_ign_ms": tau})
import pandas as pd
return pd.DataFrame(results).dropna()
# Methane ignition delay
tau_df = compute_ignition_delay("CH4", phi=1.0,
T_range=np.arange(1000, 1700, 50),
P=20*101325)
fig, ax = plt.subplots(figsize=(8, 5))
ax.semilogy(1000 / tau_df["T_K"], tau_df["tau_ign_ms"], 'bs-', ms=6)
ax.set_xlabel("1000/T (K⁻¹)")
ax.set_ylabel("Ignition delay τ (ms)")
ax.set_title("Methane Ignition Delay Time (φ=1.0, P=20 atm)")
ax.grid(True, which="both", alpha=0.3)
plt.tight_layout()
plt.savefig("ignition_delay.png", dpi=150)
plt.show()
print("Ignition delays:")
print(tau_df.to_string(index=False))
Step 3: Laminar Flame Speed
import cantera as ct
import numpy as np
import matplotlib.pyplot as plt
deflaminar_flame_speed(fuel, phi_range=None, T0=300, P=101325, mechanism="gri30.yaml"):
"""
Compute laminar burning velocity using 1D freely propagating flame.
Returns
-------
pd.DataFrame: phi, Su_cm_s (laminar flame speed in cm/s)
"""if phi_range isNone:
phi_range = np.arange(0.7, 1.41, 0.1)
results = []
for phi in phi_range:
gas = ct.Solution(mechanism)
gas.set_equivalence_ratio(phi, fuel, "air")
gas.TP = T0, P
flame = ct.FreeFlame(gas, width=0.03) # 3 cm domain
flame.set_refine_criteria(ratio=3, slope=0.06, curve=0.12)
flame.set_max_jac_age(50, 50)
flame.solve(loglevel=0, auto=True)
Su = flame.velocity[0] * 100# m/s → cm/s
results.append({"phi": phi, "Su_cm_s": Su})
print(f" φ={phi:.1f}: Su={Su:.1f} cm/s")
import pandas as pd
return pd.DataFrame(results)
print("Computing methane laminar flame speeds...")
flame_df = laminar_flame_speed("CH4", phi_range=[0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3])
fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(flame_df["phi"], flame_df["Su_cm_s"], 'ro-', ms=8, linewidth=2)
ax.set_xlabel("Equivalence ratio φ")
ax.set_ylabel("Laminar burning velocity Su (cm/s)")
ax.set_title("Methane/Air Laminar Flame Speed")
ax.set_ylim(0)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("laminar_flame_speed.png", dpi=150)
plt.show()
peak_phi = flame_df.loc[flame_df["Su_cm_s"].idxmax(), "phi"]
peak_Su = flame_df["Su_cm_s"].max()
print(f"\nPeak Su = {peak_Su:.1f} cm/s at φ = {peak_phi:.1f}")
Advanced Usage
Perfectly Stirred Reactor (PSR) and Sensitivity Analysis
import cantera as ct
import numpy as np
import matplotlib.pyplot as plt
defpsr_residence_time_sweep(fuel, phi=1.0, T0=800, P=101325,
tau_range=None, mechanism="gri30.yaml"):
"""
Sweep residence times in a PSR to find extinction.
"""if tau_range isNone:
tau_range = np.logspace(-4, 0, 30) # 0.1 ms to 1 s
gas = ct.Solution(mechanism)
results = []
for tau in tau_range:
gas.set_equivalence_ratio(phi, fuel, "air")
gas.TP = T0, P
reactor = ct.IdealGasReactor(gas)
upstream = ct.Reservoir(gas)
exhaust = ct.Reservoir(gas)
intake_mfc = ct.MassFlowController(upstream, reactor)
intake_mfc.mass_flow_rate = reactor.mass / tau
ct.PressureController(reactor, exhaust, master=intake_mfc)
net = ct.ReactorNet([reactor])
try:
net.advance_to_steady_state()
results.append({"tau_s": tau, "T_K": reactor.T, "X_CO2": reactor.thermo["CO2"].X[0]})
except Exception:
breakimport pandas as pd
df = pd.DataFrame(results)
extinction_tau = df[df["T_K"] < T0 + 100]["tau_s"].min() ifnot df.empty elseNoneprint(f"PSR extinction residence time: {extinction_tau:.4f} s"if extinction_tau else"No extinction")
return df
psr_df = psr_residence_time_sweep("CH4")
fig, ax = plt.subplots(figsize=(8, 5))
ax.semilogx(psr_df["tau_s"] * 1000, psr_df["T_K"], 'b-o', ms=5)
ax.set_xlabel("Residence time τ (ms)")
ax.set_ylabel("Reactor temperature T (K)")
ax.set_title("PSR S-curve: CH4/Air, φ=1.0")
ax.grid(True, which="both", alpha=0.3)
plt.tight_layout()
plt.savefig("psr_s_curve.png", dpi=150)
plt.show()
Troubleshooting
Error: cantera._cantera.CanteraError: mechanism file not found
Fix:
import cantera as ct
# List available built-in mechanismsimport os
print([f for f in os.listdir(ct.get_data_path()) if f.endswith(".yaml")])
# Common: gri30.yaml, h2o2.yaml, air.yaml, nDodecane_Reitz.yaml
Issue: Flame solver doesn't converge
Fix:
# Use wider domain and coarser initial grid
flame = ct.FreeFlame(gas, width=0.05)
flame.set_refine_criteria(ratio=5, slope=0.2, curve=0.3)
flame.solve(loglevel=0, auto=True)
Goodwin, D.G. et al. (2023). Cantera: An object-oriented software toolkit for chemical kinetics, thermodynamics, and transport processes. Zenodo.
Examples
Example 1: NOx Formation vs. Equivalence Ratio
# =============================================# NOx equilibrium concentration vs. phi# =============================================import cantera as ct, numpy as np, matplotlib.pyplot as plt
gas = ct.Solution("gri30.yaml")
phi_arr = np.linspace(0.6, 1.4, 30)
NOx_ppm = []
for phi in phi_arr:
gas.set_equivalence_ratio(phi, "CH4", "air")
gas.TP = 300, 101325
gas.equilibrate("HP")
NO = gas["NO"].X[0]
NO2 = gas["NO2"].X[0]
NOx_ppm.append((NO + NO2) * 1e6)
fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(phi_arr, NOx_ppm, 'r-o', ms=5)
ax.set_xlabel("Equivalence ratio φ")
ax.set_ylabel("NOx (ppm)")
ax.set_title("Equilibrium NOx vs. Equivalence Ratio (CH4/air)")
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("nox_vs_phi.png", dpi=150)
plt.show()
print(f"Peak NOx: {max(NOx_ppm):.0f} ppm at φ = {phi_arr[np.argmax(NOx_ppm)]:.2f}")
Interpreting these results: NOx peaks near stoichiometric conditions due to the high adiabatic flame temperature. Lean or rich operation reduces NOx — the basis of lean-premix combustors.
Last updated: 2026-03-17 | Maintainer: @xjtulycIssues: GitHub Issues
