# Flame Speed with Sensitivity Analysis¶

In this example we simulate a freely-propagating, adiabatic, 1-D flame and

• Calculate its laminar burning velocity
• Perform a sensitivity analysis of its kinetics

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}$.

### Import Modules¶

In [1]:
from __future__ import print_function
from __future__ import division

import cantera as ct
import numpy as np

print("Running Cantera Version: " + str(ct.__version__))

Running Cantera Version: 2.5.0a2

In [2]:
# Import plotting modules and define plotting preference
%matplotlib notebook
import matplotlib.pylab as plt

plt.rcParams['axes.labelsize'] = 14
plt.rcParams['xtick.labelsize'] = 12
plt.rcParams['ytick.labelsize'] = 12
plt.rcParams['legend.fontsize'] = 10
plt.rcParams['figure.figsize'] = (8,6)

plt.style.use('ggplot')
plt.style.use('seaborn-deep')

plt.rcParams['figure.autolayout'] = True

# Import Pandas for DataFrames
import pandas as pd


### Define the reactant conditions, gas mixture and kinetic mechanism associated with the gas¶

In [3]:
#Inlet Temperature in Kelvin and Inlet Pressure in Pascals
#In this case we are setting the inlet T and P to room temperature conditions
To = 300
Po = 101325

#Define the gas-mixutre and kinetics
#In this case, we are choosing a GRI3.0 gas
gas = ct.Solution('gri30.cti')

# Create a stoichiometric CH4/Air premixed mixture
gas.set_equivalence_ratio(1.0, 'CH4', {'O2':1.0, 'N2':3.76})
gas.TP = To, Po


### Define flame simulation conditions¶

In [4]:
# Domain width in metres
width = 0.014

# Create the flame object
flame = ct.FreeFlame(gas, width=width)

# Define tolerances for the solver
flame.set_refine_criteria(ratio=3, slope=0.1, curve=0.1)

# Define logging level
loglevel = 1


### Solve¶

In [5]:
flame.solve(loglevel=loglevel, auto=True)
Su0 = flame.u[0]
print("Flame Speed is: {:.2f} cm/s".format(Su0*100))

# Note that the variable Su0 will also be used downsteam in the sensitivity analysis

************ Solving on 8 point grid with energy equation enabled ************

..............................................................................
Attempt Newton solution of steady-state problem...    failure.
Take 10 timesteps     2.136e-05      5.452
Attempt Newton solution of steady-state problem...    failure.
Take 10 timesteps     0.0003649      4.427
Attempt Newton solution of steady-state problem...    failure.
Take 10 timesteps     2.435e-05      6.061
Attempt Newton solution of steady-state problem...    failure.
Take 10 timesteps     3.468e-05       5.63
Attempt Newton solution of steady-state problem...    failure.
Take 10 timesteps      0.001333      4.122
Attempt Newton solution of steady-state problem...    success.

Problem solved on [9] point grid(s).
Expanding domain to accomodate flame thickness. New width: 0.028 m
##############################################################################
Refining grid in flame.
New points inserted after grid points 0 1 2 3 4 5 6
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 u
##############################################################################

*********** Solving on 16 point grid with energy equation enabled ************

..............................................................................
Attempt Newton solution of steady-state problem...    failure.
Take 10 timesteps     2.136e-05      5.751
Attempt Newton solution of steady-state problem...    failure.
Take 10 timesteps     4.055e-05       5.58
Attempt Newton solution of steady-state problem...    failure.
Take 10 timesteps     2.887e-05      5.947
Attempt Newton solution of steady-state problem...    failure.
Take 10 timesteps     2.055e-05       6.11
Attempt Newton solution of steady-state problem...    failure.
Take 10 timesteps     0.0001756      5.504
Attempt Newton solution of steady-state problem...    failure.
Take 10 timesteps     1.172e-05      6.764
Attempt Newton solution of steady-state problem...    success.

Problem solved on [16] point grid(s).
Expanding domain to accomodate flame thickness. New width: 0.056 m
##############################################################################
Refining grid in flame.
New points inserted after grid points 3 4 5 6 7 8 9
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 u
##############################################################################

*********** Solving on 23 point grid with energy equation enabled ************

..............................................................................
Attempt Newton solution of steady-state problem...    failure.
Take 10 timesteps     1.424e-05      6.312
Attempt Newton solution of steady-state problem...    failure.
Take 10 timesteps     3.604e-05      5.689
Attempt Newton solution of steady-state problem...    failure.
Take 10 timesteps     1.283e-05      6.137
Attempt Newton solution of steady-state problem...    failure.
Take 10 timesteps     0.0001096      5.738
Attempt Newton solution of steady-state problem...    failure.
Take 10 timesteps     2.601e-05       6.05
Attempt Newton solution of steady-state problem...    failure.
Take 10 timesteps     6.943e-06      6.869
Attempt Newton solution of steady-state problem...    failure.
Take 10 timesteps     7.909e-05      5.845
Attempt Newton solution of steady-state problem...    failure.
Take 10 timesteps     6.256e-06      6.571
Attempt Newton solution of steady-state problem...    failure.
Take 10 timesteps     0.0001069      5.566
Attempt Newton solution of steady-state problem...    failure.
Take 10 timesteps      7.61e-05      5.798
Attempt Newton solution of steady-state problem...    failure.
Take 10 timesteps     2.709e-05      5.646
Attempt Newton solution of steady-state problem...    failure.
Take 10 timesteps     0.0006942      4.499
Attempt Newton solution of steady-state problem...    success.

Problem solved on [23] point grid(s).

..............................................................................
grid refinement disabled.

******************** Solving with grid refinement enabled ********************

..............................................................................
Attempt Newton solution of steady-state problem...    success.

