GRI Mech Senstivity Analysis

Objective: -

In this assignment, we will be writing a Python/Cantera Code to find the 10 most sensitive reactions out of all the reactions in GRI mechanism.The main parameters of this assignment are given below: -

• Write a code that takes all the reactions from the GRI mechanism and calculates 10 most sensitive reactions
• The sensitivity parameter should be with respect to the Temperature
• The program should show the reaction strings in the order of greatest sensitivity to the least sensitivity

Solution: -

Sensitivity Analysis: -

Sensitivity Analysis is a process to get the relationship between the parameters and the output of a model.

Ex: Effect of reaction 1 rate on the temperature or the species concentration.

A generic ODE system is represented as: -

`(dbbZ)/dt = F(bbZ,t;a)`

where `bbZ = (T, Y_1, Y_2,....Y_k)^T` is a vector of species and temperature with respect to time and

`a` is the net reaction rate multiplier (Perturbation Factor)

The ODE has paremeters which are represented by `a`. For example, if we have hundred reactions, then the rate constants for each of those reactions are our parameters. If we change the parameters, our ODE's function will also change respectively

To have a better understanding of Sensitivity, let's look at the following graph:

Here, the X-axis and Y-axis are represented by Time and Variable (Temperature, Concentration, etc...) respectively. Let's say that for a given reaction parameter, we get the curve c1 as a function of time. If we change the parameter of c1 with `deltap` which is `p + deltap`, we obtain the c2 curve due to the change. So, we have a change of `deltay` in variable for the c2 curve when compared to the c1 curve. If we take the ratio of change in the variable to the change in the parameter `(deltay)/(deltap)`, we obtain the sensitivity of `y` with respect to `p`.

For the generic ODE discussed , we have a sensitivity matrix which is represented as : -

`w_(j,i) = (delZ_j)/(dela_i)`

where `Z_j` represents the corresponding variable and

`a_i` represents the reaction rate constant where `i` is the reaction

So, if we have n reactions with variables k species and temperature, the size of the sensitivity matrix would be `w_(k+1, n)`

Application: -

Reaction mechanisms aren’t always made with specific parameters in mind. For example, say we have a case of high EGR (Exhaust Gas Recirculation) conditions. The mechanism may not be valid for that condition or it may not give the desired results. To get the desired results, we have to change the co-efficient of the reactions i.e. tune the reactions to our problem specific parameters. But in the case of GRI mechanism, we have around 325 reactions and we cannot change them individually which would be a highly inefficient.

The purpose of Sensitivity Analysis is that we can tune the most sensitive reaction mechanisms i.e. we can adjust the rate co-efficient so that it will give us a match for the Ignition Delay. Engineers use this method to tune the reaction mechanisms to their problem specific parameters.

Program Structure and Code: -

In this assignment, we are calculating the effect of all the reactions from the GRI mechanism file on the temperature parameter of sensitivity. We are using the Ideal Gas Constant Pressure reactor for this problem. The initial conditions for this reactor are given below:

Temperature: 1500 Kelvin

Pressure: 1 bar (101325 Pascals)

Species and their mole fractions: 'CH4':1, 'O2':2, 'N2':7.52

Initially, we are going to add the reactions to the reactor for the analysis. Then we set the tolerances for the solution and the sensitivity co-efficient. There are two solution tolerances namely Relative, Absolute for which the values are 1.0e-16 and 1.0e-15. These values constrain the solution stating that any change in the solution should be less than the given values. Similarly, we set the tolerances for sensitivity co-efficient which is 1.0e-6 for both Absolute and relative change.

We are considering the analysis duration of about 2 milliseconds which means that we start our time integration from 0 to 2 milliseconds with every 5-microsecond time step and find the sensitivity with resoect to the temperature parameter for all the reactions to store them in an 2D array (Sensitivity Matrix) with respect to the particular time step.

In the current problem, we are trying to find the most sensitive reactions considering that the sensitivity is positive and negative. After obtaining the sensitivity co-efficients, we find modulus of every co-efficient and associate it to its respective reaction. Finally, we sort the reactions in the order of greatest to least sensitivity and produce the top ten.

The code for the above program is given below: -

# Constant Pressure, Adiabatic kinetics simulation with sensitivity analysis

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

gas = ct.Solution('gri30.xml')
temp = 1500
press = 101325

n_reactions = 325 # Number of Reactions

time = np.arange(0, 2e-3, 5e-6) # Simulation Time
interval = 5e-6
n_time = len(time)

gas.TPX = temp, press, {'CH4':1, 'O2':2, 'N2':7.52}
r = ct.IdealGasConstPressureReactor(gas, name='R1')
sim = ct.ReactorNet([r])

# Enable sensitivity with respect to te rates of all the 325 reactions
# reactions (0 through 324)
for i in range(n_reactions):
	r.add_sensitivity_reaction(i)
	

# set the tolerances for the solution and the sensitivity co-efficients
sim.rtol = 1.0e-6 # Realtive change in the solution < sim.rtol
sim.atol = 1.0e-15 # Absolute change in the solution < sim.atol

sim.rtol_sensitivity = 1.0e-6 #relative change in sensitivity < sim.rtol_sense
sim.atol_sensitivity = 1.0e-6 #absolute change in sensitivity

s = np.zeros((n_time, n_reactions))
t = 0

for n in range(n_time): # Iterates for every time step

	t = t + interval
	sim.advance(t) # Increments the reactor by the actual time (not time step)

	for j in range(n_reactions):
	
		s[n,j] = (sim.sensitivity(1, j)) # Sensitivity of Temperature to 100 reactions

		print('%10.3e %10.3f %10.3f %14.6e %s' %
			(sim.time, r.T, r.thermo.P, r.thermo.u, s[n,j]))


def max_abs(array):

	b_ar = np.abs(array) # Finds the Absolute values
	b_arr_f = np.max(b_ar) # Finds the maximum value

	return b_arr_f


b_array = []
fin_sens = [] # Final Sensitivities Array

for i in range(0,n_reactions):
	
	# Stores reaction names and the absolute maximum sensitivities for the respective reaction
	
	b_array = s[:,i]
	
	fin_sens.append((sim.sensitivity_parameter_name(i), max_abs(b_array)))

dtype = [('reaction name', 'S40'), ('max sensitivity', float)]

fin_sens_temp = np.array(fin_sens, dtype=dtype)

final_sort = np.sort(fin_sens_temp, order='max sensitivity')[::-1] # Sorts the reactions in descending order
print(final_sort)

Results: -

The above snapshot shows the top ten reactions which are sensitive to temperature with their sensitivity co-efficients.

The above figure gives us a pictorial representaion of the top ten sensitive reactions for Temperature in GRI mechanism. The modulus of sensitivity co-efficient and the top reactions are represented in x and y axes respectively. These results tell us that of all the reactions in the GRI mechanism, `CH3 + O2 <=> CH3O + O` reaction has the highest temperature senitivity.

When we obtain a positive sensitivity co-efficient, the temperature increases. Due to this, the Ignition Delay decrease because the ignition is happening at a faster rate. Similarly, when we have a negative sensitivity co-efficient, the burning is slow and respectively the Ignition Delay inreases.

To increase or decrease Ignition Delay, we can refer to the particular sensitivity co-efficients and change the rate co-efficients accordingly.