In [1]:
%%file reaction_coeffs.py
"""This module contains functions that return reaction rate coefficients
MEMBERS:
========
k_const: Returns a constant reaction rate coefficient
k_arr: Returns the Arrhenius reaction rate coefficient
k_mod_arr: Returns the modified Arrhenius reaction rate coefficient
"""
import numpy as np
def k_const(k=1.0):
"""Simply returns a constant reaction rate coefficient
INPUTS:
=======
k: float, default value = 1.0
Constant reaction rate coefficient
RETURNS:
========
k: float
Constant reaction rate coefficient
EXAMPLES:
=========
>>> k_const(5.0)
5.0
"""
if k < 0:
raise ValueError("Negative reaction rate coefficients are prohibited.")
return k
def k_arr(A, E, T, R=8.314):
"""Calculates the Arrhenius reaction rate coefficient
INPUTS:
=======
A: float
Arrhenius prefactor
Must be positive
E: float
Activation energy
T: float
Temperature
Must be positive
R: float, default value = 8.314
Ideal gas constant
Must be positive
RETURNS:
========
k: float
Arrhenius reaction rate coefficient
EXAMPLES:
=========
>>> k_arr(2.0, 3.0, 100.0)
1.9927962618542914
"""
if A < 0.0:
raise ValueError("A = {0:18.16e}: Negative Arrhenius prefactor is prohibited!".format(A))
if T < 0.0:
raise ValueError("T = {0:18.16e}: Negative temperatures are prohibited!".format(T))
if R < 0.0:
raise ValueError("R = {0:18.16e}: Negative ideal gas constant is prohibited!".format(R))
return A * np.exp(-E / R / T)
def k_mod_arr(A, b, E, T, R=8.314):
"""Calculates the modified Arrhenius reaction rate coefficient
INPUTS:
=======
A: float
Arrhenius prefactor
Must be positive
b: float
Modified Arrhenius parameter
E: float
Activation energy
T: float
Temperature
Must be positive
R: float, default value = 8.314
Ideal gas constant
Must be positive
RETURNS:
========
k: float
Modified Arrhenius reaction rate coefficient
EXAMPLES:
=========
>>> k_mod_arr(2.0, -0.5, 3.0, 100.0)
0.19927962618542916
"""
if A < 0.0:
raise ValueError("A = {0:18.16e}: Negative Arrhenius prefactor is prohibited!".format(A))
if T < 0.0:
raise ValueError("T = {0:18.16e}: Negative temperatures are prohibited!".format(T))
if R < 0.0:
raise ValueError("R = {0:18.16e}: Negative ideal gas constant is prohibited!".format(R))
return A * T**b * np.exp(-E / R / T)
Pretend you're a client who needs to compute the reaction rates at three different temperatures ($T = \left\{750, 1500, 2500\right\}$) of the following system of irreversible reactions: \begin{align} 2H_{2} + O_{2} \longrightarrow 2OH + H_{2} \\ OH + HO_{2} \longrightarrow H_{2}O + O_{2} \\ H_{2}O + O_{2} \longrightarrow HO_{2} + OH \end{align}
The client also happens to know that reaction 1 is a modified Arrhenius reaction with $A_{1} = 10^{8}$, $b_{1} = 0.5$, $E_{1} = 5\times 10^{4}$, reaction 2 has a constant reaction rate parameter $k = 10^{4}$, and reaction 3 is an Arrhenius reaction with $A_{3} = 10^{7}$ and $E_{3} = 10^{4}$.
You should write a script that imports your chemkin module and returns the reaction rates of the species at each temperature of interest given the following species concentrations:
You may assume that these are elementary reactions.
In [6]:
%%file chemkin.py
import numpy as np
def progress_rate(nu_react, concs, k):
"""Returns the progress rate of a system of irreversible, elementary reactions
INPUTS:
=======
nu_react: numpy array of floats,
size: num_species X num_reactions
stoichiometric coefficients for the reaction
k: array of floats
Reaction rate coefficient for the reaction
concs: numpy array of floats
concentration of species
RETURNS:
========
omega: numpy array of floats
size: num_reactions
progress rate of each reaction
EXAMPLES:
=========
>>> progress_rate_2(np.array([[2.0, 1.0], [1.0, 0.0], [0.0, 1.0]]), np.array([2.0, 1.0, 1.0]), 10.0)
array([ 40., 20.])
"""
progress = k # Initialize progress rates with reaction rate coefficients
for jdx, rj in enumerate(progress):
if rj < 0:
raise ValueError("k = {0:18.16e}: Negative reaction rate coefficients are prohibited!".format(rj))
for idx, xi in enumerate(concs):
nu_ij = nu_react[idx,jdx]
if xi < 0.0:
raise ValueError("x{0} = {1:18.16e}: Negative concentrations are prohibited!".format(idx, xi))
if nu_ij < 0:
raise ValueError("nu_{0}{1} = {2}: Negative stoichiometric coefficients are prohibited!".format(idx, jdx, nu_ij))
progress[jdx] *= xi**nu_ij
return progress
def reaction_rate(nu_react, nu_prod, rj):
"""Returns the reaction rate of a system of irreversible, elementary reactions
INPUTS:
=======
nu_react: numpy array of floats,
size: num_species X num_reactions
stoichiometric coefficients for the reactants
nu_prod: numpy array of floats,
size: num_species X num_reactions
stoichiometric coefficients for the products
k: float, default value is 10,
Reaction rate coefficient for the reaction
concs: numpy array of floats
concentration of species
RETURNS:
========
f: numpy array of floats
size: num_species
reaction rate of each specie
EXAMPLES:
=========
>>> nu_react = np.array([[1.0, 0.0], [2.0, 0.0], [0.0, 2.0]])
>>> nu_prod = np.array([[0.0, 1.0], [0.0, 2.0], [1.0, 0.0]])
>>> r = np.array([ 40., 20.])
>>> reaction_rate(nu_react, nu_prod, r)
array([-20., -40., 0.])
"""
nu = nu_prod - nu_react
return np.dot(nu, rj)
In [1]:
%%file main.py
import numpy as np
import chemkin
import reaction_coeffs as rrc
M = 3 # Number of reaction
N = 5 # Number of species
concs = np.array([2.0, 1.0, 0.5, 1.0, 1.0])
Temps = [750.0, 1500.0, 2500.0]
A1, b1, E1 = 1.0e+08, 0.5, 5.0e+04
A2, E2 = 1.0e+07, 1.0e+04
nu_react = np.array([[2.0, 0.0, 0.0], [1.0, 0.0, 1.0], [0.0, 1.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]])
nu_prod = np.array([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [2.0, 0.0, 1.0], [0.0, 0.0, 1.0], [0.0, 1.0, 0.0]])
k = np.zeros(M)
for T in Temps:
# Get the reaction rate coefficients
k[0] = rrc.k_mod_arr(A1, b1, E1, T)
k[1] = rrc.k_const(1.0e+04)
k[2] = rrc.k_arr(A2, E2, T)
# Get progress rates for reactions
omega = chemkin.progress_rate(nu_react, concs, k)
# Get reaction rates for species
f = chemkin.reaction_rate(nu_react, nu_prod, omega)
# Print reaction rates to screen
print(f)