Problem solved on [23] point grid(s).

..............................................................................
##############################################################################
Refining grid in flame.
New points inserted after grid points 6 7 8 9 10 11 12 13
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 u
##############################################################################

..............................................................................
Attempt Newton solution of steady-state problem...    failure.
Take 10 timesteps     0.0001709      5.006
Attempt Newton solution of steady-state problem...    success.

Problem solved on [31] point grid(s).

..............................................................................
##############################################################################
Refining grid in flame.
New points inserted after grid points 8 9 10 11 12 13 14 15 16 27 28
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 u
##############################################################################

..............................................................................
Attempt Newton solution of steady-state problem...    failure.
Take 10 timesteps     7.594e-05      5.228
Attempt Newton solution of steady-state problem...    success.

Problem solved on [42] point grid(s).

..............................................................................
##############################################################################
Refining grid in flame.
New points inserted after grid points 12 13 14 15 16 17 18 19 20 21 40
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 u
##############################################################################

..............................................................................
Attempt Newton solution of steady-state problem...    success.

Problem solved on [53] point grid(s).

..............................................................................
##############################################################################
Refining grid in flame.
New points inserted after grid points 15 16 17 18 19 20 21 22 23 24 25 26 27 50
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 u
##############################################################################

..............................................................................
Attempt Newton solution of steady-state problem...    success.

Problem solved on [67] point grid(s).

..............................................................................
##############################################################################
Refining grid in flame.
New points inserted after grid points 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36
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 u
##############################################################################

..............................................................................
Attempt Newton solution of steady-state problem...    success.

Problem solved on [85] point grid(s).

..............................................................................
##############################################################################
Refining grid in flame.
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
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 u
##############################################################################

..............................................................................
Attempt Newton solution of steady-state problem...    success.

Problem solved on [113] point grid(s).

..............................................................................
##############################################################################
Refining grid in flame.
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
to resolve C C2H2 C2H3 C2H4 C2H5 C2H6 C3H8 CH CH2 CH2(S) CH2CO CH2OH CH3 CH3CHO CH3O HCCO HCO
##############################################################################

..............................................................................
Attempt Newton solution of steady-state problem...    success.

Problem solved on [142] point grid(s).

..............................................................................
no new points needed in flame
Flame Speed is: 38.22 cm/s


### Plot figures¶

Check and see if all has gone well. Plot temperature and species fractions to see

#### Temperature Plot¶

In [6]:
plt.figure()

plt.plot(flame.grid*100, flame.T, '-o')
plt.xlabel('Distance (cm)')
plt.ylabel('Temperature (K)');


#### Major species' plot¶

To plot species, we first have to identify the index of the species in the array For this, cut & paste the following lines and run in a new cell to get the index

for i, specie in enumerate(gas.species()):
print(str(i) + '. ' + str(specie))
In [7]:
# Extract concentration data
X_CH4 = flame.X[13]
X_CO2 = flame.X[15]
X_H2O = flame.X[5]

plt.figure()

plt.plot(flame.grid*100, X_CH4, '-o', label=r'$CH_{4}$')
plt.plot(flame.grid*100, X_CO2, '-s', label=r'$CO_{2}$')
plt.plot(flame.grid*100, X_H2O, '-<', label=r'$H_{2}O$')

plt.legend(loc=2)
plt.xlabel('Distance (cm)')
plt.ylabel('MoleFractions');


## Sensitivity Analysis¶

See which reactions effect the flame speed the most

In [8]:
# Create a dataframe to store sensitivity-analysis data
sensitivities = pd.DataFrame(data=[], index=gas.reaction_equations(range(gas.n_reactions)))


### Compute sensitivities¶

In [9]:
# Set the value of the perturbation
dk = 1e-2

# Create an empty column to store the sensitivities data
sensitivities["baseCase"] = ""

In [10]:
for m in range(gas.n_reactions):
gas.set_multiplier(1.0) # reset all multipliers
gas.set_multiplier(1+dk, m) # perturb reaction m

# Always force loglevel=0 for this
# Make sure the grid is not refined, otherwise it won't strictly
# be a small perturbation analysis
flame.solve(loglevel=0, refine_grid=False)

# The new flame speed
Su = flame.u[0]

sensitivities["baseCase"][m] = (Su-Su0)/(Su0*dk)

# This step is essential, otherwise the mechanism will have been altered
gas.set_multiplier(1.0)

In [11]:
sensitivities.head()

Out[11]:
baseCase
2 O + M <=> O2 + M 0.00156775
H + O + M <=> OH + M 0.00110824
H2 + O <=> H + OH 0.0252473
HO2 + O <=> O2 + OH 0.00313172
H2O2 + O <=> HO2 + OH 0.000752237

### Make plots¶

In [12]:
# Reaction mechanisms can contains thousands of elementary steps. Choose a threshold
# to see only the top few
threshold = 0.03

firstColumn = sensitivities.columns[0]

# For plotting, collect only those steps that are above the threshold
# Otherwise, the y-axis gets crowded and illegible
sensitivitiesSubset = sensitivities[sensitivities[firstColumn].abs() > threshold]
indicesMeetingThreshold = sensitivitiesSubset[firstColumn].abs().sort_values(ascending=False).index
sensitivitiesSubset.loc[indicesMeetingThreshold].plot.barh(title="Sensitivities for GRI 3.0",
legend=None)
plt.gca().invert_yaxis()

plt.rcParams.update({'axes.labelsize': 20})
plt.xlabel(r'Sensitivity: $\frac{\partial\:\ln{S_{u}}}{\partial\:\ln{k}}$');

# Uncomment the following to save the plot. A higher than usual resolution (dpi) helps
# plt.savefig('sensitivityPlot', dpi=300)