REACTION MECHANISM REDUCTION FOR THE GRI 3.0 MECHANISM WITH DIFFERENT STATE VARIABLE COMBINATIONS USING PYTHON AND CANTERA

PROBLEM STATEMENT:

To write a Reaction Reduction Mechanism code for the GRI3.0 mechanism,

code should do the following,

1. First, it 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 with respect to temperature must be run on all the reactions in the mechanism and must sort the reactions from the most sensitive to least sensitive
3. Then, 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
6. Now, same steps are to be repeated for various state values

INTRODUCTION:

A reaction mechanism is a step-by-step sequence of elementary reactions by which overall chemical change occurs. It helps to study the different stages of a chemical reaction in detail. Obtaining an accurate mechanism for a given set of conditions is very challenging. Different conditions affect the reaction path from the reactants to the final products. The reaction path includes the intermediate reactions that occur between intermediate species formed. To accurately solve the chemical kinetics involved in a reaction, using an accurate chemistry reaction mechanism is essential. This reaction mechanism can contain hundreds of reactions and thousands of species being formed. The computational cost for solving a mechanism is proportional to `N^3`, where `N` is the number of species involved in the reaction mechanism. Mechanism reduction can help accelerate this computational analysis. Mechanism reduction is the process of neglecting some reactions which do not have a considerable impact on the mechanism. The impact of a reaction is compared through its sensitivity with respect to temperature. For this study, we implement the mechanism reduction for the GRI 3.0 mechanism with different combinations of the state variables as:

Temperature: `950, 1050, 1150 K`

Pressure: `3.0, 4.0, 5.0` `b ar`

Equivalence Ratio: `0.5, 1.0, 1.5`

The steps in the process of implementation of the mechanism reduction are listed below:

• Extract the input variable data from an input file.
• Initiate gas as an instance of the Solution class from Cantera.
• Create IdealGasReactor and ReactorNet objects from the gas instance based on specified state variables.
• Simulate the 0-D reactor for the above reactor network and calculate the maximum temperature and ignition delay for each of the cases along with the reaction sensitivities to temperature.
• Sort the reactions based on the absolute values of the sensitivities for the unreduced mechanism in each case.
• Prioritize the reactions based on the sorted reactions from each case, such that the overall reaction order can be obtained based on all the cases.
• For each case, cumulatively add single reactions to the mechanism and calculate the maximum temperature and ignition delay for the reduced mechanism.
• If the reactions are not happening due to a lesser number of reactions, the ignition delay will not be calculated. Catch this exception using NameError.
• Check whether the calculated values fall within the tolerance limits of the reference values. If all the calculated values meet the tolerance criteria, store the number of reactions as the reactions for the reduced mechanism.
• Continue the calculations for all the reactions in all the cases and plot the results.

REACTION MECHANISM:

A reaction mechanism is a general process where all the reactions are listed until we reach the end of species where further reactions dont happen anymore. It gives a detail process of how the reactions undergo changes throughout until it reaches the final. Intermediate reactions are also formed during where they happen to occur in one reaction and get absorbed in further reactions as the process moves on. For large fuels like Hydrocarbon a detailed mechanism is enlisted for further reference to the end user.

MECHANISM REDUCTION:

Reducing the size of the reaction mechanism in the sense helps us understand the mechanism and work on it better. So a set of conditions are imposed that help us achieve the desired output. A large computational problem has several 1000's of reactions and output consists of 1000's of species. The cost of computing the results is also very high. So there is a requirement for a feasible solution to the problem. Therefore all the redundant reactions which do not have an impact on the mechanism are removed from our final reduced mechanisms.

METHODOLGY:

• Sensitivity Analysis is done in GRI 3.0 mechanism file and the reactions are sorted in the descending order from the order of the highest to the least sensitivity.
• Tailoring the mechanism for our requirement by providing the value of ignition delay in our calculations i.e 400 K. Only reaction which can pass through the above constraint are eligible to progress. So the reactions which have less to least sensitivities are removed from the mechanism.
• Ignition Delay and maximum temperature calculations are then performed while having two gas objects i.e. The first object has less number of reactions because it is part of the reduced mechanism and has 40 chemical reactions. The second gas object will have as usual 325 reactions in its mechanism

INPUT:

The following input is provided as a text file and is shown below:

