Reaction mechanism reduction

Aim: The main objective of this project is to write a Reaction mechanism reductio code for GRI 3.0 mech. This project is done to identify the minimum number of reactions required to obtain the maximum value for ignition delay and temperature at a given tolerance value.

The following are the step by step objective of the program:

1. First, the program should run a reference Ignition delay and Max temperature Analysis for the unreduced mechanism for a set of state values (Temperature, Pressure and Equivalence ratio) for a time duration of (5 sec)
2. Then, a sensitivity analysis concerning temperature must be run on all the reactions in the mechanism and must sort the reactions from the most sensitive to least sensitive
3. Furthermore, these sorted reactions should be added to a new mechanism file/object one by one starting from a fixed minimum number of reactions (40).
4. Now a graph must be plotted plotting the reduced mechanism's ignition delay and max temperature while having the reference ignition delay and Temperature as a label for reference (with 5% tolerance)
5. This graph must have the values of Ignition Delay and Max temperature on the Y-axis, while the number of reactions in the mechanism on the X-axis.

Additionally, the reaction mechanism code is computed for all combinations of the below-mentioned state value range:

State Values Range -

Temperature - [950 , 1050 , 1150] Kelvin

Pressure - [3 , 4 , 5] bar

Equivalence ratio - [0.5 , 1 , 1.5]

A reaction mechanism is a step by step sequence of elementary reactions by which overall chemical change occurs. It is a theoretical compilation that explains in detail what happens at every stage of a chemical reaction. Mechanism Reduction is a process where we remove redundant reactions from the mechanism concerning the problem parameters. This helps us remove reactions that do not play a major role in the overall chemical reaction. Mechanism reduction would help us customized solvers with a bunch of pertinent reactions that would play a pivotal role in a given problem parameter. The other reactions are discarded. This technique helps us reduce the simulation time by a considerable amount and simplifies our understanding of the output and would help us make targeted correction wherever necessary depending upon the requirement.

In this project, the initial temperature, pressure, and equivalence ratio values are obtained by file parsing techniques which are shown in the code.

The solution is obtained by following a simple procedure. Initially, a reference ignition delay and maximum temperature value are computed for the reaction using the all possible state value combination. Now, a sensitivity analysis for all the 325 reactions and all the state values is computed. The most sensitive reactions from this are sorted and arranged in the descending order of their sensitivities. All the other reactions are discarded. After this step, an ignition delay and maximum temperature calculation are run form a fixed set of reactions starting from 40 and ending at 325. These values are computed in their descending order. The values obtained are compared with the help of a graph and the required output is printed.

Source code:

"Text for file parsing"

Temperature     950 1050 1150
Pressure	300000 400000 500000
Phi		0.5 1.0 1.5

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

# Extracting mech_reduction date from 'mech_reduction_file' using user-defined class 'FileReader'
class FileReader:

	def __init__(self, filename):
		self.filename = filename
		self.mech_reduction_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.mech_reduction_dict[line.split()[0]] = line.split()[1:]
				
		else:
			print('No file under such name')

a = FileReader('mech_reduction.txt')
a.read()
T1 = [float(x) for x in a.mech_reduction_dict['Temperature']]
P1 = [float(x) for x in a.mech_reduction_dict['Pressure']]
equivalence_ratio = [float(x) for x in a.mech_reduction_dict['Phi']]

size = 40
S = []
pw_max = []
ind1 = []
n = 0
t_end = 5
dt = 0.001
tol_Tmax = 0.05
tol_ignd = 0.05
gas = ct.Solution('gri30.xml')
S_max = [0]*gas.n_reactions
spec = ct.Species.listFromFile('gri30.xml')
K7 = []
ignd_a1 = []
Tmax_b1 = []
size_of_mech = []
case = 1

for T in T1:
	for P in P1:
		for phi in equivalence_ratio:
			gas.TP = T, P
			gas.set_equivalence_ratio(phi, 'CH4', 'O2:1, N2:3.76')
			r = ct.IdealGasReactor(gas)
			sim = ct.ReactorNet([r])

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

			sim.rtol = 1e-6
			sim.atol = 1e-15
			sim.rtol_sensitivity = 1e-6
			sim.atol_sensitivity = 1e-6
	
			for t in np.arange(0, t_end+dt, dt):
				sim.advance(t)
				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)

			for i in range(0, gas.n_reactions):
				pw_max.append([])
				pw_max[i].append(np.abs(S_max[i]))
				pw_max[i].append(i)

			# sorting pw_max array in decreasing order of sensitivity
			S_sorted = sorted(pw_max, key = lambda x: x[0], reverse = True)
			
			# appending the sorted order of reaction indices to ind1
			for i in range(0, gas.n_reactions):
				
				ind1.append([])
				ind1[i].append(S_sorted[i][1])
	
			pw_max.clear()
			n = n + 1

