Week 9 - Senstivity Analysis Assignment

Aim: Sensitivity analysis of GRI-3.0 mechanism

Objectives:

• Determine 10 most sensitive reactions from the GRI-3.0 mechanism
• Determine sensitivity with respect to temperature
• Output the reactions in decreasing order of sensitivity

Theory:

In any reaction mechanism, multiple intermediate reactions also occurs which affect the rate of the entire mechanism. These reactions depend on the operating conditions and can be controlled to achieve desired state. In order to control the reactions, we can change parameters like temperature, pressure, species, concentration, etc. that have a direct impact on the rate of the particular reaction. This impact is studied by sensitivity analysis. Hence, we get the relationship between input and output parameters of the system. Higher the sensitivity, greater the impact of any change in the paramter to the output. Suppose, `x` is the input parameter and `y` is the output, then if `x` is incremented, and thus becomes `x + deltax` then the output is deviated by `deltay`. Then, the sensitivity is given as:

Sensitivity `= (partialy)/(partialx)`

Hence, based on this sensitivity, we can vary the input to obtain desired output.

Code:

• We need to initialize the gas instance from the gri30 mechanism fileand set initial temperature to 1500K, the pressure to 1atm and species concentration according to stoichiometric coefficients. `CH_4 + 2(O_2 + 3.76N_2) -> CO_2 + 2H_2O + 7.52N_2`.
• Create an ideal gas constant pressure reactor 'r' from the gas instance and a reactor network 'sim' from the 'r' instance. Setting tolerances for the solvers.
• Initialize a solution array to store thermodynamic state of gas at each time step (and hence, initializing a time step).
• Further the reaction is passed through each time step and thermodynamic state variables are appended in solution variable. Sensitivity of each reaction is evaluated and maximum of the result is stored.
• The above data is sorted and top ten values are displayed.

"""
Program to identify top 10 temperature sensitive reactions for the GRI 3.0 mechanism

"""
from numpy import *
from matplotlib.pyplot import *
import cantera as ct

# Initial conditions for the reaction mechanism
gas = ct.Solution('gri30.xml')
T0 = 1500
P0 = ct.one_atm
gas.TPX = T0, P0, {'CH4':1, 'O2':2, 'N2':7.52}

# Initializing the reactor and the reactor network
r = ct.IdealGasConstPressureReactor(gas)
sim = ct.ReactorNet([r])
reactions = gas.reaction_equations()

# Extracting and adding reactions from the mechanism to sensitivity analysis
n = gas.n_reactions
rid = linspace(1,n,n).astype(int)
N = 10  # Number of sensitive reactions to be displayed
for i in range(n):
	r.add_sensitivity_reaction(i)

# Setting the relative and absolute tolerances for the solver
sim.rtol = 1e-06
sim.atol = 1e-15
sim.rtol_sensitivity = 1e-06
sim.atol_sensitivity = 1e-06

# Initializing the maximum temperature sensitivities of each reaction
sT_max = [0]*n

# Initializing a solution array to store the reaction properties at each time-step
states = ct.SolutionArray(gas, extra=['t_in_ms'])

# Creating the time domain for the reactor network
t = linspace(0, 2e-3, 401)

# Outer time loop
for i in range(len(t)):
	sim.advance(t[i])
	states.append(r.thermo.state, t_in_ms=t[i]*1e3)
	# Storing the maximum sensitivities of each reaction
	for j in range(n):
		sT = sim.sensitivity('temperature',j)
		if abs(sT) >= abs(sT_max[j]):
			sT_max[j] = sT

# Sorting the reactions and reaction IDs based on absolute value of sensitivities
reactions = [r for s,r in sorted(zip(np.abs(sT_max),reactions), reverse=True)]
rid = [r for s,r in sorted(zip(np.abs(sT_max),rid), reverse=True)]
sT_max = [sT for abs_sT, sT in sorted(zip(np.abs(sT_max),sT_max), reverse=True)]

# Extracting the top 10 sensitive reactions
rplot = reactions[0:N]
rid_plot = rid[0:N]
splot = sT_max[0:N]

# Printing the results
print("The top %d temperature sensitive reactions are:" % (N))
print("-------------------------------------------------------------")
print("%8s %4s %18s %12s %12s" % ('Reaction ID','|','Reaction','|','Sensitivity'))
print("-------------------------------------------------------------")
for i in range(N):
	print("%7d %8s %26s %4s %10.4f" % (rid_plot[i],'|',rplot[i],'|',splot[i]))
print("-------------------------------------------------------------")

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

# Plotting the results
fig1, ax1 = subplots()
y_pos = arange(len(rplot))
plt.grid('both',linestyle='--',linewidth=1.2)
ax1.set_axisbelow(True)
ax1.barh(y_pos, splot, 0.5, align='center', color='black')
ax1.set_yticks(y_pos)
ax1.set_yticklabels(rplot, rotation=0, fontweight='bold')
ax1.invert_yaxis()
xticks = [-20, -15, -10, -5, 0, 5, 10, 15, 20]
ax1.set_xticklabels(xticks, rotation=0, fontsize = 12)
ax1.set_xlabel('Temperature Sensitivity', fontsize = 18, fontweight='bold')
title('10 most temperature sensitive reactions for the GRI 3.0 mechanism of CH4 combustion', fontsize=20, fontweight = 'bold')
savefig('temperature_sensitivity_ranks.png')
show()

fig2, ax2 = subplots(2,2)
subplot(2,2,1)
plot(states.t_in_ms,states.T,linewidth=2)
grid('both',linestyle='--',linewidth=0.8)
xlabel('Time [ms]',fontweight='bold')
ylabel('Temperature [K]',fontweight='bold')
subplot(2,2,2)
plot(states.t_in_ms,OH,linewidth=2)
grid('both',linestyle='--',linewidth=0.8)
xlabel('Time [ms]',fontweight='bold')
ylabel('OH Mole Fraction',fontweight='bold')
subplot(2,2,3)
plot(states.t_in_ms,H,linewidth=2)
grid('both',linestyle='--',linewidth=0.8)
xlabel('Time [ms]',fontweight='bold')
ylabel('H Mole Fraction',fontweight='bold')
subplot(2,2,4)
plot(states.t_in_ms,CH4,linewidth=2)
grid('both',linestyle='--',linewidth=0.8)
xlabel('Time [ms]',fontweight='bold')
ylabel('CH4 Mole Fraction',fontweight='bold')
tight_layout()
savefig('Sensitivities.png')
show()

Output&Result:

• The positive sensititvity indicates that an increase in temperature would increase the rate of that reaction and negative sensitivity indiactes that an increase in temperature would decrease the rate of that reaction.
• The second figure describes the ignition delay and concentration change of reaction species with respect to time.