Reaction Mechanism Reduction using Cantera

Aim: To perform Reaction Mechanism Reduction for GRI3.0 mechanism using Cantera(python)

1. Introduction:

A chemical reaction undergoes a series of intermediate reactions before reaching equilibrium. These intermediate reactions combined together form a Reaction Mechanism. Every intermediate reaction is important. Certain reactions will be sensitive, which means that because of these reactions, temperature or any other property may undergo drastic changes.

If only the highest sensitive reactions are kept and other less sensitive reactions are removed, than the computational time decreases. This concept is applied in this project for a range of temperature, pressure and equivalence ratio values.

During combustion, initially a set of reactants undergo chemical reaction and form some finite number of intermediate reactants. Due to the first reaction, temperature increases. This sets off other intermediate reactions. So, even though temperature increases after a particular reaction, it decreases eventually to support other intermediary reactions. This trend can be observed in the results section.

Mechanism reduction is performed for methane gas. The GRI 3.0 mechanism involves only 325 reactions. But complex mechanism will have thousands of intermediary reactions. In such cases, to save time reaction mechanism must be reduced to a minimum number which can still yield accurate values of AFT and ignition delay.

2. Approach:

• Required modules are imported.
• An object called file reader is written to read through, validate and split the lines in a given input text file.
• The program is written to find how reduction in reaction mechanism affects Adiabatic Flame Temperature(AFT) and ignition delay. So, 27 different scenarios are considered with varying temperature, pressure and equivalence ratios.
• In the first step, sensitivity values of all the reactions are determined and they are sorted in the decreasing order.
• These reactions are zipped, which means they are merged with the reaction indices. They are zipped in such a way, only the reaction numbers are considered for next steps and not the sensitivity values.
• The previous step determines the reaction number for just one iteration or 1 set of values of temperature, pressure and equivalence ratio. Same procedure is repeated for the remaining 26 set of values.
• To store all the above found values in a single place, a list is created, where the 1st row includes all the reaction numbers of highest sensitive reactions in all set of values and the second row includes the second highest sensitive reaction numbers and so on.
• A new mechanism is created by comparing all the values in the previous list of reaction numbers. No reaction numbers are same in this mechanism.
• This mechanism containes all 325 reactions, but in decreasing sensitivity order, taken from all combinations of temperature, pressure and equivalence ratio.
• The second step is to determine AFT and ignition delay values for both old and updated mechanisms.
• For old and new mechanisms, AFT and ignition delay are determined by creating a gas object and defining its properties, creating an Ideal gas reactor and a reactor network and at last by advancing this reactor network for a certain time step size.
• For new mechanism, AFT and ignition delays are calculated for every reaction to plot the variation.
• The size of new mechanism is determined by comparing the error between AFTs and ignition delays of old and new mechanisms. If AFT or ignition delay satisfy the tolerance criteria, than the reaction number corresponding to this tolerance is the mechanism size.
• That is the mechanism size required to accurately estimate the AFT and ignition delay.
• The minimum number of reactions required to create a new mechanism is considered to be 40.
• Finally, the values of AFT and ignition delay are plotted against the number of reactions.

3. Code:

import os
import cantera as ct
import numpy as py
import matplotlib.pyplot as plt
from tabulate import tabulate

class file_reader:
	def __init__(self,file_name):
		self.file_name = file_name
		self.input_dict = {}

	def validate_input_file(self):
		if os.path.isfile(self.file_name):
			return True
		else:
			return False

	def read(self):
		if self.validate_input_file():
			for line in open(self.file_name,'r'):
				self.input_dict[line.split()[0]] = line.split()[1:]

		else:
			print('validation error')

data = file_reader('data.txt')
data.read()
T = [float(value) for value in data.input_dict['Temperature']]
P = [float(value) for value in data.input_dict['Pressure']]
phi = [float(value) for value in data.input_dict['Equivalence_ratio']]

methane = ct.Solution('gri30.xml')

rxns = methane.n_reactions
sens_max = [0]*rxns
reaction_indices = py.arange(0,rxns)

table_index = []

for temp in T:
	
	for pressure in P:

		for eqr in phi:		
			
			methane.TPX = temp, pressure, {'CH4':1, 'O2':2/eqr, 'N2':7.52/eqr}

			reactor = ct.IdealGasReactor(methane)
			reactor_network = ct.ReactorNet([reactor])

			for a in range(0,rxns):
				reactor.add_sensitivity_reaction(a)

			reactor_network.rtol = 1e-6
			reactor_network.atol = 1e-15

			reactor_network.rtol_sensitivity = 1e-6
			reactor_network.atol_sensitivity = 1e-6

			for t in py.arange(0, 5, 1e-3):
				reactor_network.advance(t)

				for b in range(0,rxns):
					sens = reactor_network.sensitivity(1,b)

					if abs(sens) >= abs(sens_max[b]):
						sens_max[b] = abs(sens)

			reaction_number = [d for e, d in sorted(zip(sens_max,reaction_indices),reverse=True)]

			for c in range(0,rxns):
				table_index.append([])
				table_index[c].append(reaction_number[c])
			
reduced_reactions = []

for i in range(0,rxns):
	for j in range(0,26):

		if table_index[i][j] not in reduced_reactions:
			reduced_reactions.append(table_index[i][j])

