Reduction of Reaction Mechanism for Methane
Objective:
The aim of this project is to reduce the reaction mechanism size from that of the GRI 3.0 Mechanism for methane (`CH_4`) fuel with various state values to compute the maximum temperature and the ignition delay for the given tolerance criteria(1%) based on the temperature sensitivity.
Procedure:
The code is done along the following stages as shown below:
- An input.txt file is created to provide the required input data.
- Input data is read and temperature sensitivity is performed for each of the combinations of state values.
- For each of the combinations of state values, GRI 3.0 Mechanism is sorted from high sensitive to least sensitive reactions. For all the combinations a .csv file is created to store the sorted reaction numbers.
- Then for each combination of state values using the sorted reaction numbers, a new gas instance is created and reactions are passed in the descending order of sensitive reactions from 40 most sensitive to all the 325 reactions to determine the maximum temperature and ignition delay with the addition of each reaction to the mechanism. Finally, to obtain the reduced mechanism size within the given tolerance criteria.
- The same is plotted along with the reference values of the GRI 3.0 mechanism.
Solution:
Input File:
Fuel CH4
Air O2:1,N2:3.76
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
The input file is used to parse input values for the main code. All the lines are self-explanatory. End_time is the simulation time, Delta_t is the time step for simulation, Min_size is minimum mechanism size, Tol_igd and Tol_Tmax are the tolerance criteria for determining the ignition delay and maximum temperature, respectively.
Python Script to find sensitive reactions order for each combination of State values:
import os
import cantera as ct
import matplotlib.pyplot as plt
import numpy as np
# Extracting input data
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']]
t_end = [float(x) for x in a.input_dict['End_time']]
dt = [float(x) for x in a.input_dict['Delta_t']]
min_reactions = [float(x) for x in a.input_dict['Min_size']]
tol_igd = [float(x) for x in a.input_dict['Tol_igd']]
tol_Tmax = [float(x) for x in a.input_dict['Tol_Tmax']]
fuel = [x for x in a.input_dict['Fuel']]
t_size=np.size(T_val)*np.size(P_val)*np.size(Phi_val)
#--------------------------------------------------------------------------------------#
gas = ct.Solution('gri30.cti')
air ='O2:1, N2:3.76'
#listing all the equations for gas instance of Solution Class of gri30.xml
n_reactions = gas.n_reactions # recording total no. of reactions involved (i.e., 325)
sensitivity_store = np.zeros((n_reactions,t_size))
reactionNo_store = np.zeros((n_reactions,t_size))
cntr = 0;
for T in T_val:
for P in P_val:
for Phi in Phi_val:
sensitivity = [0]*n_reactions
reactions = gas.reaction_equations()
gas.TP = T, P*1e5
gas.set_equivalence_ratio(Phi, fuel[0], air)
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[0]+dt[0],dt[0]):
sim.advance(t)
for i in range (gas.n_reactions):
S = sim.sensitivity('temperature', i)
if np.abs(S)>np.abs(sensitivity[i]):
sensitivity[i] = S
reaction_number = np.arange(1,n_reactions+1) # creating an array for reaction numbers
abs_sensitivity = np.abs(sensitivity) # calculating the absolute value of sensitivity
reactions = [a for b, a in sorted (zip(abs_sensitivity,reactions), reverse = True)]
reaction_number = [a for b, a in sorted (zip(abs_sensitivity,reaction_number), reverse = True)]
sensitivity = [a for b, a in sorted (zip(abs_sensitivity,sensitivity), reverse = True)]
sensitivity_store[:,cntr]=sensitivity
reactionNo_store[:,cntr]=reaction_number
cntr = cntr+1;
np.savetxt('f_array.csv',sensitivity_store, delimiter=',', fmt = '%.15f')
np.savetxt('f_array1.csv',reactionNo_store, delimiter=',', fmt = '%d')At first, the input file is extracted and sensitivity analysis is performed for each of the combinations of state values. Then, the reactions are sorted from highest to least sensitive reactions w.r.t to temperature. All the sorted reaction numbers for each combination are stored in the f_array1.csv file. Detail explanation to determine the sensitive reactions can be obtained from the following link: <https://skill-lync.com/projects/Sensitivity-Analysis-for-Combustion-of-Methane-13366>
Python Script to determine the reduced mechanism size:
# function to determine the reduced mechanism size
def reduction_mechanism(T,P,Phi,fuel,air,dt,t_end,min_size,Tol_igd, Tol_Tmax,case,reaction_numbers_sorted):
# inputs(temperature, pressure, fuel, air, time step, final time, minimum size, tolerance criteria ignition delay, tolerance for temperature, case No., reaction numbers)
import os
import cantera as ct
import matplotlib.pyplot as plt
import numpy as np
reactions_sorted = []
ign_delay_arr = []
T_max_arr = []
mech_size = []
T_ign = T+400
gas = ct.Solution('gri30.cti')
gas.TP = T,P*1e5 # Pressure in bar
gas.set_equivalence_ratio(Phi, fuel, air)
spec = ct.Species.listFromFile('gri30.cti')
r=ct.IdealGasReactor(gas)
sim = ct.ReactorNet([r])
T_max_ref =0
check = 0
# To determine the maximum temperature and ignition delay for GRI 3.0 mechanism
for t in np.arange(0,t_end+dt,dt):
sim.advance(t)
if (gas.T > T_ign and check == 0):
ign_delay_ref = t
check+=1
if (gas.T > T_max_ref):
T_max_ref = gas.T
# To find the reduced mechanism size
const=True
for i in range (gas.n_reactions):
reactions_sorted.append(gas.reaction(int(reaction_numbers_sorted[i])-1)) # to pass reactions one by one
gas1 = ct.Solution(thermo = 'IdealGas', kinetics = 'GasKinetics', species = spec, reactions = reactions_sorted) # creating gas instance
gas1.TP = T,P*1e5
gas1.set_equivalence_ratio(Phi, fuel, air)
r1=ct.IdealGasReactor(gas1)
sim1 = ct.ReactorNet([r1])
Flag = 0
ign_delay=0
T_max=0
for t in np.arange(0,t_end+dt,dt):
sim1.advance(t)
if (gas1.T > T_ign and Flag == 0):
ign_delay = t
Flag+=1
if (gas1.T > T_max):
T_max = gas1.T
try: # appending the values of ignition delay, mechanism size, maximum temperature
if(i>=min_size):
ign_delay_arr.append(ign_delay)
T_max_arr.append(T_max)
mech_size.append(i)
except NameError: # passing the NameError
pass
# 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
if (err_T_max < Tol_Tmax and err_ign_delay < Tol_igd and const is True):
red_mech_size=i+1
T_record=T_max
ign_record=ign_delay
const = False
T_max_ref=[T_max_ref]*len(T_max_arr) #increasing the size to match the lengths to ease in plot
ign_delay_ref=[ign_delay_ref]*len(ign_delay_arr)
#Plotting Maximum Temperature
plt.subplot(2,1,1)
plt.plot(mech_size,T_max_arr, color = 'red', label = 'Calculated Max Temperature')
plt.plot(mech_size,T_max_ref, color = 'blue', linestyle = 'dashed', label = 'Reference Max Temperature= '+str(format(max(T_max_ref),'0.4f'))+' [K]')
plt.xlabel('No. of sorted reactions most sensitive to Temperature')
plt.ylabel('Maximum Temperature [K]')
plt.legend(loc = 'best')
titlename='Case : '+str(case)+'\n State Variables T : '+str(T)+' [K], P : '+str(P)+' [bar], $\phi$ : '+str(Phi)
plt.title(titlename)
#Plotting ignition delay
plt.subplot(2,1,2)
plt.plot(mech_size,ign_delay_arr, color = 'green', label = 'Calculated Ignition Delay')
plt.plot(mech_size,ign_delay_ref, color = 'black', linestyle = 'dashed', label = 'Reference Ignition Delay = '+str(format(max(ign_delay_ref),'0.4f'))+' sec')
plt.xlabel('No. of sorted reactions most sensitive to Temperature')
plt.ylabel('Ignition Delay [sec]')
plt.legend(loc = 'best')
plt.tight_layout()
manager = plt.get_current_fig_manager() # to display images in full window
manager.resize(*manager.window.maxsize())
figname='Case_No_'+str(case)+'.png' # image name to saved as
plt.savefig(figname) # saving image
plt.close()
list1 = [case,T,P,Phi,T_max_ref[0],ign_delay_ref[0],T_record,ign_record,red_mech_size] # returning the values for comparison and write in separate CSV
return list1
Here for each combination of state values using the sorted reaction numbers, a new gas instance is created and reactions are passed in the descending order of sensitive reactions from 40 most sensitive to all the 325 reactions to determine the maximum temperature and ignition delay with the addition of each reaction to the mechanism. Finally, to obtain the reduced mechanism size within the given tolerance criteria.
To ease the formation of code, this is written as a function wherein state values, tolerance criteria, and all the necessary conditions are passed which returns the list of required values viz., maximum temperature, ignition delays, and images are stored in the folder for that particular case. Detail explanation to determine the auto-ignition delay and maximum temperatures can be obtained from <https://skill-lync.com/projects/Autoignition-delay-calculation-using-Cantera-60447>
Name Error:
As we are using the minimum size of 40 reactions, at times there are chances that Cantera may come across errors due to an insufficient number of reactions to obtain the solution. So we pass this error by using the exception of NameError.
Strictly speaking in the above code we are passing reactions from 1 to 325 of the sorted reactions, but we are recording only values for more than 40 reactions. In order to avoid the error from 1 to 40 reactions due to an insufficient number of reactions in mechanism size, we use exception and put the criteria of tolerance to obtain the solution within the tolerance by passing the rest of the reactions to determine the reduced size of the mechanism for the given tolerance criteria.
Python Script to call the function of reduction mechanism and writing the results in CSV file:
import csv
# reading the f_array1.csv file where all the sorted reaction numbers are stored.
def readcsv(filename):
ifile = open(filename, "rU")
reader = csv.reader(ifile, delimiter=",")
rownum = 0
a = []
for row in reader:
a.append (row)
rownum += 1
ifile.close()
return a
reaction_numbers_sorted1 = readcsv('f_array1.csv')
# Extracting input data
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 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']]
t_end = [float(x) for x in a.input_dict['End_time']]
dt = [float(x) for x in a.input_dict['Delta_t']]
min_size = [float(x) for x in a.input_dict['Min_size']]
Tol_igd = [float(x) for x in a.input_dict['Tol_igd']]
Tol_Tmax = [float(x) for x in a.input_dict['Tol_Tmax']]
fuel = [x for x in a.input_dict['Fuel']]
air = [x for x in a.input_dict['Air']]
case=0
# writing the csv file
with open('data.csv','w',newline='') as file:
writer = csv.writer(file)
writer.writerow(["Case No.","Temperature [K]", "Pressure [bar]", "Equivalence Ratio", "Maximum Temperature [K]", "Ignition Delay [s]", "Maximum Temperature [K]", "Ignition Delay [s]","Mechanism Size"])
for T in T_val:
for P in P_val:
for Phi in Phi_val:
reaction_numbers_sorted = [ x[case] for x in reaction_numbers_sorted1]
list2=reduction_mechanism(T,P,Phi,fuel[0],air[0],dt[0],t_end[0],min_size[0],Tol_igd[0],Tol_Tmax[0],case+1,reaction_numbers_sorted)
with open('data.csv','a',newline='') as fd: # appending the data of each combination
wr = csv.writer(fd,dialect='excel')
wr.writerow(list2)
case+=1
Here, the reduction mechanism function is called by passing each combination of state values and corresponding sorted reactions. For, each call of a function, an image is stored and the list is returned which is written in a data.csv file which is used to compare the results between different state values.
Results:
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
Comparison:
| Case No. | State Values | Reference Values with GRI 3.0 Mechanism | Computed with the tolerance of 1% | |||||
| Temperature [k] | Pressure [bar] | Equivalence Ratio | Maximum Temperature [K] | Ignition Delay [s] | Maximum Temperature [K] | Ignition Delay [s] | Mechanism Size | |
| 1 | 950 | 3 | 0.5 | 2283.660196 | 0.6 | 2283.635756 | 0.6 | 163 |
| 2 | 950 | 3 | 1 | 2827.509284 | 0.77 | 2827.780958 | 0.77 | 208 |
| 3 | 950 | 3 | 1.5 | 2724.848449 | 0.94 | 2724.870796 | 0.94 | 177 |
| 4 | 950 | 4 | 0.5 | 2286.378274 | 0.446 | 2286.37934 | 0.446 | 162 |
| 5 | 950 | 4 | 1 | 2854.329177 | 0.555 | 2854.111926 | 0.555 | 163 |
| 6 | 950 | 4 | 1.5 | 2733.608694 | 0.669 | 2733.76574 | 0.669 | 167 |
| 7 | 950 | 5 | 0.5 | 2288.300249 | 0.356 | 2288.301762 | 0.356 | 161 |
| 8 | 950 | 5 | 1 | 2859.971627 | 0.432 | 2860.070123 | 0.432 | 191 |
| 9 | 950 | 5 | 1.5 | 2739.985986 | 0.515 | 2740.223174 | 0.515 | 166 |
| 10 | 1050 | 3 | 0.5 | 2356.736154 | 0.103 | 2356.738104 | 0.103 | 160 |
| 11 | 1050 | 3 | 1 | 2862.492617 | 0.142 | 2862.462873 | 0.142 | 162 |
| 12 | 1050 | 3 | 1.5 | 2783.643146 | 0.176 | 2783.813168 | 0.176 | 160 |
| 13 | 1050 | 4 | 0.5 | 2359.697707 | 0.076 | 2359.703016 | 0.076 | 164 |
| 14 | 1050 | 4 | 1 | 2878.940846 | 0.103 | 2878.913899 | 0.103 | 141 |
| 15 | 1050 | 4 | 1.5 | 2792.65345 | 0.127 | 2792.37427 | 0.127 | 130 |
| 16 | 1050 | 5 | 0.5 | 2362.087139 | 0.06 | 2362.091354 | 0.06 | 167 |
| 17 | 1050 | 5 | 1 | 2913.388054 | 0.079 | 2913.222922 | 0.079 | 171 |
| 18 | 1050 | 5 | 1.5 | 2800.223274 | 0.098 | 2800.095108 | 0.098 | 133 |
| 19 | 1150 | 3 | 0.5 | 2428.116286 | 0.023 | 2428.272013 | 0.023 | 50 |
| 20 | 1150 | 3 | 1 | 2906.89638 | 0.032 | 2906.746901 | 0.032 | 166 |
| 21 | 1150 | 3 | 1.5 | 2836.809114 | 0.04 | 2837.036113 | 0.04 | 162 |
| 22 | 1150 | 4 | 0.5 | 2430.218819 | 0.018 | 2430.42293 | 0.018 | 147 |
| 23 | 1150 | 4 | 1 | 2913.105093 | 0.024 | 2913.000581 | 0.024 | 167 |
| 24 | 1150 | 4 | 1.5 | 2849.483198 | 0.03 | 2849.677225 | 0.03 | 156 |
| 25 | 1150 | 5 | 0.5 | 2435.003904 | 0.014 | 2434.918546 | 0.014 | 69 |
| 26 | 1150 | 5 | 1 | 2928.281419 | 0.019 | 2928.294251 | 0.019 | 161 |
| 27 | 1150 | 5 | 1.5 | 2864.091912 | 0.023 | 2864.364627 | 0.023 | 197 |
From the above figures, it is clear that for each of the state value combinations, there is higher inaccuracy with less number of reactions, and the accuracy increases upon adding the reactions to the mechanism size. From our analysis, it can be found that for our state values with an approximate size of 156 reactions sorted in the order of highest sensitivity to least w.r.t temperature, we can obtain the solution with 1% accuracy in determining the ignition delay and maximum temperature.
Future Suggestions:
For cases 19 & 25, it can be observed that the reduced mechanism size is 50 & 69, respectively. These are the odd one outs in comparison to other mechanism sizes. It is due to the fact that upon adding exact those numbers of reactions for that particular combination of state values, we get the result within tolerance which is evident from the figures of cases 19 & 25. Upon adding further reactions the error goes beyond 1%. In order to avoid such kind of errors, we need to put a check such that upon adding further reactions to the mechanism size, the error for at least next 25 to 30 additional reactions lies within the tolerance criteria. Then, that particular reduced size can be taken into consideration to avoid the above discrepancy.
NB: The codes and data files can be found from the link below:
https://drive.google.com/drive/folders/1ytii4CvNHnhYQVy-O0jwGSuEySW2YrGf?usp=sharing