Reaction Mechanism Reduction Project

INTRODUCTION :-

Detailed chemical kinetic mechanisms contain up to several thousand species and several ten thousand reaction steps. Within such large mechanisms, the key processes can be pointed out by the presence of many reaction steps of only marginal importance. A first step to reducing the size of a kinetic mechanism is to identify species and reaction steps which do not need to be included in order to accurately predict the key target outputs of the model.

A chemical reaction usually occurs in steps, although it may not always be obvious to an observer. We call each step in a reaction mechanism an elementary reaction. Elementary reactions occur exactly as they are written and cannot be broken down into simpler steps. Elementary reactions add up to the overall reaction right into the mechanism. The reaction mechanism describes the sequence of elementary reactions that must occur to go from reactants to products. Hence, the computational cost for solving these mechanisms is very high of order of `N^3`where N - no. of species involved in reaction mechanism.

The very purpose of reduction of these mechanism is to identify and neglect such reactions which have negligible impact on mechanism. In this problem, reduction of GRI 3.0 mechanism is done which has 325 reactions through sensitivity analysis w.r.t temperature.

CODE METHODOLOGY :-

• Input Conditions are fed using a text file. The file parser takes in input parameters and reads in form of arrays i.e. analysis is done by nest looping. Following are values parsed from text file.

State Values Range -

Temperature - [950 , 1050 , 1150] Kelvin

Pressure - [3 , 4 , 5] bar

Equivalence ratio - [0.5 , 1 , 1.5]

• Sensitivity analysis is performed w.r.t temperature and decreasing order is sort order. Starting input conditions are with base temperaure and ignition delay for each cases.
• Out of all reactions, sorting is done until error in Ignition delay is below 5%

Steps carried out :-

1. Parse the data and create gas instance of Solution class.
2. Create IdealGasReactor and ReactorNet objects and simulate 0-D reactor to calculate Max temperature and ignition delay for each case.
3. Sorting is done based on absolute values of sensitivities (un-reduced) and priortize the reactions based on sorted reactions from each case to obtain overall reaction order.
4. If reactions aren't occuring i.e. due lesser due to no. of reactions, using NameError we catch this exception.
5. Then store these reactions in reduced mechanism if all values meet the given tolerance value and proceed with calculations for rest and plot the results.

PYTHON CODE :-

"""
Program for reduction of GRI 3.0 mechanism  based on temperature sensitivity

"""

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


# Python class to import data from file
class FileReader:

	def __init__(self, filename):
		self.filename = filename
		self.input_dict = {}

	def validate_file(self):
		return os.path.isfile(self.filename)

	def read(self):
		if self.validate_file():
			for line in open(self.filename, 'r'):
				self.input_dict[line.split()[0]] = line.split()[1:]
				
		else:
			print('File name does not exist')

a = FileReader('input.txt')
a.read()
T_val = [float(x) for x in a.input_dict['Temperature']]
P_val = [float(x) for x in a.input_dict['Pressure']]
phi_val = [float(x) for x in a.input_dict['Phi']]
for x in a.input_dict['End_time']: t_end = float(x)
for x in a.input_dict['Delta_t']: dt = float(x)
for x in a.input_dict['Min_size']: min_size = int(x)
for x in a.input_dict['Tol_Tmax']: tol_T_max = float(x)
for x in a.input_dict['Tol_igd']: tol_ign_delay = float(x)
for x in a.input_dict['Fuel']: fuel = x

#################################################

# Ordering the reaction indices in decreasing order of sensitivity to temperature for all combination of state values

S = []
S_max0 = []
R_index = []
const1 = True
n = 0


air = 'O2:1, N2:3.76'

# Finding the reference Ignition Delay and Maximum Temperature
gas = ct.Solution('gri30.cti')
S_max = [0]*gas.n_reactions

