Week 9 - Senstivity Analysis Assignment

AIM:To perform sensitivty analysis of reaction from GRI Mechanism reactions

Introduction:

A Combustion reaction contains thousands of intermediate reactions each having their individual rate constants.The reaction mechanism produces a number of observable parameters such as temperature profiles,laminar flame speed, species concentration.These depend largely on operating conditions of the reaction.

Sensitivity Analysis is an approach of examining the relationship between the input parameters and the observable output.The reaction mechanism can be tuned by changing the individual rate constants and thus the reaction mechnism is optimised accordingly.

The rate of change of output function due to change in input parameter is Sensitivity Coefficient.

It is represented as follows:

S=∂y∂p

The output function is represented by matrix Z.

Z=(T,Y1,Y2.....Yk)

Wji=∂Zj∂ai

Here a is the reaction , considering i Intermediate Reactions.

Method:

The cantera module in Python is used to perform sensitivity analysis.

A gas object is created. The reactor object containing the gas is created and added into a reactor network.The reactor network is advanced with timesteps for 2ms.The SolutionArray object in cantera is used to store the state of the system for every timestep.

Sensitivity analysis is carried out for temperature, and the most sensitive reactions and found.

The python script below generates a text file listing the most sensitive reactions and also generates plots of the same. The temperature and concentration profiles of species are also plotted. The initial Temperature and Pressure are taken as 1500K and one atmosphere respectively.

Python Code:

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

#gas object
gas=ct.Solution('gri30.xml')
T=1500
P=ct.one_atm
gas.TPX=T,P,{'CH4':1,'O2':2,'N2':7.52}

#reactor and reactor network
r=ct.IdealGasConstPressureReactor(gas,name='R1')
sim=ct.ReactorNet([r])

reac=gas.reaction_equations()
#print(reac)
n_reac=gas.n_reactions
n=10
S_max=[0]*n_reac
label=[]

for i in range(n_reac):
	r.add_sensitivity_reaction(i)
	#label.append(r.sensitivity_parameter_name(i))

#tolerences
sim.rtol=1e-6
sim.atol=1e-16

sim.rtol_sensitivity=1e-6
sim.atol_sensitivity=1e-6
states=ct.SolutionArray(gas,extra='t')
reac_num=np.arange(1,n_reac+1)
for t in np.arange(0,2e-3,5e-6):
	sim.advance(t)
	states.append(r.thermo.state,t=t*1000)


	for j in range(n_reac):
		s=sim.sensitivity('temperature',j)
		#plt.plot(t,s)
		if np.abs(s)>np.abs(S_max[j]):
			S_max[j]=s


#sorting
S_abs=np.abs(S_max)
reac=[a for b, a in sorted(zip(S_abs,reac),reverse=True)]
reac_num=[a for b,a in sorted(zip(S_abs,reac_num),reverse=True)]
S_max=[a for b, a in sorted(zip(S_abs,S_max),reverse=True)]

max_r=reac[0:n]
max_n=reac_num[0:n]
max_s=S_max[0:n]

#mole fractions
H_conc=states('H').X
OH_conc=states('OH').X
CH4_conc=states('CH4').X

#listing reactions
f=open("reaction_list.txt","w+")
f.write("The 10 most Temperture sensitive reactions are as follows:\n")
f.write("********************************************************* \n")
f.write("\n %8s|%30s|%10s\n"%('Reaction No','Reaction','Sensitivity'))
f.write("--------------------------------------------------------- \n")

for i in range(n):
	f.write("%12d|%30s|%10.5f\n"%(max_n[i],max_r[i],max_s[i]))
f.write("********************************************************* \n")	
f.close()

#plotting
plt.figure(1)
plt.barh(max_r,max_s,0.4,color='brown')
plt.xlabel('Temperature Sensitivity')
plt.grid( linestyle=':', color='black')
plt.title('Highest Temperature Sensitive Reactions in GRI30 for Stoichiometric Combustion of Methane')
plt.figure(2)
plt.subplot(2,2,2)
plt.plot(states.t,H_conc)
plt.xlabel('Time(ms)')
plt.ylabel('H Mole Fraction')
plt.grid(True)
plt.subplot(2,2,1)
plt.plot(states.t,states.T)
plt.xlabel('Time(ms)')
plt.ylabel('Temperature(K)')
plt.grid(True)
plt.subplot(2,2,3)
plt.plot(states.t,OH_conc)
plt.xlabel('Time(ms)')
plt.ylabel('OH Mole Fraction')
plt.grid(True)
plt.subplot(2,2,4)
plt.plot(states.t,CH4_conc)
plt.xlabel('Time(ms)')
plt.ylabel('CH4 Mole Fraction')
plt.grid(True)
plt.tight_layout()

plt.show()

Result:

The text file generated is given in the link below. A screenshot of the same is given below.

Reaction Index	Reaction	Sensitivity
0	H + O2 <=>O+OH	16.84801
1	CH3+H(+M)<=>Ch4(+M)	-13.61874
2	CH20+02<=>HCO+H02	13.10419
3	CH3+0H(+M)<=>CH30H(+M)	-10.15745
4	H+HCO(+M)<=>CH2O(+M)	6.78009
5	CH2+H(+M)<=>CH3(+M)	6.27081
6	CH2+OH<=>CH20+H	4.34541
7	CH3+0<=>CH20+H	4.02909
8	CO2+2H<=>CO2+H2	4.02639
9	C2H5+H<=>C2H4+H2	3.43721

The bar chart below indicated the various sensitivities when the Equivalence Ratio is one. The positive sensitivity implies that increasing the rates of these reactions increases the temperature.The negetive sensitivity implies that the increase in rate of these reactions causes the temperature to reduce.

Thus if the ignition delay time has to be reduced, the temperature needs to rise sooner which can be done by incresing the rates of the reactions with positive sensitivities and reducing the rates of reactions with negetive sensitivities.

The plots below indicate various profiles during the duration of the reaction. the first plot is the plot of Temperature against time and it shows the increase in temperature as the mixture ignites.The next three plots are of the mole fractions of species H, OH and CH4. A sudden increase can be observed in the concentrations of H and OH around the time of ignition of the mixture indicating the formation of these species.the concentration of these two then starts to reduce indicating that these are intermediate products and are transitioning into some other species. While in the case of CH4 the mole fraction can be seen to suddenly fall at auto ignition as it is the main reactant and is being converted into products.