Fuel		CH4
Temperature 950 1050 1150
Pressure	300000 400000 500000
Phi			0.5 1.0 1.5
End_time	5
Delta_t		0.001
Min_size	40
Tol_igd		0.01
Tol_Tmax	0.01

PYTHON CODE :

'''
Program to reduce reaction mechanism for the GRI 3.0 mechanism 
Author : Harshith Rajamahanti 

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

# Extracting input date from 'input_file' using user-defined class 'FileReader'
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 mixture 
air = 'O2:1, N2:3.76'

# creating gas object
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)

			# setting absolute and relative tolerances
			sim.rtol = 1e-6
			sim.atol = 1e-15
			sim.rtol_sensitivity = 1e-6
			sim.atol_sensitivity = 1e-6
		

			# 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))
			plt.plot(mech_size, ign_delay_arr, color = 'green')
			plt.axhline(y = ign_delay_ref, color = 'red', linestyle = 'dotted', label = 'Reference Ignition Delay')
			plt.xlabel('No. of reactions most sensitive to Temperature')
			plt.ylabel('Ignition Delay (s)')
			plt.legend(loc = 'best')

			plt.subplot(2,1,2)
			plt.plot(mech_size, T_max_arr, color = 'blue')
			plt.axhline(y = T_max_ref, color = 'red', linestyle = 'dotted', label = 'Reference Max Temperature')
			plt.xlabel('No. of reactions most sensitive to Temperature')
			plt.ylabel('Maximum Temperature (K)')
			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()

Result:

The plots for ignition delay and maximum temperature against the mechanism size for each state value combination are shown below;

1.

2.

3.

4.

5.

6.

7.

8.

9.

10.

11.

12.

13.

14.

15.

16.

17.

18.

19.

20.

21.

22.

23.

24.

25.

26.

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 dampening of fluctuations or oscillations can be explained as follows.

We are ordering the reactions in the decreasing order of sensitivities. And we know that the sensitivity may either be positive or negative. Higher the sensitivity, lower the ignition delay and higher the maximum temperature and vice versa. Therefore reactions being appended at the early stages have the highest sensitivity and hence greater influence on the ignition delay and maximum temperature. That is the reason why the magnitude of fluctuations is high initially. And as we append more of the less sensitive reactions, the average or the net sensitivity of the reactions gradually decreases resulting in lesser magnitudes of fluctuations and after a certain point, the graph attains a steady state.

Therefore the reduced mechanism size depends on the desired accuracy i.e., how close do we need the ignition delay and maximum temperature values to be with the ones obtained from the original unreduced GRI30 mechanism. For a tolerance of 0.01% for maximum temperature and ignition delay with the reference values, the reduced mechanism size satisfying the tolerance criteria for each state value combination was found to be as follows:

Use of Error exception

The NameError exception has been used in the code as shown in the snippet below;

The error exception can be implemented with the help of try-except statements as shown from line-205 through line-225. The try-except statements work like the if-else statements. Python first tries to execute the code under try: statement. And except: code is executed only if an error is encountered while trying to execute try: statement. The error to be caught is written alongside the except keyword as shown in the snippet.

Need for the error exception

The aim of the program is to determine the ignition delay and maximum temperature values for reduced mechanism sizes starting from through all the 325 reactions. This values can be computed only if the combustion is initiated and for combustion to be initiated there must be enough reactions in the reaction object (i.e. 35).

As shown in the code in line 192, the system temperature must rise above the ignition temperature, which is 400 more than the initial temperature in order to assign the value for the variable 'ign_delay'. Until the no. of reactions in the reaction object reaches 35, the system temperature does not surpass the ignition temperature and the if condition in the line 192 will remain false and the variable 'ign_delay' will not have been initialized. This will throw an error saying 'variable ign_delay not defined' on the very first iteration when calling the variable later to calculate errors and for printing results. This will cause the program to abort.

In order to command the program to keep running despite encountering the NameError, the try-except statements are necessary. And when the no. of reactions in the reaction object increases to 35 and more, the combustion occurs, the system temperature rises above the ignition temperature and the variable ign_delay is defined. So after 35 reactions, only the try: statement is executed.

CONCLUSION:

Thus Reaction Mechanism Reduction is performed in our GRI 3.0 file using CANTERA library, based on the input values and also the tolerances allotted for ignition delay and max temperatures as well.