Reaction Mechanism Reduction
Aim: - To write a Reaction Reduction Mechanism Code for GRI3.0 Mechanism File.
Objective: -
- To perform Sensitivity Analysis with respect to Temperature For all reactions IN the Mechanism File and Arrange them in decreasing Order.
- To run Ignition Delay reference and Max Temperature Analysis for Unreduced Mechanism for 5 sec for a set of State values of temperature, Pressure and Equivalence Ratio.
State Values: The input file Consists of - Temperature 950 1050 1150
Pressure 3 4 5
phi 0.5 1 1.5
- To reduce the Mechanism file into new reduced Mechanism for tolerance within 5% Ignition delay and Max temperature with respect to reference from above objective.
- To plot the values of Ignition Delay and Max temperature Vs no. of reactions.
Theory: -
Reaction Mechanism –
It is the step-by-step sequence of elementary reactions by which overall chemical change occurs. It is a theoretical compilation that explains in detail what happens at every stage of a chemical reaction. The reactions themselves may involve the interactions of atoms, molecules, ions, electrons, and free radicals, and they may take place in gases, liquids, or solids—or at interfaces between any of these. It consists of Elements, Species, Reaction equation and its transport Properties Coefficient and Values.
GRI 3.0 is the mechanism file is an optimized mechanism of Combustion of Natural gas having Air fuel mixture and consists of total of 325 Reaction.
Reduction of Reaction Mechanism –
It is a process to exclude the redundant reaction which have very low to no sensitivity with the parameter of choice and get the solution faster with less calculation capacity on the expense of accuracy. Accuracy is set on certain tolerance range and mechanism is reduced depending on the specific range of operating conditions we want for the fuel.
Steps/Approach: -
- Sensitivity analysis – An instance of object gas is created with is setup with initial condition form set of state Values. Sensitivity is calculated by creating a Reactor Network and calculating the sensitivity for each reaction for every set of state value. The max sensitivity is stored in an array which transform into 2-D Matrix since of 325*27(number of sets of state value), the Sensitivity is sorted and zipped with the reaction equation.
- the 2-D matric is flattened out to create an array of non-repeating reactions in decreasing order of sensitivity.
- Reference Values are calculated by conditioning the temperature and storing the max temperature and the time of this first max temperature instance.
T_ign = T + 400 K
- A loop is set inside the reference values loop to add reaction one by one and repeat the same process again and check whether the value meet the tolerance criteria and then it is printed and plotted.
- Since when the number of reactions is too low python show NameError which stops the code. This could be avoided by using try and except.
Code: -
"""
program for reduction of GRI 3.0 mechanism
on the basis of Temperature Sensitivity
"""
import os
from cantera import *
from numpy import *
from matplotlib.pyplot import *
# file reader class definition
class file_reader:
def __init__(self,file_name):
self.file_name = file_name
self.input_dict = {}
def validate_file(self):
return os.path.isfile(self.file_name)
def read(self):
if self.validate_file():
for line in open(self.file_name,'r'):
# print(line)
self.input_dict[line.split()[0]] = line.split()[1:4]
else:
print('File Does Not Exist')
a = file_reader('C:/Users/Asus/Documents/py_files/input_file.txt')
a.read()
# print(a.input_dict)
T_range = [float(x) for x in a.input_dict['Temperature']]
P_range = [float(x) for x in a.input_dict['Pressure']]
phi_range = [float(x) for x in a.input_dict['phi']]
total_t = 5
dt = 1e-3
#steps = int(total_t/dt)
gas = Solution('gri30.xml')
n = gas.n_reactions
rxn_no = linspace(1,n,n).astype(int)
ph = 0
indx = []
s_sor = []
const1 = True
#print (T_range,P_range,phi_range)
# For creating decreasing sensitivity order of reaction indicies from gri mech file
for temp in T_range:
for p in P_range:
for phi in phi_range:
print(temp,p,phi)
gas.TP = temp , p*one_atm
gas.set_equivalence_ratio(phi , 'CH4' , 'O2:1 , N2:3.76')
r = IdealGasReactor(gas)
sim = ReactorNet([r])
for i in range(n): r.add_sensitivity_reaction(i)
# set the tol. for the solution and the sensitivity coeff.
sim.rtol = 1.0e-6
sim.atol = 1.0e-15
sim.rtol_sensitivity = 1.0e-6
sim.atol_sensitivity = 1.0e-6
s = zeros(n)
s_max =[0]*n
# calculating sensitivity
for t in arange(0, total_t+dt, dt):
sim.advance(t)
for j in range(n):
s = sim.sensitivity('temperature',j)
if abs(s) >= abs(s_max[j]): # Storing Max Sensitivity
s_max[j] = s
# Storing Sensitivity Order For all Phase Combinations
for k in range(n):
s_sor.append([])
s_sor[k].append(s_max[k])
s_sor[k].append(k)
# Sorting Sensitivity Order
s_max_sor= sorted(s_sor, key = lambda x: abs(x[0]), reverse = True)
# Storing Reaction Indices in Max Sensitivity Order
for k in range(n):
if const1 is True:
indx.append([])
indx[k].append(s_max_sor[k][1])
s_sor.clear()
ph += 1
const1 = False
# storing Reaction Indicies by flattening Out the 2D indices Matrix into array
in_set=set()
in_order = []
for i in range(n):
for j in range(ph):
if indx[i][j] not in in_set:
in_order.append(indx[i][j])
in_set.add(indx[i][j])
#print(in_order)
#printing the Reaction equation With its Sensitivity in Decending Order
rxn_name = [r for s_ab, r in sorted(zip(in_order,gas.reaction_equations()),key=lambda x: x[0],reverse = True)]
for k in range(len(rxn_name)): print (rxn_no[k], 't' ,rxn_name[k] ,'t :', s_max_sor[k] ,'n')
temp_autoi = []
t_autoi = []
mech_size = []
min_size = 40
rxns = []
tol_delay = 5
tol_temp = 5
# For Ignition delay And Mechanism Reduction
for temp in T_range:
for p in P_range:
for phi in phi_range:
gas.TP = temp , p*one_atm
gas.set_equivalence_ratio(phi , 'CH4' , 'O2:1 , N2:3.76')
r = IdealGasReactor(gas)
sim = ReactorNet([r])
# Calculation Ignition Temperature For Reference
T_ign = temp + 400
T_max_ref = 0
const2 = True
# Calculating Ign Delay
for t in arange(0, total_t+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
const4 = True
i = 0
print(temp,p,phi)
# Mechanism Reduction
for i in range(0, n):
rxns.append(gas.reaction(in_order[i]))
species_names = {'CH4', 'O2','N2'}
for reaction in rxns:
species_names.update(reaction.reactants)
species_names.update(reaction.products)
spec = [gas.species(name) for name in species_names]
gas1 = Solution(thermo='IdealGas', kinetics='GasKinetics', species=spec, reactions =rxns)
gas1.TP = temp,p*one_atm
gas1.set_equivalence_ratio(phi,'CH4','O2:1, N2:3.76')
r1 = IdealGasReactor(gas1)
sim1 = ReactorNet([r1])
T_max = 0
const3 = True
for t in arange(0, total_t, 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 += 1
try:
# Storing IgnDelay And Max Temperature
if (i >= min_size):
mech_size.append(i)
temp_autoi.append(T_max)
t_autoi.append(ign_delay)
# Calculating Error
temp_err = np.abs(((T_max_ref - T_max)/T_max_ref)*100)
ign_delay_err = np.abs(((ign_delay_ref - ign_delay)/ign_delay_ref)*100)
if ((tol_delay > ign_delay_err) and (tol_temp > temp_err) and const4 is True):
print('T = %d , P = %d , Phi = %f , Mechanism Size = %d' %(temp,p,phi,i))
print('T_max_reduced = %d K, Ign_delay_reduced = %f , T_max_ref = %d K, Ign_delay_ref = %f' %(T_max,ign_delay,T_max_ref,ign_delay_ref))
print('T_max_error = %f , Ign_delay_error = %f'%(temp_err,ign_delay_err))
const4 = False
except NameError:
pass
# plotting
subplot(2,1,1)
suptitle('Ignition time & Max Temperature plot at T = %d , P = %d , Phi = %.1f ' %(temp,p,phi))
plot(mech_size,t_autoi, color='blue')
axhline(y = ign_delay_ref, color = 'black', linestyle = 'dotted', label = 'Reference Ignition Delay')
ylabel('Ignition time(sec)')
xlabel('Number of reactions')
grid('both',linestyle='--',linewidth=0.8)
legend(loc = 'best')
subplot(2,1,2)
plot(mech_size,temp_autoi, color='red')
axhline(y = T_max_ref, color = 'black', linestyle = 'dotted', label = 'Reference Max Temperature')
ylabel('Maximum Temperature (K)')
xlabel('Number of reactions')
grid('both',linestyle='--',linewidth=0.8)
legend(loc = 'best')
tight_layout()
savefig(str('Ign-time & Max Temperature plot at %d K, %d atm , %.1f equivalence ratio' %(temp,p,phi))+'.png',dpi=300,bbox_inches='tight')
show()
# resetting arrays for next state values combination
rxns.clear()
temp_autoi.clear()
mech_size.clear()
t_autoi.clear()Results: -
Sensitivity Order-
Plots for States -
1.
T_max_reduced = 2283 K, Ign_delay_reduced = 0.592000 , T_max_ref = 2283 K, Ign_delay_ref = 0.591000
T_max_error = 0.023288 , Ign_delay_error = 0.169205
2.
T_max_reduced = 2848 K, Ign_delay_reduced = 0.763000 , T_max_ref = 2830 K, Ign_delay_ref = 0.759000
T_max_error = 0.620342 , Ign_delay_error = 0.527009
3.
T_max_reduced = 2765 K, Ign_delay_reduced = 0.935000 , T_max_ref = 2727 K, Ign_delay_ref = 0.925000
T_max_error = 1.385378 , Ign_delay_error = 1.081081
4.
T_max_reduced = 2286 K, Ign_delay_reduced = 0.440000 , T_max_ref = 2286 K, Ign_delay_ref = 0.440000
T_max_error = 0.003945 , Ign_delay_error = 0.000000
5.
T_max_reduced = 2867 K, Ign_delay_reduced = 0.550000 , T_max_ref = 2849 K, Ign_delay_ref = 0.547000
T_max_error = 0.629492 , Ign_delay_error = 0.548446
6.
T_max_reduced = 2777 K, Ign_delay_reduced = 0.665000 , T_max_ref = 2733 K, Ign_delay_ref = 0.659000
T_max_error = 1.621557 , Ign_delay_error = 0.910470
7.
T_max_reduced = 2288 K, Ign_delay_reduced = 0.351000 , T_max_ref = 2287 K, Ign_delay_ref = 0.352000
T_max_error = 0.025099 , Ign_delay_error = 0.284091
8.
T_max_reduced = 2882 K, Ign_delay_reduced = 0.428000 , T_max_ref = 2864 K, Ign_delay_ref = 0.426000
T_max_error = 0.618481 , Ign_delay_error = 0.469484
9.
T_max_reduced = 2786 K, Ign_delay_reduced = 0.511000 , T_max_ref = 2741 K, Ign_delay_ref = 0.507000
T_max_error = 1.626210 , Ign_delay_error = 0.788955
10.
T_max_reduced = 2357 K, Ign_delay_reduced = 0.102000 , T_max_ref = 2356 K, Ign_delay_ref = 0.102000
T_max_error = 0.046249 , Ign_delay_error = 0.000000
11.
T_max_reduced = 2880 K, Ign_delay_reduced = 0.142000 , T_max_ref = 2862 K, Ign_delay_ref = 0.140000
T_max_error = 0.639175 , Ign_delay_error = 1.428571
12.
T_max_reduced = 2813 K, Ign_delay_reduced = 0.177000 , T_max_ref = 2782 K, Ign_delay_ref = 0.173000
T_max_error = 1.136323 , Ign_delay_error = 2.312139
13.
T_max_reduced = 2360 K, Ign_delay_reduced = 0.075000 , T_max_ref = 2359 K, Ign_delay_ref = 0.075000
T_max_error = 0.050353 , Ign_delay_error = 0.000000
14.
T_max_reduced = 2901 K, Ign_delay_reduced = 0.102000 , T_max_ref = 2881 K, Ign_delay_ref = 0.101000
T_max_error = 0.692845 , Ign_delay_error = 0.990099
15.
T_max_reduced = 2827 K, Ign_delay_reduced = 0.127000 , T_max_ref = 2793 K, Ign_delay_ref = 0.125000
T_max_error = 1.234009 , Ign_delay_error = 1.600000
16.
T_max_reduced = 2363 K, Ign_delay_reduced = 0.059000 , T_max_ref = 2362 K, Ign_delay_ref = 0.059000
T_max_error = 0.039130 , Ign_delay_error = 0.000000
17.
T_max_reduced = 2916 K, Ign_delay_reduced = 0.079000 , T_max_ref = 2898 K, Ign_delay_ref = 0.078000
T_max_error = 0.619921 , Ign_delay_error = 1.282051
18.
T_max_reduced = 2837 K, Ign_delay_reduced = 0.098000 , T_max_ref = 2803 K, Ign_delay_ref = 0.096000
T_max_error = 1.233011 , Ign_delay_error = 2.083333
19.
T_max_reduced = 2428 K, Ign_delay_reduced = 0.023000 , T_max_ref = 2427 K, Ign_delay_ref = 0.023000
T_max_error = 0.041042 , Ign_delay_error = 0.000000
20.
T_max_reduced = 2912 K, Ign_delay_reduced = 0.033000 , T_max_ref = 2892 K, Ign_delay_ref = 0.032000
T_max_error = 0.677232 , Ign_delay_error = 3.125000
21.
T_max_reduced = 2861 K, Ign_delay_reduced = 0.040000 , T_max_ref = 2841 K, Ign_delay_ref = 0.039000
T_max_error = 0.689546 , Ign_delay_error = 2.564103
22.
T_max_reduced = 2433 K, Ign_delay_reduced = 0.017000 , T_max_ref = 2432 K, Ign_delay_ref = 0.017000
T_max_error = 0.015117 , Ign_delay_error = 0.000000
23.
T_max_reduced = 2934 K, Ign_delay_reduced = 0.024000 , T_max_ref = 2912 K, Ign_delay_ref = 0.024000
T_max_error = 0.742355 , Ign_delay_error = 0.000000
24.
T_max_reduced = 2877 K, Ign_delay_reduced = 0.030000 , T_max_ref = 2851 K, Ign_delay_ref = 0.029000
T_max_error = 0.881633 , Ign_delay_error = 3.448276
25.
T_max_reduced = 2436 K, Ign_delay_reduced = 0.014000 , T_max_ref = 2434 K, Ign_delay_ref = 0.014000
T_max_error = 0.096609 , Ign_delay_error = 0.000000
26.
T_max_reduced = 2951 K, Ign_delay_reduced = 0.019000 , T_max_ref = 2928 K, Ign_delay_ref = 0.019000
T_max_error = 0.764673 , Ign_delay_error = 0.000000
27.
T_max_reduced = 2888 K, Ign_delay_reduced = 0.024000 , T_max_ref = 2859 K, Ign_delay_ref = 0.023000
T_max_error = 1.006577 , Ign_delay_error = 4.347826
Conclusion: -
1. The max error is around 1.63% for max temperature which is less than 5% error tolerence, Hence every state is well within the 5% range at 40 reaction inclusion and many even before.
2. Most of the result Stablises below 70 reaction after which the error may fall below 0.01 %.
3. Decreasing the tolerence will ramp up the calculation time but accuracy will be better say 0.1%.