Reaction mechanism reduction of GRI 30 mechanism
Objectives:
- Write a reaction mechanism reduction code in python starting from 40 to 325 most temperature sensitive reactions and calculate maximum temperature and ignition delay time for the reduced mechanism at varying values of temperature, pressure and equivalence ratio
Theory:
Reaction mechanism is a sequence of series of intermediary reactions, through which the overall chemical change occurs. A file which contains all the intermediary chemical reactions in sequence is a mechanism file. The rate of occurence of each reaction as well as the rate of change of each species is of primary interest in a combustion simulation. As the reactions progresses, the state variables of gas changes also.
A sudden increase in temperature of the gas signifies the start of combustion and the time taken for it to happen is called as ignition delay. For this case, the time taken to achieve a difference in temperature of 400 K, between current state and previous state, is taken as ignition delay, which is of interest in current project.
In a reaction mechanism file, some reactions progresses more rapidly than others with time, which implies that such reactions are more 'sensitive', which again are of interest. The sensitivity can be attributed to a reaction only with respect to a state variable like temperature, pressure etc. In this case the sensitivities are found out with respect to temperature.
When certain reactions are more sensitive than others, the effect of these reactions on the areas of interest like ignition delay, maximum temperature etc, is more. Which indicates that, for a reasonable tolerance, the system doesn't need to consider all the reactions in the mechanism file, which reduces load on the machine. The minimum number of reactions (which are sensitive corresponding to a certain state variable) that require to simulate the mechanism is to be found out.
Reducing the number of reactions to simulate the mechanism for a given tolerance, is called as reaction mechanism reduction.
Introduction:
In this project, gri30 reaction mechanism file is used, which is to be reduced by conisidering only a certain number of temperature sensitive reactions. Corresponding to temperature, the top most sensitive reactions are found out, they are sorted in the descending order of their sensitivities. The sorted reactions are then simulated and the ignition delay and maximum temperature (which are of interest here, hence temperature sensitiveness) are found out. Ignition delay and maximum temperature values are then compared with the actual values (considering all the reactions). If the error is more than the tolerance, more reactions are added. The process is repeated until the tolerance is matched.
By having a reduced reaction mechanism for simulation, significant load reduction on the machine can be achieved, as not all reactions need to be simulated.
In the current case, sensitivity of reactions corresponding to temperature are found out and the reaction reduction is done based on it.
Coding in Python:
- In this project, file parsing technique is used to read the values of temperature, pressure and equivalence ratio from a text file and are stored in lists created globally
- The logic followed is, first sensitivity analysis corresponding to temperature is carried out considering all reactions in GRI 30 reaction mechanism file (325). The sensitivity values are sorted in descending order and the corresponding reaction indices are also sorted, for all cases (27) and are stored in a list (this list contains 27 sub lists, with each sub list having length of 325)
- The reaction indices are stored in a new list in a way that top most sensitive reaction index in all 27 cases are compared stored in 'first come first served' basis, if a reaction index is repeated, it is skipped. The same is done for second most sensitive reaction index in all cases. This way the new list contains 325 reaction indices, which are the most sensitive reactions
- Now, a gas object is created considering only the first reaction from the list of most senitive reactions and try to calculate maximum temperature and ignition delay, in case it fails, the second reaction from the most sensitive reactions is added to the first reaction and the process is repeated. If it fails too, then one more reaction is added to the list
- The process is repeated for all 325 reactions in the above list, the values of maximum temperature and ignition delay time of the reduced reactions are stored in a list, also maximum temperature and ignition delay time for the corresponding case is calculated too, in the same loop, considering all 325 reactions. These values are used to compare with the corresponding value of reduced mechanism
- The maximum temperature values and corresponding number of reactions, starting from a value of 40, are stored in a list and are plotted, so a total of 27 plot will be obtained corresponding to each case
The data in the text file from which the data is parsed is,
Temperature 950 1050 1150
Pressure 3 4 5
Phi 0.5 1 1.5
Python code to reduce the size of mechanism is
import cantera as ct
import numpy as np
import matplotlib.pyplot as plt
import os
class FileReader:
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('No file under such name')
gas = ct.Solution('gri30.xml')
a = FileReader('Input.txt')
a.read()
Temp = [float(x) for x in a.input_dict['Temperature']]
Press = [float(x) for x in a.input_dict['Pressure']]
Phi = [float(x) for x in a.input_dict['Phi']]
case = 0
reaction_index = []
for i in range(len(Temp)):
for j in range(len(Press)):
for k in range(len(Phi)):
print(case)
gas.TPX = Temp[i], Press[j]*101325, {'CH4':1, 'O2':2/Phi[k], 'N2':7.52/Phi[k]}
r = ct.IdealGasConstPressureReactor(gas)
sim = ct.ReactorNet([r])
for l in range(gas.n_reactions):
r.add_sensitivity_reaction(l)
reaction_index.append([])
reaction_index[case].append(l)
sim.rtol = 1e-6
sim.atol = 1e-15
sim.rtol_sensitivity = 1e-6
sim.atol_sensitivity = 1e-6
sen = np.zeros(gas.n_reactions)
T = np.arange(0,5,1e-3)
states = ct.SolutionArray(gas, extra = ['time_in_ms','t'])
for m in T:
sim.advance(m)
for n in range(gas.n_reactions):
states.append(r.thermo.state, time_in_ms = m*1e3, t = m)
s = sim.sensitivity('temperature',n)
if np.abs(s) >= np.abs(sen[n]):
sen[n] = s #Storing maximum values
for o in range(gas.n_reactions):
for p in range(o+1,gas.n_reactions):
if np.abs(sen[o]) < np.abs(sen[p]):
temporary_1 = sen[o]
sen[o] = sen[p]
sen[p] = temporary_1
temporary_2 = reaction_index[case][o]
reaction_index[case][o] = reaction_index[case][p]
reaction_index[case][p] = temporary_2
case = case + 1
sensitive_reactions = []
reaction_counter = 1
length = len(Temp)*len(Press)*len(Phi)
#print(len(reaction_index))
counter = 0
#sensitive_reactions.append(reaction_index[0][0])
for reac_num in range(gas.n_reactions):
for case_num in range(length):
if reaction_index[case_num][counter] not in sensitive_reactions:
sensitive_reactions.append(reaction_index[case_num][counter])
counter += 1
#print(sensitive_reactions)
case_2 = 0
spec = ct.Species.listFromFile('gri30.xml')
R_min = 40
#gas2 = ct.Solution('gri30.xml')
T_max_ref = []
Ignition_ref = []
gas_ref = ct.Solution('gri30.xml')
figure_index = 1
tol = 0.05
Temp_error = 9e9
Ignition_error = 9e9
for i in range(len(Temp)):
for j in range(len(Press)):
for k in range(len(Phi)):
print('Reduction has started for case %d' %case_2)
Reactions = []
T_max_plot = []
Ignition_plot = []
Reactions_plot = []
Temp_temporary_2 = []
Ign_temporary_2 = []
minimum_reactions = []
gas_ref.TPX = Temp[i], Press[j]*101325, {'CH4':1, 'O2':2/Phi[k], 'N2':7.52/Phi[k]}
r_ref = ct.IdealGasReactor(gas_ref)
sim_ref = ct.ReactorNet ([r_ref])
states_ref = ct.SolutionArray(gas_ref)
Temp_prev = Temp[i]
Temp_ignition = Temp[i]
Temp_temporary = []
time = 0
for t in range(5000):
time = time + 1e-3
sim_ref.advance(time)
#states_ref.append(r_ref.thermo.state)
#temp_new = states_ref.T
if (gas_ref.T > Temp_prev):
Temp_prev = gas_ref.T
if ((gas_ref.T - Temp_ignition) >= 400):
Ignition_delay_ref = time
#lee = states_ref.T
#Temp_prev = lee[t]
Temp_ignition = gas_ref.T
#Temp_temporary.append(np.max(Temp_prev))
Ignition_ref.append(Ignition_delay_ref)
#Temp_max_ref = np.max(Temp_temporary)
T_max_ref.append(Temp_prev)
for reac in range(gas_ref.n_reactions):
#try:
Reactions.append(gas_ref.reaction(sensitive_reactions[reac]))
gas_reduced = ct.Solution(thermo='IdealGas',kinetics='GasKinetics',species=spec,reactions=Reactions)
gas_reduced.TPX = Temp[i],Press[j]*101325,{'CH4':1, 'O2':2/Phi[k], 'N2':7.52/Phi[k]}
r_reduced = ct.IdealGasReactor(gas_reduced)
sim_reduced = ct.ReactorNet([r_reduced])
states_reduced = ct.SolutionArray(gas_reduced)
Temp_prev_2 = Temp[i]
Temp_ignition_red = Temp[i]
time2 = 0
for ti in range(5000):
time2 = time2 + 1e-3
sim_reduced.advance(time2)
#states_reduced.append(r_reduced.thermo.state)
#temp_new2 = states_reduced.T
if (gas_reduced.T >= Temp_prev_2):
Temp_prev_2 = gas_reduced.T
if ((gas_reduced.T - Temp_ignition_red)>=400):
Ignition_delay_reduced = time2
Temp_ignition_red = gas_reduced.T
#mad = states_reduced.T
#Temp_prev_2 = mad[ti]
#Temp_temp.append(np.max(mad))
#Ign_temporary_2.append(Ignition_delay_reduced)
#Temp_temporary_2.append(np.max(Temp_temp))
try:
if (reac>=R_min):
Ignition_plot.append(Ignition_delay_reduced)
T_max_plot.append(Temp_prev_2)
Reactions_plot.append(reac)
Ignition_error = np.abs(((Ignition_delay_reduced - Ignition_delay_ref)/Ignition_delay_ref)*100)
Temp_error = np.abs(((Temp_prev_2 - Temp_prev)/Temp_prev)*100)
if(Temp_error<=tol and Ignition_error<=tol):
minimum_reactions.append(reac)
except NameError:
print('Reactions are not enough')
print('minimum_reactions for case %d = %d' %(case_2,(np.min(minimum_reactions))))
plt.figure(figure_index)
plt.subplot(2,1,1)
plt.plot(Reactions_plot,T_max_plot,color='red',label='Reference maximum Temperature - %d\nTemperature - %d\nPressure - %d\nPhi - %d' %(Temp_prev,Temp[i],Press[j],Phi[k]))
#plt.plot(Reactions_plot,T_max_ref,'--',color='yellow',label='Reference maximum Temperature - %d' %(Temp_prev))
plt.xlabel('Number of reactions')
plt.ylabel('Maximum temperature (K)')
plt.legend()
#x_pos = 0.8
#y_pos = 0.6
#plt.text(x_pos,y_pos,"Reference maximum Temperature - %d" %(Temp_prev))
#plt.title('Temperature - %d K, Pressure - %d bar, Phi - %d' %(Temp[i],Press[j],Phi[k]))
plt.subplot(2,1,2)
plt.plot(Reactions_plot,Ignition_plot,color='green',label='Reference Ignition delay - %f\nTemperature - %d\nPressure - %d\nPhi - %d' %(Ignition_delay_ref,Temp[i],Press[j],Phi[k]))
#plt.plot(Reactions_plot,Ignition_ref,'--',color='brown',label='Reference Ignition delay - %d' %(Ignition_delay_ref))
#x_pos2 = 0.8
#y_pos2 = 0.6
#plt.text(x_pos2,y_pos2,"Reference Ignition delay - %d" %(Ignition_delay_ref))
plt.xlabel('Number of reactions')
plt.ylabel('Ignition delay (sec)')
plt.legend()
plt.tight_layout()
plt.show()
case_2 = case_2+1
figure_index = figure_index + 1Result:
When the above code is executed, the result is,
0
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
Reduction has started for case 0
Reactions are not enough
Reactions are not enough
Reactions are not enough
Reactions are not enough
Reactions are not enough
Reactions are not enough
Reactions are not enough
Reactions are not enough
Reactions are not enough
Reactions are not enough
Reactions are not enough
Reactions are not enough
Reactions are not enough
Reactions are not enough
Reactions are not enough
Reactions are not enough
Reactions are not enough
Reactions are not enough
Reactions are not enough
Reactions are not enough
Reactions are not enough
Reactions are not enough
Reactions are not enough
Reactions are not enough
minimum_reactions for case 0 = 37
Reduction has started for case 1
minimum_reactions for case 1 = 155
Reduction has started for case 2
minimum_reactions for case 2 = 123
Reduction has started for case 3
minimum_reactions for case 3 = 37
Reduction has started for case 4
minimum_reactions for case 4 = 70
Reduction has started for case 5
minimum_reactions for case 5 = 108
Reduction has started for case 6
minimum_reactions for case 6 = 41
Reduction has started for case 7
minimum_reactions for case 7 = 123
Reduction has started for case 8
minimum_reactions for case 8 = 108
Reduction has started for case 9
minimum_reactions for case 9 = 35
Reduction has started for case 10
minimum_reactions for case 10 = 70
Reduction has started for case 11
minimum_reactions for case 11 = 77
Reduction has started for case 12
minimum_reactions for case 12 = 50
Reduction has started for case 13
minimum_reactions for case 13 = 77
Reduction has started for case 14
minimum_reactions for case 14 = 77
Reduction has started for case 15
minimum_reactions for case 15 = 37
Reduction has started for case 16
minimum_reactions for case 16 = 123
Reduction has started for case 17
minimum_reactions for case 17 = 75
Reduction has started for case 18
minimum_reactions for case 18 = 35
Reduction has started for case 19
minimum_reactions for case 19 = 77
Reduction has started for case 20
minimum_reactions for case 20 = 123
Reduction has started for case 21
minimum_reactions for case 21 = 38
Reduction has started for case 22
minimum_reactions for case 22 = 73
Reduction has started for case 23
minimum_reactions for case 23 = 83
Reduction has started for case 24
minimum_reactions for case 24 = 35
Reduction has started for case 25
minimum_reactions for case 25 = 73
Reduction has started for case 26
minimum_reactions for case 26 = 77
[Finished in 3256.0s]
The resultant plots are,
For case - 1:
For case - 2:
For case - 3:
For case - 4:
For case - 5:
For case - 6:
For case - 7:
For case - 8:
For case - 9:
For case - 10:
For case - 11:
For case - 12:
For case - 13:
For case - 14:
For case 15:
For case 16:
For case 17:
For case 18:
For case 19:
For case 20:
For case 21:
For case 22:
For case 23:
For case 24:
For case 25:
For case 26:
For case 27:
-
When the number of reactions are not enough, for the given time (5 seconds in this case), the combustion will not happen. In the logic used in this code, the variable 'Ignition_delay_reduced' will be created only when the temperature difference between the current state and the immediate last case is greater than or equal to 400 K (which implies combustion). So if this condition is not met within the given time (5000 milli seconds), the variable will not be created and after the loop ends, when it is tried to access, a 'NameError' shows up, stating an undeclared variable is being used.
To get past this, 'try' and 'except' method is used in which it is tried to access the 'Ignition_delay_reduced' variable and when the NameError shows up for the reason mentioned above, it is excepted by using except NameError. Now, when the reactions are not enough, the code will no crash at the error, rather a new reaction is added to the reaction list and the process is repeated until all 325 reactions are used.
Minimum number of reactions required:
For an accepted deviation of 5 percent, the data of minimum number of reactions for all 27 cases is what we have. When the minimum number of reactions is the maximum of minimum number of reactions for each case, a deviation of less than or equal to five percent for all cases from the corresponding maximum temperature and ignition delay time value when all 325 reactions are considered, will be there. In that way, the maximum of minimum number of reactions required for each case is 155.
It is also evident from the above plots that, after 155 reactions, the system has reached steady state in some cases and though it has not reached steady state in some cases, the deviation is very less (less than or equal to five percent).
So 155 top most sensitive reactions are enough to compute maximum temperature and ignition delay time of the system with atmost 5 percent devaition from the corresponding original values.
Conclusions:
- The gri30 reaction mechanism has been reduced and maximum temperature and ignition delay time of the reduced mechanism are found out, for all 27 cases. A minimum of 155 reactions are required to compute maximum temperature and ignition delay time, for a maximum deviation of five percent from the corresponding values, when all the reactions are considered