Sensitivity Analysis for Combustion of Methane

The aim of this project is to determine the 10 most sensitive reactions w.r.t temperature among the list of reactions in the GRI mechanism for the combustion.

Objectives:

1. Plot the variation of temperature and molar concentration of species (viz., O2, OH & CH4 )w.r.t to time.
2. Plot the bar graph showing the 10 most sensitive reactions w.r.t temperature.

Statement:

Reactor	Constant Pressure Chamber
Fuel	Methane
Initial Temperature	1500 K
Initial Pressure	1 atm.
Mixture Species & Molar Concentration	{\'CH4\':1, \'O2\':2, \'N2\':7.52}
Simulation Time	2 ms with an integral time step of 5 µs

Methodology:

• A Parametric code is written, where a user inserts a number of maximum sensitive reactions w.r.t to a parameter (viz., \'temperature or species \'OH\' or any species from the list) to be considered.
• A gas object is created from the instance of Solution Class and initial conditions are provided
• An ideal gas Constant Pressure reactor is created by passing the gas object.
• Then the simulation network is created by passing the reactor.
• Arrays are created to store the sensitivity of all reactions and thermodynamic results at every time step.
• The reactor network is initialized and results are stored after every integral time step interval for all the reactions.
• Data is sorted to select the 10 most sensitive reactions
• Data Visualization.

\'\'\'
constant pressure, adiabatic kinetic simulation with senstivity analysis

\'\'\'

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

gas = ct.Solution(\'gri30.xml\')    # creating gas instance
T = 1500                          # initial temperature
P = ct.one_atm                    # initial pressure

gas.TPX = T,P,{\'CH4\':1, \'O2\':2, \'N2\':7.52}        #initial state of gas

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

#List of Reactions
reactions = gas.reaction_equations()              #listing all the equations for gas instance of Solution Class of gri30.xml
n_reactions = gas.n_reactions                     # recording total no. of reactions involved (i.e., 325)

#Number of most sensitive reactions that we want to calculate
n=10                                              # reactions need to be considered
parameter = \'temperature\'                         # variable w.r.t either temperature or molar concentration of species

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

for i in range (n_reactions):
	r.add_sensitivity_reaction(i)                 # record sensitivity of each reaction

# relative and absolute tolerance for ignition delay calculation
sim.rtol = 1.0e-6                                   
sim.atol = 1.0e-15

# relative and absolute tolerance for sensitivity calculation
sim.rtol_sensitivity = 1.0e-6
sim.atol_sensitivity = 1.0e-6

states = ct.SolutionArray(gas, extra =\'t\') # initiializing the state object to store thremodynamic properties

for t in np.arange(0,2e-3,5e-6):   # t range with integral time interval
	sim.advance(t)                 # determining the solution at time \'t\'
	states.append(r.thermo.state, t=1000*t) # appends the results at every time step

	for j in range (n_reactions):   #initiating loop to check senstivity for all reactions
		s = sim.sensitivity(parameter,j) # calculating sensitivity wrt a parameter (temperature or molar concentration of species) for every reaction

		if np.abs(s)>np.abs(sensitivity[j]): #storing the most sensitive value through out the calculation for each reaction
			sensitivity[j] = s


reaction_number = np.arange(1,n_reactions+1)  # creating an array for reaction numbers
abs_sensitivity = np.abs(sensitivity)         # calculating the absolute value of sensitivity 

#first combining the list of two columns with same indices and sorting wrt to first column in descending order.
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 first \'n\' prioritized reactions, number and sensitivity

n_reactions = reactions[0:n]
n_reaction_number = reaction_number[0:n]
n_sensitivity = sensitivity [0:n]


#printing the report of first n reactions with number, reaction and sesitivity
print(\'The %d most sensitive reactions are as follows: \' %n)

for k in range (n):
	log_message = \'Reaction number = %5d || Reaction = %30s  || Senstivity = %10.6f\' %(n_reaction_number[k], n_reactions[k],n_sensitivity[k])

	print(log_message)

# reversing to plot, else maximum sensitive eqn is plotted at the bottom on bar graph
n_reactions.reverse()
n_reaction_number.reverse()
n_sensitivity.reverse()

plt.figure(1)
plt.subplot(2,2,1)
plt.plot(states.t,states.T, color= \'r\')
plt.title(\'Temp. Vs Time\')
plt.xlabel(\'Time (ms)\')
plt.ylabel(\'Temperature [K]\')

plt.subplot(2,2,2)
plt.plot(states.t,states(\'O2\').X)
plt.title(\'[O2] Vs Time\')
plt.xlabel(\'Time (ms)\')
plt.ylabel(\'[O2]\')

plt.subplot(2,2,3)
plt.plot(states.t,states(\'OH\').X)
plt.title(\'[OH] Vs Time\')
plt.xlabel(\'Time (ms)\')
plt.ylabel(\'[OH]\')

plt.subplot(2,2,4)
plt.plot(states.t,states(\'CH4\').X)
plt.title(\'[CH4] Vs Time\')
plt.xlabel(\'Time (ms)\')
plt.ylabel(\'[CH4]\')
plt.tight_layout()

plt.figure(2)
plt.barh(n_reactions, n_sensitivity, color = \'k\')
plt.title(\'{0} most sensitive reactions w.r.t [{1}] for Methane Combustion\'.format(n,parameter))

plt.xlabel(\'Senstivity\')
plt.ylabel(\'Reactions\')
plt.tight_layout()
plt.show()

Results:

Discussion:

The above bar chart represents the 10 most sensitive reactions to temperature. The reactions with +ve sensitivity imply upon increasing the reaction rate of that particular reaction increases the temperature thereby decreasing the ignition delay (i.e., ignites faster). The one with -ve sensitivity implies upon increasing the reaction rate of that particular reaction decreases the temperature thereby increasing the ignition delay (i.e., ignites slower).