Reduced Reaction Mechanism for GRI3.0

Objective: To write a code to compile a reduced reaction mechanism for the GRI mechanism under different state values.

• To 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).
• To perform a sensitivity analysis with respect to temperature on all the reactions in the mechanism and sort the reactions from the most sensitive to least sensitive.
• Addition of the sorted reactions to a new mechanism file starting from a minimum of 40 reactions.
• To plot a graph containing the reduced mechanism ignition delay and compare it with reference ignition delay for number of reactions added to the new mechanism.

Introduction:

A reaction Mechanism is a detailed process by which chemical substances are transformed from one substance to another. Each step in a mechanism is called an elementary reaction.Reaction intermediates are formed in one step and consumed in another.

Detailed chemical mechanism are available for a large number of hydrocarbon fuels.These involve thousands of reactions and thousands of species.This can take a huge computation time. Mechanism reduction reduces these mechanism to a smaller size that can lead to an accurate result as produced by the original mechanism by utilising lesser number of reactions than those present in the original mechanism.

In this process the redundant reactions are removed for a specific range of conditions.

Methodology:

The following methodology is followed to reduce the GRI Mechanism for a specific range of operating conditions that are as follows:

Temperature =[950 1050 1150] Kelvin

Pressure=[3 4 5]bar

Equivalence Ratio=[0.5 1 1.5]

The Input conditions are fed using a text file. The file parser takes in the input parameters and reads them into arrays. The analysis is performed for a combination of the input conditions by using nested looping.

The sensitivity analysis is first performed with respect to temperature and the reactions are sorted accordingly in the decreasing order of their sensitivity. A base temperature and ignition delay reference is calculated by using all the reactions.

The next step is to reduce the mechanism. Startinf with 40 reactions the sorted reactions are added one by one into the mechanism until the error between the ignition delay and the Maximum Temperature of the new reduced mechanism and the reference values is below 5%.

The python script generated a plot for each case and also a text file is generated which gives the summary of the states and the number of reactions that are added in the reduced mechanism for each set of operating condition to obtain the error below 5%.

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

class File_reader():

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

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

	def read_file(self):
		if self.validate_file():
			for x in open(self.file_name,'r'):
				self.input_dict[x.split()[0]]=x.split()[1:]
		else:
			print('No File Found!!!')


f=File_reader("input.txt")
f.read_file()
Temp=[float(T) for T in f.input_dict['Temperature']]
Pressure=[int(T) for T in f.input_dict['Pressure']]
E_R=[float(T) for T in f.input_dict['Phi']]
Temp_tol=[float(T) for T in f.input_dict['Temp_tol']]
time_tol=[float(T) for T in f.input_dict['time_tol']]


Pressure=ct.one_atm*np.array(Pressure)


t_sim=5
dt=0.001
gas=ct.Solution('gri30.xml')
spec=ct.Species.listFromFile('gri30.xml')
case=1
R_tab=[]
R_red=[]
mech_size=[]
delay=[]
temp_max=[]

for T in Temp:
	for P in Pressure:
		for phi in E_R:

			print(case)
			gas.TP=T,P


			gas.set_equivalence_ratio(phi,'CH4','O2:1,N2:3.76')

#reactor and reactor network
			r=ct.IdealGasReactor(gas,name='R1')
			sim=ct.ReactorNet([r])


			n_reac= gas.n_reactions
			S_max=[0]*n_reac

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

			#tolerences
			sim.rtol=1e-6
			sim.atol=1e-16

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

			reac_num=np.arange(0,n_reac)
			for t in np.arange(0,t_sim+dt,dt):
				sim.advance(t)
	


				for j in range(n_reac):
					s=sim.sensitivity('temperature',j)
				#plt.plot(t,s)
					if np.abs(s)>np.abs(S_max[j]):
						S_max[j]=s

			S_abs=np.abs(S_max)

			reac_num=[a for b,a in sorted(zip(S_abs,reac_num),reverse=True)]

			#print(reac_num)
			case=case+1
			# Creating Table
			for i in range(0,n_reac):
				R_tab.append([])
				R_tab[i].append(reac_num[i])
			#print(R_tab)	

#Creating rxn list
for i in range(0,n_reac):
	for j in range(0,case-1):
		if R_tab[i][j] not in R_red:
			R_red.append(R_tab[i][j])
#print(R_red)

# Ignition delay and formation of new mechanism
wr=open("Output.txt","w+")