# state values combination loop
for T in T_val:

	for P in P_val:

		for phi in phi_val:
			
			gas.TP = T, P
			gas.set_equivalence_ratio(phi, fuel, air)

			# creating a reactor and reactor net objcets
			r = ct.IdealGasReactor(gas)
			sim = ct.ReactorNet([r])

			# including all reactions from gri mech to perform sensitivity analysis
			for i in range(0, gas.n_reactions):
				r.add_sensitivity_reaction(i)

			# # set the tolerances for the solution and the sensitivity co-efficients
			sim.rtol = 1e-6 # Realtive change in the solution < sim.rtol
			sim.atol = 1e-15 # Absolute change in the solution < sim.atol

			sim.rtol_sensitivity = 1e-6 #relative change in sensitivity < sim.rtol_sense
			sim.atol_sensitivity = 1e-6 #absolute change in sensitivity =400):
		

			# solving the reaction rate odes and sensitivity odes
			for t in np.arange(0, t_end+dt, dt):

				sim.advance(t)

				# finding maximum sensitivities for all reactions
				for i in range(0, gas.n_reactions):
					
					S = sim.sensitivity(2, i)

					if np.abs(S) > np.abs(S_max[i]):
						S_max[i] = np.abs(S)

			# zipping the max sensitivities with respective reaction indices
			for i in range(0, gas.n_reactions):
				S_max0.append([])
				S_max0[i].append(np.abs(S_max[i]))
				S_max0[i].append(i)

			# sorting S_max0 array in decreasing order of sensitivity
			S_sorted = sorted(S_max0, key = lambda x: x[0], reverse = True)
			
			# appending the sorted order of reaction indices to R_index
			for i in range(0, gas.n_reactions):

				if const1 is True:
					R_index.append([])

				R_index[i].append(S_sorted[i][1])
			
			# resetting values
			const1 = False
			S_max0.clear()
			n = n + 1

#################################################

# flattening out R_index to create a single list of reactions from all state value combinations in the order of decreasing sensitivities

R_set = set()
R_ordered = []

for i in range(0, gas.n_reactions):

	for j in range(0, n):

		if R_index[i][j] not in R_set:
			R_ordered.append(R_index[i][j])
			R_set.add(R_index[i][j])

print(R_ordered,'n')
#################################################

# Determining the size of the reduced mechanism

spec = ct.Species.listFromFile('gri30.cti')
rxns = []
ign_delay_arr = []
T_max_arr = []
mech_size = []

for T in T_val:
	for P in P_val:
		for phi in phi_val:
			
			# computing reference ign_delay and T_max for each state
			T_ign = T + 400	# K

			gas.TP = T, P
			gas.set_equivalence_ratio(phi, fuel, air)

			r = ct.IdealGasReactor(gas)
			sim = ct.ReactorNet([r])
			
			const2 = True
			T_max_ref = 0

			for t in np.arange(0, t_end+dt, dt):

				sim.advance(t)

				if (gas.T > T_ign and const2 is True):
					ign_delay_ref = t
					const2 = False
				if (gas.T > T_max_ref):
					T_max_ref = gas.T

			# performing mechanism reduction for each state
			err_T_max = 100
			err_ign_delay = 100
			i = 0
			const4 = True
			
			# appending reactions one by one and comparing T_max and ign_delay results
			for i in range(0, gas.n_reactions):

				rxns.append(gas.reaction(R_ordered[i]))
				gas1 = ct.Solution(thermo = 'IdealGas', kinetics = 'GasKinetics', species = spec, reactions = rxns)
				gas1.TP = T, P
				gas1.set_equivalence_ratio(phi, fuel, air)
				r1 = ct.IdealGasReactor(gas1)
				sim1 = ct.ReactorNet([r1])
				const3 = True
				T_max = 0

				for t in np.arange(0, t+dt, dt):

					sim1.advance(t)

					if (gas1.T > T_ign and const3 is True):
						ign_delay = t
						const3 = False

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

				# incrementing reaction index
				i = i + 1

				# catching exception until the no. of reactions increase to suitable number
				try :
					
					# storing values for the plot
					if (i >= min_size) :
						ign_delay_arr.append(ign_delay)
						T_max_arr.append(T_max)
						mech_size.append(i)

					# computing accuracy
					err_T_max = np.abs(((T_max_ref - T_max)/T_max_ref)*100) # % error in maximum temperature
					err_ign_delay = np.abs(((ign_delay_ref - ign_delay)/ign_delay_ref)*100)	# % error in ignition delay

					# printing mechanism size and other results when the tolerance criteria is satisfied
					if (err_T_max < tol_T_max and err_ign_delay < tol_ign_delay and const4 is True):
						
						print('T = %g, P = %g, phi = %g, T_max_ref = %g, ign_delay_ref = %g, T_max = %g, ign_delay = %g, T_max_error = %g, Ign_delay_error = %g, red_mech_size = %g' %(T, P, phi, T_max_ref, ign_delay_ref, T_max, ign_delay, err_T_max, err_ign_delay, i))
						const4 = False

				except NameError:

					pass

			
			# plotting results
			plt.subplot(2,1,1)
			plt.title('State Values: T = %g K, P = %g bar, phi = %g' %(T, P/100000, phi),fontweight='bold',fontsize=17)
			plt.plot(mech_size, ign_delay_arr, color = 'blue')
			plt.axhline(y = ign_delay_ref, color = 'red', linestyle = 'dotted', label = 'Reference Ignition Delay')
			plt.xlabel('No. of reactions most sensitive to Temperature',fontweight='bold',fontsize=11)
			plt.ylabel('Ignition Delay [sec]',fontweight='bold',fontsize=9)
			plt.grid('both',linestyle='--',linewidth=1.2)
			plt.legend(loc = 'best')

			plt.subplot(2,1,2)
			plt.plot(mech_size, T_max_arr, color = 'red')
			plt.axhline(y = T_max_ref, color = 'red', linestyle = 'dotted', label = 'Reference Max Temperature')
			plt.xlabel('No. of reactions most sensitive to Temperature',fontweight='bold',fontsize=11)
			plt.ylabel('Maximum Temperature (K)',fontweight='bold',fontsize=9)
			plt.grid('both',linestyle='--',linewidth=1.2)
			plt.legend(loc = 'best')

			plt.tight_layout()
			plt.show()

			# resetting arrays for next state values combination
			rxns.clear()
			ign_delay_arr.clear()
			T_max_arr.clear()
			mech_size.clear()

