Week 9 - Senstivity Analysis Assignment

Aim: - To perform a Sensitivity Analysis on GRI mech file.

Objective: -

1. To perform Sensitivity Analysis on GRI mech file reactions with respect to temperature and plot 10 most sensitive reactions
2. To plot the change in temperature and Mole fraction of {CH₄, H, OH} with respect Time.

Theory: -

Sensitivity –

It is the response of a certain mathematical model/system to the change in input parameters. The study of how the uncertainty in the output of the model/system (numerical or otherwise) can be divided and allocated to different sources of uncertainty in its inputs is called Sensitivity Analysis.

 say if:          Z = f(X)    `=>          S_i = (dZ)/dx_i`

where S_i is the sensitivity coefficient of Z with respect to x_i parameter.

Coding Approach: -

1. Instantiate a gas object using Solution in Cantera for GRI mech File.
2. Temperature, pressure and Mole Fraction of mixture is assigned using gas.TPX.
3. An object r is created using gas object in IdealGasConstPressureReactor, which specify the equation to be solved for the species such as species equation and energy equation. There is other reactor as well, which is used in specific conditions prescribed for it such as IdealGasConstantVolume reactor, constant volume reactor, constant pressure reactor etc.
4. A reactor Network is setup using r as input and named sim.
5. Respective Tolerances are set.
6. Attributes of Solution Array from Cantera is grouped as an instance of states with addition of Time in milli seconds.
7. A for loop is set and value are calculated for specific Time period using Advance command.
8. For each time loop sensitivity of each reaction is calculated and max sensitivity is saved is s_max.
9. states are appended for each time loop.
10. Reactions are shorted in descending order according to s_max and both values are stored and plotted as required

Code: -

"""
    progrm to perorm senstivity analysis of first 10 reaction in GRI mech - chemkin
    with respect to temperture
"""

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


temp = 1500
pressure = ct.one_atm

gas = ct.Solution('gri30.xml')
gas.TPX = temp,pressure,{ 'CH4':1 , 'O2':2 , 'N2':7.52}
rxn_name= gas.reaction_equations()
n = gas.n_reactions
r = ct.IdealGasConstPressureReactor(gas)
sim = ct.ReactorNet([r])

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

# set the tol. for the solution and the sensitivity coeff.
sim.rtol = 1.0e-6
sim.atol = 1.0e-15 

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

states = ct.SolutionArray(gas, extra = ['t_ms'])

s = np.zeros(n)
s_max = [0]*n
# loop for time 
for t in np.arange(0,2e-3,5e-6):
 sim.advance(t)
 for j in range(n):
  s= sim.sensitivity('temperature', j)
  # for storing max sensitivity 
  if np.abs(s) >= np.abs(s_max[j]):
   s_max[j] = s 
 states.append(r.thermo.state,t_ms = 1000*t)

# sorting rections base on absolute value of max sensitivity
rxn_name = [r for s_ab, r in sorted(zip(np.abs(s_max),rxn_name),reverse = True)]
s_max = [s for s_ab , s in sorted(zip(np.abs(s_max),s_max),reverse = True)]

#Extracting rections and sensitivity for selected no. of rections to be printed
n_p = 10   # max = 325
rxn =  rxn_name[0:n_p]
st =  s_max[0:n_p]

# Extracting mole fractions of OH, H and CH4
mf_OH = states.X[:, states.species_index('OH')]
mf_H = states.X[:, states.species_index('H')]
mf_CH4 = states.X[:, states.species_index('CH4')]

# printing the top most seected reactions and their sensitivity
print('10 Most Temperature Sensitive Rection are: -')
print('s \t\t\t  s ' %('Reactions','Sensitivity'))
for k in range(n_p):
 print('%s \t  f ' %(rxn[k],st[k]))

# potting the sensitivity bar graph
plt.figure(1)
pos = np.arange(n_p)
plt.barh(pos,st,0.4, align='center', tick_label =rxn,  color=['black','brown'])
plt.gca().invert_yaxis()
plt.yticks(fontsize=6)
plt.xticks([-20,-15,-10,-5,0,5,10,15,20], fontsize=6)
plt.xlabel('Sensitivity')
plt.title('10 Most Temperature Sensitive Rection in GRI Mechanism')
plt.tight_layout()
plt.grid(color='b', ls = '-.', lw = 0.25)


# potting the chnges in temp and moe ractions
ax2=plt.figure(2)
plt.subplot(2,2,1)
plt.suptitle('Change in Parameters with time (in millisecond)')
plt.plot(states.t_ms,states.T)
plt.xlabel('Time (ms)')
plt.ylabel('Temperature (K)')
plt.grid(color='b', ls = '-.', lw = 0.25)

plt.subplot(2,2,2)
plt.plot(states.t_ms,mf_OH)
plt.xlabel('Time (ms)')
plt.ylabel('OH Mole Fraction')
plt.grid(color='b', ls = '-.', lw = 0.25)

plt.subplot(2,2,3)
plt.plot(states.t_ms,mf_H)
plt.xlabel('Time (ms)')
plt.ylabel('H Mole Fraction')
plt.grid(color='b', ls = '-.', lw = 0.25)

plt.subplot(2,2,4)
plt.plot(states.t_ms,mf_CH4)
plt.xlabel('Time (ms)')
plt.ylabel('$CH_4$ Mole Fraction')
plt.grid(color='b', ls = '-.', lw = 0.25)

plt.tight_layout()
plt.show()

Results: - 

Conclusion: -

The most sensitive reaction to temperature is ‘CH3 + O2 <=> CH3O + O’ with sensitivity coefficient 16.834.