Reaction mechanism reduction for the GRI 3.0 mechanism with different state variable combinations using Python and Cantera

INTRODUCTION:

A reaction mechanism is a step-by-step sequence of elementary reactions by which overall chemical change occurs. It helps to study the different stages of a chemical reaction in detail. Obtaining an accurate mechanism for a given set of conditions is very challenging. Different conditions affect the reaction path from the reactants to the final products. The reaction path includes the intermediate reactions that occur between intermediate species formed. To accurately solve the chemical kinetics involved in a reaction, using an accurate chemistry reaction mechanism is essential. This reaction mechanism can contain hundreds of reactions and thousands of species being formed. The computational cost for solving a mechanism is proportional to `N^3`, where `N` is the number of species involved in the reaction mechanism. Mechanism reduction can help accelerate this computational analysis. Mechanism reduction is the process of neglecting some reactions which do not have a considerable impact on the mechanism. The impact of a reaction is compared through its sensitivity with respect to temperature. For this study, we implement the mechanism reduction for the GRI 3.0 mechanism with different combinations of the state variables as:

Temperature: `950, 1050, 1150` `K`

Pressure: `3.0, 4.0, 5.0` `b ar`

Equivalence Ratio: `0.5, 1.0, 1.5`

The steps in the process of implementation of the mechanism reduction are listed below:

• Extract the input variable data from an input file.
• Initiate gas as an instance of the Solution class from Cantera.
• Create IdealGasReactor and ReactorNet objects from the gas instance based on specified state variables.
• Simulate the 0-D reactor for the above reactor network and calculate the maximum temperature and ignition delay for each of the cases along with the reaction sensitivities to temperature.
• Sort the reactions based on the absolute values of the sensitivities for the unreduced mechanism in each case.
• Prioritize the reactions based on the sorted reactions from each case, such that the overall reaction order can be obtained based on all the cases.
• For each case, cumulatively add single reactions to the mechanism and calculate the maximum temperature and ignition delay for the reduced mechanism.
• If the reactions are not happening due to a lesser number of reactions, the ignition delay will not be calculated. Catch this exception using NameError.
• Check whether the calculated values fall within the tolerance limits of the reference values. If all the calculated values meet the tolerance criteria, store the number of reactions as the reactions for the reduced mechanism.
• Continue the calculations for all the reactions in all the cases and plot the results.

IMPLEMENTING THE PROGRAM IN PYTHON USING CANTERA MODULE:

The input state variables are provided as:

Temperature		950 1050 1150
Pressure		3 4 5
Eq_ratio		0.5 1 1.5

The program contains different functions so that the program looks clearer:

"""
Program to reduce reaction mechanism for the GRI 3.0 mechanism


Program by: Shrey Shah
Date: 26-Oct-19
"""
import os
import cantera as ct
import matplotlib.pyplot as plt
import numpy as np
#=====================================================================================#
# Python class to import data from file
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]] = [float(x) for x in line.split()[1:]]
#=====================================================================================#
# Function to perform 0D reactor simulation - "ref" is boolean to perform sensitivity
# analysis or not
def zeroD_reactor(gas, T, P, phi, t, ref):
    gas.TP = T, P
    gas.set_equivalence_ratio(phi, 'CH4', 'O2:1.0, N2:3.76')
    r = ct.IdealGasReactor(gas)
    sim = ct.ReactorNet([r])
    sim.rtol = 1e-06
    sim.atol = 1e-15
    Tmax = 0
    t_ig_flag = True
    # Reactor with sensitivity analysis to return max. temp., ignition delay and 
    # sensitivities
    if ref:
        for i in range(gas.n_reactions):
            r.add_sensitivity_reaction(i)

        sim.rtol_sensitivity = 1e-06
        sim.atol_sensitivity = 1e-06

        sT_max = [0]*gas.n_reactions
        for time in t:
            sim.advance(time)
            for j in range(gas.n_reactions):
                sT = sim.sensitivity('temperature', j)
                if np.abs(sT) >= np.abs(sT_max[j]):
                    sT_max[j] = sT
            if gas.T >= Tmax:
                Tmax = gas.T
            if gas.T >= T + 400 and t_ig_flag:
                t_ig = time
                t_ig_flag = False
        return [Tmax, t_ig, sT_max]
    # Only the reactor simulation to return max. temp. and ignition delay
    else:
        for time in t:
            sim.advance(time)
            if gas.T >= Tmax:
                Tmax = gas.T
            if gas.T >= T + 400 and t_ig_flag:
                t_ig = time
                t_ig_flag = False
        return [Tmax, t_ig]