PLOTS :-

1. T = 950, P = 300000, phi = 0.5, T_max_ref = 2283.66, ign_delay_ref = 0.6, T_max = 2283.63, ign_delay = 0.6, T_max_error = 0.0011619, Ign_delay_error = 0, red_mech_size = 159

2. T = 950, P = 300000, phi = 1, T_max_ref = 2827.51, ign_delay_ref = 0.77, T_max = 2827.78, ign_delay = 0.77, T_max_error = 0.00959652, Ign_delay_error = 0, red_mech_size = 160

3. T = 950, P = 300000, phi = 1.5, T_max_ref = 2724.85, ign_delay_ref = 0.94, T_max = 2724.58, ign_delay = 0.94, T_max_error = 0.00990105, Ign_delay_error = 0, red_mech_size = 114

4. T = 950, P = 400000, phi = 0.5, T_max_ref = 2286.38, ign_delay_ref = 0.446, T_max = 2286.34, ign_delay = 0.446, T_max_error = 0.00150478, Ign_delay_error = 0, red_mech_size = 159

5. T = 950, P = 400000, phi = 1, T_max_ref = 2854.33, ign_delay_ref = 0.555, T_max = 2854.57, ign_delay = 0.555, T_max_error = 0.00851258, Ign_delay_error = 0, red_mech_size = 205

6.T = 950, P = 400000, phi = 1.5, T_max_ref = 2733.61, ign_delay_ref = 0.669, T_max = 2733.64, ign_delay = 0.669, T_max_error = 0.00126037, Ign_delay_error = 0, red_mech_size = 114

7.T = 950, P = 500000, phi = 0.5, T_max_ref = 2288.3, ign_delay_ref = 0.356, T_max = 2288.48, ign_delay = 0.356, T_max_error = 0.00805527, Ign_delay_error = 0, red_mech_size = 65

8. T = 950, P = 500000, phi = 1, T_max_ref = 2859.97, ign_delay_ref = 0.432, T_max = 2859.98, ign_delay = 0.432, T_max_error = 0.000197673, Ign_delay_error = 0, red_mech_size = 114

9. T = 950, P = 500000, phi = 1.5, T_max_ref = 2739.99, ign_delay_ref = 0.515, T_max = 2740.15, ign_delay = 0.515, T_max_error = 0.00580735, Ign_delay_error = 0, red_mech_size = 114

10. T = 1050, P = 300000, phi = 0.5, T_max_ref = 2356.74, ign_delay_ref = 0.103, T_max = 2356.85, ign_delay = 0.103, T_max_error = 0.00467677, Ign_delay_error = 0, red_mech_size = 101

11. T = 1050, P = 300000, phi = 1, T_max_ref = 2862.49, ign_delay_ref = 0.142, T_max = 2862.76, ign_delay = 0.142, T_max_error = 0.00924147, Ign_delay_error = 0, red_mech_size = 181

12. T = 1050, P = 300000, phi = 1.5, T_max_ref = 2783.64, ign_delay_ref = 0.176, T_max = 2783.48, ign_delay = 0.176, T_max_error = 0.00587908, Ign_delay_error = 0, red_mech_size = 179

13. T = 1050, P = 400000, phi = 0.5, T_max_ref = 2359.7, ign_delay_ref = 0.076, T_max = 2359.85, ign_delay = 0.076, T_max_error = 0.00663949, Ign_delay_error = 0, red_mech_size = 73

