Sensitivity Analysis for Combustion of Methane from GRI 3.0 mechanism file using Python & Cantera

Aim:

To write a code from GRI mechanism for the combustion of Methane and to calculate 10 most sensitive reactions with respect to Temperature

Objectives:

1. Plot the variation of Temperature and Molar concentration of different species w.r.t. time

2. Plot a bar graph showing most Sensitive reactions w.r.t. temperature

Problem Statement:

The reaction takes place in constant pressure Chamber

The Methane Gas is an Ideal gas

Initial Conditions are: T = 1500 K and P = 1 atm

The combustion of Methane is represented as:

`CH_4 + 2O_2 + 7.52N_2 =>`Products of combustion

The simulation is run for 2 milliseconds having a time step of 0.005 milliseconds

Methodology:

1. An object gas is created from Solution class using the mechanism in GRI 3.0 mechanism file and Initial conditions for Pressure, Temperature and Mole Fraction are provided

2. This gas object is inserted in a constant pressure Ideal reactor

3. The reactor is then inserted in a reactor network to solve for all species of combustion at an instance

4. A null matrix for Sensitivity is defined to store the Sensitivity for all the reactions

5. State array is defined to store the thermodynamic state of gas at every time instant

6. Relative & Absolute tolerances for molar concentration and Sensitivity coefficients are defined

7. Reactor network is initialized and after each time step, sensitivities of all reactions are stored, which are compared to the sensitivities from the next consecutive time step and the maximum of the two is stored in Sensitivity matrix

8. These sensitivities are then sorted from maximum to minimum along with their corresponding Reaction equations and Reaction number; the results are then plotted

9. The code is made parametric where the user input determines the number of maximum sensitivities reactions to be calculated by the program

Code:

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

