Week 9 - Senstivity Analysis Assignment

#Aim:1.) To study sensitivity of reactions involved in combustion of Methane in presence of air in Stoiochiometric ratio.(Reactions taken from GRI mechanism).

2.)To figure out 10 most sensitive reactions w.r.t Temperature.

Introduction:

For a particular reaction like combustion of Methane ,at certain initial Temperature and Pressure it has fixed ignition delay.From GRI mech file of Methane,we can see that when combustion of methane takes place ,so many intermediate reactions takes place to form different species.With the help of Sensitive analysis of these reactions,either it is with respect to Temperature or species involved in the reactions.Ignition delay time and product of the reaction involved can be controlled as per our requirement.

#Coding Approach:

1.Required packages imported.

2.Initial condition and gas composition set in gas object.

3.Reactor object 'r' defined and Reactornet object 'sim' defined.

4.All 325 reactions imported imported from GRI mech.

5.Number of sensitivity reactions to be evaluated defined.

6.Total time for combustion of reaction and time step size defined.

7.Object 'states' defined which is a Solutionarray storing all properties of gas object.

8.Then 'sim' object advanced with time having time step size 5e-6 to 2e-3 seconds.

9.With this change in time ,sensitivity of reactions with temperature is captured.

10.From the data ,10 most sensitive reactions are found w.r.t sensitivity coefficient and plotted on bar graph.

-------------------------------------------*****************************--------------------------------

"""
Constant-pressure, adiabatic kinetics simulation with sensitivity analysis
"""

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

gas = ct.Solution('gri30.xml')                          #A instance object (gas) of Solution class created
temp = 1500.0                                           #Temperature defined
pres = ct.one_atm                                       #Pressure set one_atm.

gas.TPX = temp, pres, 'CH4:1, O2:2, N2:7.52'            #Composition of gas object is defined.
r = ct.IdealGasConstPressureReactor(gas,name='R1')      #Reactor object r from class  IdealGasConstPressureReactor defined,it contains gas object

sim = ct.ReactorNet([r])                                #Object sim from class ReactorNet defined having reactor r 


# Enable sensitivity with respect to the rates of all 325 reactions (reactions 0 through 324)
n=325                                                   #Total reactions in Gri mech file.

#No of sensitive reaction to be found out
x=10

Total_time=2e-3
Time_step_size=5e-6
N=int(Total_time/Time_step_size)                    #No of iteration


for i in range(n):
    r.add_sensitivity_reaction(i)

# Set the tolerances for the solution.
sim.rtol = 1.0e-6                                   #Relative Tolerance
sim.atol = 1.0e-15                                  #Absolute Tolerance


# Set the tolerances  for the sensitivity coefficients.
sim.rtol_sensitivity = 1.0e-6                       #Relative Tolerance
sim.atol_sensitivity = 1.0e-6                       #Absolute Tolerance




reaction_name=[]                                    #List containing all the reactions
for k in range(0,n):
    sens_para=sim.sensitivity_parameter_name(k)     #Gives reaction name
    reaction_name.append(sens_para)



states = ct.SolutionArray(gas, extra=['s'])         #Solution array defined having all properties of gas object


                
t=0                                                 #Initial time set to zero
s=np.zeros((N,n))                                   #Array of zeros defined for sensitivity


for i in range(N):                                  #Outer loop for time advancing
    t=t+5e-6
    for j in range(0,n):                            #Inner loop for selecting all the reaction
        
        sim.advance(t)                              #advance() method useed to change state of object w.r.t time (t)
       
        s[i,j] =  sim.sensitivity('temperature',j)
        states.append(r.thermo.state, s=s)

        print('%10.3e %10.3f %10.3f %14.6e  %10.5f ' %
              (sim.time, r.T, r.thermo.P, r.thermo.u,s[i,j]))



s_abs=np.abs(s)                                     #Make all sensitivity values absolute in s
s_abs_max=s_abs.max(axis=0)                         #Max value of sensitivity found for each reaction
s_abs_max_asc=np.argsort(s_abs_max)[-x:]            #All sensitivity values for each reaction arranged in ascending order and last 10 values index chosen

s_10sens_reaction_name=np.array([sim.sensitivity_parameter_name(i) for i in s_abs_max_asc])     #Array storing name of 10 most sensitive reactions

s_10_value=np.zeros((x,N))                                  #Array for storing sensitivity values of 10 most sensitive reactions
s_10_value_max=np.zeros(x)                                  #Array storing max values of all 10 reactions (Positive or Negative max)


for m in range(x):
    a=s_abs_max_asc[m]
    for p in range (N):
        s_10_value[m,p]=s[p,a]



for r in range(10):
    s_10_value[r].sort()                                    #Values are arranged in ascending order
    if np.abs(s_10_value[r,0])<np.abs(s_10_value[r,-1]) :
        s_10_value_max[r]=s_10_value[r,-1]
    else:
        s_10_value_max[r]=s_10_value[r,0]




y_pos = np.arange(len(s_abs_max_asc))
plt.rcdefaults()
fig, ax = plt.subplots()
ax.barh(y_pos, s_10_value_max,align='edge')
ax.set_yticks(y_pos)
ax.set_yticklabels(s_10sens_reaction_name)
ax.invert_yaxis()                                          # labels read top-to-bottom
ax.set_xlabel('Sensitivity')
ax.set_title('Plot showing trend of sensitivity of a reaction with Temperature')
plt.tight_layout()
plt.show()

--------------------------------------------***************************------------------------------

Analysis:

Above shown is the plot of 10 most sensitive reactions involved in the combustion of Methane.So

what it indicates,right now if we do combustion using given initial condition,there would be some

value for time for ignition delay.If this ignition delay value we want to increase or decrease ,due to

lesser efficiency of combustion inside the system we are combusting it.We can adjust the ignition

delay time by increasing or decreasing the reaction rate of above reactions.If we want to increase

ignition delay, we will increase the reaction rate of the above reactions having negative sensitivity

coefficient and vice versa.