Sensitivity analysis of Methane combustion using Python and Cantera
Aim:-> To compute sensitivities of various reactions in GRI 3.0 mechanism for methane combustion.
Theory:-
As we know any combustion reaction consists of number of intermediate chemical reactions and chemical species before complete combustion takes place which affects the entire reaction process. The bunch of such reaction can be called as a reaction mechanism. Since there are many parameters involved in the chemistry of combustion, it is needed to understand the effects of each parameter on the system. One of the best technique is to analyze it by performing a sensitivity check for given chemistry mechanism.
Sensitivity analysis is nothing but getting a relationship between the sensitivities of reaction and the parameters like temperature or in mathematical term it is the study of change of dependent variable with respect to the independent variable i.e. slope of the curve. The main moto behind this sensitivity analysis process is to estimate the sensitivity for parameter which indirectly tells about ignition delay in the combustion process for intermediate reactions. In analysis the reaction with more sensitivity results less ignition delay and by analyzing behavior of intermediate reactions we can regulate the mechanism reactions to tune the ignition delay.
As the combustion process is extremely complicated to model and also requires high computation cost and time to simulate complete process accurately. In a combustion reaction there are almost 325 intermediate reactions and 53 species as per GRI 3.0 mechanism. So it will became tremendously hectic and expensive to simulate for each and every reaction. Hence to reduce the computational cost only few reactions are taken in to account while studying. To find these reactions we carry out sensitivity analysis. The reactions which are highly affect parameters has maximum sensitivity and only those reactions are considered in study.
Coding approach is provided in following steps:
- Create “gas" as an instance of the “Solution" class based on the mechanism from the “gri30.xml" file.
- Consider initial conditions as-
- Initial temperature = T = 1500k
- Initial temperature = P = 1atm
- species concentrations as per methane’s stochiometric combustion reaction
- Methane's combustion reaction is `CH_4+2(O_2 + 3.76N_2)<=>CO_2 + 2H_2O + 7.52N_2`
- Provide initial condition to reactor (TPX) using “gas.TPX” command.
- Create an ideal gas constant pressure reactor "r" from the "gas" object.
- Then create reactor network "sim" from the "r" instance.
- Consider the number of reaction to take in account from GRI 3.0 mechanism using “gas.n_reactions” command.
- Extract the reactions from "r" and add them to “reac_n” array for the sensitivity analysis. [Using for loop]
- Set the relative and absolute tolerance of the reactor solver to 10−6 and 10-15
- Set the relative and absolute tolerance of the sensitivity solver to 10−6 and 10-6
- Create a solution array called as “states” to store the thermodynamic state of "gas" at every time-step.
- Create time domain for simulation. The time domain is 0 to 2 milli-seconds and with 5 micro-second timestep size.
- Then advance time to get results at various time steps i.e. at various states using “sim.advance(time)” command. [Using for loop of time]
- Again use ‘for loop’ inside ‘time advancement loop’ to obtained sensitivities for each reactions at that time step using “sim.sensitivity('temperature',i)” command. [Using for loop of space]
- Use ‘if loop’ to get maximum sensitivity values inside ‘space – for loop’
- Append states object to store values.
- Take absolute values of sensitivities to compare most sensitive reactions.
- Use “sorted” and “zip” function appropriately to sort and find out most sensitive reactions with their sensitivities.
- Print required values like sensitivity and reactions.
- And then plot all the results.
 |