R_set = set()
R_ordered = []
for i in range(0, gas.n_reactions):
	for j in range(0, n):
		if ind1[i][j] not in R_set:
			R_ordered.append(ind1[i][j])
			R_set.add(ind1[i][j])
print(R_ordered, n)
# Determining the size of the reduced mechanism
for T in T1:
	for P in P1:
		for phi in equivalence_ratio:
			# computing reference ign_delay and T_max for each state
			T_ign = T + 400	# K

			gas.TP = T, P
			gas.set_equivalence_ratio(phi, 'CH4', 'O2:1, N2:3.76')

			r = ct.IdealGasReactor(gas)
			sim = ct.ReactorNet([r])
		
			Tmax_r = 0
			c1 = True
			for t in np.arange(0, t_end+dt, dt):

				sim.advance(t)

				if (T_ign < gas.T and c1 is True):
					ignd_r = t
					c1 = False

				if (Tmax_r < gas.T):
					Tmax_r = gas.T

			error_Tmax = 0
			error_ignd = 0
			i = 0
			c3 = True

			for i in range(0, gas.n_reactions):
				K7.append(gas.reaction(R_ordered[i]))
				gas1 = ct.Solution(thermo = 'IdealGas', kinetics = 'GasKinetics', species = spec, reactions = K7)
				gas1.TP = T, P
				gas1.set_equivalence_ratio(phi, 'CH4', 'O2:1, N2:3.76')
				r1 = ct.IdealGasReactor(gas1)
				sim1 = ct.ReactorNet([r1])
				T_max = 0
				c2 = True

				for t in np.arange(0, t+dt, dt):
					sim1.advance(t)
					if ( T_ign < gas1.T and c2 is True):
						ign_delay = t
						c2 = False

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

				# reaction index increment
				i = i + 1
	# finding exceptions until the no. of reactions increases to a suitable number
				try :
					if (i >= size) :
						ignd_a1.append(ign_delay)
						Tmax_b1.append(T_max)
						size_of_mech.append(i)

					# calculating the error
					error_Tmax = np.abs(((Tmax_r - T_max)/Tmax_r)*100) 
					error_ignd = np.abs(((ignd_r - ign_delay)/ignd_r)*100)

# printing the values of phi, temperature etc. when the tolerance value is satisfied
					if (error_Tmax < tol_Tmax and error_ignd < tol_ignd and c3 is True):
						print('T = %g, P = %g, phi = %g, Tmax_r = %g, ignd_r = %g, T_max = %g, ign_delay = %g, T_max_error = %g, Ign_delay_error = %g, red_size_of_mech = %g' %(T, P, phi, Tmax_r, ignd_r, T_max, ign_delay, error_Tmax, error_ignd, i))
						c3 = False
				except NameError:

					pass
			plt.subplot(211)
			plt.title('State Values: T = %g K, P = %g bar, phi = %g' %(T, P/100000, phi))
			plt.plot(size_of_mech, ignd_a1,color = 'brown',label= ' Ign del ref = ' + str(ignd_r) + 'sec' + '\n Temperature = ' + str(T) + 'K' + '\n Pressure = ' + str(P/100000) + 'bar' + '\n Equi Ratio = ' + str(phi))
			plt.axhline(y = ignd_r, color = 'yellow', linestyle = ':', label = 'Ref ign delay line')
			plt.xlabel('No. of temperature sensitive reactions')
			plt.ylabel('Ignition Delay (s)')
			plt.legend(loc = 'best')
			plt.subplot(212)
			plt.plot(size_of_mech, Tmax_b1, color = 'orange', label= 'max ref Temp = ' + str(Tmax_r) + ' Kelvin' + '\n Temperature = ' + str(T) + 'K' + '\n Pressure = ' + str(P/100000) + 'bar' + '\n Equi Ratio = ' + str(phi))
			plt.axhline(y = Tmax_r, color = 'red', linestyle = ':', label = 'Ref temp line')
			plt.xlabel('No. of temperature sensitive reactions')
			plt.ylabel('Max Temperature (K)')
			plt.legend(loc = 'best')
			
# To save the graph as a PNG file
			plt.savefig(' Case '+ str(case) + '.png')
			plt.tight_layout()
			plt.show()

			K7.clear()
			ignd_a1.clear()
			Tmax_b1.clear()
			size_of_mech.clear()
			case += 1

Output:

The output for the code is as follows:

T = 950, P = 300000, phi = 0.5, Tmax_r = 2283.66, ignd_r = 0.6, T_max = 2284.21, ign_delay = 0.6, T_max_error = 0.0238931, Ign_delay_error = 0, red_size_of_mech = 65