gas = ct.Solution(\'gri30.xml\')
temp = 1500 #K
Pres = ct.one_atm

gas.TPX = temp,Pres,{\'CH4\':1,\'O2\':2,\'N2\':7.52}

#Setting up the reactor and reactor network
r = ct.IdealGasConstPressureReactor(gas)
sim = ct.ReactorNet([r])

#List of reactions
reactions = gas.reaction_equations()
n_reactions = gas.n_reactions

#Number of most sensitive reactions that we want to calculate
n = 10

#Initialising null matrix to store sensitivity
sensitivity = [0]*n_reactions

#adding the reactions to sensitivity class to calculate corresponding sensitivity
for i in range (n_reactions):
	r.add_sensitivity_reaction(i)

#setting tolerances for molar concentration
sim.rtol = 1e-6
sim.atol = 1e-15

#setting tolerances for sensitivity coefficients
sim.rtol_sensitivity = 1e-6
sim.atol_sensitivity = 1e-6

#Initialising a states object to store thermodynamic properties of gas at every instance of time
states = ct.SolutionArray(gas,extra = \'t\')

#Time integration of reactor network and storing properties at each instance
for t in np.arange(0,2e-3,5e-6): # means t starts from 0 and ends at 5e-6 with time step size of 2e-3
	sim.advance(t)
	states.append(r.thermo.state, t= 1000*t)

	#Calculating the sensitivity and Maximum sensitivity wrt temp
	for j in range(n_reactions):
		s=sim.sensitivity(\'temperature\',j)  # Storing sensitivity as \'s\' for all reactions (n_reactions) at each time step

		if np.abs(s)>np.abs(sensitivity[j]): # Comparing sensitivity of each reaction obtained in new time step with previous time step and storing it in sensitivity matrix
			sensitivity[j] = s


#Sorting maximum sensitivities along with their reaction and reaction numbers

reaction_number = np.arange(1,n_reactions+1)
abs_sensitivity = np.abs(sensitivity)

reactions = [a for b, a in sorted (zip(abs_sensitivity,reactions), reverse = True)]
reaction_number = [a for b, a in sorted (zip(abs_sensitivity,reaction_number), reverse = True)]
sensitivity = [a for b, a in sorted (zip(abs_sensitivity,sensitivity), reverse = True)]

#Storing the n maximum sensitivities along with  their reactions and reaction numbers
max_reactions = reactions[0:n]
max_reaction_number = reaction_number[0:n]
max_sensitivity = sensitivity[0:n]

# Printing the n most sensitive reactions
print(\'The %d most sensitive reactions are as follows: \' % n)

for k in range(n):
	log_message = \'Reaction number = %5d | Reaction = %30s | Sensitivity = %0.6f\' %(max_reaction_number[k], max_reactions[k], max_sensitivity[k])
	print(log_message)

#Reversing the array for plotting the results
max_reactions.reverse()
max_reaction_number.reverse()
max_sensitivity.reverse()

#Plotting the results

plt.figure(1)
plt.subplot(2,2,1)
plt.plot(states.t,states.T,\'r\')
plt.title(\'Variation of Temperature w.r.t. time\')
plt.xlabel(\'Time ($ms$)\')
plt.ylabel(\'Temperature($K$)\') # \'$\' is used for italic font
plt.grid(which = \'major\')
plt.minorticks_on() # divides the adjacent values of x and y cordinates by ticking 
plt.grid(which = \'minor\', alpha = 0.2) #making grid for minor variable i.e. time; & alpha for darkness of minor grid lines

plt.subplot(2,2,2)
plt.title(\'Variation of concentration of $O2$ n w.r.t time\')
plt.plot(states.t,states.X[:, gas.species_index(\'O2\')])
plt.xlabel(\'Time ($ms$)\')
plt.ylabel(\'$O_2$ concentration\') # \'_\' is used for subscript
plt.grid(which = \'major\')
plt.minorticks_on() 
plt.grid(which = \'minor\', alpha = 0.2)

plt.subplot(2,2,3)
plt.title(\'Variation of concentration of $OH$ n w.r.t time\')
plt.plot(states.t,states.X[:, gas.species_index(\'OH\')])
plt.xlabel(\'Time ($ms$)\')
plt.ylabel(\'$OH$ concentration\')
plt.grid(which = \'major\')
plt.minorticks_on() 
plt.grid(which = \'minor\', alpha = 0.2)

plt.subplot(2,2,4)
plt.title(\'Variation of concentration of $CH4$ n w.r.t time\')
plt.plot(states.t,states.X[:, gas.species_index(\'CH4\')])
plt.xlabel(\'Time ($ms$)\')
plt.ylabel(\'$CH_4$ concentration\')
plt.grid(which = \'major\')
plt.minorticks_on() 
plt.grid(which = \'minor\', alpha = 0.2)
plt.tight_layout() # prevents overlapping of contents within a plot

plt.figure(2)
plt.barh(max_reactions,max_sensitivity, color = \'blue\')
plt.title(\'{0} most sensitive reactions w.r.t. Temperature in GRI 3.0 for Methane combustion\'.format(n))
plt.xlabel(\'Sensitivity\')
plt.ylabel(\'Reactions\')
plt.grid(which = \'major\')
plt.minorticks_on() 
plt.grid(which = \'minor\', alpha = 0.2)
plt.tight_layout()

plt.show()

Results:

10 most sensitive reactions with respect to temperature

Conclusions:

The sensitivity indicates the influence of the reaction on the output parameter such as Temperature or molar fraction of species. In the bar graph, the positive sensitivity indicates that the increase in sensitivity will increase the rate of increase of the temperature of combustion and hence will decrease the ignition delay, i.e. the fuel mixture ignites faster if the positive sensitivities reaction are used

On the other hand, if we want to increase the ignition delay time, then the sensitivities of reaction with negative sensitivity behaviour should be used

For example: By increasing the Sensitivities of second and fourth reactions (from the bar graph), the ignition delay will be increased. While increasing the sensitivities of rest of the reactions will lead to a decrease in ignition delay