Author: CAChemE.org ~ License: CC BY 4.0
In [1]:
# Loading our favorite python libraries
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from __future__ import division
plt.style.use('fivethirtyeight')
In [2]:
from IPython.html.widgets import interact
from IPython.display import clear_output, display, HTML
Ok we are ready to go!
In [3]:
def STreactor(C_A0 = 10.2, C_B0 = 0, C_C0 = 0, cat_load = 0.2):
''' Batch STR reactor model
initial conditions (concentrations of A, B and C)
should be given in [mol*L-1]
catalyst load is given in %
'''
C_A0 = C_A0 # mol*L-1
C_B0 = C_B0 # mol*L-1
C_C0 = C_C0 # mol*L-1
C_D0 = 0 # mol*L-1
C_H0 = 19.2 # mol*L-1
mtot = 200 # g of A in a semibatch reactor
# Molecular weights
MW = [84.118, 86.132, 88.148, 160, 1.008] # g mol-1
dens = 0.86 # g.mL-1 average denstiy of A
# Note: Assumes similar densities of all the compounds
# Reaction volume
VL = mtot/MW[0]/C_A0 # in L
To = 348 # K - Initial temperature (75ºC)
# Catalyst
CatalystMass = 0.3829 # in [g]
CatAvalaible = cat_load/100 # Percent
CatMW = 106.42
nCat = (CatalystMass*CatAvalaible)/CatMW
# System reaction paramiters
# Matrix with stoichiometry coefficients
# i reactions, j components (RxS)
# j components in this order: 1 2 3 4 5
# A B C D H
alfa = np.array([
[-1, 1, 0, 0, -1],
[0, -1, 1, 0, -1],
[-1, 0, 1, 0, -2],
[-2, 0, 0, 1, -2],
], dtype=float)
Nreactions, Ncomponents = alfa.shape
def f(y,t):
''' y is a vector with the concentration
t is the indpendt variable (time)
'''
C = y[:] # mol*L-1
# If you don't understand this
# you can uncomment the following line and
# run this part of the code line by line
#import pdb; pdb.set_trace()
# k=ko*exp(-Ea/(R*T));
# Kinetic parameters
kp1 = 99627.28 # L.molCat-1.s-1
kp2 = 1884.84 # L.molCat-1.s-1
kp3 = 915.49 # L.molCat-1.s-1
kp4 = 479.29 # L.molCat-1.s-1
K_A = 17.50 # L.mol-1
K_B = 1.97 # L.mol-1
K_H = 0.04 # non-dimensional
Kest = 9.13 # L.mol-1
#(1xR) en mol.molCat-1.s-1
r1 = kp1*C[0] / (1+K_A*C[0]+np.sqrt(K_H)+K_B*C[1]+Kest*(C[2]+C[3]))**2
r2 = kp2*C[1] / (1+K_A*C[0]+np.sqrt(K_H)+K_B*C[1]+Kest*(C[2]+C[3]))**2
r3 = kp3*C[0] / (1+K_A*C[0]+np.sqrt(K_H)+K_B*C[1]+Kest*(C[2]+C[3]))**2
r4 = kp4*C[0] / (1+K_A*C[0]+np.sqrt(K_H)+K_B*C[1]+Kest*(C[2]+C[3]))**2
# D.Eqs
dCAdt = nCat/VL*(-r1-r3-2*r4)
dCBdt = nCat/VL*(r1-r2)
dCCdt = nCat/VL*(r2+r3)
dCDdt = nCat/VL*r4
dCHdt = 0
return [dCAdt, dCBdt, dCCdt, dCDdt, dCHdt]
Ci0 = [C_A0, C_B0, C_C0, C_D0, C_H0] # mol*L-1
t = np.linspace(0, 400.*60, 100000) # time grid
# solving the DEs
soln = odeint(f, Ci0, t)
C_A0 = soln[:, 0] # mol*L-1
C_B0 = soln[:, 1] # mol*L-1
C_C0 = soln[:, 2] # mol*L-1
C_D0 = soln[:, 3] # mol*L-1
C_H20 = soln[:, 4] # mol*L-1
t = t/60 # conversion to min
# Plotting the results
fig, ax1 = plt.subplots(figsize=(8,6))
ax1.plot(t, C_A0, t, C_B0, t, C_C0)
ax1.set_xlabel('time (min)')
ax1.set_ylabel('$C_{i}$ (mol/L)')
ax1.legend(['$C_A$','$C_B$','$C_C$'])
plt.title('Heterogeneous STR: A -> B -> C')
plt.suptitle('(CAChemE.org)')
ax1.set_xlim([0,400])
ax1.set_ylim([0,12])
#ax2 = ax1.twinx()
#ax2.plot(t, C_D0, color='y', ylim=(0,0.20) )
#ax2.set_ylabel('Concentration [mol/L] of D')
plt.show()
In [4]:
interact(STreactor,
C_A0=(0.1,20,0.1), C_B0=(0.1,20,0.1),
C_C0=(0.1,20,0.1), cat_load=(0.01,2.01,0.01) )
In [4]: