Week 9 - Senstivity Analysis Assignment

"""
python code to conduct sensitivity analysis

"""

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

gas = ct.Solution('gri30.xml')

#applying initial variables
T = 1500
P = 101325
species = {'CH4':1,'O2':2,'N2':7.52}

gas.TPX = T,P,species
r = ct.IdealGasConstPressureReactor(gas,name = 'R1')    #
sim = ct.ReactorNet([r])

# enable sensitivity with respect to the rates of first 10
# reaction (reaction 0 through 9)
reaction_n = []
num_reactions = gas.n_reactions
for i in range(num_reactions):
	r.add_sensitivity_reaction(i)
	reaction_n.append(gas.reaction_equations([i]))

# set the tolerance for the solution and for the sensitivity coefficients relative and absolute tolerance
sim.rtol = 1e-6   # relative change in solution should be less than this number
sim.atol = 1e-15  # absolute change in solution should be less than this number

sim.rtol_sensitivity = 1e-6   # relative change in sensitivity should be less than this number
sim.atol_sensitivity = 1e-6   # absolute change in sensitivity should be less than this number

states = ct.SolutionArray(gas,extra = ['t_ms'])

sensitivity_temperature = [0]*num_reactions

for t in np.arange(0, 2e-3, 5e-6):
	sim.advance(t)
	
	for j in range(num_reactions):
		gradient = sim.sensitivity('temperature',j)                 

		if (abs (gradient) >= abs(sensitivity_temperature[j])):
			sensitivity_temperature[j] = gradient
	states.append(r.thermo.state, t_ms = 1000*t)

# Sorting the reactions
s_abs = np.abs(sensitivity_temperature)
reactions = [p for q, p in sorted (zip(s_abs,reaction_n), reverse = True)]
sensitivity_n = [p for q, p in sorted (zip(s_abs,sensitivity_temperature), reverse = True)]
max_reaction = reactions[0:10]
max_sensitivity = sensitivity_n[0:10]
for i in range(10):
	print('reaction' + str(max_reaction[i]) + 'sensitivity is ' + str(max_sensitivity[i]))  

#print('%10.3e %10.3f %10.3f %14.6e %10.3f %10.3f' %(sim.time,r.T, r.thermo.P, r.thermo.u, s2, s3))



plt.subplot(2,2,1)
plt.plot(states.t_ms,states.T)
plt.xlabel('Time(ms)')
plt.ylabel('Temperature(k)')
plt.title('Plot showing the Temperature against Time')
plt.subplot(2,2,2)
plt.plot(states.t_ms,states('OH').X)
plt.xlabel('Time(ms)')
plt.ylabel('OH Mole Fraction')
plt.title('Plot showing the OH concentration against Time')
plt.subplot(2,2,3)
plt.plot(states.t_ms,states('H').X)
plt.xlabel('Time(ms)')
plt.ylabel('H Mole Fraction')
plt.title('Plot showing the H concentration against Time')
plt.subplot(2,2,4)
plt.plot(states.t_ms,states('CH4').X)
plt.xlabel('Time(ms)')
plt.ylabel('CH4 Mole Fraction')
plt.title('Plot showing the CH4 concentration against Time')
plt.tight_layout()


plt.figure(2)
plt.barh(range(10),max_sensitivity,0.6,color = 'green')   #h stands for horizontalwise
plt.yticks(range(10),max_reaction)
plt.xlabel('Sensitivity')
plt.ylabel('Reactions')
plt.title('Top 10 Temperature Sensitive Reaction of grimech 3.0')
plt.gca().invert_yaxis()
plt.grid()
plt.tight_layout()
plt.show()

AIM :

Sensitivity analysis of grimech 3.0

OBJECTIVE :

• To write a code that takes all the reactions from the GRE mechanism and calculates 10 most sensitive reactions.
• To plot a graph of sensitivities vs the reactions.

SENSITIVITY ANALYSIS :

In my previous challenge i have learned about the Ignition delay and the factors which affect the ignition delay which are temperature,pressure and fuel composition. We know that in any reaction mechanism, there will be an involement of elementary reactions,the rate of these reaction effects the entire reaction mechanism. The rate of reactions are completely depends on the temperature,pressure and fuel composition. If we change any of these parameters,the entire reaction rates will be effected. So the change in the output depends on the input parameters is calculated by performing Sensitivity analysis.

CODING APPROACH :