case=1
error_T=0
error_ign=0
n_red=40
for T in Temp:
	for P in Pressure:
		for phi in E_R:
			T_ign=T+400
			Red_Mech=[]
			flag2=1
			mech_size.clear()
			delay.clear()
			temp_max.clear()
			
			print(case)

			for i in range(0,n_reac):
				Red_Mech.append(gas.reaction(R_red[i]))
				T_ref=T
				T_max=0

				flag=1
				flag1=1

				gas_new=ct.Solution(thermo='IdealGas',kinetics='GasKinetics',species=spec,reactions=Red_Mech)
				gas_new.TP=T,P
				gas_new.set_equivalence_ratio(phi,'CH4','O2:1,N2:3.76')
				gas.TP=T,P
				gas.set_equivalence_ratio(phi,'CH4','O2:1,N2:3.76')
				r1=ct.IdealGasReactor(gas_new)
				r2=ct.IdealGasReactor(gas)
				sim1=ct.ReactorNet([r1])
				sim2=ct.ReactorNet([r2])
				

				for t in np.arange(0,t_sim+0.001,0.001):
					sim2.advance(t)
					if gas.T>=T_ref:
						T_ref=gas.T
					if gas.T>=T_ign and flag:
						t_ref=t
						flag=0
					

					sim1.advance(t)
					if gas_new.T>=T_max:
						T_max=gas_new.T
					if gas_new.T>=T_ign and flag1:
						t_delay=t
						flag1=0
					#elif gas_new.T=n_red):
						
						delay.append(t_delay)
						temp_max.append(T_max)
						mech_size.append(i)

					error_ign=np.abs(((t_delay-t_ref )/t_ref)*100)
					
					error_T=np.abs(((T_ref-T_max)/T_ref)*100)
					
					if (error_ign<time_tol) and (error_T<Temp_tol) and flag2:
						req_eqns=i
						output1='\nCase : '+str(case)+'\n State Variables T : '+str(T)+'[K] P : '+str(P/(ct.one_atm))+'[bar] Phi : '+str(phi)+'\n Minimum Number of Reactions to be added for 5 percent error :'+ str(req_eqns)+'\n'
						output2='Reference Ignition Delay:'+str(t_ref)+'[K] Reference Maximum Temperature'+str(T_ref)+'\n'
						wr.write(output1)
						wr.write(output2)

						flag2=0


				except NameError:
					
					pass
				
			t_ref=[t_ref]*len(delay)
			T_ref=[T_ref]*len(temp_max)		
			plt.figure(case)
			plt.subplot(2,1,1)
			ign_delay_label='Reference Ignition Delay'+str(t_ref[0])+'(ms)'
			plt.plot(mech_size,delay,color='green',label='IgnitionDelay for reduced Mechanism')
			plt.plot(mech_size,t_ref,linestyle=':',color='purple',label=ign_delay_label)
			plt.grid(True)
			plt.xlabel('Number of Reactions')
			plt.ylabel('Ignition Delay (ms)')
			plt.legend(loc='best')
			titlename='Case : '+str(case)+'\n State Variables T : '+str(T)+'[K] P : '+str(P/(ct.one_atm))+'[bar] Phi : '+str(phi)
			plt.title(titlename)
			plt.subplot(2,1,2)
			ref_temp_label='Reference Maximum Temperature'+str(T_ref[0])+'[K]'
			plt.plot(mech_size,temp_max,color='red',label='Maximum Temperature for reduced Mechanism')
			plt.plot(mech_size,T_ref,linestyle=':',color='blue',label=ref_temp_label)
			plt.grid(True)
			plt.xlabel('Number of Reactions')
			plt.ylabel('Maximum Temperature(K)')
			plt.legend(loc='best')
			plt.tight_layout()
			figname='Case_No_'+str(case)+'.png'
			
			plt.savefig(figname)
			plt.clf
			
			case=case+1	

wr.close()			

A NameError arises initially when the number of reactions added in the reaction mechanism are low and are insufficient to cause an ignition hence the ignition delay time cannot be recorded.

This error is fixed using the try and except command.

Result:

The output file generated can be viewed using the link below.

https://drive.google.com/open?id=1qeiwoLpmharTjDO69G0KrL00MRq8kNqk

The following images contain the plots of the ignition delay and the maximum temperature for each case.

From the plots above it can be inferred that an average 150 reactions are needed to be added to be able to predict the Maximum Temperature are the ignition delay accurately. It can be seen that in some cases the Maximum temperature is predicted with a smaller amount of reactions. However this is not stable.The reactions are added in the decreasing order of their sensitivities and hence the fluctuations are higher when the initial reactions are added and the flustuations die down as the later reactions are added as their sensitivites are low and the effects are negligible.