print(reduced_reactions)

spec = ct.Species.listFromFile('gri30.xml')

reactions = []
ignition_dealy_array = []
T_max_array = []
mechanism_size = []

min_reactions = 40
tol = 0.05

for temp1 in T:

	for pressure1 in P:

		for eqr1 in phi:

			Tig = temp1 + 400

			methane.TPX = temp1, pressure1, {'CH4':1, 'O2':2/eqr1, 'N2':7.52/eqr1}

			reactor = ct.IdealGasReactor(methane)
			reactor_network = ct.ReactorNet([reactor])

			constant1 = True
			T_max_ref = 0

			for time1 in py.arange(0, 5, 1e-3):

				reactor_network.advance(time1)

				if (methane.T > Tig and constant1 is True):
					ignition_delay_ref = time1
					constant1 = False 

				if (methane.T > T_max_ref):
					T_max_ref = methane.T

			error_T_max = 100
			error_ignition_delay = 100
			i = 0
			constant2 = True

			for i in range(0,rxns):

				reactions.append(methane.reaction(reduced_reactions[i]))

				methane1 = ct.Solution(thermo = 'IdealGas', kinetics = 'GasKinetics', species = spec, reactions = reactions)
				methane1.TPX = temp1, pressure1, {'CH4':1, 'O2':2/eqr1, 'N2':7.52/eqr1}

				reactor1 = ct.IdealGasReactor(methane1)
				reactor_network1 = ct.ReactorNet([reactor1])

				constant3 = True
				T_max = 0

				for time2 in py.arange(0,5,1e-3):

					reactor_network1.advance(time2)

					if (methane1.T > Tig and constant3 is True):
						ignition_delay = time2
						constant3 = False

					if (methane1.T > T_max):
						T_max = methane1.T

				i = i +1

				try:


					if (i >= min_reactions):
						ignition_dealy_array.append(ignition_delay)
						T_max_array.append(T_max)
						mechanism_size.append(i)

					error_T_max = py.abs(((T_max_ref - T_max)/T_max_ref)*100)
					error_ignition_delay = py.abs(((ignition_delay_ref - ignition_delay)/ignition_delay_ref)*100)

					if (error_T_max < tol and error_ignition_delay < tol and constant2 is True):

						print(tabulate([['T','%gk'], ['P', '%gpa'], ['Equivalence_ratio', '%g'], ['T_max', '%gk'], ['Ignition_Delay', '%gms'], ['T_max_ref', '%gk'],['Ignition Dealy Ref','%gms'],['Mechanism size', '%g']], headers = ['State Variables', 'Values'], tablefmt='orgtbl')%(temp1, pressure1, eqr1,T_max, ignition_delay,T_max_ref, ignition_delay_ref, i))
						constant2 = False

				except NameError:

					pass

			plt.subplot(2,1,1)
			plt.plot(mechanism_size, ignition_dealy_array)
			plt.axhline(y = ignition_delay_ref, linestyle = 'dotted', color = 'black', label = 'Ignition Delay reference value')
			plt.ylabel('Ignition Delay in ms')
			plt.xlabel('Number of reactions')
			plt.title('Variation of Ignition Delay wrt most sensitive reactions\n for state variable values: T = %g k, P = %g pa, phi = %g' %(temp1,pressure1,eqr1))

			plt.subplot(2,1,2)
			plt.plot(mechanism_size, T_max_array, color = 'red')
			plt.axhline(y = T_max_ref, linestyle = 'dotted', color = 'black')
			plt.xlabel('Number of reactions')
			plt.ylabel('Maximum Temperature in k')
			plt.title('Variation of AFT wrt most sensitive reactions\n for state variable values: T = %g k, P = %g pa, phi = %g' %(temp1,pressure1,eqr1))

			plt.tight_layout()
			plt.show()

			reactions.clear()
			ignition_dealy_array.clear()
			T_max_array.clear()
			mechanism_size.clear()

4. Results:

• The above table indicates the final values of AFT, ignition delay and mechanism size for old and new mechanisms.

4.1 Results for Temperature of 950k with all 3 combinations of pressure and equivalence ratio

4.2 Results for Temperature of 1050k with all 3 combinations of pressure and equivalence ratio

4.3 Results for Temperature of 1150k with all 3 combinations of pressure and equivalence ratio

• After a reaction has happened, it releases some amount of energy. This energy causes the temperature to increase and this can be seen in the above graphs, where after a certain reaction number temperature rises exponentially.
• This same energy is also consumed by the next reaction. So, the temperature increase or decrease after the second reaction is dictated by the energy remaining after consumption and energy released after reaction. If the sum total of both of these energies is higher than the energy before the second reaction, than the temperature increases further. If not temperature decreases.
• This trend of increasing and decreasing of temperature can be observed in all of the above graphs.
• Once the reactions consume and release the same amount energy, the reaction reaches equilibrium and there will be lesser variations in temperature.
• As the temperature dips after reaching its peak and if it doesnt change after its dip, then the reaction number corresponding to the dip is the minimum reactions required to determine AFT and ignition delay.
• Also, with increase in the initial temperature, reaction mechanism size decreases.
• For ignition delay, once temperature increases exponentially, combustion or ignition begins. This time value for ignition doesnt change throughout the reaction.