Sensitivity Analysis of CH4 air combustion   (Cantera, Skill Lync, Hard)

Aim:

To carry out sensitivity analysis and find 10 most sensitive reactions affecting temperature of methane-air combustion.

Introduction:

The real life combustion processes have multiple intermediate reactions and species, which increase the total number of reactions in mechanism. Simulating combustion using these detailed mechanisms required high computational power. So, to reduce the size of mechanism only those reactions are taken in to account which highly affect the parameter under study. To find these reactions, we carry out sensitivity analysis.

Sensitivity of a reaction means how much does a small change in reaction rate affect the overall state parameter, like temperature. It can be understood by an example in figure below:

A small change in reaction A led to a big change in y, whereas a small change in reaction B led to a small change. This means that y variable is more sensitive to reaction A than reaction B. So, while reducing the reaction file, the reactions like B can be neglected.

Sensitivity is mathematically defined as:

`S = (dy)/(dp)`

In this report, the sensitivity is calculated using cantera's inbuilt function which returns sensitivity of given parameter with respect to a reaction. Cantera uses normalized sensiivity which it defines as:

Data relating to sensitivity and cantera code:

1) Initial temperature = 1500 K

2) Initial pressure = 1 atm

3) Stoichiometric mixture of CH4 and air.

4) Absolute Tolerance = 1e-16 (Tolerance between two iteration in a paticular timestep, used by cantera's solver while simulating)

5) Relative Tolerance = 1e-6 (Tolerance between two timestep solutions)

6) Absolute tolerance of sensitivity = 1e-6 (Absolute tolerance calcuated for sensitivities)

7) Relative Tolerance of sensitivity = 1e-6 (Relative tolerance calcuated for sensitivities)

How does the code work?

1) Mentioning initial temperature, pressure and mole fractions.

2) Passing 'gas' object to Ideal gas constant pressure reactor, which is passed onto object 'r'

3) Passing object 'r' to reactor network to object 'sim'

4) Saving all reactions equations as strings in list 'rxn'

5) Saving the total number of reactions in mechanism file in 'n_rxn'

6) Asking Cantera to also compute sensitivity of reaction 'i' while running simulation.

7) Defining the absolute and relative tolerances.

8) Defining the timestep and total time for simulation to run for.

9) Advancing the solution till 'time'

10) Saving the thermodynamic state of reactor object 'r'

11) Saving the sensitivity of reaction only if its is maximum than its previous sensitivity value.

Each reaction reaches a maximum sensitivity and then it falls down. The array 'max_s' stores these max values.

12) Creating a zip data variable which contains absolute and real values of sensitivities of each reaction and their index numbers.

This is done in order to sort the data keeping the values of index numbers with them. The sorting variable is kept at first place in each of the member of zip.

13) Creating a tuple of zip in order to access the data. The members of tuple are also a tuple which contain 3 variables.

14) Sorting the tuple using first variable of each member.

15) Saving top 'n' reactions

16) Printing the reactions.

17) Reversing the reactions in order to plot bar graphs from highest to lowest sensitivites.

18) Plotting the data.

Code:

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

gas = ct.Solution('gri30.xml')
temp_T = 1500		#K
pres_P = ct.one_atm

gas.TPX = temp_T,pres_P, {'CH4':1,'O2':2,'N2':7.52}
r = ct.IdealGasConstPressureReactor(gas)
sim = ct.ReactorNet([r])

rxn = gas.reaction_equations()		# list of all reactions in GRI file. Each of reaction represented as string
n_rxn = gas.n_reactions				# number of reactions in GRI file. Equal to 'len(rxn)'

for i in range(n_rxn):
	r.add_sensitivity_reaction(i)	# Asking Cantera to also compute sensitivity of reaction at index 'i'

max_s = [0]*n_rxn					# Array to store maximum sensitivity of each reaction in whole time
max_s_abs = [0]*n_rxn				# Array to store maximum absolute sensitivity of each reaction in whole time

sim.rtol=1e-6						# The relative error tolerance used while integrating the reactor equations
sim.atol=1e-16						# The absolute error tolerance used while integrating the reactor equations

sim.rtol_sensitivity=1e-6			# The relative error tolerance for sensitivity analysis
sim.atol_sensitivity=1e-6			# The absolute error tolerance for sensitivity analysis
states = ct.SolutionArray(gas,extra = 't')		# Storing Solution array in states object

tot_time = 2e-3				# Total time till which simulation must run
dt = 5e-6					# Time step used