T = 950, P = 300000, phi = 1, Tmax_r = 2827.51, ignd_r = 0.77, T_max = 2827.22, ign_delay = 0.77, T_max_error = 0.0103219, Ign_delay_error = 0, red_size_of_mech = 114

T = 950, P = 300000, phi = 1.5, Tmax_r = 2724.85, ignd_r = 0.94, T_max = 2724.58, ign_delay = 0.94, T_max_error = 0.00990105, Ign_delay_error = 0, red_size_of_mech = 114

T = 950, P = 400000, phi = 0.5, Tmax_r = 2286.38, ignd_r = 0.446, T_max = 2286.7, ign_delay = 0.446, T_max_error = 0.0140882, Ign_delay_error = 0, red_size_of_mech = 65

T = 950, P = 400000, phi = 1, Tmax_r = 2854.33, ignd_r = 0.555, T_max = 2853.25, ign_delay = 0.555, T_max_error = 0.0378236, Ign_delay_error = 0, red_size_of_mech = 130

T = 950, P = 400000, phi = 1.5, Tmax_r = 2733.61, ignd_r = 0.669, T_max = 2733.64, ign_delay = 0.669, T_max_error = 0.00126037, Ign_delay_error = 0, red_size_of_mech = 114

T = 950, P = 500000, phi = 0.5, Tmax_r = 2288.3, ignd_r = 0.356, T_max = 2288.48, ign_delay = 0.356, T_max_error = 0.00805527, Ign_delay_error = 0, red_size_of_mech = 65

T = 950, P = 500000, phi = 1, Tmax_r = 2859.97, ignd_r = 0.432, T_max = 2859.98, ign_delay = 0.432, T_max_error = 0.000197673, Ign_delay_error = 0, red_size_of_mech = 114

T = 950, P = 500000, phi = 1.5, Tmax_r = 2739.99, ignd_r = 0.515, T_max = 2740.15, ign_delay = 0.515, T_max_error = 0.00580735, Ign_delay_error = 0, red_size_of_mech = 114

T = 1050, P = 300000, phi = 0.5, Tmax_r = 2356.74, ignd_r = 0.103, T_max = 2357.5, ign_delay = 0.103, T_max_error = 0.0324359, Ign_delay_error = 0, red_size_of_mech = 65

T = 1050, P = 300000, phi = 1, Tmax_r = 2862.49, ignd_r = 0.142, T_max = 2861.18, ign_delay = 0.142, T_max_error = 0.0459318, Ign_delay_error = 0, red_size_of_mech = 130

T = 1050, P = 300000, phi = 1.5, Tmax_r = 2783.64, ignd_r = 0.176, T_max = 2782.97, ign_delay = 0.176, T_max_error = 0.0243533, Ign_delay_error = 0, red_size_of_mech = 163

T = 1050, P = 400000, phi = 0.5, Tmax_r = 2359.7, ignd_r = 0.076, T_max = 2359.85, ign_delay = 0.076, T_max_error = 0.00663949, Ign_delay_error = 0, red_size_of_mech = 73

T = 1050, P = 400000, phi = 1, Tmax_r = 2878.94, ignd_r = 0.103, T_max = 2878.92, ign_delay = 0.103, T_max_error = 0.000798105, Ign_delay_error = 0, red_size_of_mech = 105

T = 1050, P = 400000, phi = 1.5, Tmax_r = 2792.65, ignd_r = 0.127, T_max = 2792.38, ign_delay = 0.127, T_max_error = 0.00980892, Ign_delay_error = 0, red_size_of_mech = 106

T = 1050, P = 500000, phi = 0.5, Tmax_r = 2362.09, ignd_r = 0.06, T_max = 2363.14, ign_delay = 0.06, T_max_error = 0.0447686, Ign_delay_error = 0, red_size_of_mech = 89

T = 1050, P = 500000, phi = 1, Tmax_r = 2913.39, ignd_r = 0.079, T_max = 2913.83, ign_delay = 0.079, T_max_error = 0.0150559, Ign_delay_error = 0, red_size_of_mech = 159

T = 1050, P = 500000, phi = 1.5, Tmax_r = 2800.22, ignd_r = 0.098, T_max = 2800.1, ign_delay = 0.098, T_max_error = 0.00437602, Ign_delay_error = 0, red_size_of_mech = 106

T = 1150, P = 300000, phi = 0.5, Tmax_r = 2428.12, ignd_r = 0.023, T_max = 2428.66, ign_delay = 0.023, T_max_error = 0.0225851, Ign_delay_error = 0, red_size_of_mech = 65

