Sensitivity Analysis of Stoichiometric Combustion of Methane to find most temperature sensitive reactions using GRI 3 0 Database

Objective: Using the GRI 3.0 database, we are writing a Python/Cantera Code to find the 10 most sensitive reactions from the Stoichiometric Combustion of Methane. We are analyzing the sensitivity of the reactions with respect to the Temperature for the given initial conditions. Performing the code that shows the reaction strings in the order of greatest sensitivity to the least sensitivity and it should be parametric with the number of reactions. Simply, by changing the \'n\' value in the program we need to get \'n\' reactions.

Introduction:

• A chemical mechanism is a theoretical conjecture that tries to describe in detail what takes place at each stage of an overall chemical reaction.
• The detailed steps of a reaction are not observable in most cases.
• In a reaction or in combustion contains thousands of the intermediate reactions. These reactions have their own individual rate constants.
• There are number of observable parameters that a reaction can produce are temperature, pressure, mole fraction, mass fraction, flame speed, and species concentration, etc.
• Sensitivity analysis is one that can increase understanding of the relationships between input and output variables in a system or model.

Sensitivity analysis:

• Sensitivity analysis is the study of how the uncertainty in the output of a mathematical model or system (numerical or otherwise) can be divided and allocated to different sources of uncertainty in its inputs.

• The above plot is a function of \"Y = f(X,p)\".
• Here \'X\' is the input and \'p\' is the independent value that can forward the function in the direction of \'X\'.
• Here we have two sets of plots. It represents how the \'Y\' behaves when \'p\' is changed from \'p\' to \'p+deltap\'.
• When we change the \'p\' value, \'Y\' value affected. The affected value is \'delat y\'.

Sensitivity `S = (deltay)/(deltap)`

• The rate of change of output function due to the change in input parameter is the sensitivity coefficient.
• We can adjust the value of \'p\' in order to obtain the desired value of \'Y\' on the basis of obtained sensitivity value.
• Same as above, in the reaction or in combustion there are some parameters that can actually control the entire mechanism.
• Every reaction in the combustion has sensitivity. Here we are calculating the sensitivity for the Stoichiometric Combustion of Methane with the air and nitrogen.
• Printing ten highly sensitive reactions with respect to the temperature in the manner of greatest sensitivity to the least sensitivity.

Code:

• Taking the gas or combustion reactants from the \'gri30.xml\' file.
• Initial conditions are temperature = 1500k, and pressure p = 1atm
• The reaction taken is CH4+2(O2+3.76N2) <=> CO2+2H2O+7.52N2
• A solver \'IdealGasConstPressureReactor\' is taken in the instance of \'r\'. And \'ReactoNet\' in the instance of \'sim\'.
• Here we maintain two tolerance values.
• Relative tolerance is a limit for the rate of change of variable from each iteration. which is denoted as `(deltax)/x`. It tells us the variable is converged or not.
• In chemical Kinetics, relative tolerance is not enough. So, we need to add absolute tolerance. Absolute tolerance deals with only `deltax`. If the output variable does not exceed the tolerance limit, then it is converged.
• We need to use both the tolerances because radicals which have got very low concentrations. Some species will show a difficult pattern to converge. Here we need to use both tolerances to converge the solution.
• The simulation time we are running is 0 to 2 milliseconds with the time step of 5e-06.
• The \'sim.advance\' is used to advances the solution in time or to perform time integration.
• Initializing the solution array to store the thermodynamic state of \'gas \' at each timestep.
• Then \'sorted\', the sorted() function returns a sorted list from the items in an iterable. The sorted() function sorts the elements of a given iterable in a specific order (either ascending or descending).
• Finally, plotting for the top ten highly sensitive reactions to the temperature. With that also plotting the species and temperature behavior.

\"\"\"
 writing a Python/Cantera Code to find the 10 most sensitive reactions
 In Stoichiometric Combustion of Methane 
 Reaction: CH4+2(O2+3.76N2) <=> CO2+2H2O+7.52N2
\"\"\"
import numpy as np
import cantera as ct
import matplotlib.pyplot as plt