Python code for sensitivity analysis:-
'''
Code to compute sensitivities of reactions in GRI 3 mechanisms
'''
import cantera as ct
import numpy as np
import matplotlib.pyplot as plt
#Creating gas object with solution class
gas = ct.Solution('gri30.xml')
# Reactor initial conditions
T = 1500
P = ct.one_atm
Species_dict = {'CH4':1,'O2':2,'N2':7.52}
# Providing initial values to gas object
gas.TPX = T, P, Species_dict
# Consider constant pressure ideal gas reactor
r = ct.IdealGasConstPressureReactor(gas,name='R1')
# Preapare reactornetwork
sim = ct.ReactorNet([r])
# Number of gas reaction to be consider
n = gas.n_reactions
#help(sim)
# Reaction array
reac_n = []
# Maximum sensitivity array
s_max =[0]*n
# Loop for appending reaction
for i in range(n):
# Adding sensitivity reactions
r.add_sensitivity_reaction(i)
reac_n.append(gas.reaction_equations([i]))
# Providing tolerences for numerical calculations
sim.rtol = 1e-6
sim.atol = 1e-15
sim.rtol_sensitivity = 1e-6
sim.atol_sensitivity = 1e-6
# Creating solution array for storing each state data
states = ct.SolutionArray(gas, extra = ['t_ms'])
# Time data
t_init = 0
t_end = 2e-3
dt = 5e-6
# Time array
time = np.arange(t_init,t_end,dt)
# Time loop
for t in time:
# Advancing time
sim.advance(t)
# Loop for various reactions sensitivity
for i in range(n):
s = sim.sensitivity('temperature',i)
if np.abs(s_max[i]) <= np.abs(s):
s_max[i] = s
states.append(r.thermo.state, t_ms = t*1e3)
# Loop for printing all reaction with their sensitivities
for i in range(n):
print('For reaction '+str(reac_n[i])+' sensitivity is '+str(s_max[i]))
# Number of reactions required to show in plot
N = 10
# Absolute values of sensitivities
s_abs = np.abs(s_max)
# Sorting logic using zip function
reac_n = [p for q, p in sorted(zip(s_abs,reac_n), reverse = True)]
s_max = [p for q, p in sorted(zip(s_abs,s_max), reverse = True)]
# Array to store top N sensitive reactions
max_reactions = reac_n[0:N]
# Array to store top N sensitivities
max_sensitivity = s_max[0:N]
# Loop to print top N sensitive reactions with their sensitivities
print('nTop '+str(N)+' maximum sensitive reactions in GRI3 mechanism in first '+str(n)+' reactions are as follows:')
for i in range(N):
print('For reaction '+str(max_reactions[i])+' sensitivity is '+str(max_sensitivity[i]))
# Plotting
plt.figure(1)
plt.subplot(2,3,1)
plt.plot(states.t_ms,states.T,color='red',linewidth=2)
plt.grid('both',linestyle='--',linewidth=0.8)
plt.xlabel('time [$ms$]',fontweight='bold')
plt.ylabel('Temperature [$k$]',fontweight='bold')
plt.subplot(2,3,2)
plt.plot(states.t_ms,states('OH').X,color='blue',linewidth=2)
plt.grid('both',linestyle='--',linewidth=0.8)
plt.xlabel('time [$ms$]',fontweight='bold')
plt.ylabel('OH mole fraction',fontweight='bold')
plt.subplot(2,3,3)
plt.plot(states.t_ms,states('H').X,color='green',linewidth=2)
plt.grid('both',linestyle='--',linewidth=0.8)
plt.xlabel('time [$ms$]',fontweight='bold')
plt.ylabel('H mole fraction',fontweight='bold')
plt.subplot(2,3,4)
plt.plot(states.t_ms,states('CH4').X,color='orange',linewidth=2)
plt.grid('both',linestyle='--',linewidth=0.8)
plt.xlabel('time [$ms$]',fontweight='bold')
plt.ylabel('CH4 mole fraction',fontweight='bold')
plt.subplot(2,3,5)
plt.plot(states.t_ms,states('N2').X,color='lavender',linewidth=2)
plt.grid('both',linestyle='--',linewidth=0.8)
plt.xlabel('time [$ms$]',fontweight='bold')
plt.ylabel('N2 mole fraction',fontweight='bold')
plt.subplot(2,3,6)
plt.plot(states.t_ms,states('CO2').X,color='purple',linewidth=2)
plt.grid('both',linestyle='--',linewidth=0.8)
plt.xlabel('time [$ms$]',fontweight='bold')
plt.ylabel('CO2 mole fraction',fontweight='bold')
plt.suptitle("Sensitivities of various species w.r.t. time",fontweight='bold')
plt.tight_layout()
plt.figure(2)
plt.gca().invert_yaxis()
plt.barh(range(N), max_sensitivity,color=['red', 'blue', 'purple', 'green', 'lavender'],height=0.5)
plt.yticks(range(N),max_reactions,fontsize=8)
plt.grid('both',linestyle='--',linewidth=0.5)
plt.ylabel('Chemical reactions',fontweight='bold')
plt.xlabel('Sensitivity',fontweight='bold')
plt.suptitle('Top '+str(N)+' reactions with higher sensitivities',fontweight='bold')
plt.tight_layout()
plt.show()
Output plots for above program:-
 |
| Fig a: Sensitivities for various species |
 |
| Fig b: Most sensitive reactions in methane combustion |
Output for above program:-
Top 10 maximum sensitive reactions in GRI3 mechanism in first 325 reactions are as follows:
For reaction ['CH3 + O2 <=> CH3O + O'] sensitivity is 16.84749748246271
For reaction ['2 CH3 (+M) <=> C2H6 (+M)'] sensitivity is -13.618280538390644
For reaction ['H + O2 <=> O + OH'] sensitivity is 13.103707537648402
For reaction ['CH4 + H <=> CH3 + H2'] sensitivity is -10.157244261797375
For reaction ['CH3 + O2 <=> CH2O + OH'] sensitivity is 6.779897254994966
For reaction ['CH3 + HO2 <=> CH3O + OH'] sensitivity is 6.27045632052444
For reaction ['2 CH3 <=> C2H5 + H'] sensitivity is 4.3453112093858275
For reaction ['CH2O + CH3 <=> CH4 + HCO'] sensitivity is 4.028984416275573
For reaction ['CH2O + O2 <=> HCO + HO2'] sensitivity is 4.026263736663865
For reaction ['CH3 + HO2 <=> CH4 + O2'] sensitivity is 3.437058351914913
The sensitivity analysis is essential in optimizing the performance of IC engines, especially the CI engines where the combustion process is initiated by the auto-ignition of charge. In such types of combustion process, the auto-ignition delay needs to precisely controlled with the movement of the piston to obtain optimum efficiency. Here, the most sensitive reactions are determined and the reaction rates are altered to get the desired ignition delay.
The sensitivity indicates the influence of the reaction on the output parameter such as temperature or molar fraction of species. The positive sensitivity indicates that increase in temperature which results in ignition delay. And on other side to increase ignition delay, the reactions with negative sensitivity are considered.
Considering above bar graph (fig b), by increasing the sensitivities of second and fourth reactions, the ignition delay will be increased. While increasing the sensitivities of rest of the reactions will lead to a decrease in ignition delay
Conclusion:-
By observing above results and plots it can be concluded that most sensitive reactions are sorted and for further analysis.
Above provided code is generated for computing sensitivities and it is parameterized for following two factors-
- 'N' which is number of reactions with maximum sensitivities to show on bar graph
- 'n' which is number of reactions to be consider in simulation from GRI mechanism
Reference:-