Detailed reaction mechanism reduction project

A reaction mechanism is something that dictates which products are formed when two or more reactants react under different conditions. It also dictates the rate of change of concentration of species and heat of combustion.

When a reaction takes place, a huge number of intermediary elementary reactions takes place which is accurately captured by detailed reaction mechanisms, for this reason a detailed reaction mechanism can used in variety of applications like premixed, partially premixed, diffusion combustion modelling. But using a detailed reaction mechanism is expensive computationally and takes a large amount of time as the reaction mechanism can be quite stiff i.e., when one reaction takes place at a faster rate and another one takes place at a slower rate, to obtain the accurate results of the slow rate reaction the simulation should be carried at a much slower pace this can increase the time taken for the simulation to complete significantly. In order to avoid this, a reduced reaction mechanism is used. The disadvantages of using a reduced reaction mechanism is that the mechanism can only be used in specific application but the mechanism is generally quite faster than the detailed one and less stiff so computationally less taxing.

In this project, a detailed reaction mechanism (GRI 3.0) is taken and converted to a reduced mechanism which can be used to obtain ignition time and max temperature at a specific condition (temperature, pressure and equivalence ratio) with the results varying a maximum of 1% with significantly lower reactions. The process is carried out using Cantera via Python.

The following steps were taken to reduce the detailed reaction mechanism:

1. First, to get a reference Ignition delay and Max temperature analysis of unreduced mechanism for a set of state values (Temperature, Pressure and Equivalence ratio) is carried out.
2. A sensitivity analysis is carried out with respect to temperature for all the reactions in the detailed mechanism and the reactions are sorted from the most sensitive to least sensitive (absolute value of sensitivity is used for sorting)
3. These sorted reactions are then added to a new mechanism one by one starting from a fixed minimum number of reactions (40).
4. A graph is plotted showing the reduced mechanism's ignition delay and max temperature against the number of reactions in the reduced mechanism. The reference ignition delay and Temperature are mentioned in the label.
5. Now, the same steps are repeated various state values such that every possible combinations of the following state values are covered:

Temperature - [950 , 1050 , 1150] Kelvin,

Pressure - [3 , 4 , 5] bar

Equivalence ratio - [0.5 , 1 , 1.5].

So, in total 27 possible combinations. For the all the analysis, the input is given via a txt file. The program is written in such a way that the text file is read and coverted into lists which are then used in the program and the computation is carried out for 5s with time steps of 1e-2s. If the starting number of reaction was significantly lower than 40, error might have occured during the calculation of ignition time which can be bypassed using the try and exception method. The starting number of reactions is taken as 40 for the project which ran without causing any error so no such methods were adopted here.

The link for the code is given below:

https://drive.google.com/file/d/1EnS2kVC6dl_77OF6EnNhvkgbsyNdHAAY/view?usp=sharing

And the code itself is pasted below:

Temperature 950 1050 1150
Pressure 300000 400000 500000
Phi 0.5 1 1.5

import os
import cantera as ct 
import numpy as np
import matplotlib.pyplot as plt
import itertools

class Input_file_reader():
 def __init__(self, file_name):
  self.file_name = file_name
  self.input_dict = {}

 def validate_input_file(self):
  if os.path.isfile(self.file_name):
   return True
  else:
   return False

 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')

