cs207 Team 10 Final Project: Chemical Kinetics library
Team Members: Hidenori Tanaka, Jiachen Song and Xiangru Shu
Introduction
This program is a chemical kinetics library, which could be easily installed by users and used for various applications.
Basic usage includes:
- Constant reaction rate coefficient
- Arrhenius reaction rate coefficient
- Modified Arrhenius reaction rate coefficient
- Backward reaction rate coefficient
- Progress rate for irreversible and reversible reactions
- Reaction rate for species
Reaction rate coefficients
The library is able to compute 4 different kinds of reaction rate coefficient: constant reaction rate coefficient, arrhenius reaction rate coefficient, modified Arrhenius reaction rate coefficient and backward reaction rate coefficient as specified below.
Each variable stands for A: Arrhenius prefactor, b: Modified Arrhenius parameter, E: Activation Energy, T: Temperature, and R: Ideal gas constant.
For an elementary reaction, we have
where is the equilibrium coefficient for reaction . The final expression for the equilibrium coefficient is,
where and is the pressure of the reactor (take it to be Pa). is the entropy change of reaction and the enthalpy change of reaction .
Progress rate for irreversible and reversible reactions
We used the principle of mass action to obtain the progress rate of each reaction.
Basically, the progress rate of a reaction is proportional to the concentrations of the reactants. Thus, the forward progress rate is:
where is the forward reaction rate coefficient.
For irreverisible reactions, the total progress rate is equal to the forward progress rate. But in reality, it is often the case that the products can react and produce the reactants. This is called a reversible reaction.
For reversible reaction, the total progress rate is:
where is the backward reaction rate coefficient.
Reaction rate for species
The reaction rate for each species was given as a linear combination of reaction progress rates.
That is,
Installation
All the classes and functions necessary to run the library were stored in chemkin_g10. All the tests were stored in tests folder.
To install program as library (make sure you have pip installed for python3.5+)
pip install chemkin_g10
Two samples were provided in the samples directory: irreversible.py and reversible.py.
To run our samples, go to the samples directory, and run
python irreversible.py
python reversible.py
To run our test program, simply run
pytest
Basic Usage and Examples
Our library includes four separate modules: chemkin
, computation
, db
and thermo
.
chemkin module
chemkin
is a module that contains a ReactionSystem
class, a Reaction
class and a Simulator
class.
-
ReactionSystem
class:buildFromList
buildFromXml
getProgressRate
getReactionRate
parse
-
Reaction
class:updateCoeff
updateReaction
-
Simulator
class:solveODE
check_equilibrium
equilibrium_graph
plot_specie
plot_specie_all
plot_reaction_all
computation module
computation
module contains all necessary functions to compute the coefficients and rates.
It includes functions:
k_const
k_arr
k_mod_arr
progress_rate
reaction_rate
equilibrium_constant
thermo module
thermo
includes functions: H_over_RT
, S_over_R
, backward_coeffs
db module
db
contains the functions to read NASA polynomials from the database.
Examples
A typical workflow starts from initializing a ReactionSystem
object. We need to set up all the needed variables: T
(temperature), R
(universal gas constant) and concs
(the concentration of each species, and the order should be same with the one in the input file). Here's an example:
from chemkin_g10 import chemkin as ck
import numpy as np
T = 1500
R = 8.314
concs = np.array([1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0])
rsystem = ck.ReactionSystem(T, R, "tests/data/db/nasa.sqlite") # your path to db file
Then, we feed an input .xml
file (see example1 or example2 for the format) and a numpy array of concerntrations of each species to the system by calling buildFromXml
function:
rsystem.buildFromXml("tests/data/xml/rxns_reversible.xml", concs) # your path to xml file
We then can easily retrieve the progress rate and reaction rate by:
print("Progress rate: \n", rsystem.getProgressRate(), "\n")
print("Reaction rate: \n", rsystem.getReactionRate(), "\n")
If we want to get a summary of each reaction and also the whole system, we can do this by:
print("System info: \n", rsystem, "\n")
Extensibility
-
Reversible/Non-elementary reactions
We store these types of information in the metadata of each reaction, for example:
rsystem.reactionList[0].reactMeta >> {'reversible': 'no', 'type': 'Elementary', 'id': 'reaction01'}
-
Reaction rate coefficients not discussed in class
We can easily add a new method of computing reaction rate coefficient to our
computation
module. -
Other extensions
We can update the property (metadata) of a single reaction or its reaction rate coefficient (re-compute the coefficient), with any valid parameters. Here's an example of changing the coefficient of the 3rd reaction to be the same as the 1st reaction:
print(rsystem.reactionList[2].k, rsystem.reactionList[0].k) # using the parameters of the first equation rsystem.reactionList[2].updateCoeff(type="modifiedArrhenius", A=100000000.0, b=0.5, E=50000.0) # both ks are same now print(rsystem.reactionList[2].k, rsystem.reactionList[0].k)
4484938.47318 70279405.1912 70279405.1912 70279405.1912
Then rebuild the system with the new variables:
rsystem.buildFromList(rsystem.reactionList, ["H","O","OH","H","H2O","O2","HO2","H2O2"], np.array([1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]))
New Feature
We implemented a new class Simulator
under the chemkin
module. The class includes an ODE solver function, two equilibrium check functions and three plot functions.
To initialize a Simulator
object, users should pass in a ReactionSystem
object, an ending time for simulation (start at 0), a number of integration time steps(default = 100), a time scale(default = 1e9, s to ns) and an equilibrium threshold(default = 1e-05).
A new example:
from chemkin_g10 import simulator as sim
from chemkin_g10 import chemkin as ck
import numpy as np
T = 900
R = 8.314
concs = np.array([0.5, 0, 0, 2, 0, 1, 0, 0])
rsystem = ck.ReactionSystem(T, R, "tests/data/db/nasa.sqlite")
rsystem.buildFromXml("tests/data/xml/rxns_reversible.xml", concs)
# Initialize a new simulation from time 0 to 0.1/1e9
simulation = sim.Simulator(rsystem, 0.1, numSample=100, timeScale=1e9, eqThreshold=1e-05)
ODE Solver
Based on Le Chatelier's principle, changing the concentration of a chemical will have an effect to the reaction equilibrium. Thus, the reaction rate, extent, and yield of products will be altered corresponding to the impact on the system.
For irreversible reactions, users may want to know the progress of the reaction by checking the concentration for each reactant and product. For reversible reactions, users may also want to check the concentration for each specie, in order to determine whether the reaction has reached the equilibrium.
Thus, we added a new feature to compute the concentration for each specie at each time stamp between time interval (0,t), and t will be entered by users.
To implement this feature, we added a new function solveODE
under the Simulator
class. It has a nested function fun
, which is to compute reaction rates. Then we used odeint
function from library scipy
to integrate reaction rates and then get concentrations for each specie at each time stamp.
The formula of reaction rate for each specie is specified below:
To use this function, users can simply call:
simulation.solveODE()
Users can also get all the concentrations at each time stamp by calling the following methods:
print(simulation.yout)
...
[ 1.62865785e+00 9.55229994e-02 9.03020740e-03 2.06975115e-01
1.22418085e+00 3.35632969e-01 4.49199685e-15 1.58667651e-20]
[ 1.62865785e+00 9.55229994e-02 9.03020740e-03 2.06975115e-01
1.22418085e+00 3.35632969e-01 4.49199685e-15 1.58667651e-20]
[ 1.62865785e+00 9.55229994e-02 9.03020740e-03 2.06975115e-01
1.22418085e+00 3.35632969e-01 4.49199685e-15 1.58667651e-20]]
The last array of concentrations is the final concentrations for each specie at different timestamps.
Equilibrium
In chemistry, chemical equilibrium is the state in which both reactants and products are present in concentrations which have no further tendency to change with time. Usually, this state results when the forward reaction proceeds at the same rate as the reverse reaction.
To check equilibrium, we developed two methods. The first method is based on the definition of equilibrium and compare reaction quotient with equilibrium constants.
For example, a reversible reaction:
The equation for reaction quotient is written by multiplying the concentrations for the species of the products and dividing by the concentrations of the reactants. If any component in the reaction has a coefficient, indicated above with lower case letters, the concentration is raised to the power of the coefficient. The above equation is therefore:
A comparison of reaction quotient(Q) with equilibrium constant(K) indicates which way the reaction shifts and which side of the reaction is favored:
- If Q>K, then the reaction favors the reactants.
- If Q<K, then the reaction favors the products.
- If Q=K, then the reaction is already at equilibrium. There is no tendency to form more reactants or more products at this point. No side is favored and no shift occurs.
First, we added a new function equilibrium_constant
under the computation.py
module to compute the equilibrium constant(k) for each reaction.
Then, in the solveODE
function, after calling odeint
to integrate concentrations, we computed the reaction quotients for each reaction, and compare the reaction quotients with the equilibrium constants by computing the following percentage:
If the percentage is less than the equilibrium threshold, then we consider the reaction has reached equilibrium.
For each reaction, we found the first time stamp that the reaction reaches equilibrium and stored them into an array for future use.
Users can simply call the following methods to see the equilibrium time stamp for each reaction:
print(simulation.eq_point)
[3.4343434343434346e-11, 6.7676767676767674e-11, 6.8686868686868691e-11, 6.2626262626262625e-11, 6.6666666666666669e-11, 5.858585858585858e-11, 5.858585858585858e-11, 6.4646464646464647e-11, 5.4545454545454549e-11, 6.5656565656565652e-11, 5.6565656565656564e-11]
Users may also want to check equilibrium with arbitrary time stamp, so we added a new function check_equilibrium
which takes an index and a time stamp as parameters. By specifying index and time stamp, users can check equilibrium for a specific reaction at a specific time stamp. The function will simply compare this time stamp with the first equilibrium time stamp for that reaction and determine whether the reaction has reached equilibrium.
An example of checking equilibrium for reaction 5 at time 5e-11(5e-11 is less than 5.858585858585858e-11, so it has not reached equilibrium):
print(simulation.check_equilibrium(5, 5e-11))
>> False
Furthermore, we came up with another method to check equilibrium, which was developed based on slopes of concentration plots.
If the largest concentration among chemical species at time t is C, then the characteristic slope of the c(t) curves can be calculated as C/t. We judge the system to be in equilibrium if all the slopes of the concentrations at the last two time steps are less than the critical slope value "1e-7*C/t". The choice of "1e-7" is our definition of equilibrium and the number could be changed to another small number.
Example:
print(simulation.equilibrium_graph())
>> True
Basic Plot (Non-Interactive)
It is always helpful for users to see how concentrations change over time graphically. Therefore, we added plot functions to visualize concentration change and we used matplotlib
library to plot graphs.
To plot concentrations for the entire reaction system, we added function plot_specie_all
, which will plot the concentrations for all the species in the system over time.
simulation.plot_specie_all()
To plot concentration for an individual specie, we added function plot_specie
, which will take a species as parameter to specify which specie to plot.
simulation.plot_specie("H2O")
We also added a function plot_reaction_all
to plot (reaction quotients - equilibrium constants)/ equilibrium constants.
simulation.plot_reaction_all()
Web Visualization (Interactive)
(Updated: Website now deployed on AWS elastic beanstalk)
We've also implemented a web visualization. With the help of javascript, the trend of how reactions/species change during simulation is much more straightforward. In addition, it also shows the detailed information about each reaction/species. It still has some bugs as the odeint
library are inconsistent with same parameters. To access the visualization of current simulation, simply run
simulation.visualize() # the website will pop up
To get a direct view without building your ReactionSystem
and Simulator
, here's the link.
Check out the demo [here] (https://giphy.com/gifs/l3mZ5WJKdimnC1zqM/fullscreen)