Week 9 - Senstivity Analysis Assignment

Sensitivity Analysis is used to identify the parameters that affet a computed solution the most. It has proven invaluable in combustion studies, where the number of parameters is typically large. The main aim of the following code is to calculate the ignition delay of methane using a sensitivity analysis which will show us which reaction are more sensitive towards the ignition delay.

Aim:- To write a code that takes all the reaction from the GRI mechanism and calculate 10 most sesitivity reactions. The sensitivity parameters should be with respect to temperature.

Code:-

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

gas = ct.Solution('gri30.xml')
temp=1500
press=ct.one_atm


gas.TPX=temp,press, {'CH4':1,'O2':2,'N2':7.52}
r=ct.IdealGasConstPressureReactor(gas)
sim=ct.ReactorNet([r])

n=gas.n_reactions
N=10

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

sim.rtol=1.0e-6
sim.atol=1.0e-15

sim.rtol_sensitivity=1.0e-6
sim.atol_sensitivity=1.0e-6

sT_m=[0]*n
states=ct.SolutionArray(gas,extra=['time_in_ms','t'])

for t in np.arrange(0,2e-3,5e-6):
	sim.advance(t)

	for j in range(n):
      sT = sim.sensitivity('temperature', j)
      states.append(r.thermo.state,time_in_ms=t*1000,t=t)

      if np.abs(sT)>=np.abs(st_m[j]):
      	sT_m[j]=sT

s_plot=[]
s=np.zeros((N))

print("10 Most sensitivity reactions")
print("%8s  %4s %3s %4s " % ('Reaction ID','|', 'sensitivity','|','Reaction'))


for m in range (N):
  Maximum = np.Max(sT_m)
  Maximum=np.min(sT_m)

  if np.abs(Maximum)>np.abs(minimum):
    s[m]=Maximum
    Reaction_index=sT_m.index(Maximum)
    s_plot.append(sim.sensitivity_parameter_name(Reaction_index))
    print('%7d %8s %10.4f %4s %26s'%(Reaction_index, '|',minimum,'|',s_plot[m]))
    sT_m.remove(minimum)

  else:
  	s[m]=minimum
    Reaction_index=sT_m.index(minimum)
    s_plot.append(sim.sensitivity_parameter_name(Reaction_index))
    print('%7d %8s %10.4f %4s %26s'%(Reaction_index, '|',minimum,'|',s_plot[m]))
    sT_m.remove(minimum


for k in range (N):
	s_plot[k]=s_plot[k].replace('IdealGasConstPressureReactor_1:','')


reactions=list(reversed(s_plot))
sensitivity=list(reversed(s))


plt.figure(1)
plt.subplot(2,2,1)
plt.plot(states.t,states.T,color='black')
plt.xlabel('Simulation Time[ms]')
plt.ylabel('temperature')
plt.subplot(2,2,2)
plt.plot(states.t,states('OH').X,color='black')
plt.xlabel('Simulation Time[ms]')
plt.ylabel('OH mole fraction')
plt.subplot(2,2,3)
plt.plot(states.t,states('O2').X,color='black')
plt.xlabel('Simulation Time[ms]')
plt.ylabel('O2 mole fraction')
plt.subplot(2,2,4)
plt.plot(states.t,states('CH4').X,color='black')
plt.xlabel('Simulation Time[ms]')
plt.ylabel('CH4 mole fraction')
plt.tight_layout()

plt.figure(2)
plt.barh(reactions,sensitivity,color='black')
plt.xlabel('sensiivity value')
plt.title('10 Most sensitivity reactions')
plt.show()

Results:

In the SA with respective to temperature the reaction CH3 + O2 `iff`CH3+ O has the most positive sensitivity value reaction with respective to tempearutre which mean that increasing the rate of this reaction or any positive value we would be seeing a decreases in the ignition delay.

On the other hand , the reaction CH3+H2O2 `iff`CH4+HO2 has the most negative sensitive value , which in thw other words means , ignition delay is high coz of low temp

The first plot shows the variation of temp with respective to time which can be used to find the combustion start and calculate the ignition delay time . And the other plots shows the change in the concentration of the species OH, O2, CH4