a = Input_file_reader('D:/Cantera/Mechanism reduction/input.txt')
a.read()
T_input = [int(x) for x in a.input_dict['Temperature']]
P_input = [int(x) for x in a.input_dict['Pressure']]
Phi_input = [float(x) for x in a.input_dict['Phi']]
combi = 0
for b in T_input:
 for c in P_input:
  for d in Phi_input:
   T = b
   P = c
   Phi = d 
   gas = ct.Solution('gri30.xml')
   t_autoig = 1
   steps = int((5)/(1e-2))
   sensitivity = [[0 for b in range(500)] for v in range(325)]
   sensitivity_max = []
   sensitivity_max_arranged = []
   reactionsmas = []
   reactions = []
   t_master_autoig = []
   T_master = []
   rnlist = []
   gas.TP = b,c
   gas.set_equivalence_ratio(d, 'CH4', 'O2:1.0, N2:3.76')

   solution = ct.SolutionArray(gas,extra=['time'])
   r = ct.IdealGasReactor(gas)
   s = ct.ReactorNet([r])
   t_bool = True
   t = 0

   for i in range(325):
    r.add_sensitivity_reaction(i)

   s.rtol = 1e-6
   s.atol = 1e-15

   s.rtol_sensitivity = 1e-6
   s.atol_sensitivity = 1e-6

   for n in range(steps):
    
    t = t + 1e-2
    s.advance(t)
    
    for j in range(325):
     sensitivity[j][n] = s.sensitivity(1,j)

    solution.append(r.thermo.state,time=t*1e3)

    if (solution.T[n]-T) >= 400:
     if t_bool==True:
      t_autoig=solution.time[n]
      t_bool=False

   sensitivityar = np.array(sensitivity)
   for y in sensitivity:
    sensitivity_max.append(np.min(y) if np.max(np.abs(y)) == abs(np.min(y)) else np.max(y))

   sorted_d = sorted(zip(gas.reactions(),sensitivity_max), key=lambda x: abs(x[1]), reverse=True)

   # number of reactions present rn
   rn = 40
   #condition for while loop ending
   cond = 1 

   while(cond==1):
    reactionsmas,sensitivity_max_arranged = zip(*sorted_d) 
    reactions = reactionsmas[:rn]
    #Including the minimum species
    species_names = {'CH4', 'O2','N2'}
    for reaction in reactions:
     species_names.update(reaction.reactants)
     species_names.update(reaction.products)

    # Get the species objects
    species = [gas.species(name) for name in species_names]

    # create the new reduced mechanism
    gas2 = ct.Solution(thermo='IdealGas', kinetics='GasKinetics', species=species, reactions=reactions)
    t_autoig2 = 1
    t_bool2 = True
    t2 = 0
    gas2.TP = b,c
    gas2.set_equivalence_ratio(d, 'CH4', 'O2:1.0, N2:3.76')
    solution2 = ct.SolutionArray(gas2,extra=['time'])
    r2 = ct.IdealGasReactor(gas2)
    s2 = ct.ReactorNet([r2])
    for o in range(steps):
    
     t2 = t2 + 1e-2
     s2.advance(t2)
     solution2.append(r2.thermo.state,time=t2*1e3)
     if (solution2.T[o]-T) >= 400:
      if t_bool2==True:
       t_autoig2=solution2.time[o]
       t_bool2=False

    t_master_autoig.append(t_autoig2)
    T_master.append(max(solution2.T))
    if((abs(t_autoig2-t_autoig)/t_autoig)*100 < 0 xss=removed xss=removed label='Ignition time plot for reduced mechanism. The Ignition time for detailed mechanism is %i ms' plt.axhline(y=t_autoig, xmin=0, xmax=1, linestyle ="--" label='Ignition time with detailed mechanism (%i ms) ' xss=removed plt.show(block=False) dpi=300,bbox_inches label='Max Temperature plot for reduced mechanism. The Max Temperature for detailed mechanism %i K' plt.axhline(y=max(solution.T), xmin=0, xmax=1, linestyle ="--" label='Max Temperature with detailed mechanism (%i K) ' xss=removed plt.show(block=False) dpi=300,bbox_inches xss=removed>

Results:

For 950K initial temperature, 3bar initial pressure and 0.5 equivalence ratio:

For 950K initial temperature, 3bar initial pressure and 1 equivalence ratio:

For 950K initial temperature, 3bar initial pressure and 1.5 equivalence ratio:

For 950K initial temperature, 4bar initial pressure and 0.5 equivalence ratio:

For 950K initial temperature, 4bar initial pressure and 1 equivalence ratio:

For 950K initial temperature, 4bar initial pressure and 1.5 equivalence ratio:

For 950K initial temperature, 5bar initial pressure and 0.5 equivalence ratio:

For 950K initial temperature, 5bar initial pressure and 1 equivalence ratio:

For 950K initial temperature, 5bar initial pressure and 1.5 equivalence ratio:

For 1050K initial temperature, 3bar initial pressure and 0.5 equivalence ratio:

For 1050K initial temperature, 3bar initial pressure and 1 equivalence ratio:

For 1050K initial temperature, 3bar initial pressure and 1.5 equivalence ratio:

For 1050K initial temperature, 4bar initial pressure and 0.5 equivalence ratio:

For 1050K initial temperature, 4bar initial pressure and 1 equivalence ratio:

For 1050K initial temperature, 4bar initial pressure and 1.5 equivalence ratio:

For 1050K initial temperature, 5bar initial pressure and 0.5 equivalence ratio:

For 1050K initial temperature, 5bar initial pressure and 1.0 equivalence ratio:

For 1050K initial temperature, 5bar initial pressure and 1.5 equivalence ratio:

For 1150K initial temperature, 3bar initial pressure and 0.5 equivalence ratio:

For 1150K initial temperature, 3bar initial pressure and 1.0 equivalence ratio:

For 1150K initial temperature, 3bar initial pressure and 1.5 equivalence ratio:

For 1150K initial temperature, 4bar initial pressure and 0.5 equivalence ratio:

For 1150K initial temperature, 4bar initial pressure and 1.0 equivalence ratio:

For 1150K initial temperature, 4bar initial pressure and 1.5 equivalence ratio:

For 1150K initial temperature, 5bar initial pressure and 0.5 equivalence ratio:

For 1150K initial temperature, 5bar initial pressure and 1.0 equivalence ratio:

For 1150K initial temperature, 5bar initial pressure and 1.5 equivalence ratio:

The above results are tabulated below for better comparisons: (Results with 1% variation from the detailed mechanism)

No.	Initial Temperature (K)	Initial Pressure (Pa)	Phi	Ignition time (ms)	Max temperature (K)	No. of reactions in mechanism
1	950	300000	0.5	600	2283	131
2	950	300000	1	770	2827	256
3	950	300000	1.5	940	2724	189
4	950	400000	0.5	450	2283	149
5	950	400000	1	560	2845	137
6	950	400000	1.5	670	2732	155
7	950	500000	0.5	360	2285	115
8	950	500000	1	440	2859	138
9	950	500000	1.5	520	2738	122
10	1050	300000	0.5	110	2348	200
11	1050	300000	1	150	2859	188
12	1050	300000	1.5	180	2781	307
13	1050	400000	0.5	80	2354	244
14	1050	400000	1	110	2878	154
15	1050	400000	1.5	130	2792	104
16	1050	500000	0.5	60	2362	145
17	1050	500000	1	80	2894	116
18	1050	500000	1.5	100	2800	120
19	1150	300000	0.5	30	2413	200
20	1150	300000	1	40	2890	180
21	1150	300000	1.5	40	2836	209
22	1150	400000	0.5	20	2424	178
23	1150	400000	1	30	2911	130
24	1150	400000	1.5	30	2849	104
25	1150	500000	0.5	20	2418	281
26	1150	500000	1	20	2927	155
27	1150	500000	1.5	30	2859	44

 As you can see from the above results, using a reduced mechanism we can obtain results with 1% tolerance with less number of reactions than the detailed reaction mechanism. In some of the above ignition time graphs, flat lines can be observed, this is beacuse with significantly less number of reactions the reaction has reached the 1% tolerance for the ignition delay whereas one can observe the max temperature plot not showing the same flat lines, this is because the temperature has not achieved the 1% tolerance and therefore dragging the simulation till the required tolerance for the max temperature is also achieved. Further, to get a more detailed graph a time step of 1e-3s can be used. This can significantly reduce the computational cost and reduce the time taken for the simulation to complete by reducing the stiffness of the solver(as lesser number of reactions are present). This project successfully demonstrates the advantages of using a reduced reaction mechanisms.