applying our knowledge to solve Chem Reactor problems
We now know ODEs
It is time everything together with other cool python modules to solve engineering problems.
If you don't know how to solve these problems it's ok, start small and make your way up.
The easiest example we could do would be to analyze the chemical reaction kinetics in the simplest case:
A -> B
We assume it is a first order reaction, and we use the design equation from the Batch reactor:
So that means that $-r_A = k C_A$.
Let's assume $k = 10^{-3} 1/s$ and $C_A(t=0) = 10$ moles
We can expect the half-time of A to be around 10000 seconds, which is ~3 hours.
Ok, let's see that through Python.
We want to track the concentration, which we will call $y_A$.
$dy_A/dt = - k C_A = -k y_A$
An important tool we need to use here is the ODEINT module:
from scipy.integrate import odeint
In this module if you have an equation $\frac{dy}{dt} = f(y,t)$ you need:
Some additional commands that you may need in your script:
#lets import the modules # plotting everything inline
# %matplotlib inline
import matplotlib.pyplot as plt # plotting modules
plt.style.use('presentation') # just have in your script for prettier plotting
import numpy as np # our matlab-like module
from scipy.integrate import odeint # integration of ODEs so you don't have to write your finite difference yourself
To plot your data use these commands:
plt.plot(times, y_A_exact,'k--', label='analytical')
plt.xlabel('$t [s]$')
plt.ylabel('$C_A [moles]$')
plt.legend()
plt.savefig('odeint.pdf')
plt.show()
#let's import the modules # plotting everything inline
%matplotlib inline
import matplotlib.pyplot as plt # plotting modules
#plt.style.use('presentation') # just have in your script for prettier plotting
import numpy as np # our matlab-like module
from scipy.integrate import odeint # integration of ODEs so you don't have to write your finite difference yourself
#reaction order
def f(y, t):
return -k*y
k = 0.001 # L/(mol*sec)
time_start = 0 #s
time_finish = 10000 #s
N_points = 50
times = np.linspace(time_start,time_finish,N_points)
y_A0 = 10. # moles
y_A_calc = odeint(f, y_A0, times)
plt.plot(times, y_A_calc, 'ro-', label='integrated')
y_A_exact = y_A0*np.exp(-k*times)
plt.plot(times, y_A_exact,'k--', label='analytical')
plt.xlabel('$t [s]$')
plt.ylabel('$C_A [moles]$')
plt.legend()
plt.savefig('odeint.pdf')
plt.show()
Now let's do the same calculation as case number 1 (of this problem), but this time we will use the conversion factor $X_A$ as our variable:
$X_A = \dfrac{N_{A0} - N_A}{N_{A0}} = \dfrac{C_{A0} - C_A}{C_{A0}}$ if the volume $V=const$
So the design equation looks like:
$\dfrac{d X_A}{dt} = (-r_A/N_{A0}) V $
%reset
%matplotlib inline
import matplotlib.pyplot as plt # plotting modules
plt.style.use('ggplot') # just have in your script for prettier plotting
import numpy as np # our matlab-like module
from scipy.integrate import odeint # integration of ODEs so you don't have to write your finite difference yourself
plt.clf()
k=0.001
y_A0 = 10.
time_start = 0 #s
time_finish = 10000 #s
N_points = 50
times = np.linspace(time_start,time_finish,N_points)
y_A_exact = y_A0*np.exp(-k*times)
x_A = ((y_A0 - y_A_exact)/y_A0)
# print(y_A_exact)
# print(x_A)
plt.plot(times, x_A, 'ro-')
plt.xlabel('$t [s]$')
plt.ylabel('$X_A$')
# plt.savefig('conversion factor.pdf')
plt.show()
2A ->B
We assume that it is a first order reaction and we use the design equation from the Batch reactor:
So that means that $-r_A = k [C_A]^2$.
Let's assume $k = 10^{-3} 1/s$ and $C_A(t=0) = 10$ moles
So we can expect the half-time of A to be around 10000 seconds which is ~3 hours.
Ok let's see that through Python.
We want to track the concentration which we will call $y_A$.
$dy_A/dt = - k C_A = -k y^2_A$
The exact solution here would be
$ 1/y_A = kt + C$ $y_A(0) = 10$ => $C = 1/10$
$1/y_A = kt + 1/y_A_0$
So now it is easier to track $1/y_A$ rather than the function itself
# %matplotlib inline
%reset
import matplotlib.pyplot as plt # plotting modules
#plt.style.use('presentation') # just have in your script for prettier plotting
import numpy as np # our matlab-like module
from scipy.integrate import odeint # integration of ODEs so you don't have to write your finite difference yourself
#reaction order
def f(y, t):
return -k*y**2
k = 0.001 # 1/sec
time_start = 0 #s
time_finish = 10 #s
N_points = 50
times = np.linspace(time_start,time_finish,N_points)
y_A0 = 10. # moles
y_A_calc = odeint(f, y_A0, times)
plt.plot(times, 1./y_A_calc, 'ro-', label='integrated')
inverse_y_A_exact = (1/y_A0 + k*times)
plt.plot(times, inverse_y_A_exact,'k--', label='analytical')
plt.xlabel('$t [s]$')
plt.ylabel('$1/C_A [moles]$')
plt.legend()
# plt.savefig('odeint.pdf')
plt.show()
Having had solved ODEs, now we are going to solve "Simple" algebraic equations.
Given a continuously stirred tank reactor with a volume of $66 m^3$ where the reaction $A \to B$ occurs, at a rate of $−r_A= k C^2_A$ ($k=3 L/mol/h$), with an entering molar flow of $F_{A0} = 5 mol/h$ and a volumetric flowrate of $\upsilon = 10 L/h$, what is the exit concentration of A?
From a mole balance we know that at steady state $0=F_{A0}−F_A+ V \cdot r_A$. That equation simply states the sum of the molar flow of A in in minus the molar flow of A out plus the molar rate A is generated is equal to zero at steady state. This is directly the equation we need to solve. We need the following relationship:
#Copyright (C) 2013 by John Kitchin
from scipy.optimize import fsolve
Fa0 = 5.0
v0 = 10.
V = 66000.0 # reactor volume L^3
k = 3.0 # rate constant L/mol/h
def func(Ca):
"Mole balance for a CSTR. Solve this equation for func(Ca)=0"
Fa = v0 * Ca # exit molar flow of A
ra = -k * Ca**2 # rate of reaction of A L/mol/h
return Fa0 - Fa + V * ra
# CA guess that 90 % is reacted away
CA_guess = 0.1 * Fa0 / v0
CA_sol, = fsolve(func, CA_guess)
print('The exit concentration is {0} mol/L'.format(CA_sol))
# Now let's solve the same problem using Python's
# built-in functionality through a Python package odeint
# Step 1:
# Import modules into Python
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
plt.style.use('fivethirtyeight')
# integration of ODEs so you don't have
# to write your finite difference yourself
from scipy.integrate import odeint
# Step 2:
# define the required parameters and functions
def f(y, t, k):
dydt = -k*t
return dydt
# define initial values, arrays, intervals
t_0 = 0.
t_f = 10.
npoints = 10
t_array = np.linspace(t_0, t_f, npoints + 1)
y_0 = 1.
k = 2.
h = t_array[1] - t_array[0]
y_num = np.zeros_like(t_array)
# let's also define a numerical solution
t_array_analyt = np.linspace(t_0, t_f, 1000 + 1)
y_analyt = -k*t_array_analyt**2/2. + y_0
# here the syntax is
# module odeint(right hand side function,
# initial value of y,
# array for the time values
# additional parameters like k in a tuple (k,))
# if we needed two parameters like k1, k2
# then args=(k1, k2)
y_num = odeint(f, y_0, t_array, args=(k,))
plt.xlabel('$t, s$')
plt.ylabel('$y$')
plt.plot(t_array_analyt, y_analyt, 'k--', lw=3,label='analytical')
plt.plot(t_array, y_num, 'go',alpha=0.6, label='numerical')
plt.legend()
Let's compare a set number of intervals vs. 1001 intervals from the analytical, at a specific time. For the sake of simplicity, I chose the last point. And because typically, error increases as we increase the number of steps, since every step is dependant on the last.
When we have 5 intervals the difference in $y$ is ~4.13922407e-08
When we have 10 intervals the difference in $y$ is ~2.98024077e-08
When we have 20 intervals the difference in $y$ is ~7.4505806e-09
When we have 30 intervals the difference in $y$ is ~3.31137073e-09
When we have 100 intervals the difference in $y$ is ~2.98030045e-10
In these cases, we have an overestimate, which gets more accurate as we increase the number of intervals.
Although it is harder to see visually, most approximations are fairly accurate.
#https://en.wikipedia.org/wiki/Euler_method
#https://rosettacode.org/wiki/Euler_method
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
plt.style.use('fivethirtyeight')
from scipy.integrate import odeint
def f(y, t, k):
dydt = -k*t
return dydt
t_0 = 0.
t_f = 10
npoints0 = 5
t_array0 = np.linspace(t_0, t_f, npoints0 + 1)
npoints1 = 10
t_array1 = np.linspace(t_0, t_f, npoints1 + 1)
npoints2 = 20
t_array2 = np.linspace(t_0, t_f, npoints2 + 1)
npoints3 = 30
t_array3 = np.linspace(t_0, t_f, npoints3 + 1)
npoints4 = 100
t_array4 = np.linspace(t_0, t_f, npoints4 + 1)
y_0 = 1.
k = 2.
h = t_array1[1] - t_array1[0]
t_array_analyt = np.linspace(t_0, t_f, 1000 + 1)
y_analyt = -k*t_array_analyt**2/2. + y_0
y_num0 = np.zeros_like(t_array0)
y_num0 = odeint(f, y_0, t_array0, args=(k,))
y_num1 = np.zeros_like(t_array1)
y_num1 = odeint(f, y_0, t_array1, args=(k,))
y_num2 = np.zeros_like(t_array2)
y_num2 = odeint(f, y_0, t_array2, args=(k,))
y_num3 = np.zeros_like(t_array3)
y_num3 = odeint(f, y_0, t_array3, args=(k,))
y_num4 = np.zeros_like(t_array4)
y_num4 = odeint(f, y_0, t_array4, args=(k,))
plt.xlabel('$t, s$')
plt.ylabel('$y$')
plt.plot(t_array_analyt, y_analyt, 'k--', lw=3,label='analytical')
plt.plot(t_array0, y_num0, 'go',alpha=0.4, label='numerical0')
plt.plot(t_array1, y_num1, 'bo',alpha=0.4, label='numerical1')
plt.plot(t_array2, y_num2, 'co',alpha=0.4, label='numerical2')
plt.plot(t_array3, y_num3, 'mo',alpha=0.4, label='numerical3')
plt.plot(t_array4, y_num4, 'yo',alpha=0.4, label='numerical3')
plt.legend()
print(y_num0[5] - y_analyt[1000])
print(y_num1[10] - y_analyt[1000])
print(y_num2[20] - y_analyt[1000])
print(y_num3[30] - y_analyt[1000])
print(y_num4[100] - y_analyt[1000])
This problem should be familiar, but we will throw in a twist.
When $A -> B$ as simple reaction with the rate $-r_A = k C_A$ and the solution for A could be found in this way: $\dfrac{d C_A}{dt} = -k C_A(t)$ $C_{A0} = 10$ mole/L
import matplotlib.pyplot as plt # plotting modules
%matplotlib inline
plt.style.use('ggplot') # just have in your script for prettier plotting
# instead of 'ggplot' you can use 'presentation'
import numpy as np # our matlab-like module
from scipy.integrate import odeint # integration of ODEs so you don't have to write your finite difference yourself
#step 1 definition of the function:
def f(C_A, t):
dcadt = - k*C_A
return dcadt
k = 1
C_A0 = 10.
time_start = 0.
time_finish = 10000.
N_points = 100
time_array = np.linspace(0, 10., N_points)
C_num = odeint(f,C_A0, time_array)
plt.plot(time_array, C_num, 'r-')
plt.plot(time_array, C_num, 'go', alpha=0.8, fillstyle='none') # alpha - 80% transparency
# fillstyle - none - means no filling of the circles
plt.xlabel('t, s')
plt.ylabel('$C_A(t)$, mole/L')
plt.show()
Now the same thing for B $ -(r_B) = -(-r_A) = - k C_A$ so the design equation will have an additional equation:
$\dfrac{d C_B}{dt} = k C_A(t)$ $C_{B0} = 0$ mol/L
import matplotlib.pyplot as plt # plotting modules
%matplotlib inline
plt.style.use('ggplot') # just have in your script for prettier plotting
# instead of 'ggplot' you can use 'presentation'
import numpy as np # our matlab-like module
from scipy.integrate import odeint # integration of ODEs so you don't have to write your finite difference yourself
#step 1 definition of the function:
def f(C_list, t):
C_A = C_list[0]
C_B = C_list[1]
dcadt = - k*C_A
dcbdt = k*C_A
return [dcadt, dcbdt]
k = 2
C_A0 = 10.
C_B0 = 0.
time_start = 0.
time_finish = 10.
N_points = 100
time_array = np.linspace(0, 10., N_points)
C0_list = [1, 0.]
C_num_list = odeint(f, C0_list, time_array)
print(C_num_list.shape)
C_A_num = C_num_list[:,0]
C_B_num = C_num_list[:,1]
plt.plot(time_array, C_A_num, 'r--', fillstyle='none', label='$A$')
plt.plot(time_array, C_B_num, 'g--', fillstyle='none', label='$B$')
plt.xlabel('t, s')
plt.ylabel('$C(t)$, mole/L')
plt.legend()
plt.show()
To see precisely what the values were from the integration, we can use
print(C_num_list[:,0])
C_num_list[:,0]
This one is more complicated:
A system of ODES
$ A + B \to C$, $k_1 = 1$ L/(mol s)
$ B + C \to D$, $k_1 = 1.5$ L/(mol s)
with initial concentrations $C_{A0} = 1$ mol/L $C_{B0} = 1$ mol/L $C_{C0} = 0$ mol/L $C_{D0} = 0$ mol/L
r1 = k1C_AC_B r2 = k2C_BC_C
dcAdt = -r1 dcBdt = -r1 - r2 dcCdt = r1 - r2 dcDdt = r2
import matplotlib.pyplot as plt # plotting modules
%matplotlib inline
plt.style.use('ggplot') # just have in your script for prettier plotting
# instead of 'ggplot' you can use 'presentation'
import numpy as np # our matlab-like module
from scipy.integrate import odeint # integration of ODEs so you don't have to write your finite difference yourself
#step 1 definition of the function:
def f(C_list, t):
k1 = 1.
k2 = 1.5
C_A = C_list[0]
C_B = C_list[1]
C_C = C_list[2]
r1 = k1*C_A*C_B
r2 = k2*C_B*C_C
dcAdt = -r1
dcBdt = -r1 - r2
dcCdt = r1 - r2
dcDdt = r2
return [dcAdt,dcBdt,dcCdt,dcDdt]
k1 = 1.
k2 = 1.5
C0_list = [1.0, 1.0, 0.0, 0.0]
time_start = 0.
time_finish = 3.
N_points = 100
time_array = np.linspace(0, 3., N_points)
C_num_list = odeint(f, C0_list, time_array)
print(C_num_list.shape)
C_A_num = C_num_list[:,0]
C_B_num = C_num_list[:,1]
C_C_num = C_num_list[:,2]
C_D_num = C_num_list[:,3]
plt.plot(time_array, C_A_num, 'r-', label='$A$')
plt.plot(time_array, C_B_num, 'g-', label='$B$')
plt.plot(time_array, C_C_num, 'm--', label='$C$')
plt.plot(time_array, C_D_num, 'k--', label='$D$')
plt.xlabel('t, s')
plt.ylabel('$C(t)$, mole/L')
plt.legend()
plt.show()
Normal butane, C4H10, is to be isomerized to isobutane in a plug-flow reactor, because isobutane is worth more. This elementary reversible reaction is to be carried out adiabatically in the liquid phase under high pressure using trace amounts of a liquid catalyst which gives a specific reaction rate of 31.1/hr at 360K. The feed enters at 330K.
Additional information:
a) Calculate the PFR volume necessary to process 100,000 gal/day (163 kmol/h) at 70% conversion of a mixture of 90 mol% n-butane and 10 mol% i-pentane, which is considered an inert.
b) Plot and analyze $X, X_e, T$, and $-r_A$ down the length of the reactor.
You may use this code in your program. You just need to understand and plot it.
DHrx = -6900. # J/mol
Ea = 65700. # J/mol
R = 8.314 # J/mol/K
T0 = 330 # K
FA0 = 0.9 * 163e3 # mol/hr
CA0 = 9.3e3 # mol/m3
def Kc(T):
"Equilibrium constant, as a function of temperature"
return 3.03 * np.exp(( -1* DHrx / R) * (1./T - 1./(273.15+60)) )
assert Kc(273.15+60) == 3.03
assert Kc(273.15+65) < 3.03 # le Chatellier's principle
def k(T):
"Specific reaction rate in 1 hour, as a function of temperature"
return 31.1 * np.exp((-1*Ea/R) * (1./T - 1./360.))
assert k(360) == 31.1
assert k(365) > 31.1
def Teb(X):
"Temperature from energy balance, as a function of conversion."
Cps = np.array([141, 141, 161])
Thetas = np.array([1., 0., 0.1/0.9])
return T0 - DHrx * X / sum(Cps*Thetas)
assert Teb(0) == T0
assert Teb(0.5) > T0 # Exothermic
def dXdV(X,V):
"""
dX/dV in a plug flow reactor is -rA/ FA0
"""
T = Teb(X)
CA = CA0*(1.-X)
CB = CA0*X
rA = -k(T)*(CA - CB/Kc(T))
return -rA/FA0
Xeq = Kc(Teb(X))/ (1+Kc(Teb(X)))
rate = k(Teb(X))*(CA0*(1.-X) - CA0*X/Kc(Teb(X)))
#For Levenspiel plot you may use this:
plt.plot(X, FA0/rate)
plt.ylim(0,10)
plt.title("Levenspiel Plot")
plt.xlabel("X")
plt.ylabel("$\\frac{F_{A0}}{-r_A}$")
plt.plot((Xeq[-1],Xeq[-1]),(0,FA0/rate[-1]), ':')
plt.show()
import matplotlib.pyplot as plt # plotting modules
%matplotlib inline
plt.style.use('ggplot') # just have in your script for prettier plotting
# instead of 'ggplot' you can use 'presentation'
import numpy as np # our matlab-like module
import scipy
from scipy.integrate import odeint
DHrx = -6900. # J/mol
Ea = 65700. # J/mol
R = 8.314 # J/mol/K
T0 = 330 # K
FA0 = 0.9 * 163e3 # mol/hr
CA0 = 9.3e3 # mol/m3
def Kc(T):
"Equilibrium constant, as a function of temperature"
return 3.03 * np.exp(( -1* DHrx / R) * (1./T - 1./(273.15+60)) )
assert Kc(273.15+60) == 3.03
assert Kc(273.15+65) < 3.03 # le Chatellier's principle
def k(T):
"Specific reaction rate in 1/hour, as a function of temperature"
return 31.1 * np.exp((-1*Ea/R) * (1./T - 1./360.))
assert k(360) == 31.1
assert k(365) > 31.1
def Teb(X):
"Temperature from energy balance, as a function of conversion."
Cps = np.array([141, 141, 161])
Thetas = np.array([1., 0., 0.1/0.9])
return T0 - DHrx * X / sum(Cps*Thetas)
assert Teb(0) == T0
assert Teb(0.5) > T0 # Exothermic
def dXdV(X,V):
"""
dX/dV in a plug flow reactor is -rA/ FA0
"""
T = Teb(X)
CA = CA0*(1.-X)
CB = CA0*X
rA = -k(T)*(CA - CB/Kc(T))
return -rA/FA0
volumes = np.linspace(0,5,501)
X = scipy.integrate.odeint(dXdV, 0., volumes)
plt.plot(volumes,X,label="Conversion")
Xeq = Kc(Teb(X))/ (1+Kc(Teb(X)))
plt.plot(volumes, Xeq,label="Equilibrium conversion")
plt.legend(loc="best")
plt.xlabel("Volume")
plt.ylabel("Conversion X")
plt.show()
rate = k(Teb(X))*(CA0*(1.-X) - CA0*X/Kc(Teb(X)))
plt.plot(volumes, rate, label="Rate")
plt.xlabel("Volume")
plt.ylabel("Rate")
plt.show()
for volume, conversion in zip(volumes, X):
if conversion>0.7:
break
print(volume)
print(conversion)
plt.plot(X, FA0/rate)
plt.ylim(0,10)
plt.title("Levenspiel Plot")
plt.xlabel("X")
plt.ylabel("$\\frac{F_{A0}}{-r_A}$")
plt.plot((Xeq[-1],Xeq[-1]),(0,FA0/rate[-1]), ':')
plt.show()