14. T = 1050, P = 400000, phi = 1, T_max_ref = 2878.94, ign_delay_ref = 0.103, T_max = 2878.92, ign_delay = 0.103, T_max_error = 0.000798105, Ign_delay_error = 0, red_mech_size = 105

15. T = 1050, P = 400000, phi = 1.5, T_max_ref = 2792.65, ign_delay_ref = 0.127, T_max = 2792.38, ign_delay = 0.127, T_max_error = 0.00980892, Ign_delay_error = 0, red_mech_size = 106

16. T = 1050, P = 500000, phi = 0.5, T_max_ref = 2362.09, ign_delay_ref = 0.06, T_max = 2362.13, ign_delay = 0.06, T_max_error = 0.00197306, Ign_delay_error = 0, red_mech_size = 159

17. T = 1050, P = 500000, phi = 1, T_max_ref = 2913.39, ign_delay_ref = 0.079, T_max = 2913.44, ign_delay = 0.079, T_max_error = 0.00192204, Ign_delay_error = 0, red_mech_size = 179

18. T = 1050, P = 500000, phi = 1.5, T_max_ref = 2800.22, ign_delay_ref = 0.098, T_max = 2800.1, ign_delay = 0.098, T_max_error = 0.00437602, Ign_delay_error = 0, red_mech_size = 106

19. T = 1150, P = 300000, phi = 0.5, T_max_ref = 2428.12, ign_delay_ref = 0.023, T_max = 2428.28, ign_delay = 0.023, T_max_error = 0.00670668, Ign_delay_error = 0, red_mech_size = 101

20. T = 1150, P = 300000, phi = 1, T_max_ref = 2906.9, ign_delay_ref = 0.032, T_max = 2907, ign_delay = 0.032, T_max_error = 0.00358204, Ign_delay_error = 0, red_mech_size = 179

21. T = 1150, P = 300000, phi = 1.5, T_max_ref = 2836.81, ign_delay_ref = 0.04, T_max = 2837.05, ign_delay = 0.04, T_max_error = 0.00845644, Ign_delay_error = 0, red_mech_size = 172

22. T = 1150, P = 400000, phi = 0.5, T_max_ref = 2430.22, ign_delay_ref = 0.018, T_max = 2430.39, ign_delay = 0.018, T_max_error = 0.00724013, Ign_delay_error = 0, red_mech_size = 159

23. T = 1150, P = 400000, phi = 1, T_max_ref = 2913.11, ign_delay_ref = 0.024, T_max = 2913.02, ign_delay = 0.024, T_max_error = 0.00308688, Ign_delay_error = 0, red_mech_size = 163

24. T = 1150, P = 400000, phi = 1.5, T_max_ref = 2849.48, ign_delay_ref = 0.03, T_max = 2849.47, ign_delay = 0.03, T_max_error = 0.000541003, Ign_delay_error = 0, red_mech_size = 106

25. T = 1150, P = 500000, phi = 0.5, T_max_ref = 2435, ign_delay_ref = 0.014, T_max = 2435.08, ign_delay = 0.014, T_max_error = 0.00302839, Ign_delay_error = 0, red_mech_size = 159

26. T = 1150, P = 500000, phi = 1, T_max_ref = 2928.28, ign_delay_ref = 0.019, T_max = 2928.55, ign_delay = 0.019, T_max_error = 0.00909904, Ign_delay_error = 0, red_mech_size = 172

27. T = 1150, P = 500000, phi = 1.5, T_max_ref = 2864.09, ign_delay_ref = 0.023, T_max = 2863.87, ign_delay = 0.023, T_max_error = 0.00788527, Ign_delay_error = 0, red_mech_size = 172

INFERENCE :-

1. Initially both graphs for different combination of state values show fluctuations and dampen to approach steady values as reactions are added. This is due to the fact that the sorting order chosen is decreasing and higher the sensitivity have great impact on ignition delay and max temperature.
2. There is a need for error exception as ignition delay is 400K + `T_i`and until a certain set of reactions are passed the variable ignition_delay in our program will not be initiated i.e. insufficient to cause an ignition hence ignition delay time cannot be recorded. Minimum no. of reactions was observed in range of 35-40 i.e. to prevent any reactions having no combustion effect and was preferably set to 40.
3. It can be observed from the graphs that in some cases Maximum temperature and ignition delay values can be predicted from few no. of reactions. However a minimum of 160 reactions are required in the mechanism to predict accurate ignition delay and max temperature results