Reaction Reduction mechanism for GRI 3.0 mechanism using Python and Cantera

Objective: Performing simulation on Reaction Reduction Mechanism for the GRI 3.0 mechanism

1. Run a reference ignition delay and maximum temperature analysis for the unlowered function for a time duration
2. Consequently, a temperature sensitivity analysis should be performed on all system reactions and the reactions should be sorted by the most sensitive to least sensitive
3. These sorted reactions should be added to a new mechanism file one by one starting from a fixed number of reactions (40)
4. Now a graph will display the reduced ignition delay and maximum temperature of the system while the reference ignition delay and temperature are marked as a reference label.
5. The value of Ignition delay and Maximum temperature on the Y-axis must be present, while the number of reactions on the X-axis must be recorded in the mechanism. The code is written for other state value combinations automatically

The combinations of the state values are:

Temperature = [950, 1050, 1150] Kelvin

Pressure = [3, 4, 5] bar

Equivalence ratio = [0.5, 1, 1.5]

Input.txt file:

Fuel CH4

Temperature [950 1050 1150]

Pressure [3 4 5]

Phi [0.5 1 1.5]

End_time 5

Delta_t 0.001

Min_size 40

Tol_igd 0.01

Tol_Tmax 0.01

Detailed chemical mechanism:

1. The chemical mechanism is theoretical speculation with tries at every stage open overall chemical reaction to describe in detail
2. Step by step reaction mechanism is an important reaction sequence that enables us in detail various stages of a chemical reaction a complete chemical change occurs
3. The inter media reactions are chemical species, typically unstable and short lived that are not reactants or product of the overall chemical reaction but are temporary product or reactants in the reaction steps of the mechanism
4. In order to resolve the exactness opti chemical Kinetic involved in this reaction accurate chemical reaction mechanism is needed
5. It contains many hydrocarbons, of course. this means that thousand and hundreds of species are involved
6. It potentially costs great deal of processing time. Computational time (N^3) of species actually engaged in process responses will be cube cycle (N^3)
7. Now comes the methodology of the reduction of the reaction process, which actually reduces computer time
8. The reduction mechanism technique is, it ignores some reactions that don’t have considerable effect on the mechanism

Procedure:

1. The numbers as input to the problem are extracted from the input file with the FileReader class and read the class to read the file
2. The gas objet is given as an instance of Solution from the Cantera
3. IdealGasReactor is built for combustion of gas object under a particular state values for Temperature, pressure and equivalence ratio
4. This is called a Zero-D reactor in the above simulation and determined maximum temperature and time of ignition for each combination of states
5. The unreduced processes in each reaction are ordered on the basis of the absolute sensitivities with respect to temperature
6. After the responses are sorted for each event, single reactions are added to the mechanism and the maximum temperature and ignition delay are measured for the reduced mechanism
7. The ignition delay cannot be calculated if the reactions are not happening due to less no. of reactions. To avoid this error called ‘NameError’, we are using the ‘try’ and ‘except’ classes.
8. The error tolerance are set to test if measured values are obtained within reference tolerance limit
9. When the value is within the limit, the number of reactions are recorded as reactions for the process of reduction. Calculations is done for all reactions for all combinations of state values and the results are computed

Python code:

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

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 doesnt 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_Tmax=float(x)
for x in a.input_dict['Tol_igd']:tol_ign_delay=float(x)
for x in a.input_dict['Fuel']:fuel=x

