Sensitivity analysis of the GRI 3.0 mechanism to find the most temperature-sensitive reactions in the methane combustion mechanism

INTRODUCTION:

In any particular reaction mechanism, there are multiple intermediate reactions happening which affect the rate of the entire mechanism. These reactions depend on different operating conditions and can be controlled to alter the mechanism to the desired state. 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. 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.

Consider the following figure:

The figure represents a function `Y = f(X, p)` where `X` is input and `p` is an independent parameter. The two plots represent the variation of `Y` when the value of `p` is changed from `p` to `p + delta p`. This change in the parameter affects the value of `Y` by a certain value `delta y`. Then the sensitivity of `Y` with respect to the parameter `p` is defined as:

Sensitivity`={delta y}/{delta p}`

Based on this sensitivity value, we can adjust the value of `p` in order to obtain the desired value of `Y`. Similarly, for a chemical reaction mechanism, there are several independent parameters that can be used to control the entire mechanism. In this project, we will study the sensitivity of the mechanism to temperature for each reaction in the GRI 3.0 mechanism and identify the `10` most temperature-sensitive reactions.

IMPLEMENTATION IN PYTHON USING CANTERA MODULE:

The following steps are taken to implement the sensitivity analysis in Python using Cantera module:

• Initialize "gas" as an instance of the "Solution" class based on the mechanism from the "gri30.xml" file.
• 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) rarr 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.
• Extract the reactions from "r" and add them to the sensitivity analysis.
• Set the relative and absolute tolerance of the reactor solver to `10^{-6}` and `10^{-15}` respectively.
• Set the relative and absolute tolerance for the sensitivity solver to `10^{-6}` each.
• Initialize a solution array to store the thermodynamic state of "gas" at each time-step.
• Create a time-domain from `0` to `2` `ms` with a time-step of `0.005` `ms`.
• 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.
• 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.
• Sort the reactions and reaction IDs based on the absolute value of their temperature sensitivity and extract the top `10` values from each.
• Print and plot the results.

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


Program by: Shrey Shah
Date: 14-Oct-19
"""
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='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')
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()

The list of the top `10` temperature-sensitive reactions for the GRI 3.0 mechanism are shown below:

The plot below shows the bar graph of the sensitivities of the above `10` reactions. A positive sensitivity indicates that an increase in the temperature would increase the rate of that reaction and a negative sensitivity indicates that an increase in temperature would decrease the rate of that reaction.

The figure below shows the variation of temperature, `OH` mole fraction, `H` mole fraction and `CH_4` mole fraction as the mechanism progresses. The spike in each of their values indicates the ignition delay of the mechanism. We can see that the mole fractions of the products `(OH` and `H)` increase after the auto-ignition and the mole fraction of the reactant `(CH_4)` decreases after the auto-ignition till it reaches a mole fraction of `0`.