Initially the program is written to find the sensitivity analysis for all the reaction in methane (CH4) combustion using the sourse code from the video. The number of reaction happening is found using the reaction object from gri mechanism. Then all reaction are needed to be sorted with higher to lower sensitivity. The source coding for sorting in decending order for the sensitivity is obtained from python coding blogs and sorted the sesitivity. Then first 10 n is seperated and appended in a new variable. Then they are plotted. There are 325 reaction happening. The program is written in a manner where n is variable of an int which will show the top n(int as input ) and plot the graph for them.

CODE PROCEDURE :

• Importing all the necessary modules to python like cantera,numpy,matplotlib.
• Creating a gas object which is instance of the solution class, so to instaniate the solution class we basically give the 'Gri30' mechanism.
• Give value of temperature is 1500k and pressure value is 1 atm.
• Gives Species Dict for the mole fractions of the reacting species as {'CH4' :1, 'O2': 2 , 'N2' : 7.52}
• Create sim instance for the ReactorNet class using r reactor
• Enable sensitivity for each reaction, add tolerance values for the solution as well as the sensitivity coefficients.
• What does absolute tolerance mean ,so when we doing Newton Raphson techniques we basically have to say that for each iteration the solution should not be higher than the value so that's called as the absolute tolerance. For each succesive Newton Rhapson iteration the species concentration should not change by more than 1e-15 ,if it does than it goes to iterate again.Now, there are certain species for example when a chemical reaction take place we talk about stiff system w.r.t chemical reaction stiff system is one where certain reactions take place very fast and certain reaction take place very slow.The reaction take place very slow the species part is spliting into those reactions might not change a lot, there absolute concentration might be 1e-17 like that, after one new iteration the concentration may be 2e- 17 so that's why we have 2 types of tolerance. Absolute tolerance is basically going to measures exactly how much is changed by and other is relative tolerance is basically when species concentration changes a lot like really quickly. Its basically (current value - old value) / old value.
• So in a similar way in our sensitive ODE we have these two tolerance relative and absolute and for sensitivity analysis the tolerance values are higher than our species ODE.
• Now i am going to create 'states'object now what is the state object well, that is going to store information. It is basically a data structure that helps us store the reactor state(temp,pressure,species concentration everything about the reactor at a given timestep) can be stored in our 'states' object .To use this we just need to call 'ct.SolutionArray' and i need to instantiate by a gas object and we create some extra parameters like time.
• Then we will start time integartion from start time to end time and we are advancing their network by '5e-6' at any given time. first i need to have 'advance' first and new state of reactor and based upon new state to the old state which is the initial condition and doing a sensitivity analysis.
• Create a time loop with a time step of 5e-6 sec and run time as 2e-3 sec. use advance method to to simulate the reaction at each time step store the sensitivities for all the reaction with respect to the temperature in an array and also store the maximum value of sensitivity for each reaction in another array. finally store the state for each iteration.
• Now,using the sorted zip method we can arrange all the reaction in decreasing order of their sensitivity in a given array.
• Printing the top sensitive values N(in my case N = 10) reaction on the output window.
• Plotting various points for the reaction 'time' vs 'parameters' like temperature and concentration(OH,CH4,H etc) .
• In addition to it, For 2nd figure i was plot number of top 10 sensitive reaction. for this i was plot bar graph with 'barh' command in order to make the graph horizontal.

OBSERVATION :

• By observing the above plots, it makes sense as the temperature remains almost constant equal to initial temperature upto 1ms and after the ignition take place, it raises immediately and reaches adiabatic flame temperature after sometime.
• By looking the OH plot, it also remian same for some time and raises immediately. it is an important radicle to initiate the ignition.
• By looking the CH4 , the concentration is almost equal to 1 at initial stage and simulation takes place drops which make sense likewise.

CONCLUSION :

• I have done the analysis for the sensitivity of the reaction with respect to the temperature and we have find out 10 most sensitive reactions for the auto ignition of methane.
• The code that i have written here is very parametric or u can say that we can easily change the parameter and do the analysis for the sensitivity with respect to that parameter.
• Also we can change the number for the highest sensitive reactions required.
• The positive sensitivity tells us the increase in the initial temperature will increase the rate of reaction of those reactions and the negative sensitivity is of decrease the reaction rate. If the reactions need to take place rapidly, by tuning the positive parameters to increase the reaction rates.