Senstivity Analysis Assignment

Intro

Chemical reactions such as CH4 + 2O2 + 7.52 N2- > CO2 + 2H2O + 7.52N2 are not a single reaction mechanism. Instead, the reactants have to go through many intermediary reactions for the products to form.

For the reaction equation shown above, i.e; the combustion of Methane, there are almost 300+ intermediatory reactions, with each reaction having its own reaction rate and reaction constant. Though the chemical composition and conditions formed at equilibrium are always the same, no matter the path followed, but the rate at which the chemical composition and condition are achieved is highly dependent on the reaction rates of the intermediatory reactions.

It should also be kept in mind that ignition also happens during one of these intermediatory reactions and thus ignition delay also depends on the reaction rates of intermediatory reactions.

In this project, we are going to find the intermediatory reactions which are most contributing in nature to the rise/fall of temperature, in turn, speeding/delaying the ignition process.

We use Cantera Library to import the GRI30 Mechanism. It is a detailed kinetic mechanism for methane/air combustion including nitric species for NOx predictions with 53 species and 325 reactions.

Code and its explanation

Importing Libraries

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

- Cantera for Chemical Reaction Mechanisms and its inbuilt classes which are used to calculate sensitivity and other properties.

- Numpy for Array operations and functions like linspace etc.

- Matplotlib.pyplot for plotting of parameters.

Importing Gri30 Mechanism and defining initial conditions

gas=ct.Solution('gri30.yaml')
T=1500
P=ct.one_atm
species={'CH4':1,'O2':2,'N2':7.52}
gas.TPX=T,P,species

- It is ideal to use .yaml format mechanism files instead of .xaml, as the latter is going to be discarded in the upcoming versions of Cantera, rendering it useless.

- Define initial conditions using .TPX

Defining and Solving ODE's

r=ct.IdealGasConstPressureReactor(gas)
sim=ct.ReactorNet([r])
states=ct.SolutionArray(gas,extra=['t','sens'])

- ODE's are defined using the Reactor class. In this case, we are using an IdealGasConstPressureReactor class.

- ReactorNet class is used to solve different Reactor classes simultaneously. ReactorNet class must be initialized even if there is a single Reactor class.

- The solution obtained by solving the ODE's are saved in the SolutionArray class

Reactions in GRE30

#Total number of reactions in gas
TR=gas.n_reactions
reactions = gas.reaction_equations()

- n_reactions gives the total number of reactions present in the GRE30mechanism

- reaction_equations gives the equation of the reaction present in the GRE30mechanism

Input

n=int(input('How many most sensitve reactions to plot ? : '))

- Asking the user to feed in the input for "How many sensitive reactions to plot"

Adding reactions for sensitivity study and defining tolerances

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

sim.rtol=1e-6
sim.atol=1e-15
sim.rtol_sensitivity=1e-6
sim.atol_sensitivity=1e-6

- atol - Absolute tolerance = absolute (a-b)

- rtol - Relative tolerance = absolute(A-B)/min(A,B).

Relative tolerances are relative to the solution value, whereas absolute tolerances are the maximum absolute error allowed in a solution.

Advancing solution in time

time=1e-3

for t in range(0,500):
    sim.advance(time)    
    s=[]
    
    for j in range(TR):
        s.append(sim.sensitivity('temperature',j))
        
    states.append(r.thermo.state,t=time,sens=s)
    time=time+1e-3

- Advancing using the advance method

- Appending the sensitivities to s

- Appending the s, temperature, and time to states

Finding the most sensitive equations

F=[]
temp={}

for k in range(TR):
    abc=np.transpose(states.sens)[k]
    
    if max(abc)>=max(abs(abc)):
        ijk=max(abc)
    else:
        ijk=max(abs(abc))
        ijk=-ijk

    F.append(ijk)
    temp[k]=F[k]

- States is a matrix of size 500*325 i.e, number of time steps * number of reactions

- Transpose the matrix to make sensitivities at each time step for a reaction into rows. i.e; Row 1 = sensitivity of reaction 1 at each timestep

- if else loop: Finding the extreme value for each reaction from all timesteps and saving them in variable ijk

- Append those values to F (F has both negative and positive sensitivity values)

- Create a dictionary called temp, with keys as reaction numbers and values as extreme of sensitivity values for each reaction obtained from F

Searching for the most sensitive reactions

temp_raw=temp
temp_raw_values=list(temp.values())
temp_pos={key: abs(val) for key, val in temp.items()}
sorted_temp_pos=dict(sorted(temp_pos.items(),key=lambda x:x[1],reverse=True))
rxn_no=list(sorted_temp_pos.keys())[0:n]

- Dummy of temp is created as temp_raw

- Making all the values of sensitivities positive and sorting them in descending order

- Storing the reaction number of top 'n' sensitive reactions in variable rxn_no

sens1=[]
for i in range(len(rxn_no)):
    sens1.append(temp_raw_values[rxn_no[i]])

equation=[]
for i in range(len(rxn_no)):
    equation.append('Reaction ' + str(rxn_no[i]) + ' : ' + gas.reaction_equations()[int(rxn_no[i])])

- Getting the sensitivity values from temp raw using their reaction numbers from rxn_no variable

- Reaction equations of the most sensitive reactions stored in the equation variable

Printing output

rxn_no=[str(i) for i in rxn_no]
print('The most sensitive reactions are : ')
[print(equation[i]) for i in range(len(equation))]

Plotting

plt.title('Most sensitive ' + str(n) + ' reactions to Temperature')
plt.barh(equation,sens1)
plt.gca().invert_yaxis()
plt.xlabel('Sensitivity')
plt.ylabel('Reaction')
plt.show()

Output

Inference/Conclusion

- The reactions with sensitivities in the positive range help in reducing the ignition delay, thus increasing the combustion efficiency

- The reactions with sensitivities in the negative range increase the ignition delay, thus reducing the combustion efficiency

- Increasing the rate of reaction number 154: CH3 + O2 <==> CH3O + O has a significant effect on reducing the ignition delay.

To read more about ignition delay and the effect of initial conditions on ignition delay, check this project - Auto ignition using Cantera