SENSITIVITY ANALYSIS TO DETERMINE MOST TEMPERATURE SENSITIVE REACTIONS FOR METHANE COMBUSTION IN THE GRI 3 0 REACTION MECHANISM

INTRODUCTION

1. In any particular reaction, there are multiple intermediate reactions happening which affect the rate of the entire mechanism.
2. These reactions depend on different operating conditions and can be controlled to alter the mechanism to the desired state.
3. In order to control the reactions, we can change different conditions like temperature, pressure, species concentration, etc. that have a direct impact on the rate of the particular reaction.
4. This impact can be studied with the help of a sensitivity analysis. Sensitivity analysis is the process of finding a relationship between the input parameters and the output of a system. The higher the sensitivity, the greater the impact of any change in the parameter to the output.

MATHEMATICAL EXPLANATION OF SENSITIVITY

NOTE: Image from MIT OCW

Consider `Y = f(X,p)`

X = Input parapemet

p = independent parameter

The two plots represent the variation of `Y` when the value of `p` is changed from `p` to `p + deltap`. This change in the parameter affects the value of `Y` by a certain value . Then the sensitivity of with respect to the parameter `p` is defined as:

`Sensitivity = (deltay)/(deltap)`

We can adjust the value of p to obtain the desired output Y. Similarly, for a chemical reaction mechanism, there are several independent parameters that can be used to control the entire mechanism.

PROCEDURE FOR THE COMPUTATIONAL EXPERIMENT

1. Initialize \"gas\" as an instance of the \"Solution\" class based on the mechanism from the \"gri30.xml\" file.
2. Set the initial temperature to 1500 K, the pressure to 1 atm and species concentration according to the stoichiometric reaction `CH_4+2(O_2+3.76N_2 )→CO_2+2H_2 O+7.52N_2`.
3. Create an ideal gas constant pressure reactor \"r\" from the \"gas\" instance and a reactor network \"sim\" from the \"r\" instance.
4. Extract the reactions from \"r\" and add them to the sensitivity analysis.
5. Set the relative and absolute tolerance of the reactor solver to 10^-6 and 10^-15 respectively.
6. Set the relative and absolute tolerance for the sensitivity solver to 10^-6 each.
7. Initialize a solution array to store the thermodynamic state of \"gas\" at each time-step.
8. Create a time-domain from 0 to 2 ms with a time-step of 0.005 ms.
9. Further the reaction through each time-step by using the \"self.advance(t)\" function for the reactor network and append the thermodynamic state variables in the solution array.
10. Calculate the sensitivity of each reaction with respect to temperature using the \"self.sensitivity(\'temperature\', reactionID)\" and store only the maximum sensitivity values for each reaction.
11. Sort the reactions and reaction IDs based on the absolute value of their temperature sensitivity and extract the top 10 values from each.
12. Print and plot the results.

PYTHON CODE

\"\"\"
Program to identify top 10 temperature sensitive reactions for the GRI 3.0 mechanism
Program by: Shouvik Bandopadhyay
\"\"\"

import numpy as np
import matplotlib.pyplot as plt
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 = np.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 = np.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 np.abs(sT) >= np.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 = plt.subplots()
y_pos = np.arange(len(rplot))
plt.grid(\'both\',linestyle=\'--\',linewidth=1.2)
ax1.set_axisbelow(True)
ax1.barh(y_pos, splot, 0.5, align=\'center\', color=\'green\')
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\')
plt.title(\'10 most temperature sensitive reactions for the GRI 3.0 mechanism of CH4 combustion\', fontsize=20, fontweight = \'bold\')
figManager = plt.get_current_fig_manager()
figManager.window.state(\'zoomed\')

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

RESULTS

1. TOP 10 MOST SENSITIVE REACTIONS

2. TEMPERATURE SENSITIVITY BAR PLOT

Positive Sensitivity: Increase in temperature denotes increase in the rate of reaction.

Negative Sensitivity: Increase in temperature denotes decrease in the rate of reaction.

3. TEMPORAL VARIATION OF PARAMETERS

• Spike in the values in each of the above graph indicates ignition delay of the reaction mechanism.
• Mole fraction of the products increase after auto-ignition while mole fraction of the reactants decrease after auto ignition till the mole fraction becomes 0.

CONCLUSION

Sensitivity Analysis has been performed to successfully find the 10 most sensitive reactions in methane combustion.