Week 9 - Senstivity Analysis Assignment

SENSITIVITY ANALYSIS OF GRI MECH3.0

Objective

• Write a code that takes all the reactions from the GRI mechanism and calculates 10 most sensitive reactions
• The sensitivity parameter should be with respect to the Temperature (It should find the sensitivity with respect to the temperature not 'OH' as shown in the video)
• Consider the initial conditions as shown in the video
• program should show the reaction strings in the order of greatest sensitivity to the least sensitivity

Introduction

It is a process to get a relationship between parameters and the output of a model.Combustion process have multiple intermediary reactions and species.which increases the total number of reaction in mechanism.Simulating combustion using a detailed mechanism takes a lot of computational power and time.So to reduce the size of mechanism,only those reactions are taken in account which highly affect the parameter under study.To find the reaction sensitivity analysis is carried out.

Sensitivity of a reaction means how much does a small change in reaction rate affect the overall state parameter.

-ve sensitive coefficients:

Increasing the reaction rate will increase ignition delay and vice-versa

+ve sensitive coefficients:

Increasing reaction rate will reduce ignition delay and vice-versa

• By figuring out the sensitive reaction we can tune the mechanism for our need
• Even build our own mechanism
• Identify which reactions controls the ignition delay

example:

reference:https://www.researchgate.net/publication/245352598_Ignition_of_Lean_Methane-Based_Fuel_Blends_at_Gas_Turbine_Pressures

Code:

for all 325 reactions in GRI mech 3.0

#importing Cantera,Matplotlib,time 
import numpy as np
import cantera as ct
import matplotlib.pyplot as plt

#total number of reactions you want to check sensitivity
n=325

#no.of reaction you want to plot
N=10


#Calling the GRI Mech 3.0 ,creating a gas object
gas=ct.Solution('gri30.xml')


#Setting Temperature,Pressure,Molefraction or concentrations
gas.TPX=1500,1*101325,{'CH4':1,'O2':2,'N2':7.52}

#create a reactor object which takes gas object
r=ct.IdealGasConstPressureReactor(gas)
  #create a reactor net object
sim=ct.ReactorNet([r])

#loop foradd sensitivity
for i in range(n):
    r.add_sensitivity_reaction(i)
    

sim.rtol=1e-6 #Relative Tolerance
sim.atol=1e-15 #absolute tolerance

sim.rtol_sensitivity=1e-6 #relative tolerance for sensitivity
sim.atol_sensitivity=1e-6 # absolute tolarance for sensitivity

#create an empty array
s_max=[0]*n

splot=[]#maximum temperature sensitivities of each rxn
states=ct.SolutionArray(gas,extra=['time_in_ms','t'])


#loop for calculating sensitivty for the n rxn
for t in np.arange(0,2e-3,5e-6):
    sim.advance(t)
    
    #storing the max sensitivites of each rxn
    for j in range(n):
        s=sim.sensitivity('temperature',j)
        states.append(r.thermo.state,time_in_ms=t*1000,t=t)
        
        #condition applied for sensitive rxns
        if np.abs(s) >= np.abs(s_max[j]):
            s_max[j]=s
    
#storing sensitive values
s=np.zeros((N))


#loop for extracting top 10 reation
for z in range(N):
    maxi=np.max(s_max)
    mini=np.min(s_max)
    
    if np.abs(maxi)>np.abs(mini):
        s[z]=maxi
        reaction_index=s_max.index(maxi)
        splot.append(sim.sensitivity_parameter_name(reaction_index))
        s_max.remove(maxi)
                     
    else:
        s[z]=mini
        reaction_index=s_max.index(mini)
        splot.append(sim.sensitivity_parameter_name(reaction_index))
        s_max.remove(mini)
 
  
 #reversing the plot for preference   
react=list(reversed(splot))
sensitivity=list(reversed(s))

#plot
plt.figure(1)
plt.title('Top 10 reaction from 325 reaction from GRI 3.0')
plt.xlabel('sensitivity')
plt.barh(react,sensitivity,color='g')


plt.figure(2)
plt.plot(states.t,states.T)
plt.title('temp v/s ig delay')
plt.xlabel('ig delay in ms')
plt.ylabel('temperature')
        
plt.figure(3)
    
plt.plot(states.t,states('oh').Y, color='b')
plt.xlabel('time in ms')
plt.ylabel('molefraction of OH')

plt.figure(4)
    
plt.plot(states.t,states('O').Y, color='r')
plt.xlabel('time in ms')
plt.ylabel('molefraction of O')
    
plt.figure(5)

plt.plot(states.t,states('ch4').Y, color='g')
plt.xlabel('time in ms')
plt.ylabel('molefraction of CH4')

        
plt.show()    

Results:

CH3 + O2 <=> CH3O + O	18.634545339295258
CH3 + H2O2 <=> CH4 + HO2	-14.407046644685927
H + O2 <=> O + OH	13.192167401983914
CH3 + H (+M) <=> CH4 (+M)	-10.70490452590123
CH2(S) + CO2 <=> CH2O + CO	7.56074121958509
CH2 + HO2 <=> CH2O + OH	6.220250346300111
CH2(S) + CO2 <=> CH2O + CO	4.572536959129847
C2H6 + CH2(S) <=> C2H5 + CH3	4.362578578510714
CH2O + O2 <=> HCO + HO2	4.355616300367177
2 HO2 <=> H2O2 + O2	4.1166238776298005

here comparing the data we can see that" CH3 + O2 <=> CH3O + O"----18.63 is one of the most sensitive reaction.

tuning these reactions

"CH3 + O2 <=> CH3O + O" or "CH3 + H2O2 <=> CH4 + HO2" can change the ignition delay parameter.

Suppose "CH3 + O2 <=> CH3O + O" with sensitivity +18.63,if we are increasing the reaction rate ,this will reduce ignition delay and vice-versa.

"CH3 + H2O2 <=> CH4 + HO2" with sensitivity -14.40,reducing the rate of this reaction will decrease the ignition delay and vice-versa

the figures below shows the variation of temperature,OH,O,CH4 molefractions of the mechanism and the spikes indicating the ignition delay

References:

https://www.researchgate.net/publication/245352598_Ignition_of_Lean_Methane-Based_Fuel_Blends_at_Gas_Turbine_Pressures

http://www.cerfacs.fr/cantera/