#=====================================================================================#
# Function to perform sorting of the reactions based on all combinations of state
# variables
def group_sorting(rxn):
    Rid_sorted = []
    for j in range(len(rxn[0])):
        for i in range(len(rxn)):
            if rxn[i][j] not in Rid_sorted:
                Rid_sorted.append(rxn[i][j])
    return Rid_sorted

#=====================================================================================#
# Main code
#=====================================================================================#
a = FileReader('input.in')
a.read()
T_in = [x for x in a.input_dict['Temperature']]
P_in = [1e05*x for x in a.input_dict['Pressure']]
phi_in = [x for x in a.input_dict['Eq_ratio']]
case_dict = {}
case = 1

# Establishing the cases based on combinations of state variables
for T in T_in:
    for P in P_in:
        for phi in phi_in:
            case_dict.update({case:[T, P, phi]})
            case += 1

# Initializing the gas, reactor and reactor network
gas = ct.Solution('gri30.xml')
rxn = [gas.reaction_equations() for i in range(len(case_dict))]
rid = [np.linspace(1,gas.n_reactions,gas.n_reactions).astype(int) for i in range(len(case_dict))]
spec = ct.Species.listFromFile('gri30.xml')
rxn_dict = {}
for i in range(1,gas.n_reactions+1):
    rxn_dict.update({i: gas.reaction_equations()[i-1]})

# Calculating the reference values for each case along with sensitivity analysis
t = np.linspace(0, 5, 5001)
T_ref = []
t_ig_ref = []
sT = [[0 for i in range(gas.n_reactions)] for j in range(len(case_dict))]
for case in case_dict:
    print('Calculating reference values for case '+str(case)+'...', sep=' ', end='', flush=True)
    T, P, phi = case_dict[case]
    Tmax, t_ig, sT_max = zeroD_reactor(gas, T, P, phi, t, True)
    T_ref.append(Tmax)
    t_ig_ref.append(t_ig)
    sT[case-1] = sT_max
    rid[case-1] = [r for s,r in sorted(zip(np.abs(sT_max),rid[case-1]), reverse=True)]
    rxn[case-1] = [r for s,r in sorted(zip(np.abs(sT_max),rxn[case-1]), reverse=True)]
    print(' [Done]')

# Sorting the reactions by sensitivities for each case
Rid_sorted = group_sorting(rid)
Reactions = [rxn_dict[i] for i in Rid_sorted]

# Calculating the max. temp. and ignition delay by cumulatively adding reactions to the
# mechanism for each case
T_case = [[0 for i in range(gas.n_reactions)] for j in range(len(case_dict))]
t_ig_case = [[0 for i in range(gas.n_reactions)] for j in range(len(case_dict))]
T_tol = 0.05
t_ig_tol = 0.05
min_red = 40
red_flag = True
mech_red = [0 for i in range(len(case_dict))]
ct_add = 1
rxns = []
for case in case_dict:
    for r in Rid_sorted:
        rxns.append(gas.reaction(r-1))
        gas1 = ct.Solution(thermo='IdealGas', kinetics='GasKinetics', species=spec, reactions=rxns)
        T, P, phi = case_dict[case]
        print('Case: '+str(case) + ', Reactions: '+ str(len(rxns)))
        # Catching exceptions in case no reaction occurs and ignition delay is not
        # calculated
        try:
            Tmax, t_ig = zeroD_reactor(gas1, T, P, phi, t, False)
            T_case[case-1][len(rxns)-1] = Tmax
            t_ig_case[case-1][len(rxns)-1] = t_ig
            T_err = np.abs(Tmax - T_ref[case-1])*100/T_ref[case-1]
            t_ig_err = np.abs(t_ig - t_ig_ref[case-1])*100/t_ig_ref[case-1]
            if T_err <=T_tol and t_ig_err <=t_ig_tol and ct_add >=min_red and red_flag:
                mech_red[case-1] = ct_add
                red_flag = False
        except NameError:
            pass
        ct_add += 1
    ct_add = 1
    rxns = []
    red_flag = True

