Week 9 - Senstivity Analysis Assignment

Aim:- The aim of this project is to find 10 most sensitive reactions affecting temperature of methane air combustion.

Introduction:- Chemical reactions comprise a number of intermediary reactions and species. Attempting to simulate the detailed intermediary steps is computationally expensive. An alternative, to this approach would be to select only those reactions which impact the parameter under study.To compute these reactions we use sensitivity analysis. By sensitivity of a reaction, we mean how a small change in reaction rate can affect the overall state parameter.

Mathematically it is defined as 

`S=-(dy)/(dp)`

In this project, sensitivity is computed using inbuilt function in Cantera.

Details about sensitivity and cantera code

i) Initial Temperature=1500L

ii) Initial Pressure = 1 atm

iii) Mixture- Stoichiometric Mixture of CH4 and air

iv) Absolute tolerance= 1e-16

vi) Relative Tolerance=1e-6

vii)Absolute tolerance of sensitivities= 1e-6

viii) Relative Tolerance of sensitivities= 1e-6

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

gas=ct.Solution('gri30.xml')

# Defining number of reaction and initial conditions
N=10
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])

# Extracting the reaction equationd and their numbers
reactions=gas.reaction_equations()
n=gas.n_reactions
r_id=np.linspace(1,n)

sim.rtol=1e-6
sim.atol=1e-15
sim.rtol_sensitivity=sim.rtol
sim.atol_sensitivity=sim.rtol

#No of sensitive reactions to be displayed
Tf=[0]*n

#Solution array that stores reaction properties at each time step
states=ct.SolutionArray(gas,extra=['t_in_ms'])

#Adding all the reaction from mechanism file to sensitivity analysis
for i in range(n):
    r.add_sensitivity_reaction(i)
    
    
#Outer time loop
for i in np.arange(0,2e-3,5e-6):
    sim.advance(i)
    states.append(r.thermo.state,t_in_ms=i*1e3)
    for j in range(n):
        T=sim.sensitivity('temperature',j)
        if abs(T)>=abs(Tf[j]):
            Tf[j]=T
#Sorting the reaction,reaction ids and maximum sensitivity using for loop in list comprehension            
reactions=[r for s, r in sorted(zip(np.abs(Tf),reactions),reverse=True)]
rid=[r for s, r in sorted(zip(np.abs(Tf),r_id),reverse=True)]
Tf=[T for abs_T, T in sorted(zip(np.abs(Tf),Tf),reverse=True)]

#Assigning 10 sensitive value of reactions,ids and sensitivity
r_id_plt=r_id[0:N]
splt=Tf[0:N]
rplt=reactions[0:N]

#Extracting the Mole Fractions
for i in range(N):
    print(r_id_plt[i],rplt[i],splt[i])
OH=states.X[:,states.species_index('OH')]
H=states.X[:,states.species_index('H')]
CH4=states.X[:,states.species_index('CH4')]
CH3=states.X[:,states.species_index('CH3')]
O=states.X[:,states.species_index('O')]
HO2=states.X[:,states.species_index('HO2')]
CH3O=states.X[:,states.species_index('CH3O')]

plt.style.use('seaborn-dark-palette')
plt.figure(1)
plt.grid('TRUE')
plt.barh(rplt,splt,height=0.6)
plt.gca().invert_yaxis()
plt.xlabel("Sensitivity")
plt.title("Temperature Sensitive Reactions : GRI Mechanism")

plt.tight_layout()
plt.figure(2)
plt.subplot(2,4,1)
plt.plot(states.t_in_ms,states.T)
plt.xlabel("Time in Miliseconds")
plt.title("Temperature in Kelvin")
plt.tight_layout()
plt.grid('TRUE')
plt.subplot(2,4,2)
plt.plot(states.t_in_ms,OH)
plt.xlabel("Time in Miliseconds")
plt.title("OH Mole Fraction")
plt.tight_layout()
plt.grid('TRUE')
plt.subplot(2,4,3)
plt.plot(states.t_in_ms,H)
plt.xlabel("Time in Miliseconds")
plt.title("H Mole Fraction")
plt.tight_layout()
plt.grid('TRUE')
plt.subplot(2,4,4)
plt.plot(states.t_in_ms,CH4)
plt.xlabel("Time in Miliseconds")
plt.title("CH4 Mole Fraction")
plt.tight_layout()
plt.grid('TRUE')
plt.subplot(2,4,5)
plt.plot(states.t_in_ms,CH3)
plt.xlabel("Time in Miliseconds")
plt.title("CH3 Mole Fraction")
plt.tight_layout()
plt.grid('TRUE')
plt.subplot(2,4,6)
plt.plot(states.t_in_ms,O)
plt.xlabel("Time in Miliseconds")
plt.title("O Mole Fraction")
plt.tight_layout()
plt.grid('TRUE')
plt.subplot(2,4,7)
plt.plot(states.t_in_ms,HO2)
plt.xlabel("Time in Miliseconds")
plt.title("HO2 Mole Fraction")
plt.tight_layout()
plt.grid('TRUE')
plt.subplot(2,4,8)
plt.plot(states.t_in_ms,CH3O)
plt.xlabel("Time in Miliseconds")
plt.title("CH3O Mole Fraction")
plt.tight_layout()
plt.grid('TRUE')

plt.show()

Results

1) The top 10 sensitive reactions are:

-----------------------------------------------------------------
Reaction ID    |           Reaction            |  Sensitivity
-----------------------------------------------------------------
    1.0        |      CH3 + O2 <=> CH3O + O    |    16.8347
   7.61        |   2 CH3 (+M) <=> C2H6 (+M)    |   -13.6080
  14.22        |          H + O2 <=> O + OH    |    13.0937
  20.84        |       CH4 + H <=> CH3 + H2    |   -10.1495
  27.45        |     CH3 + O2 <=> CH2O + OH    |     6.7748
  34.06        |    CH3 + HO2 <=> CH3O + OH    |     6.2657
  40.67        |         2 CH3 <=> C2H5 + H    |     4.3420
  47.29        |   CH2O + CH3 <=> CH4 + HCO    |     4.0259
   53.9        |    CH2O + O2 <=> HCO + HO2    |     4.0232
  60.51        |     CH3 + HO2 <=> CH4 + O2    |     3.4344
---------------------------------------------------------------

A positive sensitivity means that increase in rate of reaction would increase temperature and negative sensitivity means increase in rate of reacton would decrease temperature.

The above image shows mole fractions of various species. Temperature rises from 1500 K to 2775K at 1.15 miliseconds. Auti ignition is achieved.CH4 concentration reduces to 0 which means all fuel is spent. Other species also show a change in concentration which indicated that reaction progresses at 1500 K.

Conclusion:- Sensitivity analysis is carried out at top 10 reactions are found that affect temperature