S=[]
S_max0=[]
R_index=[]
const1=True
n=0
air='O2:1,N2:3.76'
gas=ct.Solution('gri30.cti')
S_max=[0]*gas.n_reactions

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)

			#Setting up the reactor and reactor network
			r = ct.IdealGasReactor(gas)
			sim = ct.ReactorNet([r])

			#adding the reactions to sensitivity class to calculate corresponding sensitivity
			for i in range (0, gas.n_reactions):
				r.add_sensitivity_reaction(i)

			#setting tolerances for molar concentration
			sim.rtol = 1e-6
			sim.atol = 1e-15

			#setting tolerances for sensitivity coefficients
			sim.rtol_sensitivity = 1e-6
			sim.atol_sensitivity = 1e-6

			#Time integration of reactor network and storing properties at each instance
			for t in np.arange(0,t_end+dt,dt): # means t starts from 0 and ends at 2e-3 with time step size of 5e-6
				sim.advance(t)

				#Calculating the sensitivity and Maximum sensitivity wrt temp
				for j in range(0, gas.n_reactions):
					S=sim.sensitivity(2,j)  # Storing sensitivity as 's' for all reactions (n_reactions) at each time step

					if np.abs(S)>np.abs(S_max[j]): # Comparing sensitivity of each reaction obtained in new time step with previous time step and storing it in sensitivity matrix
						 S_max[j]=np.abs(S)

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

			S_sorted = sorted(S_max0, key=lambda x: x[0], reverse = True)


			for i in range (0, gas.n_reactions):
				if const1 is True:
					R_index.append([])
				R_index[i].append(S_sorted[i][1]) 			

			const1 = False
			n=n+1
			S_max0.clear()
		
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')

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:
			T_ign=T+400
			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

			err_T_max=100
			err_ign_delay=100
			i=0
			const4=True

			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_end+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	
				i=i+1

				try:
					if (i>=min_size):
						ign_delay_arr.append(ign_delay)
						T_max_arr.append(T_max)
						mech_size.append(i)

					err_T_max=np.abs(((T_max_ref-T_max)/T_max_ref)*100)
					err_ign_delay=	np.abs(((ign_delay_ref-ign_delay)/ign_delay_ref)*100)

					if (err_T_max

Results:

1. Temperature = 950 K,Pressure= 3 bar,Equivalence ratio= 0.5

2. Temperature = 950 K,Pressure= 3 bar,Equivalence ratio= 1

3. Temperature = 950 K,Pressure= 3 bar,Equivalence ratio= 1.5

4. Temperature = 950 K,Pressure= 4 bar,Equivalence ratio= 0.5

5. Temperature = 950 K,Pressure= 4 bar,Equivalence ratio= 1

6. Temperature = 950 K,Pressure= 4 bar,Equivalence ratio= 1.5

7.Temperature = 950 K,Pressure= 5 bar,Equivalence ratio= 0.5

8.Temperature = 950 K,Pressure= 5 bar,Equivalence ratio= 1

9. Temperature = 950 K,Pressure= 5 bar,Equivalence ratio= 1.5

10. Temperature = 1050 K,Pressure= 3 bar,Equivalence ratio= 0.5

11. Temperature = 1050 K,Pressure= 3 bar,Equivalence ratio= 1

12. Temperature = 1050 K,Pressure= 3 bar,Equivalence ratio= 1.5

13. Temperature = 1050 K,Pressure= 4 bar,Equivalence ratio= 0.5

14. Temperature = 1050 K,Pressure= 4 bar,Equivalence ratio= 1

15. Temperature = 1050 K,Pressure= 4 bar,Equivalence ratio= 1.5

16. Temperature = 1050 K,Pressure= 5 bar,Equivalence ratio= 0.5

17. Temperature = 1050 K,Pressure= 5 bar,Equivalence ratio= 1

18. Temperature = 1050 K,Pressure= 5 bar,Equivalence ratio= 1.5

19. Temperature = 1150 K,Pressure= 3 bar,Equivalence ratio= 0.5

20. Temperature = 1150 K,Pressure= 3 bar,Equivalence ratio= 1

21. Temperature = 1150 K,Pressure= 3 bar,Equivalence ratio= 1.5

22. Temperature = 1150 K,Pressure= 4 bar,Equivalence ratio= 0.5

23. Temperature = 1150 K,Pressure= 4 bar,Equivalence ratio= 1

24. Temperature = 1150 K,Pressure= 4 bar,Equivalence ratio= 1.5

25. Temperature = 1150 K,Pressure= 5 bar,Equivalence ratio= 0.5

26. Temperature = 1150 K,Pressure= 5 bar,Equivalence ratio= 1

27. Temperature = 1150 K,Pressure= 5 bar,Equivalence ratio= 1.5

Conclusions:

We can see that the initial fluctuations in Maximum Temperature and Ignition Delay is large since the most sensitive reactions are used initially

After a specific number of reactions, the Maximum Temperature and Ignition delay plots converge towards their Reference values, which is obvious. Using this, we can determine the minimum number of most sensitive reactions withina reaction mechanism which gives same results as that of a reference values (obtained from unreduced reactio mechanism)