Preparing reduced chemical reaction mechanism
Aim:-> To sort reactions present in gri.30 mechanism as per sensitivity then create reduced mechanism using the sorted reactions and finally compare ignition delay and maximum temperature (Adiabatic Flame Temperature) of reduced mechanism with original gri.30 mechanism.
Theory:-
Reaction mechanism:
A balanced equation for a chemical reaction indicates what is reacting and what is produced, but it reveals nothing about how the reaction actually takes place. The reaction mechanism (or reaction path) is the process, or pathway, by which a reaction occurs. So in chemistry, it is the sequence of individual steps, or elementary reactions, by which reactants are converted into products during the course of a reaction is called the reaction mechanism. Usually, a chemical reaction occurs in steps, although it may not always be visible to an observer. While studying reaction mechanism, we call each step in a reaction mechanism an elementary reaction. Elementary reactions occur exactly as they are written and cannot be broken down into further simpler steps. And elementary reactions add up to the overall reaction. The elementary reaction possess reaction intermediates. So, reaction intermediates are chemical species, often unstable and short-lived, which are not reactants or products of the overall chemical reaction, but are temporary products and / or reactants in the mechanism's reaction steps. Most of the times reaction intermediates are free radicals or ions.
Need of reduced reaction mechanism:
As we know that the reaction mechanism describes the sequence of elementary reactions that must occur to go from reactants to products. And the simple chemical reaction mechanism contains many hydrocarbons, means it involves thousands of species and reactions. Since the computational time required to solve any reacting mechanism is very high i.e. in order of N3 of species actually engaged in process. So, it is expensive to solve such reacting mechanism which contains high number of species and reactions. To optimize the problem, scientists found a method called ‘reduction of reaction mechanism’, which reduces computational time and cost. In this reduction mechanism technique, we ignore some reactions that don't have shown any proper effect on the mechanism i.e. less sensitive reaction.
Coding procedure:-
- Created a class with some functions which can use to read text file and then extract data from text file.
- For each state values combination, a sensitivity analysis was performed for each reaction and then reactions were ordered in the descending order of sensitivities.
- The most sensitive reactions from all state values combination were appended one by one into a new reaction array and by eliminating any repeated reactions.
- The reference values for maximum temperature and ignition delay were computed for the original unreduced gri30 mechanism for each state value combination.
- Then consider reduced mechanism now.
- The maximum temperature and ignition delay values were computed for the new reduced mechanism.
- Hence after getting values for reduced mechanism and unreduced mechanism compare them using graph.
- For each of the state values combination, a graph was plotted for ignition delay and maximum temperature along with the reference values against the reduced mechanism size starting from the most sensitive 40 reactions to all the 325 reactions arranged in descending order of sensitivities.
Python program to create .txt file:-
'''
Program to write in text file
'''
text = ["Temperature 950 1050 1150n"
"Pressure 300000 400000 500000n"
"phi 0.5 1.0 1.5"]
# Writing to a file
a = open('reduce_react_mech.txt', 'w')
a.writelines(text)
a.close()
Here the above provided small code is compiled to create text file which consists of state variables like temperature, pressure and equivalence ratio with values given below.
`"Temperature"=[950, 1050, 1150] k`
`"Pressure"=[3, 4, 5]Bar`
`phi = [0.5, 1.0, 1.5]`
And this values of state variables are parsed using read function and used as initial condition for program.
Text file output:-
Temperature 950 1050 1150
Pressure 300000 400000 500000
phi 0.5 1.0 1.5
Python source program:-
'''
Program to create reduced reaction mechanism
'''
# Importing libraries
import numpy as np
import cantera as ct
import matplotlib.pyplot as plt
import os
# ~Creating user-defined class 'File_Reader' read data from file
class File_Reader:
def __init__(self, filename):
self.filename = filename
self.input_dict = {}
# Function to validate file is present
def validate_file(self):
return os.path.isfile(self.filename)
# Function to read file
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('No file under such name')
# Executing File_Reader class
a = File_Reader('reduce_react_mech.txt')
# Read file
a.read()
# Extracting values from input dictionary present in text file
## For temperature values
T_val = (a.input_dict['Temperature']) # in kelvin
## For pressure value
P_val = (a.input_dict['Pressure']) # in Pascal
## For equivalence ratio value
phi_val = (a.input_dict['phi'])
#----------------------------------------------------------------------------------------------------------------#
# ~Defining initial boundary conditions
# Creating gas object for gri.30 mechanism file
gas = ct.Solution('gri30.cti')
# Consider fuel as methane
fuel = 'CH4:1'
# Consider oxidizer
oxidizer = {'O2':0.21, 'N2':0.79}
# Case number at initial
case = 0
# Number of reactions to be consider from gri 30 mechanism
n = gas.n_reactions
min_reac = 40
# Time data
t_init = 0 # Initial time in second
t_end = 5 # Total simulation time in second
dt = 1e-3 # Time step in second
# Time array
time = np.arange(t_init,t_end,dt)
# Array for reduuced reaction index
R_index = []
const1 = True
#----------------------------------------------------------------------------------------------------------------#
# ~To calculate sensitivity
# Loop to vary temperature
for T in (T_val):
# Loop to vary pressure
for P in (P_val):
# Loop to vary equivalence ratio
for ɸ in (phi_val):
#print((T, P, ɸ))
# Reaction array
reac_n = []
# Reaction index array
index = []
# Maximum sensitivity array
s_max =[0]*n
# Introducing initial conditions
gas.TP = int(T), int(P)
gas.set_equivalence_ratio(float(ɸ), fuel, oxidizer)
# Reactor
r = ct.IdealGasReactor(gas)
# Preparing reactor network
sim = ct.ReactorNet([r])
# Loop for appending reaction
for i in range(n):
# Adding sensitivity reactions
r.add_sensitivity_reaction(i)
reac_n.append(gas.reaction_equations([i]))
index.append(i)
# Providing tolerences for numerical calculations
sim.rtol = 1e-6
sim.atol = 1e-15
sim.rtol_sensitivity = 1e-6
sim.atol_sensitivity = 1e-6
# Creating solution array for storing each state data
states = ct.SolutionArray(gas, extra = ['t_ms'])
# Time loop
for t in time:
# Advancing time
sim.advance(t)
# Loop for various reactions sensitivity
for i in range(n):
s = sim.sensitivity(2,i)
if np.abs(s_max[i]) <= np.abs(s):
s_max[i] = s
states.append(r.thermo.state, t_ms = t*1e3)
# Taking absolutevalue for sorting
s_abs = np.abs(s_max)
# Sorting reaction in descending order as per sensitivity
s_max = [p for q, p in sorted(zip(s_abs,s_max), reverse = True)]
reac_n = [p for q, p in sorted(zip(s_abs,reac_n), reverse = True)]
index = [p for q, p in sorted(zip(s_abs,index), reverse = True)]
# Loop for printing all reaction with their sensitivities and reaction index
for i in range(n):
print(str(index[i])+'.For reaction '+str(reac_n[i])+' sensitivity is '+str(s_max[i]))
# appending the sorted order of reaction indices to R_index
for i in range(0, n):
if const1 is True:
R_index.append([])
R_index[i].append(index[i])
# resetting values
const1 = False
# Upgrading case number and printing it
case += 1
print(case)
#----------------------------------------------------------------------------------------------------------------#
# ~Logic to create reduced mechanism
R_set = set()
R_order = []
# Loop for varying reactions
for i in range (0, n):
# Loop to vary state values
for j in range (0, case):
# Loop to neglect repeated reaction index
if R_index[i][j] not in R_set:
R_order.append(R_index[i][j])
R_set.add(R_index[i][j])
# Printing reaction order with number of total cases
print(R_order, case)
#----------------------------------------------------------------------------------------------------------------#
# ~To calculate AFT and ignition delay and comparing with reduced mechanism
# Extracting species from gri.30 mechanism file
spec = ct.Species.listFromFile('gri30.cti')
# Setting tolerence values
Tol_AFT = 5 # For AFT in %
Tol_IGD = 5 # For ignition delay in %
# Case number
case_n = 0
# Loop to vary temperature
for T in (T_val):
# Loop to vary pressure
for P in (P_val):
# Loop to vary equivalence ratio
for ɸ in (phi_val):
# Introducing initial conditions
gas.TP = int(T), int(P)
gas.set_equivalence_ratio(float(ɸ), fuel, oxidizer)
# Reactor
r = ct.IdealGasReactor(gas)
# Preparing reactor network
sim = ct.ReactorNet([r])
# Array to store reduced reaaction mechanism
Red_reac = []
# Array to store index of reduced reaction mechanism
Mech_size = []
# Array to store ignition delay of reduced reaction mechanism
Ig_delay_rm = []
# Array to store AFT of reduced reaction mechanism
Max_temp_rm = []
loop_counter_1 = 0
# Initialize values
Auto_ig_time = 0
AFT = 0
# Time loop for advancing
for t in time:
# Advancing time
sim.advance(t)
# To calculate Auto ignition delay at (T+400)k
if ((gas.T >= (int(T) + 400)) and loop_counter_1 == 0):
Auto_ig_time = t
loop_counter_1 = 1
# To calculate AFT
if (AFT < gas.T):
AFT = gas.T
# For reaction index
i = 0
Const2 = True
# Initializing error value
Error_aft = 0 # For AFT
Error_igd = 0 # For ignition delay
# Loop to add reactions
for i in range (0, n):
# Array consisting reduced reaction mechanism
Red_reac.append(gas.reaction(R_order[i]))
# Creating gas_rm object of solution class for reduced reaction mechanism
gas_rm = ct.Solution(thermo='IdealGas',kinetics='GasKinetics',species=spec,reactions=Red_reac)
# Introducing initial values
gas_rm.TP = int(T), int(P)
gas_rm.set_equivalence_ratio(float(ɸ), fuel, oxidizer)
# Reactor
r_rm = ct.IdealGasReactor(gas_rm)
# Preparing reactor network for reduced reaction mechanism
sim_rm = ct.ReactorNet([r_rm])
loop_counter_2 = 0
# Initializing values
AFT_rm = 0
Auto_ig_time_rm = 0
# Loop for time advancement
for t in time:
# Advancing time
sim_rm.advance(t)
# To compute auto ignition delay for reduced reaction mechanism
if ((gas_rm.T >= (int(T) + 400)) and loop_counter_2 == 0):
Auto_ig_time_rm = t
loop_counter_2 = 1
# To compute AFT for reduced reaction mechanism
if (AFT_rm < gas_rm.T):
AFT_rm = gas_rm.T
# Incrementing reaction index
i += 1
# Using try and except functions as per requirements
try:
if (i >= min_reac):
Ig_delay_rm.append(Auto_ig_time_rm)
Max_temp_rm.append(AFT_rm)
Mech_size.append(i)
# Error computing
Error_aft = np.abs(((AFT - AFT_rm)/AFT)*100)
Error_igd = np.abs(((Auto_ig_time - Auto_ig_time_rm)/Auto_ig_time)*100)
# Comparing errors
if (Error_aft < Tol_AFT and Error_igd < Tol_IGD and Const2 == True):
text_to_print=('T = '+str(int(T))+'k, P = '+str(int(P)/100000)+'Bar, ɸ = '+str(float(ɸ))+
', AFT = '+str("{:.4f}".format(AFT))+'k, Ig. delay = '+str("{:.4f}".format(Auto_ig_time))+
'S, AFT_rm = '+str("{:.4f}".format(AFT_rm))+'k, Ig. delay_rm = '+
str("{:.4f}".format(Auto_ig_time_rm))+'S,n Error in AFT = '+str("{:.4f}".format(Error_aft))+
'%, Error in delay = '+str("{:.4f}".format(Error_igd))+'%, Reduced mechanism size = '+str(i))
print(text_to_print)
Const2 = False
except NameError:
pass
# Case increment
case_n += 1
# Plotting
plt.figure(case_n)
## Sublot for raction v/s ignition delay
plt.subplot(2,1,1)
plt.plot(Mech_size, Ig_delay_rm, color='orange')
plt.axhline(y=Auto_ig_time,linestyle='--', color='blue')
plt.legend(['Computed ignition delay','Reference delay='+str("{:.4f}".format(Auto_ig_time))+' S'], loc = 0)
plt.ylabel('Ignition delay $[S]$', fontweight='bold')
plt.xlabel('Number of reactions', fontweight='bold')
plt.grid('both',linestyle='--',linewidth = 0.8)
## Subplot for reaction v/s AFT
plt.subplot(2,1,2)
plt.plot(Mech_size, Max_temp_rm, color='green')
plt.axhline(y=AFT, linestyle='--', color='red')
plt.legend(['Computed AFT','Reference temp='+str("{:.4f}".format(AFT))+'k'], loc = 0)
plt.ylabel('AFT $[k]$', fontweight='bold')
plt.xlabel('Number of reactions', fontweight='bold')
plt.grid('both',linestyle='--',linewidth = 0.8)
plt.suptitle('Case:'+str(case_n)+'nT = '+str(int(T))+'k, P = '+str(int(P)/100000)+' Bar & ɸ = '+str(float(ɸ)))
plt.tight_layout()
## Saving plots in .png format
plot_name = 'Case_no'+str(case_n)+'.png'
plt.savefig(plot_name)
print('Case '+str(case_n)+' is completed')
plt.show()
The reduced reaction mechanism as per gri.30 mechanism index:-
[118, 157, 31, 155, 154, 160, 52, 115, 169, 97, 35, 120, 37, 100, 56, 84, 117, 114, 162, 156, 33, 112, 32, 10, 167, 158,
166, 164, 86, 14, 57, 51, 45, 44, 174, 119, 96, 73, 94, 83, 34, 26, 77, 3, 161, 87, 177, 98, 103, 104, 111, 47, 144, 88,
163, 141, 55, 286, 67, 283, 289, 311, 68, 9, 46, 43, 165, 293, 143, 284, 178, 306, 172, 290, 24, 303, 17, 74, 305, 159,
239, 4, 134, 60, 2, 287, 181, 149, 179, 138, 198, 95, 294, 30, 207, 126, 124, 18, 285, 197, 65, 75, 70, 64, 146, 82,
80, 125, 11, 205, 53, 168, 238, 66, 184, 147, 20, 102, 85, 113, 121, 91, 282, 99, 196, 323, 191, 54, 313, 89, 42, 16,
25, 79, 48, 193, 183, 139, 182, 116, 62, 211, 137, 153, 90, 189, 6, 150, 29, 292, 131, 203, 314, 279, 76, 22, 324, 204,
61, 92, 317, 145, 213, 12, 304, 194, 170, 13, 192, 135, 27, 315, 130, 214, 309, 152, 310, 322, 129, 72, 312, 308, 175, 316,
28, 136, 180, 298, 108, 93, 299, 200, 173, 63, 40, 307, 297, 190, 300, 38, 5, 50, 301, 71, 128, 148, 23, 202, 59, 321,
106, 295, 110, 0, 123, 1, 320, 151, 319, 49, 318, 210, 78, 227, 105, 288, 291, 212, 240, 101, 209, 19, 58, 81, 132, 186,
171, 140, 228, 15, 206, 39, 215, 107, 201, 296, 21, 127, 302, 109, 7, 208, 188, 122, 268, 221, 230, 258, 133, 185, 237, 8,
241, 275, 231, 195, 199, 264, 223, 256, 263, 244, 259, 260, 187, 273, 219, 249, 242, 248, 276, 255, 245, 274, 266, 225, 277,
233, 176, 247, 222, 254, 217, 236, 226, 243, 265, 232, 218, 252, 69, 253, 216, 272, 41, 257, 251, 250, 267, 246, 281, 271,
235, 234, 229, 220, 270, 261, 269, 224, 278, 262, 280, 142, 36] 27Output data:-
 |
