SENSITIVITY ANALYSIS IN THE STOICHIOMETRIC COMBUSTION OF METHANE WITH RESPECT TO TEMPERATURE

                                          SENSITIVITY ANALYSIS IN THE STOICHIOMETRIC COMBUSTION OF METHANE CONCERNING TEMPERATURE

l. OBJECTIVE

1. Write a code to calculate the 10 most temperature-sensitive reactions in the GRI-Mech 3.0 mechanism file.

2. Plot the variation of temperature and mole fraction of different species concerning time.

3. Plot a bar graph to show the 10 most temperature-sensitive reactions in the GRI-Mech 3.0 mechanism file.

ll. INTRODUCTION

There are several mechanisms involved during the combustion of any gas. These mechanisms together formulate the final states of the gas. Every reaction mechanism has its own sensitivity, which can determine the output parameters of the combustion process.

Sensitivity analysis of any reaction helps us to determine the dependency of parameters on the reaction rate of a particular reaction mechanism. This can help us determine whether an increase or decrease in the reaction rate of a particular mechanism has a positive or negative impact on the output parameter.

In this project, we shall be determining the sensitivity of temperature concerning the reaction rates in the GRI-Mech 3.0 mechanism file.

lll. PROBLEM STATEMENT

1. The combustion of Methane takes place in a constant pressure chamber.

2. The methane gas is an ideal gas having T = 1500 K and P = 1 atm.

3. The equation representing the stoichiometric combustion of methane is as follows - 

• `CH_4+2(O_2+3.76N_2)rarrCO_2+2H_2O+7.52N_2`

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

lV METHODOLOGY

1. A 'gas ' object is created as an instance of the 'Solution' class based on the GRI-Mech 3.0 mechanism file

2. The temperature, pressure, and the number of moles of the reactants are defined.

3. A constant pressure reactor object 'r' from the "gas"  instance is created and is placed in a rector network "sim".

4. All the reactions involved in the process are stored and the number of reactions whose sensitivity is determined is defined.

5. A certain number of most sensitive reactions is predetermined. This number is specified by the user. In this case, we shall consider the 10 most sensitive reactions.

6. A null matrix is initialized to store the sensitivity of all the reactions.

7. All the relations whose sensitivity needs to be determined are added to the reactor.

8. The relative and absolute tolerance of the solution and the sensitivity coefficients are defined.

9. The state object is initialized to store the properties of the gas at different instances of time.

10. Time integration is performed on the reactor network and the states of the gas are stored. The sensitivity of the reactions is computed w.r.t the temperature and is stored in an array.

11. Sort the reaction number, reactions, and sensitivity based on decreasing order of absolute sensitivity.

12. Store the reaction numbers, reactions, and sensitivity of the 'n' most sensitive reactions to display and plot.

13. Printing the 'n' most sensitive reactions along with their reaction number and sensitivity.

14. The graph showing the variation of temperature and mole fractions of different species involved in the combustion process w.r.t. time is plotted. A bar graph is also plotted showing the 'n' most sensitive reactions in descending order.

V. CODE

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

# Creating a gas object using instance of solution class having gri mechanism
gas = ct.Solution('gri30.xml')

# Defining the Temperature, Pressure and Mole fractions
T = 1500 # K
P = ct.one_atm
gas.TPX = T, P, {'CH4': 1, 'O2': 2, 'N2': 7.52}

# Setting up the reactor and the Reactor Network

r = ct.IdealGasConstPressureReactor(gas, name = 'R1')
sim = ct.ReactorNet([r])

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

# User Input - The total number of most sensitive reactions to be calculated 
n = 10

# Initialising a Null Matrix to store sensitivity
sensitivity = [0]* n_reactions

# Adding the reactions to the sensitivity class to calculate the corresponding sensitivity

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

# Set the tolerances of the solution 
sim.rtol = 1e-6
sim.atol = 1e-15

# Set the tolerance of the sensitivity coefficients
sim.rtol_sensitivity = 1e-6
sim.atol_sensitivity = 1e-6


# Initializing a 'States' object to store the properties of the gas at different instances of time
states = ct.SolutionArray(gas, extra = 't')

# Time Integration on the reactor Network and storing the properties at each state
for t in np.arange(0, 2e-3, 5e-6):
	sim.advance(t)
	states.append(r.thermo.state, t= 1000*t)

	# Calculating the sensitivity and maximum sensitivity with respect to Temperature	

	for j in range(n_reactions):
		s = sim.sensitivity('temperature', j)

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


# Sorting the maximum sensitivity along with their reactions 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 reaction 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 = %5f' % (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.plot(states.t, states.T)
plt.title('Variation of Temperature w.r.t. Time')
plt.xlabel('Time [ms]')
plt.ylabel('Temperature [K]')
plt.grid(which='major')
plt.minorticks_on()
plt.grid(which='minor', alpha=0.2)
plt.tight_layout()

plt.figure(2)
plt.subplot(2, 2, 1)
plt.plot(states.t, states('O2').X)
plt.title('Variation of Mole Fraction of O2 w.r.t. Time')
plt.xlabel('Time [ms]')
plt.ylabel('O2 Mole Fraction')
plt.grid(which='major')
plt.minorticks_on()
plt.grid(which='minor', alpha=0.2)

