Week 9 - Senstivity Analysis Assignment

Aim: To analyze the sensitivity of gri30 mechanism of all reactions. The sensitivity variable is Temperature.

"""
Sensitivity analysis code for Methane with respect to temperature
"""

import numpy as np
import cantera as ct

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

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


# enable sensitivity with respect to the rates for all reactions 
for i in range(gas.n_reactions):
	r.add_sensitivity_reaction(i)

# Setting the tolerance for the solution 
sim.rtol = 1e-6 # relative change in solution should be less than this value (relative => Xold - Xnew/Xold)
sim.atol = 1e-15 # absolute change in solution should be less than this value (abs0lute => Xold -X new)

# setting the tolerance for the sensitivity coefficients
sim.rtol_sensitivity = 1e-6 # relative change should be less than this value
sim.atol_sensitivity = 1e-6 # absolute change should be less than this value

# empty array to store sensitivity for each reaction 
s =np.zeros(gas.n_reactions) 
# empty array to store the maximum sensitivity of each reaction
sens_max = np.zeros(gas.n_reactions)

# empty array to store the reactions  
reaction_list = []
# array to store the id of each reaction
reaction_id = np.linspace(1,325,325).astype(int)

# appending the reactions in the reaction_list array
for i in range(gas.n_reactions):
	reaction_list.append(sim.sensitivity_parameter_name(i))

# to store the state at each time step  
states = ct.SolutionArray(gas, extra = ['t','sens'])

# time loop with time step 5e-6 and total time 2e-3
for t in np.arange(0,2e-3,5e-6):
	# simulating using the advance method
	sim.advance(t)
	# to find the sensitivity with respect to temperature for each reaction
	for j in range(gas.n_reactions):
		s[j] = sim.sensitivity(r.component_index('temperature'),j)
		# to update the maximum value of the sensitivity for each reaction
		if np.abs(sens_max[j]) < np.abs(s[j]):
			sens_max[j] = s[j]

	# appending the state at each time step 
	states.append(r.thermo.state,t = 1000*t,sens = s)
	

# using the sorted zip method 
# to store the reaction in decreasing order
reac_n = [a for b, a in sorted(zip(np.abs(sens_max),reaction_list), reverse = True)]
# to store the sensitivity in decreasing order
s_max = [a for b, a in sorted(zip(np.abs(sens_max),sens_max), reverse = True)]

# top N sensitive reactions 
N=10

# Array to store top N sensitive reactions
max_reactions = reac_n[0:N]
#Array to store top N sensitivities
max_sensitivity = s_max[0:N]


# Loop to print top N sensitive reactions with their sensitivities
print(str(N)+' most sensitive reactions in GRI3.0 mechanism in '+str(gas.n_reactions)+' reactions are as follows:')
for i in range(N):
	print('For reaction ' +str(max_reactions[i])+' sensitivity is '+str(max_sensitivity[i]))


# plotting
import matplotlib.pyplot as plt
plt.figure(1)
plt.subplot(3,2,1)
plt.plot(states.t,states.T)
plt.xlabel('time(ms)')
plt.ylabel('temperature(K)')
plt.subplot(3,2,2)
plt.plot(states.t, states('OH').X)
plt.xlabel('time(ms)')
plt.ylabel('OH mole fraction')
plt.subplot(3,2,3)
plt.plot(states.t, states('H').X)
plt.xlabel('time(ms)')
plt.ylabel('H mole fraction')
plt.subplot(3,2,4)
plt.plot(states.t, states('O2').X)
plt.xlabel('time(ms)')
plt.ylabel('O2 mole fraction')
plt.subplot(3,2,5)
plt.plot(states.t,states('CH3').X)
plt.xlabel('time(ms)')
plt.ylabel('CH3 mole fraction')
plt.subplot(3,2,6)
plt.plot(states.t,states('CH4').X)
plt.xlabel('time(ms)')
plt.ylabel('CH4 mole fraction')
plt.tight_layout()

plt.figure(2)
plt.plot(states.t,states.sens)
plt.xlabel('time(ms)')
plt.ylabel('Temperature sensitivity')


plt.figure(3)
plt.gca().invert_yaxis()
plt.barh(range(N), max_sensitivity)
plt.yticks(range(N),max_reactions)
plt.grid()
plt.ylabel('Chemical Reactions')
plt.xlabel('Temperature Sensitivity')
plt.suptitle('Top '+str(N)+' reactions with higher sensitivities')
plt.tight_layout()

plt.show()

All the reactions of the mechanism are included with the aid of 'gas.n_reactions' this method helped to include all the elementary reactions.

For solving the sensitivity of all the reactions, all the gas object reactions are added with the aid of add_sensitivity_reactions() methods.

As mentioned in the lecture the relative and absolute tolerances are required for chemical reactions sensitivity solution.

Empty arrays are initialized for storing necessary values for plotting and printing.

The states object is created and the state is integrated from 0 to 2e-3 seconds. The dt value is selected by Cantera itself.

Results:

The sensitivity of Temperature for the reaction CH3 + O2 <--> CH3O + O is the highest. If we reduce the temperature the sensitivity will reduce that will increase the ignition delay.

Hence sensitivity analysis is done and sensitivity is achieved for temperature variable.