{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Flame Speed with Sensitivity Analysis\n",
"In this example we simulate a freely-propagating, adiabatic, 1-D flame, calculate its laminar burning velocity and perform a sensitivity analysis of its kinetics.\n",
"\n",
"The figure below illustrates the setup, in a flame-fixed co-ordinate system. The reactants enter with density $\\rho_{u}$, temperature $T_{u}$ and speed $S_{u}$. The products exit the flame at speed $S_{b}$, density $\\rho_{b}$ and temperature $T_{b}$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Import Modules"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Running Cantera Version: 3.0.0\n"
]
}
],
"source": [
"import cantera as ct\n",
"import numpy as np\n",
"import pandas as pd\n",
"\n",
"print(f\"Running Cantera Version: {ct.__version__}\")"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"# Import plotting modules and define plotting preference\n",
"%config InlineBackend.figure_formats = [\"svg\"]\n",
"%matplotlib inline\n",
"import matplotlib.pylab as plt\n",
"\n",
"plt.rcParams[\"axes.labelsize\"] = 14\n",
"plt.rcParams[\"xtick.labelsize\"] = 12\n",
"plt.rcParams[\"ytick.labelsize\"] = 12\n",
"plt.rcParams[\"legend.fontsize\"] = 10\n",
"plt.rcParams[\"figure.figsize\"] = (8, 6)\n",
"plt.rcParams[\"figure.dpi\"] = 120\n",
"\n",
"# Get the best of both ggplot and seaborn\n",
"plt.style.use(\"ggplot\")\n",
"plt.style.use(\"seaborn-v0_8-deep\")\n",
"\n",
"plt.rcParams[\"figure.autolayout\"] = True"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Define the reactant conditions, gas mixture and kinetic mechanism associated with the gas"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"# Inlet Temperature in Kelvin and Inlet Pressure in Pascals\n",
"# In this case we are setting the inlet T and P to room temperature conditions\n",
"To = 300\n",
"Po = 101325\n",
"\n",
"# Define the gas-mixutre and kinetics\n",
"# In this case, we are choosing a GRI3.0 gas\n",
"gas = ct.Solution(\"gri30.yaml\")\n",
"\n",
"# Create a stoichiometric CH4/Air premixed mixture\n",
"gas.set_equivalence_ratio(1.0, \"CH4\", {\"O2\": 1.0, \"N2\": 3.76})\n",
"gas.TP = To, Po"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Define flame simulation conditions"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"# Domain width in metres\n",
"width = 0.014\n",
"\n",
"# Create the flame object\n",
"flame = ct.FreeFlame(gas, width=width)\n",
"\n",
"# Define tolerances for the solver\n",
"flame.set_refine_criteria(ratio=3, slope=0.1, curve=0.1)\n",
"\n",
"# Define logging level\n",
"loglevel = 1"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Solve"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"scrolled": true,
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"************ Solving on 8 point grid with energy equation enabled ************\n",
"\n",
"..............................................................................\n",
"Attempt Newton solution of steady-state problem... failure. \n",
"Take 10 timesteps 2.136e-05 5.455\n",
"Attempt Newton solution of steady-state problem... failure. \n",
"Take 10 timesteps 0.0003649 4.425\n",
"Attempt Newton solution of steady-state problem... failure. \n",
"Take 10 timesteps 3.653e-05 5.86\n",
"Attempt Newton solution of steady-state problem... failure. \n",
"Take 10 timesteps 1.734e-05 6.005\n",
"Attempt Newton solution of steady-state problem... failure. \n",
"Take 10 timesteps 0.0006666 4.488\n",
"Attempt Newton solution of steady-state problem... success.\n",
"\n",
"Problem solved on [9] point grid(s).\n",
"Expanding domain to accommodate flame thickness. New width: 0.028 m\n",
"##############################################################################\n",
"Refining grid in flame.\n",
" New points inserted after grid points 0 1 2 3 4 5 6 \n",
" to resolve C C2H2 C2H3 C2H4 C2H5 C2H6 C3H7 C3H8 CH CH2 CH2(S) CH2CHO CH2CO CH2O CH2OH CH3 CH3CHO CH3O CH3OH CH4 CO CO2 H H2 H2O H2O2 HCCO HCCOH HCN HCNO HCO HNCO HO2 N N2 N2O NCO NH NO NO2 O O2 OH T velocity \n",
"##############################################################################\n",
"\n",
"*********** Solving on 16 point grid with energy equation enabled ************\n",
"\n",
"..............................................................................\n",
"Attempt Newton solution of steady-state problem... failure. \n",
"Take 10 timesteps 2.136e-05 5.75\n",
"Attempt Newton solution of steady-state problem... failure. \n",
"Take 10 timesteps 4.055e-05 5.579\n",
"Attempt Newton solution of steady-state problem... failure. \n",
"Take 10 timesteps 5.773e-05 6.027\n",
"Attempt Newton solution of steady-state problem... failure. \n",
"Take 10 timesteps 2.74e-05 5.924\n",
"Attempt Newton solution of steady-state problem... failure. \n",
"Take 10 timesteps 7.802e-05 5.655\n",
"Attempt Newton solution of steady-state problem... failure. \n",
"Take 10 timesteps 1.852e-05 5.771\n",
"Attempt Newton solution of steady-state problem... success.\n",
"\n",
"Problem solved on [16] point grid(s).\n",
"Expanding domain to accommodate flame thickness. New width: 0.056 m\n",
"##############################################################################\n",
"Refining grid in flame.\n",
" New points inserted after grid points 3 4 5 6 7 8 9 \n",
" to resolve C C2H2 C2H3 C2H4 C2H5 C2H6 C3H7 C3H8 CH CH2 CH2(S) CH2CHO CH2CO CH2O CH2OH CH3 CH3CHO CH3O CH3OH CH4 CO CO2 H H2 H2O H2O2 HCCO HCCOH HCN HCNO HCO HNCO HO2 N N2 N2O NCO NH NH2 NO NO2 O O2 OH T velocity \n",
"##############################################################################\n",
"\n",
"*********** Solving on 23 point grid with energy equation enabled ************\n",
"\n",
"..............................................................................\n",
"Attempt Newton solution of steady-state problem... failure. \n",
"Take 10 timesteps 1.424e-05 6.312\n",
"Attempt Newton solution of steady-state problem... failure. \n",
"Take 10 timesteps 3.604e-05 5.689\n",
"Attempt Newton solution of steady-state problem... failure. \n",
"Take 10 timesteps 9.622e-06 6.227\n",
"Attempt Newton solution of steady-state problem... failure. \n",
"Take 10 timesteps 0.0001644 5.518\n",
"Attempt Newton solution of steady-state problem... failure. \n",
"Take 10 timesteps 1.3e-05 6.218\n",
"Attempt Newton solution of steady-state problem... failure. \n",
"Take 10 timesteps 2.083e-05 6.329\n",
"Attempt Newton solution of steady-state problem... failure. \n",
"Take 10 timesteps 5.932e-05 5.87\n",
"Attempt Newton solution of steady-state problem... failure. \n",
"Take 10 timesteps 1.056e-05 6.871\n",
"Attempt Newton solution of steady-state problem... failure. \n",
"Take 10 timesteps 0.0001203 5.583\n",
"Attempt Newton solution of steady-state problem... failure. \n",
"Take 10 timesteps 8.561e-05 5.809\n",
"Attempt Newton solution of steady-state problem... failure. \n",
"Take 10 timesteps 6.095e-05 5.255\n",
"Attempt Newton solution of steady-state problem... failure. \n",
"Take 10 timesteps 0.001562 4.092\n",
"Attempt Newton solution of steady-state problem... success.\n",
"\n",
"Problem solved on [23] point grid(s).\n",
"\n",
"..............................................................................\n",
"grid refinement disabled.\n",
"\n",
"******************** Solving with grid refinement enabled ********************\n",
"\n",
"..............................................................................\n",
"Attempt Newton solution of steady-state problem... success.\n",
"\n",
"Problem solved on [23] point grid(s).\n",
"\n",
"..............................................................................\n",
"##############################################################################\n",
"Refining grid in flame.\n",
" New points inserted after grid points 6 7 8 9 10 11 12 13 \n",
" to resolve C C2H2 C2H3 C2H4 C2H5 C2H6 C3H7 C3H8 CH CH2 CH2(S) CH2CHO CH2CO CH2O CH2OH CH3 CH3CHO CH3O CH3OH CH4 CO CO2 H H2 H2O H2O2 HCCO HCCOH HCN HCNO HCO HNCO HO2 N N2 N2O NCO NH NH2 NO NO2 O O2 OH T velocity \n",
"##############################################################################\n",
"\n",
"..............................................................................\n",
"Attempt Newton solution of steady-state problem... failure. \n",
"Take 10 timesteps 0.0001709 5.009\n",
"Attempt Newton solution of steady-state problem... success.\n",
"\n",
"Problem solved on [31] point grid(s).\n",
"\n",
"..............................................................................\n",
"##############################################################################\n",
"Refining grid in flame.\n",
" New points inserted after grid points 8 9 10 11 12 13 14 15 16 27 28 \n",
" to resolve C C2H C2H2 C2H3 C2H4 C2H5 C2H6 C3H7 C3H8 CH CH2 CH2(S) CH2CHO CH2CO CH2O CH2OH CH3 CH3CHO CH3O CH3OH CH4 CO CO2 H H2 H2O H2O2 HCCO HCCOH HCN HCNO HCO HNCO HO2 N N2 N2O NCO NH NH2 NH3 NO NO2 O O2 OH T velocity \n",
"##############################################################################\n",
"\n",
"..............................................................................\n",
"Attempt Newton solution of steady-state problem... failure. \n",
"Take 10 timesteps 7.594e-05 5.214\n",
"Attempt Newton solution of steady-state problem... success.\n",
"\n",
"Problem solved on [42] point grid(s).\n",
"\n",
"..............................................................................\n",
"##############################################################################\n",
"Refining grid in flame.\n",
" New points inserted after grid points 12 13 14 15 16 17 18 19 20 21 40 \n",
" to resolve C C2H C2H2 C2H3 C2H4 C2H5 C2H6 C3H7 C3H8 CH CH2 CH2(S) CH2CHO CH2CO CH2O CH2OH CH3 CH3CHO CH3O CH3OH CH4 CO CO2 H H2 H2O H2O2 HCCO HCCOH HCN HCNO HCO HNCO HO2 N N2 N2O NCO NH NH2 NH3 NO NO2 O O2 OH T point 40 velocity \n",
"##############################################################################\n",
"\n",
"..............................................................................\n",
"Attempt Newton solution of steady-state problem... success.\n",
"\n",
"Problem solved on [53] point grid(s).\n",
"\n",
"..............................................................................\n",
"##############################################################################\n",
"Refining grid in flame.\n",
" New points inserted after grid points 15 16 17 18 19 20 21 22 23 24 25 26 27 50 \n",
" to resolve C C2H C2H2 C2H3 C2H4 C2H5 C2H6 C3H7 C3H8 CH CH2 CH2(S) CH2CHO CH2CO CH2O CH2OH CH3 CH3CHO CH3O CH3OH CH4 CO CO2 H H2 H2O H2O2 HCCO HCCOH HCN HCNO HCO HNCO HO2 N N2 N2O NCO NO NO2 O O2 OH T velocity \n",
"##############################################################################\n",
"\n",
"..............................................................................\n",
"Attempt Newton solution of steady-state problem... success.\n",
"\n",
"Problem solved on [67] point grid(s).\n",
"\n",
"..............................................................................\n",
"##############################################################################\n",
"Refining grid in flame.\n",
" New points inserted after grid points 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 \n",
" to resolve C C2H C2H2 C2H3 C2H4 C2H5 C2H6 C3H7 C3H8 CH CH2 CH2(S) CH2CHO CH2CO CH2O CH2OH CH3 CH3CHO CH3O CH3OH CH4 CO CO2 H H2 H2O H2O2 HCCO HCCOH HCN HCNO HCO HNCO HO2 N N2 NO NO2 O O2 OH T velocity \n",
"##############################################################################\n",
"\n",
"..............................................................................\n",
"Attempt Newton solution of steady-state problem... success.\n",
"\n",
"Problem solved on [85] point grid(s).\n",
"\n",
"..............................................................................\n",
"##############################################################################\n",
"Refining grid in flame.\n",
" New points inserted after grid points 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 \n",
" to resolve C C2H C2H2 C2H3 C2H4 C2H5 C2H6 C3H7 C3H8 CH CH2 CH2(S) CH2CHO CH2CO CH2O CH2OH CH3 CH3CHO CH3O CH3OH CH4 CO H H2 H2O H2O2 HCCO HCCOH HCN HCO HNCO HO2 N2 NO2 O O2 OH T velocity \n",
"##############################################################################\n",
"\n",
"..............................................................................\n",
"Attempt Newton solution of steady-state problem... success.\n",
"\n",
"Problem solved on [113] point grid(s).\n",
"\n",
"..............................................................................\n",
"##############################################################################\n",
"Refining grid in flame.\n",
" New points inserted after grid points 34 35 36 37 38 39 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 \n",
" to resolve C C2H2 C2H3 C2H4 C2H5 C2H6 C3H8 CH CH2 CH2(S) CH2CO CH2OH CH3 CH3CHO CH3O HCCO HCO \n",
"##############################################################################\n",
"\n",
"..............................................................................\n",
"Attempt Newton solution of steady-state problem... success.\n",
"\n",
"Problem solved on [142] point grid(s).\n",
"\n",
"..............................................................................\n",
"no new points needed in flame\n",
"Flame Speed is: 38.25 cm/s\n"
]
}
],
"source": [
"flame.solve(loglevel=loglevel, auto=True)\n",
"Su0 = flame.velocity[0]\n",
"print(f\"Flame Speed is: {Su0 * 100:.2f} cm/s\")\n",
"\n",
"# Note that the variable Su0 will also be used downstream in the sensitivity analysis"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Plot figures\n",
"\n",
"Check and see if all has gone well. Plot temperature and species fractions to see\n",
"\n",
"#### Temperature Plot"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.figure()\n",
"\n",
"plt.plot(flame.grid * 100, flame.T, \"-o\")\n",
"plt.xlabel(\"Distance (cm)\")\n",
"plt.ylabel(\"Temperature (K)\");"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Major species' plot"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To plot species, we first have to identify the index of the species in the array.\n",
"For this, cut & paste the following lines and run in a new cell to get the index\n",
"\n",
"```python\n",
"for i, specie in enumerate(gas.species()):\n",
" print(f\"{i}. {specie}\")\n",
"```"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Extract concentration data\n",
"X_CH4 = flame.X[13]\n",
"X_CO2 = flame.X[15]\n",
"X_H2O = flame.X[5]\n",
"\n",
"plt.figure()\n",
"\n",
"plt.plot(flame.grid * 100, X_CH4, \"-o\", label=\"$CH_{4}$\")\n",
"plt.plot(flame.grid * 100, X_CO2, \"-s\", label=\"$CO_{2}$\")\n",
"plt.plot(flame.grid * 100, X_H2O, \"-<\", label=\"$H_{2}O$\")\n",
"\n",
"plt.legend(loc=2)\n",
"plt.xlabel(\"Distance (cm)\")\n",
"plt.ylabel(\"MoleFractions\");"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Sensitivity Analysis\n",
"\n",
"See which reactions effect the flame speed the most"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"# Create a dataframe to store sensitivity-analysis data\n",
"sensitivities = pd.DataFrame(index=gas.reaction_equations(), columns=[\"base_case\"])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Compute sensitivities"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"# Set the value of the perturbation\n",
"dk = 1e-2"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"for m in range(gas.n_reactions):\n",
" gas.set_multiplier(1.0) # reset all multipliers\n",
" gas.set_multiplier(1 + dk, m) # perturb reaction m\n",
"\n",
" # Always force loglevel=0 for this\n",
" # Make sure the grid is not refined, otherwise it won't strictly\n",
" # be a small perturbation analysis\n",
" # Turn auto-mode off since the flame has already been solved\n",
" flame.solve(loglevel=0, refine_grid=False, auto=False)\n",
"\n",
" # The new flame speed\n",
" Su = flame.velocity[0]\n",
"\n",
" sensitivities.iloc[m, 0] = (Su - Su0) / (Su0 * dk)\n",
"\n",
"# This step is essential, otherwise the mechanism will have been altered\n",
"gas.set_multiplier(1.0)"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
base_case
\n",
"
\n",
" \n",
" \n",
"
\n",
"
2 O + M <=> O2 + M
\n",
"
0.001548
\n",
"
\n",
"
\n",
"
H + O + M <=> OH + M
\n",
"
0.00109
\n",
"
\n",
"
\n",
"
H2 + O <=> H + OH
\n",
"
0.025247
\n",
"
\n",
"
\n",
"
HO2 + O <=> O2 + OH
\n",
"
0.003123
\n",
"
\n",
"
\n",
"
H2O2 + O <=> HO2 + OH
\n",
"
0.000735
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" base_case\n",
"2 O + M <=> O2 + M 0.001548\n",
"H + O + M <=> OH + M 0.00109\n",
"H2 + O <=> H + OH 0.025247\n",
"HO2 + O <=> O2 + OH 0.003123\n",
"H2O2 + O <=> HO2 + OH 0.000735"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sensitivities.head()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Make plots"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
"