plt.subplot(2, 2, 2)
plt.plot(states.t, states('OH').X)
plt.title('Variation of Mole Fraction of OH w.r.t. Time')
plt.xlabel('Time [ms]')
plt.ylabel('OH Mole Fraction')
plt.grid(which='major')
plt.minorticks_on()
plt.grid(which='minor', alpha=0.2)

plt.subplot(2, 2, 3)
plt.plot(states.t, states('H').X)
plt.title('Variation of Mole Fraction of H w.r.t. Time')
plt.xlabel('Time [ms]')
plt.ylabel('H Mole Fraction')
plt.grid(which='major')
plt.minorticks_on()
plt.grid(which='minor', alpha=0.2)

plt.subplot(2, 2, 4)
plt.plot(states.t, states('CH4').X)
plt.title('Variation of Mole Fraction of CH4 w.r.t. Time')
plt.xlabel('Time [ms]')
plt.ylabel('CH4 Mole Fraction')
plt.grid(which='major')
plt.minorticks_on()
plt.grid(which='minor', alpha=0.2)
plt.tight_layout()

plt.figure(3)
plt.subplot(2, 2, 1)
plt.plot(states.t, states('O').X)
plt.title('Variation of Mole Fraction of O w.r.t. Time')
plt.xlabel('Time [ms]')
plt.ylabel('O Mole Fraction')
plt.grid(which='major')
plt.minorticks_on()
plt.grid(which='minor', alpha=0.2)

plt.subplot(2, 2, 2)
plt.plot(states.t, states('H2').X)
plt.title('Variation of Mole Fraction of H2 w.r.t. Time')
plt.xlabel('Time [ms]')
plt.ylabel('H2 Mole Fraction')
plt.grid(which='major')
plt.minorticks_on()
plt.grid(which='minor', alpha=0.2)

plt.subplot(2, 2, 3)
plt.plot(states.t, states('CO2').X)
plt.title('Variation of Mole Fraction of CO2 w.r.t. Time')
plt.xlabel('Time [ms]')
plt.ylabel('CO2 Mole Fraction')
plt.grid(which='major')
plt.minorticks_on()
plt.grid(which='minor', alpha=0.2)

plt.subplot(2, 2, 4)
plt.plot(states.t, states('N2').X)
plt.title('Variation of Mole Fraction of N2 w.r.t. Time')
plt.xlabel('Time [ms]')
plt.ylabel('N2 Mole Fraction')
plt.grid(which='major')
plt.minorticks_on()
plt.grid(which='minor', alpha=0.2)
plt.tight_layout()

plt.figure(4)
plt.barh(max_reactions, max_sensitivity, color = 'maroon')
plt.title('{0} Most Sensitive Reactions  w.r.t. Temperature'. 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()

Vl. OUTPUTS

The 10 most sensitive reactions are as follows -
Reaction Number =   155 | Reaction =          CH3 + O2 <=> CH3O + O | Sensitivity = 16.834746
Reaction Number =   158 | Reaction =       2 CH3 (+M) <=> C2H6 (+M) | Sensitivity = -13.607978
Reaction Number =    38 | Reaction =              H + O2 <=> O + OH | Sensitivity = 13.093667
Reaction Number =    53 | Reaction =           CH4 + H <=> CH3 + H2 | Sensitivity = -10.149532
Reaction Number =   156 | Reaction =         CH3 + O2 <=> CH2O + OH | Sensitivity = 6.774795
Reaction Number =   119 | Reaction =        CH3 + HO2 <=> CH3O + OH | Sensitivity = 6.265732
Reaction Number =   159 | Reaction =             2 CH3 <=> C2H5 + H | Sensitivity = 4.341987
Reaction Number =   161 | Reaction =       CH2O + CH3 <=> CH4 + HCO | Sensitivity = 4.025934
Reaction Number =    32 | Reaction =        CH2O + O2 <=> HCO + HO2 | Sensitivity = 4.023221
Reaction Number =   118 | Reaction =         CH3 + HO2 <=> CH4 + O2 | Sensitivity = 3.434430
[Finished in 33.5s]

Figure 1 -

Figure 2 -

Figure 3 -

Figure 4 -

Vll. RESULTS

1. Figure 1 shows that there is a sudden increase in the temperature of the system at a particular instance of time called ignition delay.

2. FIgure 2 and 3 shows the mole fraction and behavior of different species at different instances of time. These plots help us to understand the consumption and production of different species in the entire combustion process.

3. Figure 4 shows the 10 most sensitive reactions in GRI-Mech 3.0 mechanism file. We can see that some of the reactions have positive sensitivity and some of the reactions have a negative sensitivity

4. If the reaction has a positive sensitivity, then increasing the reaction rate causes a faster rise in the temperature. Similarly, if the reaction has a negative sensitivity, then increasing the reaction rate causes a slower rise in the temperature.

5. Faster temperature rise is required to decrease the ignition delay. This can be achieved by increasing the reaction rate of the mechanism having positive sensitivity and decreasing the reaction rates of mechanisms having negative sensitivity.

Vlll. CONCLUSION

The 10 most temperature-sensitive reactions in the GRI-Mech 3.0 mechanism file are computed and a bar graph is plotted. Also, the variation of temperature and mole fraction of different species with respect to time is plotted.