Sensitivity Analysis on Combustion of methane using Cantera

INTRODUCTION:

Sensitivity Analysis is a process to get a relationship between input parameters and output of the model. Eg. Effect of a particular species concentration of the adiabatic flame temperature.

In computational combustion, sensitivity has a separate system of ODE which has to be solved for obtaining the relationship. Generally, sensitivity is given by :

`(dZ)/dt = F(Z,t,a)` where Z is the vector of temperature and species whos size is (K+1) where K denotes the number of species. 't' is the time and 'a' is the reaction rate constant.

The sensitivity ODE system has its own tolerance values. The two types of tolerances specified are absolute and relative tolerances. The absolution tolerance is used for the Newton Raphson convergence and the relative tolerance is used for species which have a very slow/fast rate of change of parameters.

A Positive value of sensitivity means the particular parameter is has a positive coefficient effect on the component and vice versa. Eg. A particular reaction might favour temperature of the state. Higher the rate of this reaction, higher the temperature obtained in that state and lesser is the ignition delay.

In this project, the sensitivity analysis of tempertaure is done on the reactions of methane from the GRI reaction mechanism file. The GRI reaction file has a total of 325 reactions listed. Out of it, the top 10 sensitive reactions to temperature is determined in this project. Additionally, the concentration of OH, O2 and CH4 is plotted against the simulation time.

Cantera solves the system of ODE's for the sensitivity analysis. The reactor considered is the Ideal Gas Constant Pressure Reactor.

PYTHON CODE [USING CANTERA] FOR SENSITIVITY ANALYSIS:

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

gas = ct.Solution('gri30.xml')
temp = 1500
pres = ct.one_atm

gas.TPX = temp, pres, 'CH4:1, O2:2, N2:7.52'
r = ct.IdealGasConstPressureReactor(gas, name = 'R1')
sim = ct.ReactorNet([r])


# enable sensitivity wrt rates of 1st 10
#reactions ( reactions 0 thru 9)
no_rxn = gas.n_reactions

for i in range(no_rxn):
	r.add_sensitivity_reaction(i)


#set tolerancces for species loop  [ NR]
sim.rtol = 1e-6
sim.atol = 1e-15
sense = final = []

#set tolerancces for sensitivity loop
sim.rtol_sensitivity  = 1e-6
sim.atol_sensitivity = 1e-6

rxn = []
maxi_1 = []
maxi = [0]*no_rxn
states = ct.SolutionArray(gas, extra = ['t', 's2'])

for t in np.arange(0, 2e-3, 5e-6):	
	sim.advance(t)
	
	for j in range (no_rxn):
		s = sim.sensitivity('temperature',j)	#sensitivity of OH to reaction 2
		states.append(r.thermo.state, t = 1000*t, s2 = s)

		if(np.abs(s) >= np.abs(maxi[j])):
			
			maxi[j] = s

#print(maxi)
maxi_1 = np.copy(maxi)
maxi_1.sort()
left_ptr = 0
right_ptr = len(maxi_1) - 1
ctr = 0

while (ctr < 10):
	if np.abs(maxi_1[left_ptr]) > np.abs(maxi_1[right_ptr]):
		final.append(maxi_1[left_ptr])
		left_ptr += 1
		ctr  = 	ctr + 1

	else:
		final.append(maxi_1[right_ptr])
		right_ptr = right_ptr - 1
		ctr  = ctr + 1
			

ctr = 0
final.reverse()
rxn.reverse()

while (ctr < 10):
	index = maxi.index(final[ctr])
	rxn.append(sim.sensitivity_parameter_name(index))
	ctr += 1

print(final)
	
plt.figure(1)
plt.barh(rxn,final)
plt.xlabel('Temp sensitivity')
plt.ylabel('Reactions')
plt.tight_layout()

plt.figure(2)
plt.subplot(2,2,1)
plt.plot(states.t, states.T)
plt.xlabel('Simulation Time')
plt.ylabel('Temperature')

plt.subplot(2,2,2)
plt.plot(states.t, states('OH').X)
plt.xlabel('Simulation Time')
plt.ylabel('Conc of OH')

plt.subplot(2,2,3)
plt.plot(states.t, states('O2').X)
plt.xlabel('Simulation Time')
plt.ylabel('Conc of O2')

plt.subplot(2,2,4)
plt.plot(states.t, states('CH4').X)
plt.xlabel('Simulation Time')
plt.ylabel('Conc of CH4')

plt.show()	

RESULTS:

A) SENSITIVITY PLOT

B) OTHER PLOTS:

INFERENCE:

1. From the sensitivity plot, it can be clearly seen that the reaction with highest positive sensitivity is `CHunderset(3)() + Ounderset(2)() iff CHunderset(3)()O + O` and the reaction with highest negative sensitivity is `2CHunderset(3)() + M iff Cunderset(2)()Hunderset(6) + M`.
2. This means if the rate of the first reaction is increased, the state temperture increases and hence we can obtain lesser ignition delay. Also, if we increase the rate of the second reaction, the state temperature decreases and hence the ignition delay will increase.
3. The other eight reactions affect the ignition delay as per their magnitude and their sign.
4. The other plots show the variation of temperature, concentration of OH, O2 and CH4 wrt simulation time from which the ignition delay can be calculated.
5. We can see that ignition takes place and there is a spike in temperature [ (T_init +400) K] . The ignition delay can be calculated to be around 1.12 ms.
6. During ignition, we can see that there is a spike in concentration of OH, decrease in concentration of O2 and CH4. This is because OH is the product of combustion is formed during combustion and O2 as well as CH4 are consumed during combustion.