Heterogeneous Stirred Tank Reactor (with Python!)

Author: CAChemE.org ~ License: CC BY 4.0

In [1]:

%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]

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

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),