# Plotting the results
count = np.arange(40,325,step=1)
for case in case_dict:
    fig, ax = plt.subplots()
    plt.subplot(2,1,1)
    title = 'Case '+str(case)+': T = '+str(case_dict[case][0])+' K, P = '+str(format(1e-05*case_dict[case][1],'0.1f'))+' bar, $\phi$ = '+str(case_dict[case][2])
    plt.title(title,fontweight='bold', fontsize = 16)
    plt.plot(count,t_ig_case[case-1][40:],color='blue',linewidth=2, label='Calculated Ignition Delay')
    plt.hlines(y=t_ig_ref[case-1],xmin=30,xmax=350,linewidth=2.5, linestyle='--',color='black', label='Reference Value = '+str(format(t_ig_ref[case-1],'0.3f')) + ' sec')
    plt.grid('both',linestyle='--',linewidth=1.2)
    plt.legend(fontsize=11)
    frame1 = plt.gca()
    frame1.axes.xaxis.set_ticklabels([])
    plt.ylabel('Ignition Delay [sec]',fontsize=10.5,fontweight='bold')
    plt.subplot(2,1,2)
    plt.plot(count,T_case[case-1][40:],color='red',linewidth=2, label='Calculated Max. Temperature')
    plt.hlines(y=T_ref[case-1],xmin=30,xmax=350,linewidth=2.5 ,linestyle='--',color='black', label='Reference Value = '+str(format(T_ref[case-1],'0.3f')) +' K')
    plt.grid('both',linestyle='--',linewidth=1.2)
    plt.legend(fontsize=11)
    plt.xlabel('No. of Reactions',fontsize=12,fontweight='bold')
    plt.ylabel('Max. Temperature [K]',fontsize=10.5,fontweight='bold')
    plt.tight_layout()
    # pic_caption = 'Case '+str(case)+'.png'
    # plt.savefig(pic_caption, format='png', dpi=600)
plt.show()

`{:(ul{text{Case No.}},,ul{text{Temp. [K]}},,ul{text{Pressure [bar]}},,ul{text{Eq. Ratio}},,ul{text{Reduced # of Reactions}}),(1,,950,,3,,0.5,,41),(2,,950,,3,,1,,112),(3,,950,,3,,1.5,,112),(4,,950,,4,,0.5,,40),(5,,950,,4,,1,,132),(6,,950,,4,,1.5,,112),(7,,950,,5,,0.5,,40),(8,,950,,5,,1,,74),(9,,950,,5,,1.5,,119),(10,,1050,,3,,0.5,,42),(11,,1050,,3,,1,,72),(12,,1050,,3,,1.5,,119),(13,,1050,,4,,0.5,,49),(14,,1050,,4,,1,,71),(15,,1050,,4,,1.5,,76),(16,,1050,,5,,0.5,,73),(17,,1050,,5,,1,,183),(18,,1050,,5,,1.5,,89),(19,,1150,,3,,0.5,,61),(20,,1150,,3,,1,,183),(21,,1150,,3,,1.5,,76),(22,,1150,,4,,0.5,,183),(23,,1150,,4,,1,,71),(24,,1150,,4,,1.5,,76),(25,,1150,,5,,0.5,,72),(26,,1150,,5,,1,,71),(27,,1150,,5,,1.5,,124):}`

We can see that because of the reaction sorting based on all the cases, the initial reactions seem to reach the reference value for the unreduced mechanism. The minimum number of reactions to be included in the reduced mechanism is set to `40` to avoid very small number of reactions appearing in the reduced mechanism.