Week 9 - Senstivity Analysis Assignment

Aim

To write a code to import all the reactions from the GRI mechanism and calculate the 10 most sensitive reactions with respect to temperature parameter

1. Introduction

Sensitivity analysis is the study performed to determine the influence of a reaction parameter such as species concentration, pressure, temperature etc. on the output results. Since a combustion reaction has several intermediary species involved, the sensitivity index needs to be calculated for each reaction. Based on this study, it can be concluded that higher is the sensitivity of the intermediary reactions, more is the influence on the output parameters. In other words, to observe controlled change in the output results, the reactions with the highest sensitivity index has to be controlled with respect to the reaction parameter.

The sensitivity analysis plays a crucial role in the study of ignition delay of a combustion process. To reduce the ignition delay, the first step is to determine the sensitivity coefficient of all intermediary reaction with respect to the reaction parameter. Once this has been determined, the most sensitive reactions are observed and the reaction parameter is changed to achieve the desired results. In this study, we use Cantera module to determine the sensitivity coefficients of the intermediary reactions of CH4 combusion of GRI mechanism. This result is used to determine the ten most sensitive reactions with respect to reaction parameter temperature. The concentrations of the intermediary species and the variation of temperature with respect to time are plotted to comment on the ignition delay.

2. Program algorithm and code

The following steps are followed to write the program for conducting the sensitivity analysis

1. Import numpy and cantera module
2. Create ‘gas’ object of Solution Method
3. Initialize the temperature, pressure and molar concentrations
4. Create ‘r’ object of IdealGasConstPressureReactor Method
5. Create ‘sim’ object of Reactor Network
6. Enable sensitivity with respect to reactions to 1 to N
7. Define add_sensitivity_reaction(i) for each reaction
8. Define tolerances for solution and sensitivity coefficients (absolute and relative)
9. Create ‘states’ object of Solution Array by passing the ‘gas’ object
10. Advance ‘sim’ object and define sim.sensitivity for species of reaction and store the maximum value in the variable s_max
11. Using sorting, the values of s_max are sorted and the top 10 reaction values are stored in sen_order
12. Using plot commands, a horizontal bar graph for Reaction vs Sensitivity coefficient is plotted
13. Using plot commands, subplots for variation of reaction parameters with respect to time is plotted.

The program is written as follows:

""" 

Program to perform sensitivity analysis on methane combustion

Written by

Vignesh V

"""

import numpy as np 
import cantera as ct 
import matplotlib.pyplot as plt

#object of Solution method
gas= ct.Solution("gri30.xml")
Temp=1500 
Pressure=ct.one_atm
gas.TPX = Temp,Pressure,{'CH4':1,'O2':2,'N2':7.52}

#create object of Ideal Gas Reactor
r = ct.IdealGasConstPressureReactor(gas, name='R1')

#create object of Reactor Network
sim = ct.ReactorNet([r])

#number of reactions in GRI mechanism
n=gas.n_reactions

#enable sensitivity wrt reactions 1 to n
for i in range(n):
	r.add_sensitivity_reaction(i)


#set relative and absolute tolerances for solution and sensitivity analysis
sim.rtol = 1e-06
sim.atol = 1e-15

sim.rtol_sensitivity= 1e-06
sim.atol_sensitivity= 1e-06

states = ct.SolutionArray(gas, extra =['t'])

s_max = [0]*n

for t in np.arange(0,2e-3,5e-6):
	sim.advance(t)
	for i in range(n):
		s = sim.sensitivity('temperature',i)
		if(abs(s)>abs(s_max[i])):
			s_max[i] = s
	states.append(r.thermo.state, t=1000*t)

#sorting top 10 sensitive reactions

N=10

sen = [0]*n
sen_order= [0]*N
react_number= [0]*N

for i in range(N):
	for j in range(n):
			if(abs(np.min(s_max))>np.max(s_max)):
				sen_order[i] = np.min(s_max)
			else:
				sen_order[i] = np.max(s_max)	
			react_number[i] = s_max.index(sen_order[i])
	s_max.remove(sen_order[i])			

sens_react = [0]*N

print('Sensitivity index for the 10 most sensitive reactions')
for i in range(N):
	sens_react[i]= sim.sensitivity_parameter_name(react_number[i])
	print('s=',sen_order[i],'for reaction',sens_react[i])		

plt.figure(1)
plt.barh(sens_react,sen_order)
plt.xlabel('Sensitivity')
plt.title('Sensitivity index of 10 most sensitive reaction')

plt.figure(2)
plt.subplot(2,2,1)
plt.plot(states.t,states.T)
plt.xlabel('Time (in ms)')
plt.ylabel('Temperature (in K)')

plt.subplot(2,2,2)
plt.plot(states.t,states('CH4').X)
plt.xlabel('Time (in ms)')
plt.ylabel('Mole fraction of CH4')


plt.subplot(2,2,3)
plt.plot(states.t,states('H').X)
plt.xlabel('Time (in ms)')
plt.ylabel('Mole fraction of H')


plt.subplot(2,2,4)
plt.plot(states.t,states('OH').X)
plt.xlabel('Time (in ms)')
plt.ylabel('Mole fraction of OH')
plt.show()

3. Results

The output of the program shows the sensitivity coefficient for the ten most sensitive intermediary reactions.

The same is plotted using a bar graph as follows:

It can be seen from the figure that the most sensitive reaction is CH3 + O2 ó CH3O + O with a sensitivity coefficient of 16.84. Another observation is the presence of positive and negative values of sensitivity coefficients.

The positive coefficient of sensitivity indicates that with a rise in the reaction parameter- temperature, the rate of the reaction will increase and contribute to faster net combustion reaction. This will result in decrease of the ignition delay. Whereas a negative sensitivity coefficient indicates slower sensitivity to temperature rise and will result in increase of ignition delay.

The ignition delay can be determined by plotting the reaction parameter variation with respect to the time advancement. The above figure shows the variation of temperature and molar concentrations of CH4, OH and H species with respect to time. It can be concluded that the ignition delay for the given reaction parameters is 1 ms.