Week 9 - Senstivity Analysis Assignment

Finding top 10 most sensitive reactions during combustion of methane

Explanation:

Step1:Set initial conditions,Temperature=1500K, pressure=1 atm

Step2: Define reactor, reactor network and add sensivity to all the reactions, 325 reactions are given in the gri mech3.0 page for combustion of methane. This can also be found by using gas.n_reactions command .To find sensitivity of species with respect to temperature we write sim.sensitivity(1,reaction no.)

Step 3: Running a time loop with final time as 2e-3 and dt as 5e-6 (arbitarily taking small values). Running a for loop inside this loop, for each reaction and comparing the value of sensitivity with earlier time step's max value to find the maximum sensitivity for each reaction from all the time steps .

step4:sens_max is the array here which keeps all the max values of sensitivity each reaction, sabs finds the absolute of these sensitivity values which helps in finding the order of sensitivities(high to low). Running a loop to find the max value and removing that max value from the loop after that to find the next maximum.

Step5 : change n value just to change the no. of reactions that are to be plotted.

The code:

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



	
gas=ct.Solution('gri30.xml')
temp=1500
pressure=ct.one_atm
gas.TPX=temp,pressure,'CH4:1,O2:2,N2:7.52'
r=ct.IdealGasConstPressureReactor(gas,name='R1')
sim=ct.ReactorNet([r])
nor=gas.n_reactions
print(nor)

for i in range(nor):
		r.add_sensitivity_reaction(i)
sim.rtol=1e-6
sim.atol=1e-15
sim.rtol_sensitivity=1e-6
sim.atol_sensitivity=1e-6
states=ct.SolutionArray(gas,extra=['t','s2'])
#no. of reactions=20

sens_max=[0]*nor	
for t in np.arange(0,2e-3,5e-6):
	sim.advance(t)
	for j in range(nor):
		s2=sim.sensitivity(1,j)
		if (abs(s2)>abs(sens_max[j])):
			sens_max[j]=s2
			sim.sensitivity_parameter_name(j)

			
s3=sens_max
print(sens_max)


#print(len(s3))
sabs=[abs(elem) for elem in s3 ]
#print(sabs)
c=[]
reaction_all=[]
#n=no. of reactions that are to be plotted
n=10
for m in range(n):
	B=sabs.index(max(sabs))
	reaction=(sim.sensitivity_parameter_name(B))
	reaction_all.append(reaction)
	print(sabs[B])
	print(s3[B])
	c.append(sens_max[B])
	sabs.remove(sabs[B])
	s3.remove(s3[B])
print(c)
print('interval')
print(reaction_all)
	
nor_list=list(range(1,n+1))
fig=plt.figure(figsize=(10,50))
plt.barh(reaction_all,c)

plt.ylabel('sensitivity')
plt.show()

Observation: y-axis shows the reactions and x axis shows the sensitivities for these reactions.

Table of sensitivities:

Reactions	Sensitivity
CH3 + O2 <=> CH3O + O	16.84749748
CH3 + H2O2 <=> CH4 + HO2	-13.61828054
H + O2 <=> O + OH	13.10370754
CH3 + H (+M) <=> CH4 (+M)	-10.15724426
CH2(S) + CO2 <=> CH2O + CO	6.779897255
CH2 + HO2 <=> CH2O + OH	6.270456321
CH2(S) + CO2 <=> CH2O + CO	4.345311209
C2H6 + CH2(S) <=> C2H5 + CH3	4.028984416
CH2O + O2 <=> HCO + HO2	4.026263737
2 HO2 <=> H2O2 + O2	3.437058352