for time in np.arange(0,tot_time,dt):
	sim.advance(time)								# Asking cantera to integrate the solution from '0' till 'time'
	states.append(r.thermo.state,t = time*1000)		# Saving thermodynamic state of reactor 'r'

	for j in range(n_rxn):
		s = sim.sensitivity('temperature',j)		# Saving sensitivity of sensitivity parameter 'j' (i.e. rxn[j]) on component(temperature)

		if np.abs(s) > np.abs(max_s[j]):			# Checking and saving maximum sensitivity in 'tot_time' of each reaction 'j'
			max_s[j] = s
			max_s_abs[j] = abs(s)

sen_zip = zip(max_s_abs,max_s,range(n_rxn))				# Creating a zip datatype which holds 3 variables stated.
sen_tuple = tuple(sen_zip)								# Converting zip to tuple inorder to access the data inside the zip
sen_tuple_sorted = sorted(sen_tuple, reverse = True)	# Sorting the tuple based on first element of each member tuple i.e. max abs sensitivity

top_n = 10						# Top N reactions which affect the temperature
top_rxn_index = [0]*top_n 		# Array to save index of reaction
top_rxn = [0]*top_n 			# Array to save reaction equation as string
top_sen = [0]*top_n 			# Array to save sensitivity of reaction with sign
for i in range(top_n):
	top_rxn_index[i] = int(sen_tuple_sorted[i][2])
	top_rxn[i] = rxn[top_rxn_index[i]]
	top_sen[i] = float(sen_tuple_sorted[i][1])

print("The top %d temperature sensitive reactions are:" % (top_n))
print("------------------------------------------------------------------")
print("%10s %s %18s %12s %12s" % ('Reaction Index','|','Reaction','|','Sensitivity'))
print("------------------------------------------------------------------")
for i in range(top_n):
	print("%7d %8s %26s %4s %10.4f" % (top_rxn_index[i],'|',top_rxn[i],'|',top_sen[i]))
print("------------------------------------------------------------------")

top_rxn_index.reverse()   	# Reversing data inorder to plot bar graphs from maximum sensitivity to minimum sensitivity
top_rxn.reverse()
top_sen.reverse()

plt.figure(1)
plt.barh(range(top_n),top_sen,0.75,align = 'center')
plt.yticks(range(10),top_rxn)
plt.xlabel('Temperature Sensitivity', fontweight = 'bold')
plt.title('{0} Most Temparture Sensitive Reactions'.format(top_n))
plt.grid(which='major')
plt.minorticks_on()
plt.grid(which='minor', alpha=0.2)

fig1,a1 = plt.subplots(2,4)
a1[0,0].plot(states.t, states.T)
a1[0,0].set_title('Temperature vs time')
a1[0,0].set_xlabel('Time (ms)')
a1[0,0].set_ylabel('Temperature (K)')
a1[0,0].grid('both')

a1[0,1].plot(states.t,states('CH4').X)
a1[0,1].set_title('CH4 mole fraction')
a1[0,1].set_xlabel('Time (ms)')
a1[0,1].grid('both')

a1[0,2].plot(states.t,states('O2').X)
a1[0,2].set_title('O2 mole fraction')
a1[0,2].set_xlabel('Time (ms)')
a1[0,2].grid('both')

a1[0,3].plot(states.t,states('N2').X)
a1[0,3].set_title('N2 mole fraction')
a1[0,3].set_xlabel('Time (ms)')
a1[0,3].grid('both')

a1[1,0].plot(states.t,states('CO2').X)
a1[1,0].set_title('CO2 mole fraction')
a1[1,0].set_xlabel('Time (ms)')
a1[1,0].grid('both')

a1[1,1].plot(states.t,states('H2O').X)
a1[1,1].set_title('H2O mole fraction')
a1[1,1].set_xlabel('Time (ms)')
a1[1,1].grid('both')

a1[1,2].plot(states.t,states('OH').X)
a1[1,2].set_title('OH mole fraction')
a1[1,2].set_xlabel('Time (ms)')
a1[1,2].grid('both')

a1[1,3].plot(states.t,states('CO').X)
a1[1,3].set_title('CO mole fraction')
a1[1,3].set_xlabel('Time (ms)')
a1[1,3].grid('both')
plt.show()

Results:

The above image lists the top 10 reactions to which the temperature is sensitive. A positive sensitivity means that increase in rate of reaction would increase temperature and negative sensitivity means that increase in rate of reaction would decrease temperature.

This bar graph shows a pictorial representation of the sensitivity of temperature due to reactions.

The above image shows mole fractions of various species. The Temperature is seen rising from 1500 K to 2775 K at 1.15 milliseconds. This shows that mixture autoignited. The CH4 concentration was seen lowering down to 0, indicating that all fuel was spent. The N2 mole fraction is seen to lower down and then rising up very little, indicating formation of NOx. Other species also show change in their mole fractions.

Conclusion:

The sensitivity analysis is carried out using cantera and top 10 reactions are found which affect temperature.