T = 1150, P = 300000, phi = 1, Tmax_r = 2906.9, ignd_r = 0.032, T_max = 2907.56, ign_delay = 0.032, T_max_error = 0.0226661, Ign_delay_error = 0, red_size_of_mech = 159

T = 1150, P = 300000, phi = 1.5, Tmax_r = 2836.81, ignd_r = 0.04, T_max = 2836.37, ign_delay = 0.04, T_max_error = 0.0153953, Ign_delay_error = 0, red_size_of_mech = 106

T = 1150, P = 400000, phi = 0.5, Tmax_r = 2430.22, ignd_r = 0.018, T_max = 2430.39, ign_delay = 0.018, T_max_error = 0.00724013, Ign_delay_error = 0, red_size_of_mech = 159

T = 1150, P = 400000, phi = 1, Tmax_r = 2913.11, ignd_r = 0.024, T_max = 2911.8, ign_delay = 0.024, T_max_error = 0.044659, Ign_delay_error = 0, red_size_of_mech = 105

T = 1150, P = 400000, phi = 1.5, Tmax_r = 2849.48, ignd_r = 0.03, T_max = 2849.47, ign_delay = 0.03, T_max_error = 0.000541003, Ign_delay_error = 0, red_size_of_mech = 106

T = 1150, P = 500000, phi = 0.5, Tmax_r = 2435, ignd_r = 0.014, T_max = 2435.08, ign_delay = 0.014, T_max_error = 0.00302839, Ign_delay_error = 0, red_size_of_mech = 159

T = 1150, P = 500000, phi = 1, Tmax_r = 2928.28, ignd_r = 0.019, T_max = 2927.91, ign_delay = 0.019, T_max_error = 0.0126998, Ign_delay_error = 0, red_size_of_mech = 105

T = 1150, P = 500000, phi = 1.5, Tmax_r = 2864.09, ignd_r = 0.023, T_max = 2863.87, ign_delay = 0.023, T_max_error = 0.00788527, Ign_delay_error = 0, red_size_of_mech = 172
[Finished in 3648.2s]

Case 1:

Case 2:

Case 3:

Case 4:

Case 5:

Case 6:

Case 7:

Case 8:

Case 9:

Case 10:

Case 11:

Case 12:

Case 13:

Case 14:

Case 15:

Case 16:

Case 16:

Case 17:

Case 18:

Case 19:

Case 20:

Case 21:

Case 22:

Case 23:

Case 24:

Case 25:

Case 26:

Case 27:

Inference:

All the above graphs for different combinations of state values show a similar trend. Both ignition delay and maximum temperature graphs show fluctuations of high magnitude initially. As more and more reactions are added, the oscillations dampen and the graphs approach towards the reference values. This happens because we have added reactions in such a way that the most sensitive reactions are stored first. Therefore, the values for the most sensitive reactions show a greater degree of fluctuation during the start and as more and more reactions are added the fluctuations of the curve reduce and they approach the reference values.

From the above graph, it can be concluded that a minimum of 170 reactions is needed to predict both ignition delay and maximum temperature values accurately. This because the fluctuations are very minuscule or are completely non-existent in all the cases at this point.

"""TRY AND EXCEPT RULE"""

				try :
					if (i >= size) :
						ignd_a1.append(ign_delay)
						Tmax_b1.append(T_max)
						size_of_mech.append(i)

					# calculating the error
					error_Tmax = np.abs(((Tmax_r - T_max)/Tmax_r)*100) 
					error_ignd = np.abs(((ignd_r - ign_delay)/ignd_r)*100)

					# printing the values of phi, temperature etc. when the tolerance value is satisfied
					if (error_Tmax < tol_Tmax and error_ignd < tol_ignd and c3 is True):
						print('T = %g, P = %g, phi = %g, Tmax_r = %g, ignd_r = %g, T_max = %g, ign_delay = %g, T_max_error = %g, Ign_delay_error = %g, red_size_of_mech = %g' %(T, P, phi, Tmax_r, ignd_r, T_max, ign_delay, error_Tmax, error_ignd, i))
						c3 = False
				except NameError:

					pass

During the simulation of the code, a name error arises for ignition delay. This error occurs because, during the start of the simulation, there aren't enough reactions to initiate the combustion process and hence a name error gets displayed but this does not mean that the logic behind the code is incorrect. Therefore, to overcome this problem, the try and exception statements are used where the try command lets the program run and when and the exception rule checks to see if an error stated under the exception rule occurs while running the program. If it does then, the exception rule lets the program run in spite of the error encountered. Eventually, when enough reactions have taken place, the name error automatically goes away but it is during the start of the simulation that this error occurs and the try and except rule helps us overcome that.

Therefore, a reaction reduction mechanism for the GRI 3.0 mechanism had been successfully carried out using the Cantera module in Python.