Senstivity Analysis - GRI Mechanism

Introduction:

Sensitivity analysis given an overview of magnitude of variation of an independent variable over dependent variable. It is the study that helps in investigation of factors that has the potential to alter the given outcome.

Since our study is focussed on the chemical reactions, the reactants and their initial temperature and concentration are the input parameters and raise in temperature is the desired output parameter.

Objective:

Present study focusses on finding out the 10 most reactions that can affect our desired raise in temperature of products. 

The study will focus on making the program parametric for N no. of reactions.

Simulation:

Python has been used as the general programing language, using cantera library for combustion reactions formulation, below is the code prepared for the given objective.

The code:

"""
The given program uses GRI 3.0 mechanism for formulating combustion calculation.


"""

# Importing Modules
import cantera as ct
import matplotlib.pyplot as plt
import numpy as np
import warnings
warnings.filterwarnings("ignore")

# Initial conditions for the reaction mechanism
# Initializing the gas variable as a solution object
gas = ct.Solution('gri30.xml')
# Setting temperature variable in Kelvin of the gas object 
T0 = 1500
# Setting pressure variable of the gas object to be one atmospheric pressure
P0 = ct.one_atm
gas.TPX = T0, P0, {'CH4':1, 'O2':2, 'N2':7.52}

# Initializing the reactor and the reactor network
r = ct.IdealGasConstPressureReactor(gas)
# Defining Reactor net
sim = ct.ReactorNet([r])
# Saving the reaction equations
reactions = gas.reaction_equations()

# Getting the total no. of reactions
n = len(reactions)


# Given no. of sensitive reactions to be displayed 
N = 10  
for i in range(n):
	r.add_sensitivity_reaction(i)

# Setting the relative and absolute tolerances for the solver
sim.rtol = 1e-06
sim.atol = 1e-15
sim.rtol_sensitivity = 1e-06
sim.atol_sensitivity = 1e-06

# Initializing the maximum temperature sensitivities of each reaction
sT_max = [0]*n

# Initializing a solution array to store the reaction properties at each time-step
states = ct.SolutionArray(gas, extra=['t_in_ms'])

# Creating the time series for the reactor network
t = np.linspace(0, 2e-3, 401)

# Outer time loop
for i in range(len(t)):
	sim.advance(t[i])
	states.append(r.thermo.state, t_in_ms=t[i]*1e3)
	# Storing the maximum sensitivities of each reaction
	for j in range(n):
		sT = sim.sensitivity('temperature',j)
		if abs(sT) >= abs(sT_max[j]):
			sT_max[j] = sT

# Sorting the reactions and reaction IDs based on absolute value of sensitivities
reactions = [r for s,r in sorted(zip(np.abs(sT_max),reactions), reverse=True)]
sT_max = [sT for abs_sT, sT in sorted(zip(np.abs(sT_max),sT_max), reverse=True)]


def sensitivity_n(N):
    # Extracting the top 10 sensitive reactions
    # Storing plot variables in an array
    reaction_plot = reactions[0:N]
    splot = sT_max[0:N]

    for i in range(N):
        print("Reaction ",str(i+1)," => ",reaction_plot[i],'Sensitivity Value =>',splot[i])


    # Extracting the mole fractions of OH, H and CH4
    OH = states.X[:, states.species_index('OH')]
    H = states.X[:, states.species_index('H')]
    CH4 = states.X[:, states.species_index('CH4')]

    # Setting indices for x axis and y-axis
    xticks = [-20, -15, -10, -5, 0, 5, 10, 15, 20]
    y_pos = np.arange(len(reaction_plot))

    # Plotting the Temperature Sensitivity plot

    fig1, ax1 = plt.subplots()

    # A bar chart 
    ax1.barh(y_pos, splot, 0.5, align='center', color='blue')
    ax1.set_yticks(y_pos)
    # Setting chemical reactions as the y labels
    ax1.set_yticklabels(reaction_plot, rotation=0, fontweight='bold')
    ax1.invert_yaxis()
    # Setting the ticks below the axis
    ax1.set_axisbelow(True)
    ax1.set_xticklabels(xticks, rotation=0, fontsize = 12)
    ax1.set_xlabel('Temperature Sensitivity', fontsize = 18, fontweight='bold')
    plt.title('Sensitivity Plot n for mechanism of CH4 combustion n 10 Most Temperature Sensative Reactions', fontsize=20, fontweight = 'bold')
    fig1.set_size_inches(18.5, 10.5, forward=True)
    plt.show(block=False)

    # Plotting the mole fraction as a function of reaction time
    #
    fig2, ax2 = plt.subplots(2,2)
    # Plotting concentration of CH4 with time in ms
    plt.subplot(2,2,1)
    plt.plot(states.t_in_ms,CH4,linewidth=2)
    plt.grid('both',linestyle='--',linewidth=0.8)
    plt.xlabel('Time [ms]',fontweight='bold')
    plt.ylabel(' CH4',fontweight='bold')

    # Plotting concentration of OH with time in ms
    plt.subplot(2,2,2)
    plt.plot(states.t_in_ms,OH,linewidth=2)
    plt.grid('both',linestyle='--',linewidth=0.8)
    plt.xlabel('Time [ms]',fontweight='bold')
    plt.ylabel('OH',fontweight='bold')

    # Plotting concentration of H with time in ms
    plt.subplot(2,2,3)
    plt.plot(states.t_in_ms,H,linewidth=2)
    plt.grid('both',linestyle='--',linewidth=0.8)
    plt.xlabel('Time [ms]',fontweight='bold')
    plt.ylabel('H',fontweight='bold')

    # Plotting Temperature with time in ms
    plt.subplot(2,2,4)
    plt.plot(states.t_in_ms,states.T,linewidth=2)
    plt.grid('both',linestyle='--',linewidth=0.8)
    plt.xlabel('Time [ms]',fontweight='bold')
    plt.ylabel('Temperature [K]',fontweight='bold')

    plt.suptitle("Plotting the Molar Fraction n of different Species and Temperature n with Reaction-Time")
    # For seperating the plot for visibility
    plt.tight_layout()
    fig2.subplots_adjust(top=0.88)

    plt.show()






if __name__ == "__main__":
    # N is the no. of sensative reactions to be studied
    # the value of N defaults to 10 if none input is given. By pressing enter, 
    # the code can be fully executed with default value of N = 10
    N = int(input("Enter the no. of Temperature Sensative reactions for analysis: ") or "10")
    sensitivity_n(N)

Sensitivity Plot

Having temperature as the desired output, the most sensative chemical reaction is  CH3 + O2 <=> CH3O + O. Thus, if we somehow manage to increase the rate of this chemical reaction, we can expect a relatively high order of magnitude in rise of temperature of products.

Concentration plot

As it is evident from the plot, after some igition delay, the temperature rises that signifies combustion taking place. Same can be seen with the decreasing concetration of fuel (ch4) and increasing products (H and OH).