# taking the gas object from solution class of gri30 file
gas = ct.Solution(\'gri30.xml\')

# providing the initial temperature and pressure values
T = 1500
P = 101325

# providing temperature, pressure, and mole fraction values to the solver
gas.TPX = T, P, {\'ch4\':1, \'o2\':2, \'n2\':7.52}

# Governing equations to be solved
r = ct.IdealGasConstPressureReactor(gas, name = \'R1\')
sim = ct.ReactorNet([r])

reac = gas.reaction_equations()

# total number of reactions in the combusion
n_reac = gas.n_reactions
# reaction to be printed
n = 10
S_max = [0]*n_reac


# providing the reactions to be added to sensititvity class for the Sensitivity
for i in range (n_reac):
	r.add_sensitivity_reaction(i)

# relative tolerance for each variable in the matrix
sim.rtol = 1e-6
# absolutr tolerance for each variable in the matrix
sim.atol = 1e-15

# relative tolerance for the sensitivity solution
sim.rtol_sensitivity = 1e-6
# absolute tolerance for the sensitivity solution
sim.atol_sensitivity = 1e-6

# here we are calling the states array to store all the states of the gas properties
states = ct.SolutionArray(gas, extra = \'t\')
# calling total number of reactions
reac_num = np.arange(1,n_reac+1)

# time integration
for t in np.arange(0, 2e-3, 5e-6):
	# defining the advance class, that integrates the function with respect to time
	sim.advance(t)
	states.append(r.thermo.state, t = 1000*t)

	# storong of sensitivity values wrt temperature
	for j in range(n_reac):
		s = sim.sensitivity(\'temperature\', j) # Here j means independent reaction
		# equalising the maximum values to normal, for the printing of first 10 most sensitive reactions
		if np.abs(s)>np.abs(S_max[j]):
			S_max[j] = s


S_abs = np.abs(S_max)
# To get the reactions printed in higher to lower values we need to use the sorted() function
# zip function, compare and prints the relative values in the two alternate files
# a for b, if a and b are two variables of a 2*2 matrix, then it prints whole variables in the matrix as a 1*n matrix
# Source: Stack over flow, github, and w3schools.com
reac = [a for b, a in sorted(zip(S_abs,reac),reverse=True)]
reac_num = [a for b, a in sorted(zip(S_abs,reac_num),reverse=True)]
S_max = [a for b, a in sorted(zip(S_abs,S_max),reverse=True)]

# setting up arrays for the purpouse of plotting
max_r = reac[0:n]
max_n = reac_num[0:n]
max_s = S_max[0:n]

# extracting the behaviour of mole fraction of H, OH, and CH4
H_conc = states(\'H\').X
OH_conc = states(\'OH\').X
CH4_conc = states(\'CH4\').X

plt.figure(1)
plt.grid(\'on\')
plt.barh(max_r,max_s,color=\'blue\')
plt.gca().invert_yaxis()
plt.xlabel(\'Sensitivity\')
plt.ylabel(\'Reactions\')

plt.title(\'Highest Temperature Sensitive Reactions in GRI30 for Stoichiometric Combustion of Methane\')

plt.figure(2)
plt.subplot(2,2,1)
plt.plot(states.t, states.T)	
plt.xlabel(\'Time [ms]\')
plt.ylabel(\'Temperature [k]\')
plt.grid(\'on\')

plt.subplot(2,2,2)
plt.plot(states.t, states(\'oh\').X)	
plt.xlabel(\'Time [ms]\')
plt.ylabel(\'OH mole fraction\')
plt.grid(\'on\')

plt.subplot(2,2,3)
plt.plot(states.t, states(\'h\').X)	
plt.xlabel(\'Time [ms]\')
plt.ylabel(\'H mole fraction\')
plt.grid(\'on\')

plt.subplot(2,2,4)
plt.plot(states.t, states(\'ch4\').X)	
plt.xlabel(\'Time [ms]\')
plt.ylabel(\'CH4 mole fraction\')
plt.tight_layout()
plt.grid(\'on\')

plt.show()

Plots:

• The above plot shows the sensitivities of the top ten reactions with respect to the temperature.
• Here if we observe there are positive and negative sensitivities that are obtained in the reaction.
• The positive sensitivity explains that if we increase the temperature the reaction rate of that particular reactants would be increased.
• The negative sensitivity explains that if we increase the temperature the reaction rate of that particular reactants would be decreased.

• The above figure shows us how the temperature and mole fractions of species have behaved after the reaction.
• The temperature we are given is 1500k and at the end of the reaction it rises up to 2750k.
• If we observe the graph, there is a spike in every plot and it indicated the ignition delay. The ignition delay observed is nearly 1.2 ms.
• The mole fractions of OH and H are getting spiked and stabled after some time, where CH4 species are continuously decreasing and after the ignition, it completely becomes zero.
• At the time between ignition starts and ignition ends, there is a rapid decrease in CH4 species. It means the consumption of carbon atoms is happening and becomes zero after the ignition completes.