|
 |
|
 |
|
 |
|
 |
|
 |
|
 |
|
 |
|
 |
|
 |
|
 |
|
 |
|
 |
|
 |
|
 |
|
 |
|
 |
|
 |
|
 |
|
 |
|
 |
|
 |
|
 |
|
 |
|
 |
|
 |
|
 |
|
In above data its seen that each plot behaves differently for different state variables. Every image provided above include 2 sub plots, 1st one shows the plots for Ignition Delay v/s number of reactions in reduced mechanism and 2nd subplot shows Maximum Temperature v/s number of reactions in reduced mechanism. For both of the subplots it can be observed that values are getting fluctuated and then becomes stable.
Earlier we have created reduced reaction mechanism which consists intermediate reactions and sorted in descending order on basis of sensitivity, thus the most sensitive reactions are at initial which gives most fluctuated values.
While plotting we did not consider first 40 reaction in account and plotted graphs after 40th reaction to observe at which point the ignition delay and maximum temperature matches with reference values. Which indirectly tells us about number of reactions needed and number of reactions can ignore to get same results as unreduced mechanism.
Conclusion:-
By observing above output data and plots it can be concluded that on an average last 100 reactions from reduced reaction mechanism (Those reactions which possess minimum sensitivity value) can be ignored during simulation, in order to optimize computational cost and time.