Sensitivity Analysis of GRI Mechanism.

Sensitivity Analysis:

Ignition delay is sensitive to the reaction mechanism, i.e some reaction is more sensitive to ignition delay and some are less sensitive. Hence, \'sensitivity analysis\' plays a vital role in tune the reaction mechanism.

What is the sensitivity analysis and why?

• Sensitivity analysis is the process to get a relationship between the parameters and output of a model (eg: effect of reaction 1 rate on temperature or species concentration).
• It is used to identify the most sensitive reactions and tune the reaction mechanism(i.e, to adjust the rate coefficient ).
• For these sensitive reactions tune the rate coefficient so that will match for ignition delay.

Generally, we are concerned with the temperature rather than other parameters to measure sensitivity because the ignition delay is more sensitive to the temperature when compared to other parameters.

The generic ODE system is represented as:

`dZ/dt=f(Z,t,a)`

where a is the net reactant rate multiplier (Perturbation factor)

Z is the vector of species and temperature and it is given by

`Z= 〖(T,y_1,y_2………y_k)〗^T`

Where T is the temperature

`y_k` is the mass fraction of the species k

The sensitivity matrix is given by the

`w_(j,i)=(∂Z_j)/(∂a_i )`

where `a_i` is the reaction rate constant

The size of the sensitivity matrix is \'k+1,n\'. here k represents species and n represents reaction.

Objective:

The objective of this challenge is to write a code to find 10 most sensitive reactions in the GRI mechanism using python and cantera.

Code:

• Create an instant gas object in solution class of cantera and define temperature and pressure attributes.
• create a reactant instant which uses the \'IdealGasConstPressureReactor\' governing equation and then create instant for the net reactor which takes the r as the input.
• Add all the sensitivity reactions into the reactant instant object which present in the \'gri30 mechanism\'. we have 325 reactions and 45 species.
• Define absolute tolerance (`∆x`) and relative tolerance (`∆x⁄x`). Here, we need to check both relative and absolute tolerances for every variable because these tolerances are converging criteria.
• Define the absolute tolerance (`∆x`) and relative tolerance (`∆x⁄x`) for the sensitivity reactions.
• Create an instant object called states of solution array class which takes gas as an argument. This \'SolutionArrray\' describes the state of the reactor object.
• The time loop has been defined and considered it would run for 2 milliseconds with a step size of 5e-6.
• Calculate the sensitivity of the reaction at every time step by advancing through \'t\'. Find the maximum sensitivity value for each reaction ( through the length of timestep), and append this maximum values to the \'max_s\' and reactions to \'reaction\'.
• Then, create a tuple and couple both the arrays and the sort the tuple based on the sensitivity value that obtained at each reaction by using \'sorted\' command.
• From the sorted list append the top ten most sensitivity reaction to the \'s3\' and create a table for the top ten sensitivity reactions
• Plot the graph between the reaction and the sensitivity.

import numpy as np 
import matplotlib.pyplot as plt 
import cantera as ct 
# create a gas instant obeject from the solution class
gas = ct.Solution(\'gri30.xml\')

T=1500 #k
P = ct.one_atm
# set temperature pressure and molefraction
gas.TPX = T,P,{\'CH4\':1,\'O2\':2,\'N2\':7.52}

# create a reactant instant which uses the IdealGasConstPressureReactor governing equation
r = ct.IdealGasConstPressureReactor(gas)
# create instant for the net reactor which takes the r as the input and basically r is an array 
sim = ct.ReactorNet([r])

for i in range(325):
	r.add_sensitivity_reaction(i) # adding the all the sensitivity reactions into the analysis
 
sim.rtol = 1e-6 # realtive change sensitivity should be less than the realtive tolerance
sim.atol = 1e-15 # absolute change in the tolerance
# set tolerance for the sensitive spieces because these are also the stiff system 
sim.rtol_sensitivity = 1e-6 # these are for the sensitive reactions
sim.atol_sensitivity = 1e-6 # the values are higher than the species value

n=0 # for integer values
T=np.arange(0,2e-3,5e-6)
states=ct.SolutionArray(gas,extra=[\'t\'])

s=np.zeros((325,len(T)))
\"\"\"
So, here we are calculating the sensitivity value at the each time step for every reaction
\"\"\"
for t in T:
	sim.advance(t)
	# j represents reactions.
	# 1 reprents the temperature.
	for j in range(325):
		s[j,n] =sim.sensitivity(1,j)
	n=n+1
	states.append(r.thermo.state, t=t*1000)

#print(s)

reaction = []
max_s=[]
c=0
for i in range(325):
	m=np.max(abs(s[i,:]))
	max_s.append(m)
	c=c+1
	reaction.append(c)
	#print(max_s)
	#print(\'maximum value is :\',max_s)

m_list = tuple(zip(reaction,max_s)) # couple the two arrays
s1=sorted(m_list,key=lambda m_list:-m_list[1],reverse = False) # sort the m_list based on sensitivity 
#print(s1)

a=0
s3=[]
for i in range(0,10):
	s2=s1[i]
	s3.append(s2)
	a=a+1
	# Creating the table for the top 10 sensitivity reactions
	if i==0:
		print(\"s.no\\t(reaction_number,Sensitivity)\\n\")
	print(i+1,\'\\t\\t\',s1[i])

#print(s3)
# plotting
plt.subplot(221)
plt.plot(states.t,states.T)
plt.xlabel(\'Time(ms)\')
plt.ylabel(\'Temperature(k)\')
plt.subplot(222)
plt.plot(states.t,states(\'OH\').X)
plt.xlabel(\'Time(ms)\')
plt.ylabel(\'OH mole fraction\')
plt.subplot(223)
plt.plot(states.t,states(\'H\').X)
plt.xlabel(\'Time(ms)\')
plt.ylabel(\'H mole fraction\')
plt.subplot(224)
plt.plot(states.t,states(\'CH4\').X)
plt.xlabel(\'Time(ms)\')
plt.ylabel(\'CH4 mole fraction\')
plt.tight_layout()

plt.figure(2)
plt.plot(s[0,:])
plt.xlabel(\'Time step\')
plt.ylabel(\'Sensitivity\')
plt.title(\'Sensitivity at each Timestep\')
plt.figure(3)
plt.barh(*zip(*s3),linewidth=10)
plt.xlabel(\'Temperature Sensitivity\')
plt.ylabel(\'Reactions\')
plt.show()

Result:

Topmost sensitive reaction n the GRI mechanism file ( which contains 325 reactions).

Ignition Delay and Molar concentration Plot of Methane gas:

The above plots describe the ignition delay of the methane gas and the concentration of the H, OH, CH4 moles in the methane mechanism.

The plot of sensitivity at each Time Step:

The above graph illustrates the sensitivity of each reaction at every time step (Sensitivity of the reaction 1 at each value of T(time step)).

The plot of 10 most sensitive Reactions:

The above graph illustrates the sensitivity of the 10 most sensitive reactions out of 325 reactions in the GRI30 mechanism file. Hence, we can tune this ten reaction coefficient to increase or decrease the ignition delay.