CSTRCOM: Isothermal Reactor with Complex Reaction
A lab report on partial fulfillment of the course Modeling and Simulation in
Chemical Engineering (CHEG 305)
Submitted by:
Bishal Lamichhane
Roll number: 16
Third Year/Sixth semester
Department of Chemical Science and Engineering, Kathmandu University
Submitted to:
Dr. Kundan Lal Shrestha
Department of Environmental Science and Engineering, Kathmandu University
Jan 2022
Table of Contents
1. INTRODUCTION ..................................................................................................................................... 1
2. METHOD .................................................................................................................................................... 2
i. Problem definition ........................................................................................................................... 2
ii. Understanding the CSTR model .................................................................................................. 2
iii. Process modeling ......................................................................................................................... 2
iv. Software selection ............................................................................................................................ 2
v. Solution of model equations......................................................................................................... 2
vi. Model validation................................................................................................................................ 2
vii. Simulation study........................................................................................................................... 2
3. OBSERVATIONS...................................................................................................................................... 3
4. RESULT AND DISCUSSION ................................................................................................................. 3
1.1 Concentration against Time ................................................................................................ 3
2.1 Parametric runs for varying time residence ................................................................ 3
3.1 Parametric runs for varying feed concentration ........................................................ 3
4.1 Parametric runs for varying rate constant ................................................................... 3
5. APPENDIX ................................................................................................................................................. 4
6. REFRENCES .............................................................................................................................................. 5
1. INTRODUCTION
Continuous Stirred Tank Reactor is an open reactor with continuous agitation where
the reaction components are entered and exited at the same time (Fig 1). Here, a series
of complex reactions are occurred inside the CSTR in isothermal environment (Ingham
et al., 2007).
𝐴 + 𝐵 ⟶ 𝑋, 𝑟𝑒𝑎𝑐𝑡𝑖𝑜𝑛 𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡 𝑘1
𝐵+𝑋 ⟶ 𝑌, 𝑟𝑒𝑎𝑐𝑡𝑖𝑜𝑛 𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡 𝑘2
𝐵 + 𝑌 ⟶ 𝑍, 𝑟𝑒𝑎𝑐𝑡𝑖𝑜𝑛 𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡 𝑘3
F, CA0, CB0
F, CA, CB, CX, CY, CZ
Fig 1: Stirred tank with series-parallel reaction
1. Compare the steady-state values of A, B, X, Y, Z with calculated values obtained by
solving the steady-state component balance equations for this problem.
2. Study the effect of varying residence time t (by changing feed rate), feed
concentration, and rate constants on reactor performance. This can be done using
“Parametric Runs” to obtain plots of the steady state concentrations as final values
versus the corresponding change in parameter.
The unsteady state mass balance equation for each component is given as
Unsteady Mass Balance equation
𝑅𝑎𝑡𝑒 𝑜𝑓 𝑅𝑎𝑡𝑒 𝑜𝑓 𝑅𝑎𝑡𝑒 𝑜𝑓 𝑅𝑎𝑡𝑒 𝑜𝑓
( )− ( )+( )= ( )
𝑓𝑙𝑜𝑤 𝑖𝑛 𝑓𝑙𝑜𝑤 𝑜𝑢𝑡 𝐺𝑒𝑛𝑒𝑟𝑎𝑡𝑖𝑜𝑛 𝐴𝑐𝑐𝑢𝑚𝑢𝑙𝑎𝑡𝑖𝑜𝑛
𝑑𝑁𝐴
𝐹𝐶𝐴0 − 𝐹𝐶𝐵0 + 𝑟𝐴 𝑉 = (1)
𝑑𝑡
𝑑𝐶𝐴 𝑉
𝐹(𝐶𝐴0 − 𝐶𝐴 ) + 𝑟𝐴 𝑉 = (2)
𝑑𝑡
𝐹 𝑑𝐶𝐴
(𝐶 − 𝐶𝐴 ) + 𝑟𝐴 = (3)
𝑉 𝐴0 𝑑𝑡
1 𝑑𝐶𝐴
(𝐶 − 𝐶𝐴 ) + 𝑟𝐴 = (4)
𝜏 𝐴0 𝑑𝑡
𝐹 1 𝑚3 /𝑠 1
∵ = , =
𝑉 𝜏 𝑚3 𝑠
Rate laws,
𝑟𝐴 = −𝑘1 𝐶𝐴 𝐶𝐵 (5)
𝑟𝐵 = −𝑘1 𝐶𝐴 𝐶𝐵 − 𝑘2 𝐶𝐵 𝐶𝑋 − 𝑘3 𝐶𝐵 𝐶𝑌 (6)
Similarly for X, Y and Z
𝑟𝑋 = 𝑘1 𝐶𝐴 𝐶𝐵 − 𝑘2 𝐶𝐵 𝐶𝑋
𝑟𝑌 = 𝑘2 𝐶𝐵 𝐶𝑋 − 𝑘3 𝐶𝐵 𝐶𝑌
𝑟𝑍 = 𝑘3 𝐶𝐵 𝐶𝑌
Mass Balance for each component,
𝑑𝐶𝐴 1
= (𝐶𝐴0 − 𝐶𝐴 ) − 𝑘1 𝐶𝐴 𝐶𝐵 (7)
𝑑𝑡 𝜏
𝑑𝐶𝐵 1
= (𝐶𝐵0 − 𝐶𝐵 ) − 𝑘1 𝐶𝐴 𝐶𝐵 − 𝑘2 𝐶𝐵 𝐶𝑋 − 𝑘3 𝐶𝐵 𝐶𝑌 (8)
𝑑𝑡 𝜏
𝑑𝐶𝑋 1
= (𝐶𝑋0 − 𝐶𝑋 ) + 𝑘1 𝐶𝐴 𝐶𝐵 − 𝑘2 𝐶𝐵 𝐶𝑋 (9)
𝑑𝑡 𝜏
𝑑𝐶𝑌 1
= (𝐶𝑌0 − 𝐶𝑌 ) + 𝑘2 𝐶𝐵 𝐶𝑋 − 𝑘3 𝐶𝐵 𝐶𝑌 (10)
𝑑𝑡 𝜏
𝑑𝐶𝑍 1
= (𝐶𝑍0 − 𝐶𝑍 ) + 𝑘3 𝐶𝐵 𝐶𝑌 (11)
𝑑𝑡 𝜏
Where,
𝐶𝐴0 , 𝐶𝐵0 , 𝐶𝑋0 , 𝐶𝑌0, 𝐶𝑍0 = Initial Feed concentration (kmol/m3)
𝐶𝐴 , 𝐶𝐵 , 𝐶𝑋 , 𝐶𝑌 , 𝐶𝑍 = Concentration in reactor (kmol/m3)
F = Flow rate (m3 /s)
k = Reaction rate constant (m3 /mol.s)
n = Reaction order
t = Time (s)
V= Volume (m3)
Parametric runs
Residence time, feed concentration, and rate constants are changed and the
corresponding change in steady state value of concentration is observed.
1
2. METHOD
i. Problem definition
The goal of this modeling was to analyze the reaction progression in an isothermal
CSTR.
ii. Understanding the CSTR model
The model was obtained by combining mole balance, rate law, and stoichiometry of a
CSTR. It relates rate of change of concentration of each component with specific rate
constant, concentration and residence time.
iii. Process modeling
Equations (7), (8), (9), (10) and (11) represent the CSTR model. The process of coding
the model in python is as follows:
1. Import all the necessary libraries.
2. Define a function for the Ordinary differential equations (ODEs). It should
return the ODE as an array as well as set the arrays of differential equation in
this function.
3. Set the values of initial moles per unit volume for each component.
4. Solve the ODE and store it in an array.
The differential equations are solved using python. The codes written in python are
explained below.
Libraries and functions used in the program:
numpy To access and manipulate arrays.
scipy For mathematical operation such as
integration and derivative.
linspace To divide the points.
x.T For transpose purpose
dxdt Function to store the components
concentration as an array
Infodict Information dictionary
Several plotting functions were used in the program. Some of them are:
1. plt.plot = for plotting the graph
2. plt.legend = for providing the legend in the graph
3. plt.xlabel = for providing the name of contains in the x-axis
4. plt.ylabel = for providing the name of contains in the y-axis
5. plt.show = displays the graph
iv. Software selection
Proper software should be chosen based on the problem's requirements and the
process model. Python was used as the primary programming language due to its large
collection of modelling and visualization related packages such as numpy,
matplotlib,scipy etc. Furthermore, Jupyter notebook was used as the text editor of
choice due to its simple and easy to visualize the layout.
v. Solution of model equations
Functions odeint and matplotlib.pyplot were used to solve ordinary differential
equations and plot graphs.
vi. Model validation
The nature of graphs showed that reactant concentration was decreasing overtime and
product concentration was increasing overtime. Therefore, the model was validated.
vii. Simulation study
The model was used to generate graphs for time versus concentration. Parametric runs
were also done.
2
3. OBSERVATIONS
• Variation of concentration with time
• Study of varying time residence, feed concentration and rate constant.
4. RESULT AND DISCUSSION
1.1 Concentration against Time
i.Fig 1.1: Concentration against Time for each component
From the figure 1.1 we can observe that at the beginning of the reaction, the
concentrations of A and B decrease while the concentration of Z increases. This
demonstrates that a reaction occurs. After some time, however, all of the
concentrations reach a steady state value. Furthermore, because X and Y are
intermediate products, their concentrations have not altered.
The final concentration for components A, B, X, Y and Z are seen to be 0.28, 0.2, 0, 0 and
0.13 respectively.
Whereas, the exact value obtains by running code 0.264507909, 0.204896844,
0.00259443057, 0.00618425365, 0.126713405 respectively.
2.1 Parametric runs for varying time residence
Fig 2.1: Effect of varying time residence on different components
Here we can observe from figure 2.1 that the concentration A is decreased as tau is
increased. This is because as residence time was increased, the extent of conversion of
reactants also increased and hence the extent of formation of product also increases
simultaneously. It can be interpreted that the conversion of A into the product is
increase as tau increased. Same response can be observed on varying tau with the
concentration of B. Similarly, the concentration of Z seems to be increasing as tau is
increased. Whereas, the concentration of X and Y are unchanged.
3.1 Parametric runs for varying feed concentration
Fig 3.1: Effect of varying feed concentrations
Similarly, as we can see from the figure 3.1 above on increasing the values of initial
molar concentration i.e., moles per unit volume (A0 and B0), the concentration of
formation of product is higher.
4.1 Parametric runs for varying rate constant
Fig 4.1: Effect of varying specific rate constants
As seen in the figure 4.1 above, increasing the specific rate constants has increased the
concentration of product and decreased the concentration of reactants for the same
point of time.
3
5. APPENDIX
#Source code
The initial given conditions are (Oberndorfer, 2001)
A0 = 0.4 #initial concentration of A, Kmols/m^3
B0 = 0.6 #initial concentration of B, Kmols/m^3
X0 = 0 #initial concentration of X, Kmols/m^3
Y0 = 0 #initial concentration of Y, Kmols/m^3
Z0 = 0 #initial concentration of Z, Kmols/m^3
K1 = 5 * 10**(-4) #Rate Constant for reaction 1st, m^3/Kmol.s
K2 = 5 * 10**(-2) #Rate Constant for reaction 2nd, m^3/Kmol.s
K3 = 2 * 10**(-2) #Rate Constant for reaction 3rd, m^3/Kmol.s
tau = 5000 #Residence time, s
Cint = 50
Plotting this value in python and solving for results:
#importing the requires libraries
from numpy import *
import matplotlib.pyplot as plt
from scipy import integrate,optimize
def dxdt(x):
Ca, Cb, Cx, Cy, Cz= x
eq1=((A0-Ca)/tau)-(k1*Ca*Cb)
eq2=((B0-Cb)/tau)-(k1*Ca*Cb)-(k2*Cb*Cx)-
(k3*Cb*Cy) eq3=((X0-Cx)/tau)+(k1*Ca*Cb)-
(k2*Cb*Cx)
eq4=((Y0-Cy)/tau)+(k2*Cb*Cx)-
(k3*Cb*Cy) eq5=((Z0-
Cz)/tau)+(k3*Cb*Cy)
return [eq1, eq2, eq3, eq4, eq5]
Ca, Cb, Cx, Cy, Cz = optimize.fsolve(dxdt, (0.35,0.2,0.15,0,0))
print(Ca, Cb, Cx, Cy, Cz)
In [3]:
Ca, Cb, Cx, Cy, Cz = optimize.fsolve(dxdt, (0.35,0.2,0.15,0,0))print(Ca, Cb, Cx, Cy,
Cz)
def dxdt(x,t):
return array([(A0-x[0])/tau-k1*x[0]*x[1],
(B0-x[1])/tau-k1*x[0]*x[1]-k2*x[1]*x[2]-k3*x[1]*x[3],
(X0-x[2])/tau+k1*x[0]*x[1]-k2*x[1]*x[2],
(Y0-x[3])/tau+k2*x[1]*x[2]-k3*x[1]*x[3], (Z0-
x[4])/tau+k3*x[1]*x[3]
])
In [4]:
t=linspace(0,Tfin,cint)
x, infodict = integrate.odeint(dxdt, x0, t, full_output=True)
print(infodict['message'])
ca,cb,cx,cy,cz=x.T
plt.plot(t, ca, 'r', label ='ca')
plt.plot(t, cb, 'b', label ='cb')
plt.plot(t, cx, 'g', label ='cx')
plt.plot(t, cy, 'c', label ='cy')
plt.plot(t, cz, 'y', label ='cz')
plt.legend(['CA','CB','CX','CY','CZ'], loc='best')plt.xlabel("Time")
plt.ylabel("Concentration(Kmol/m^3)") plt.title("Concentration
against Time")
Integration successful.
Out[4]: Text(0.5, 1.0, 'Concentration against Time')
In [5]:
for tau in range(500,10500,2000):
x, infodict = integrate.odeint(dxdt, x0, t, full_output=True)ca,cb,cx,cy,cz=x.T
plt.plot(t,ca,'r',label='CA, tau={}s'.format(tau))
plt.plot(t,cb,'b',label='CB, tau={}s'.format(tau))
plt.plot(t,cx,'g',label='CX, tau={}s'.format(tau))
plt.plot(t,cy,'c',label='CY, tau={}s'.format(tau))
plt.plot(t,cz,'y',label='CZ, tau={}s'.format(tau))plt.title("Varying
time residence")
plt.xlabel("Time(sec)")
plt.ylabel("Concentration of components(Kmol/m^3)")
plt.legend(bbox_to_anchor=(-0.15,1))
Out[5]: <matplotlib.legend.Legend at 0x24f060810a0>
In [6]:
def range_with_floats(a, b, x):
while b > a:
yield a a
=a+x
for A0 in range_with_floats(0.2, 1.0, 0.5):
x0=array([A0,B0,X0,Y0,Z0])
x, infodict = integrate.odeint(dxdt, x0, t, full_output=True)ca,cb,cx,cy,cz=x.T
plt.plot(t,ca,'r',label='CA, A0={}'.format(A0))
plt.plot(t,cb,'b',label='CB, A0={}'.format(A0))
plt.plot(t,cx,'g',label='CX, A0={}'.format(A0))
plt.plot(t,cy,'c',label='CY, A0={}'.format(A0))
plt.plot(t,cz,'y',label='CZ, A0={}'.format(A0))
plt.legend(bbox_to_anchor=(-0.15,1))
plt.xlabel("Time(sec)")
plt.ylabel("Concentration(Kmol/m^3)")
plt.title("Varying Feed Concentration A0")
Out[6]: Text(0.5, 1.0, 'Varying Feed Concentration A0')
In [7]:
def range_with_floats(a, b, x):
while b > a:
yield a a
=a+x
for B0 in range_with_floats(0.2, 1.0, 0.5):
x0=array([A0,B0,X0,Y0,Z0])
x, infodict = integrate.odeint(dxdt, x0, t, full_output=True)ca,cb,cx,cy,cz=x.T
plt.plot(t,ca,'r',label='CA, B0={}'.format(B0))
plt.plot(t,cb,'b',label='CB, B0={}'.format(B0))
plt.plot(t,cx,'g',label='CX, B0={}'.format(B0))
plt.plot(t,cy,'c',label='CY, B0={}'.format(B0))
plt.plot(t,cz,'y',label='CZ, B0={}'.format(B0))
plt.legend(bbox_to_anchor=(-0.15,1))
plt.xlabel("Time(sec)")
plt.ylabel("Concentration(Kmol/m^3)")
plt.title("Varying Feed Concentration B0")
Out[7]: Text(0.5, 1.0, 'Varying Feed Concentration B0')
In [8]:
for k1 in range_with_floats(0.0003,0.0006,0.00015):
x, infodict = integrate.odeint(dxdt, x0, t, full_output=True)ca,cb,cx,cy,cz=x.T
plt.plot(t,ca,'r',label='CA, k1={}'.format(k1))
plt.plot(t,cb,'b',label='CB, k1={}'.format(k1))
plt.plot(t,cx,'g',label='CX, k1={}'.format(k1))
plt.plot(t,cy,'c',label='CY, k1={}'.format(k1))
plt.plot(t,cz,'y',label='CZ, k1={}'.format(k1))
plt.legend(bbox_to_anchor=(-0.15,1))
plt.xlabel("Time")
plt.ylabel("Concentration(Kmol/m^3)")
plt.title("Change in Rate constant K1")
4
Out[8]: Text(0.5, 1.0, 'Change in Rate constant K1'
for k2 in range_with_floats(0.0003,0.0006,0.00015):
Text(0.5, 1.0, 'Change in Rate constant K2')
Out[9]:
for k3x,ininfodict = integrate.odeint(dxdt, x0, t, full_output=True)ca,cb,cx,cy,cz=x.T
range_with_floats(0.0003,0.0006,0.00015):
plt.plot(t,ca,'r',label='CA, k2={}'.format(k2))
x, infodict = integrate.odeint(dxdt, x0, t, full_output=True)ca,cb,cx,cy,cz=x.T
Inplt.plot(t,cb,'b',label='CB,
[10]: k2={}'.format(k2))
plt.plot(t,ca,'r',label='CA, k3={}'.format(k3))
plt.plot(t,cx,'g',label='CX, k2={}'.format(k2))
plt.plot(t,cb,'b',label='CB, k3={}'.format(k3))
plt.plot(t,cy,'c',label='CY, k2={}'.format(k2))
plt.plot(t,cx,'g',label='CX, k3={}'.format(k3))
plt.plot(t,cz,'y',label='CZ,
plt.plot(t,cy,'c',label='CY, k2={}'.format(k2))
k3={}'.format(k3))
plt.legend(bbox_to_anchor=(-0.15,1))
plt.plot(t,cz,'y',label='CZ, k3={}'.format(k3))
plt.xlabel("Time")
plt.legend(bbox_to_anchor=(-0.15,1))
plt.ylabel("Concentration(Kmol/m^3)")
plt.xlabel("Time")
plt.title("Change in Rate constant K2")
plt.ylabel("Concentration(Kmol/m^3)")
plt.title("Change in Rate constant K3")
Text(0.5, 1.0, 'Change in Rate constant K3')
[10]:
6. REFRENCES
• Ingham, J., Dunn, I. J., Heinzle, E., Přenosil, J. E., & Snape, J. B. (2007). Chemical
Engineering Dynamics: An Introduction to Modelling and Computer Simulation:
Third Edition. In Chemical Engineering Dynamics: An Introduction to Modelling
and Computer Simulation: Third Edition.
https://doi.org/10.1002/9783527614219
• Oberndorfer, C. (2001). Chemical engineering dynamics. In Zeitschrift fuer
Metallkunde/Materials Research and Advanced Techniques (Vol. 92, Issue 2).
https://doi.org/10.1002/9783527616015
