Week 9 - Senstivity Analysis Assignment

OBJECTIVE:

• To write a code that takes all the reactions from the GRI mechanism and calculates 10 most sensitive reactions

• To plot a graph of sensitivities vs the reactions.

THEORY:

SENSITIVITY:

Sensitivity analysis is a process to get relationship between the parameters and the output of a model.

FOR EXAMPLE: Effect of reaction rate on the temperature or species concentration.

A generic ODE system is represented as,

dZdt=F(Z,t,a)dZdt=F(Z,t,a)

Z=(T,Y1,Y2,….YK)TZ=(T,Y1,Y2,….YK)T

Where, Z is a vector of species and temperature

First order sensitivity,

Wi,j=∂Zj∂aiWi,j=∂Zj∂ai

Where, j is the independent variable.

i is the reaction.

a is the net reaction rate multiplier (perturbation factor)

CODE APPROACH:

• we need to initialize the gas instance from the gri30 mechanism file and set the initial temperature to 1500k and pressure 1 atm and species concentration for CH4, O2 and N2 at stoichiometric condition.
• The ideal gas constant pressure reactor ‘r’ from the gas instance and a rectornetwork ‘sim’ from the ‘r’ instance.
• The tolerance are given for the solution to converge.
• A solution array is initialised to store thermodynamic state of gas at each time step.
• The thermodynamic state variables are appended to solution variable
• Sensitivity of each reaction is evaluated and the results are stored.
• The top ten high sensitive values of the gri mechanism reactions are sorted and displayed.

CODE:

"""
	program to calculate sensitivity
"""

#importing modules
import matplotlib.pyplot as plt
import cantera as ct
import numpy as np

gas = ct.Solution('gri30.cti')
#applying the initial variables
T = 1500
P = ct.one_atm
gas.TPX = T,P,{'CH4':1,'O2':2,'N2':7.52}

r = ct.IdealGasConstPressureReactor(gas)
sim = ct.ReactorNet([r])
#extracting the reactions equations and their numbers
reactions = gas.reaction_equations()
n = gas.n_reactions
#number of reactions to be ised for sensitive calculation
N = 10
#adding the reactions from mechanism to sensitivity analysis
for i in range(n):
	r.add_sensitivity_reaction(i)

#relative and absolute tolerance
sim.rtol = 1e-6
sim.atol = 1e-15

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

#number of sensitive reactions to be displayed
sT_max = [0]*n

#solution array that stores the reaction at each time step
states = ct.SolutionArray(gas,extra='time')
reac_num = np.arange(1,n+1)

#outer time loop
for t in np.arange(0,2e-3,5e-6):
	sim.advance(t)
	states.append(r.thermo.state,time = t*1e3)
	for j in range(n):
		sT = sim.sensitivity('temperature',j)
		if abs(sT)>=abs(sT_max[j]):
			sT_max[j] = sT
s_abs = np.abs(sT_max)

#sorting the reactions
reactions = [a for b, a in sorted(zip(s_abs,reactions),reverse = True)]
reac_num = [a for b, a in sorted(zip(s_abs,reac_num),reverse = True)]
sT_max = [a for b, a in sorted(zip(s_abs,sT_max),reverse = True)]

#assigning the 10 sensitive values
max_r = reactions[0:N]
max_n = reac_num[0:N]
max_s = sT_max[0:N]

#table creation
print("The Top %d Most Temperature Sensitive Reactions :"%(N))
print("------------------------------------------------------------------")
print("%1s %8s %4s %18s %12s %12s %4s"%('|','RACTION ID','|','REACTION','|','SENSITIVITY','|'))
print("------------------------------------------------------------------")
for i in range(N):
	print("%7d %10s %26s %4s %10.4f"%(max_n[i],'|',max_r[i],'|',sT_max[i]))
print("------------------------------------------------------------------")

#plotting
#time vs temperature
plt.subplot(2,2,1)
plt.plot(states.time,states.T)
plt.xlabel('time in ms')
plt.ylabel('temperature k')

#time vs OH mole fraction
plt.subplot(2,2,2)
plt.plot(states.time,states('OH').X)
plt.xlabel('time in ms')
plt.ylabel('OH mole fraction')

#time vs H mole fraction
plt.subplot(2,2,3)
plt.plot(states.time,states('H').X)
plt.xlabel('time in ms')
plt.ylabel('H mole fraction')

#time vs CH4 mole fraction
plt.subplot(2,2,4)
plt.plot(states.time,states('CH4').X)
plt.xlabel('time in ms')
plt.ylabel('CH4 mole fraction')
plt.tight_layout()

#bar graph showing 10 sensitive reactions
plt.figure(2)
plt.barh(max_r,max_s,0.4,color='brown')
plt.xlabel('sensitivity')
plt.ylabel('Reactions')
plt.title('Top 10 Temperature Sensitive Reactions of grimech 3.0')
plt.gca().invert_yaxis()
plt.grid()
plt.tight_layout()
plt.show()

RESULTS:

The output shows the top 10 temperature sensitive reactions:

BAR GRAPH:

From the bar graph we can see that the positive sensitivity indicates that with rise in reaction parameters i.e, temperature the rate of reaction will increase and leads to faster combustion. This will decrease the ignition delay.

Whereas the negative sensitivity indicates that with rise in reaction parameters i.e, temperature the rate of reaction will decrease and leads to slower combustion. This will increase the ignition delay.

The above plot shows the variation of temperature and molar concentration of OH, H, CH4 with respect to time.