Reaction Mechanism Reduction Project
PROBLEM STATEMENT:
Write a Reaction Reduction Mechanism code for the GRI3.0 mechanism,
Your code should do the following,
- 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)
- 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
- 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).
- 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)
- 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
- Now, repeat the same steps for various state values
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:
Temperature 950 1050 1150
Pressure 3 4 5
Phi 0.5 1 1.5CODE:
import os
import numpy as np
import cantera as ct
import sys
import matplotlib.pyplot as plt
# Part 1 - To allow the import of data from file
class Input_file_reader():
def __init__(self, file_name):
self.file_name = file_name
self.input_dict = {}
def validate_input_file(self):
return os.path.isfile(self.file_name)
def read(self):
if self.validate_input_file():
for line in open(self.file_name, 'r'):
self.input_dict[line.split()[0]] = line.split()[1:]
else:
print ("Validation Failed")
# Part 2 - Format in which the parameters are written for accepting the data text file
a = Input_file_reader("input.txt")
a.read()
T_val = [float(x) for x in a.input_dict['Temperature']] # Enter the 3 values of Temperature
P_val = [float(x) for x in a.input_dict['Pressure']] # Enter the 3 values of Pressure
phi_val = [float(x) for x in a.input_dict['Phi']] # Enter the 3 values of Equivalence Ratio
# Part 3 - Various Input Parameters for our Mechanism
dt = 0.001 # Time Step
t_end = 5 # seconds (End Time for the simulation)
Tmax_tol = 0.05 # Tolerance values for Maximum Temperature
IgDelay_tol = 0.05 # Tolerance values for Ignition Delay
Min_Size = 40 # The minumum size of the mechanism
# Part 4 - Preparing the reaction mechanism for sensitivity analysis
gas=ct.Solution('gri30.cti')
spec=ct.Species.listFromFile('gri30.xml')
case = 1
Reac_tab = [] # Reaction Table
Reac_red = [] # Reaction Reduction
P_val=ct.one_atm*np.array(P_val)
for T in T_val:
for P in P_val:
for phi in phi_val:
print(case)
gas.TPX = T, P, {'CH4':1,'O2':2/phi,'N2':7.52/phi}
r = ct.IdealGasConstPressureReactor(gas, name='R1')
sim = ct.ReactorNet([r])
num_reac = gas.n_reactions
sen_max = [0] * num_reac
for i in range(0,num_reac):
r.add_sensitivity_reaction(i)
sim_rtol = 1.0e-6
sim_atol = 1.0e-15
sim.rtol_sensitivity = 1.0e-6
sim.atol_sensitivity = 1.0e-6
total_reac = np.arange(0,num_reac)
for t in np.arange(0,t_end+dt,dt):
sim.advance(t)
for j in range(num_reac):
s = sim.sensitivity('temperature', j)
if np.abs(s)>np.abs(sen_max[j]):
sen_max[j] = s
sen_abs = np.abs(sen_max)
total_reac = [a for b, a in sorted(zip(sen_abs,total_reac), reverse = True)]
case=case+1
for i in range(0,num_reac):
Reac_tab.append([])
Reac_tab[i].append(total_reac[i])
print(Reac_tab)
for i in range(0,num_reac):
for j in range(0, case-1):
if Reac_tab[i][j] not in Reac_red:
Reac_red.append(Reac_tab[i][j])
print(Reac_red)
q=open("Output.txt","w+")
case = 1
T_error = 1e6
IgDelay_error = 1e6
Mech_Size = []
Time_Delay = []
Temperature_Max = []
for T in T_val:
for P in P_val:
for phi in phi_val:
T_ign = T + 400
Red_Mech = []
const4 = True
Mech_Size.clear()
Time_Delay.clear()
Temperature_Max.clear()
print(case)
for i in range(0,num_reac):
Red_Mech.append(gas.reaction(Reac_red[i]))
Tref = T
Tmax = 0
const2 = True
const3 = True
gas1=ct.Solution(thermo='IdealGas',kinetics='GasKinetics',species=spec,reactions=Red_Mech)
gas1.TPX = T, P, {'CH4':1,'O2':2/phi,'N2':7.52/phi}
gas2=ct.Solution('gri30.cti')
gas2.TPX = T, P, {'CH4':1,'O2':2/phi,'N2':7.52/phi}
r1 = ct.IdealGasReactor(gas1)
r2 = ct.IdealGasReactor(gas2)
sim1 = ct.ReactorNet([r1])
sim2 = ct.ReactorNet([r2])
for t in np.arange(0,t_end+dt,dt):
sim2.advance(t)
if (gas2.T >= Tref):
Tref = gas2.T
if (gas2.T >= T_ign) and const2 is True:
tref = t
const2 = False
sim1.advance(t)
if (gas1.T >= Tmax):
Tmax = gas1.T
if (gas1.T >= T_ign) and const3 is True:
tdelay = t
const3 = False
# Try and Except Rule
try:
if (i >= Min_Size):
Mech_Size.append(i)
Time_Delay.append(tdelay)
Temperature_Max.append(Tmax)
IgDelay_error = np.abs( ( ( tdelay - tref ) / tref ) * 100 )
T_error = np.abs( ( ( Tref - Tmax ) / Tref ) * 100 )
if ( IgDelay_error < IgDelay_tol ) and (T_error < Tmax_tol) and const4 is True:
Req_Eq = i
Result1='Case:'+str(case)+'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_Eq)+'n'
Result2='Reference Ignition Delay:'+str(tref)+'[K] Reference Maximum Temperature:'+str(Tref)+'n'
q.write(Result1)
q.write(Result2)
const4 = False
except NameError:
pass
tref = [tref]*len(Time_Delay)
Tref = [Tref]*len(Temperature_Max)
plt.figure(case)
plt.subplot(2,1,1)
label_igdelay='Reference Ignition Delay:'+str(tref[0])+'(ms)'
plt.plot(Mech_Size,Time_Delay,color='green',label='Ignition Delay for reduced Mechanism')
plt.plot(Mech_Size,tref,linestyle=':',color='purple',label=label_igdelay)
plt.grid(True)
plt.xlabel('Number of Reactions')
plt.ylabel('Ignition Delay (ms)')
plt.legend(loc='best')
titlename='Case : '+str(case)+' State Variables T : '+str(T)+'[K] P : '+str(P/(ct.one_atm))+'[bar] Phi : '+str(phi)
plt.title(titlename)
plt.subplot(2,1,2)
label_reftemp='Reference Maximum Temperature:'+str(Tref[0])+'[K]'
plt.plot(Mech_Size,Temperature_Max,color='red',label='Maximum Temperature for reduced Mechanism')
plt.plot(Mech_Size,Tref,linestyle=':',color='blue',label=label_reftemp)
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
q.close()OUTPUT:
Based on our manual interpretation we tried to verify the minimum reaction that are required to achieve the ignition delay and maximum Temperature and the following are tabulated below:
| PRESSURE(bar) | EQUIVALENCE RATIO(phi) | MINIMUM REACTIONS REQUIRED TO PREDICT IGNITION DELAY | MINIMUM REACTIONS REQUIRED TO PREDICT MAXIMUM TEMPERATURE |
| 3 | 0.5 | 65 | 190 |
| 3 | 1 | 70 | 190 |
| 3 | 1.5 | 105 | 190 |
| 4 | 0.5 | 60 | 190 |
| 4 | 1 | 68 | 190 |
| 4 | 1.5 | 110 | 190 |
| 5 | 0.5 | 77 | 190 |
| 5 | 1 | 63 | 195 |
| 5 | 1.5 | 105 | 190 |
| 3 | 0.5 | 85 | 190 |
| 3 | 1 | 72 | 190 |
| 3 | 1.5 | 110 | 190 |
| 4 | 0.5 | 75 | 190 |
| 4 | 1 | 60 | 190 |
| 4 | 1.5 | 70 | 190 |
| 5 | 0.5 | 63 | 195 |
| 5 | 1 | 65 | 190 |
| 5 | 1.5 | 110 | 260 |
| 3 | 0.5 | 60 | 195 |
| 3 | 1 | 60 | 190 |
| 3 | 1.5 | 75 | 260 |
| 4 | 0.5 | 64 | 190 |
| 4 | 1 | 30 | 78 |
| 4 | 1.5 | 61 | 265 |
| 5 | 0.5 | 40 | 190 |
| 5 | 1 | 40 | 80 |
| 5 | 1.5 | 45 | 240 |
From the table we see that a minimum of 60 on an average will be required for us to predict the Ignition Delay and close to 190 reactions are required to predict the Maximum Temperature from the mechanism based on our requirements.
# Try and Except Rule
try:
if (i >= Min_Size):
Mech_Size.append(i)
Time_Delay.append(tdelay)
Temperature_Max.append(Tmax)
IgDelay_error = np.abs( ( ( tdelay - tref ) / tref ) * 100 )
T_error = np.abs( ( ( Tref - Tmax ) / Tref ) * 100 )
if ( IgDelay_error < IgDelay_tol ) and (T_error < Tmax_tol) and const4 is True:
Req_Eq = i
Result1='Case:'+str(case)+'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_Eq)+'n'
Result2='Reference Ignition Delay:'+str(tref)+'[K] Reference Maximum Temperature:'+str(Tref)+'n'
q.write(Result1)
q.write(Result2)
const4 = False
except NameError:
passDuring the initial stages there are less number of reactions due to which it becomes impossible to initiate the combustion reaction from our GRI 3.0 mechanism. Eventhough the code is very apt for our simulation there occurs a problem of NameError because of the limited number of reactions in our mechanism. So we except that error and run our simulation code again and again. This error will vanish when we have limited number of reactions in place to perform the ignition delay and the Maximum Temperature calculations.
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.