ChEn-3170: Computational Methods in Chemical Engineering Fall 2018 UMass Lowell; Prof. V. F. de Almeida 05Nov2018
$ \newcommand{\Amtrx}{\boldsymbol{\mathsf{A}}} \newcommand{\Bmtrx}{\boldsymbol{\mathsf{B}}} \newcommand{\Mmtrx}{\boldsymbol{\mathsf{M}}} \newcommand{\Imtrx}{\boldsymbol{\mathsf{I}}} \newcommand{\Pmtrx}{\boldsymbol{\mathsf{P}}} \newcommand{\Lmtrx}{\boldsymbol{\mathsf{L}}} \newcommand{\Umtrx}{\boldsymbol{\mathsf{U}}} \newcommand{\Smtrx}{\boldsymbol{\mathsf{S}}} \newcommand{\xvec}{\boldsymbol{\mathsf{x}}} \newcommand{\uvar}{\boldsymbol{u}} \newcommand{\fvar}{\boldsymbol{f}} \newcommand{\avec}{\boldsymbol{\mathsf{a}}} \newcommand{\bvec}{\boldsymbol{\mathsf{b}}} \newcommand{\cvec}{\boldsymbol{\mathsf{c}}} \newcommand{\rvec}{\boldsymbol{\mathsf{r}}} \newcommand{\mvec}{\boldsymbol{\mathsf{m}}} \newcommand{\gvec}{\boldsymbol{\mathsf{g}}} \newcommand{\zerovec}{\boldsymbol{\mathsf{0}}} \newcommand{\norm}[1]{\bigl\lVert{#1}\bigr\rVert} \newcommand{\transpose}[1]{{#1}^\top} \DeclareMathOperator{\rank}{rank} $
This problem is extracted from Example 13.4 in the book: Essentials of Chemical Reaction Engineering, Fogler, 2$^\text{nd}$ edition, 2018.
The series reactions 2A->B->3C are catalysed by H2SO4. All reactions are first order in the reactant concentration. However the first reaction is exothermic and the second, endothermic. The reaction is to be carried out in a semibatch reaction with heat exchange. Pure A enters the the batch reactor at a constant volumetric flow rate and temperature. Initially, there is a heel of liquid in the reactor which contains A and catalyst. The reaction rates are independent of the catalyst concentration. The initial temperature of the reactor is lower than the inflow stream. Plot the species concentrations and reactor temperature inside the reactor as a function of time.
Next, an interactive simulation shows the continuous effect of a change in coolant temperature and inflow molarity of the A species.
Name | Parameter | Value | Unit | |
---|---|---|---|---|
volumetric flow rate | $q$ | 240 | L/h | |
inflow molarity of A | $c_\text{Ain}$ | 4 | mol/L | |
inflow molarity of A | $c_\text{Bin}$ | 0 | mol/L | |
inflow molarity of A | $c_\text{Cin}$ | 0 | mol/L | |
inflow temperature | $T_\text{in}$ | 305 | K | |
intial volume | $V_0$ | 100 | L | |
mass density | $\rho$ | 1000 | g/L | |
molar heat capacity of A | $C_{pA}$ | 30 | cal/(mol K) | |
molar heat capacity of B | $C_{pB}$ | 60 | cal/(mol K) | |
molar heat capacity of C | $C_{pC}$ | 20 | cal/(mol K) | |
heat capacity of H2SO4 | $C_{pS}$ | 35 | cal/(mol K) | |
heat of reaction (rxn 1) | $\Delta H_{\text{R},1}$ | -6500 | cal/mol | |
heat of reaction (rxn 2) | $\Delta H_{\text{R},2}$ | 8000 | cal/mol | |
Arrhenius frequency (rxn 1) | $k_{0,1}$ | $1.25$ @ 320 K | $\text{h}^{-1}$ | |
Arrhenius actv. energy (rxn 1) | $E_{\text{a},1}$ | 9500 | cal/mol | |
Arrhenius frequency (rxn 2) | $k_{0,2}$ | $0.08$ @ 300 K | $\text{h}^{-1}$ | |
Arrhenius actv. energy (rxn 2) | $E_{\text{a},2}$ | 7000 | cal/mol | |
Ideal gas constant | $R$ | 1.987 | cal/(mol K) | |
coolant heat transfer coeff. | $UA$ | 35000 | cal/(h K) | |
coolant temperature | $T_\text{c}$ | 298 | K | |
initial reactor molarity | $c_\text{A}(0)$ | 1 | mol/L | |
initial reactor molarity | $c_\text{B}(0)$ | 0 | mol/L | |
initial reactor molarity | $c_\text{C}(0)$ | 0 | mol/L | |
initial reactor molarity | $c_\text{S}(0)$ | 1 | mol/L | |
initial reactor temperature | $T(0)$ | 290 | K |
Note the exothermic reaction (rxn 1) 2A → B, and the endothermic reaction (rxn 2) B → 3C.
'''Parameters'''
params = dict()
params['q_flow'] = 240.0 # L/h
params['c_a_in'] = 4.0 # mol/L
params['c_b_in'] = 0.0 # mol/L
params['c_c_in'] = 0.0 # mol/L
params['temp_in'] = 305.0 # K
params['volume_0'] = 100.0 # L
params['rho'] = 1000.0 # g/L
params['heat_capacity_a'] = 30 # cal/mol/K
params['heat_capacity_b'] = 60 # cal/mol/K
params['heat_capacity_c'] = 20 # cal/mol/K
params['heat_capacity_s'] = 35 # cal/mol/K
params['enthalpy_rxn_1'] = - 6.5e3 # cal/mol
params['enthalpy_rxn_2'] = 8.0e3 # cal/mol
params['k_0_1'] = 1.25 # 1/h
params['energy_a_1'] = 9.5e3 # cal/mol
params['temp_ref_1'] = 320 # K
params['k_0_2'] = 0.08 # 1/h
params['energy_a_2'] = 7.0e3 # cal/mol
params['temp_ref_2'] = 300 # K
params['r_gas_cte'] = 1.987
params['u_a'] = 3.5e4 # J/min/K
params['temp_c'] = 298 # K
params['c_a_0'] = 1.0 # mol/L
params['c_b_0'] = 0.0 # mol/L
params['c_c_0'] = 0.0 # mol/L
params['c_s_0'] = 1.0 # mol/L
params['temp_0'] = 290.0 # K
The fluid is assumed incompressible and ideal thus
\begin{equation*} V(t) = V_0 + q\,t . \end{equation*}where $\tau = \frac{V}{q}$ is the flow residence time in the reactor, $g_\text{A}(t)$ is the species A production rate; similarly for species B, $g_\text{B}(t)$ and C, $g_\text{C}(t)$. Note that these equations are the classical relaxation process with a generation term where the relaxation time is the flow residence time $\tau$.
The species production rates are obtained from the stoichiometric relation
\begin{equation*} \Smtrx^\top\,\rvec = \gvec \end{equation*}for this problem it gives
\begin{equation*} \begin{pmatrix} g_\text{A} \\ g_\text{B} \\ g_\text{C} \end{pmatrix} = \begin{pmatrix} -2 & 1 & 0 \\ 0 & -1 & 3 \end{pmatrix}^\top \, \begin{pmatrix} k_1\,c_\text{A} \\ k_2\,c_\text{B} \end{pmatrix} = \begin{pmatrix} -2\,k_1\,c_\text{A} \\ \,k_1\,c_\text{A} - k_2\,c_\text{B} \\ 3\,k_2\,c_\text{B} \end{pmatrix}, \end{equation*}where the reaction rate constant is expressed in the Arrhenius form $k_i = k_{0,i}\,e^{\frac{-E}{R} \bigl( \frac{1}{T} - \frac{1}{T_{\text{R}i}} \bigr) } \quad\ \text{for}\ \quad i=1,2$.
'''Stoichiometric matrix'''
# (2A → B)
# (B → 3C)
import numpy as np
stoic_mtrx = np.array([[-1.0, 0.5, 0.0],
[ 0.0,-1.0, 3.0]]) # stoichiometric matrix
params['stoic_mtrx'] = stoic_mtrx
where $r_1 = k\,c_\text{A}$ and $r_2 = k\,c_\text{B}$ are the reaction rates, and the heat transfer model for the heat exchanger is
\begin{equation*} \dot{Q} = - U\,A\bigl(T-T_\text{c}\bigr) . \end{equation*}The number of moles of the catalyst is $n_\text{S} = V_0\,c_\text{S}(0)$ and it is assumed to be a constant, the heat capacity of the mixture is
\begin{equation*} C(t) = \frac{n_\text{S}\,C_\text{pS}}{V(t)} +\sum\limits_{j=\text{A},\text{B},\text{C}} c_j\,C_{\text{p}j} . \end{equation*}There exists four equations (below) and four unknowns, namely the molarity of three species, $c_j(t)$, and the temperature of the reactor, $T(t)$. Therefore there are as many unknowns as there are equations and the problem is potentially solvable give the initial conditions for the unknowns, and all the values of the parameters in foregoing the table.
A vector notation for the foregoing system of equations greatly improves the generality of the derived computer code. First, the usage of the stoichiometric matrix is instrumental for the case with multiple chemical reactions, where the reaction rates can be related to the species production rates that appear as sources in the species mass balances
Towards this goal let us define
\begin{equation*} \frac{d\uvar}{dt} = \fvar( \uvar, t ) \end{equation*}where $\uvar(t) = (u_1,u_2,u_3,u_4)$ and we assign
\begin{align*} u_1(t)&=c_\text{A}(t),\\ u_2(t)&=c_\text{B}(t),\\ u_3(t)&=c_\text{C}(t),\\ u_4(t)&=T . \end{align*}Also for $\fvar(\uvar,t) = \bigl(f_1(\uvar,t), f_2(\uvar,t), f_3(\uvar,t), f_4(\uvar,t)\bigr)$ we assign
\begin{align*} f_1 & = - \frac{1}{\tau}\,(u_1-c_\text{Ain}) + g_\text{A}, \\ f_2 & = - \frac{1}{\tau}\,(u_2-c_\text{Bin}) + g_\text{B}, \\ f_3 & = - \frac{1}{\tau}\,(u_3-c_\text{Bin}) + g_\text{C}, \\ f_4 & = - \frac{1}{\tau}\,(u_4-T_\text{in}) \sum\limits_{j=\text{A},\text{B},\text{C}} \frac{c_{j\text{in}}\,C_{\text{p}j}}{C(t)} - \sum\limits_{i=1}^{2} \frac{\Delta H_{\text{R},i}\,r_i}{C(t)} + \frac{\dot{Q}}{V(t)\,C(t)} , \end{align*}where $C(t)$ is a function of $\uvar(t)$.
Finally, the initial conditions given are as follows
\begin{align*} u_1(0)&=c_\text{A}(0),\\ u_2(0)&=c_\text{B}(0),\\ u_3(0)&=c_\text{C}(0),\\ u_4(0)&=T(0). \end{align*}'''ODE function'''
def f_vec( t, u_vec, params ):
c_a = u_vec[0]
c_b = u_vec[1]
c_c = u_vec[2]
temp = u_vec[3]
assert c_a >= 0.0
assert c_b >= 0.0
assert c_c >= 0.0
assert temp >= 0.0
q_flow = params['q_flow']
c_a_in = params['c_a_in']
c_b_in = params['c_b_in']
c_c_in = params['c_c_in']
temp_in = params['temp_in']
volume_0 = params['volume_0']
rho = params['rho']
heat_capacity_a = params['heat_capacity_a']
heat_capacity_b = params['heat_capacity_b']
heat_capacity_c = params['heat_capacity_c']
heat_capacity_s = params['heat_capacity_s']
enthalpy_rxn_1 = params['enthalpy_rxn_1']
enthalpy_rxn_2 = params['enthalpy_rxn_2']
k_0_1 = params['k_0_1']
energy_a_1 = params['energy_a_1']
temp_ref_1 = params['temp_ref_1']
k_0_2 = params['k_0_2']
energy_a_2 = params['energy_a_2']
temp_ref_2 = params['temp_ref_2']
r_gas_cte = params['r_gas_cte']
u_a = params['u_a']
temp_c = params['temp_c']
c_s_0 = params['c_s_0']
stoic_mtrx = params['stoic_mtrx']
volume = volume_0 + q_flow * t
tau = volume / q_flow
import math
k_1 = k_0_1 * math.exp( -energy_a_1/r_gas_cte*(1/temp-1/temp_ref_1) )
k_2 = k_0_2 * math.exp( -energy_a_2/r_gas_cte*(1/temp-1/temp_ref_2) )
import numpy as np
r_vec = np.zeros(2)
r_vec[0] = k_1 * c_a # reaction rate vector
r_vec[1] = k_2 * c_b
g_vec = stoic_mtrx.transpose() @ r_vec # species production rates S^T r = g
f_0 = - 1/tau*(c_a - c_a_in) + g_vec[0]
f_1 = - 1/tau*(c_b - c_b_in) + g_vec[1]
f_2 = - 1/tau*(c_c - c_c_in) + g_vec[2]
q_dot = - u_a * (temp - temp_c)
sum_c_j_in_cp_j = c_a_in * heat_capacity_a + c_b_in * heat_capacity_b + c_c_in * heat_capacity_c
sum_c_j_cp_j = c_a * heat_capacity_a + c_b * heat_capacity_b + c_c * heat_capacity_c
n_s = volume_0 * c_s_0
c_t = n_s * heat_capacity_s / volume + sum_c_j_cp_j
f_3 = - 1/tau*(temp - temp_in) * sum_c_j_in_cp_j / c_t \
- (enthalpy_rxn_1 * r_vec[0] + enthalpy_rxn_2 * r_vec[1]) / c_t \
+ q_dot/volume/c_t
return np.array([f_0, f_1, f_2, f_3])
'''Create the CSTR run function'''
def run_cstr(f_vec, time_stamps, params):
from scipy.integrate import odeint
max_n_steps_per_time_step = 100 # max number of nonlinear algebraic solver iterations per time step
c_a_0 = params['c_a_0']
c_b_0 = params['c_b_0']
c_c_0 = params['c_c_0']
temp_0 = params['temp_0']
u_vec_0 = np.zeros(4,dtype=np.float64)
u_vec_0[0] = c_a_0
u_vec_0[1] = c_b_0
u_vec_0[2] = c_c_0
u_vec_0[3] = temp_0
(u_vec_history, info_dict) = odeint( f_vec, u_vec_0, time_stamps,
args=( params, ),
rtol=1e-7, atol=1e-7, mxstep=max_n_steps_per_time_step,
full_output=True, tfirst=True )
assert info_dict['message']=='Integration successful.',\
'Fatal: scipy.integrate.odeint failed %r'%info_dict['message']
return u_vec_history
'''Evolve the CSTR'''
tau_0 = params['volume_0'] / params['q_flow']
time_final = 4 * tau_0 # number of residence flow times to evolve
n_time_steps = 200 # number of solution values in time
import numpy as np
time_stamps = np.linspace(0.0, time_final, num=n_time_steps)
u_vec_history = run_cstr( f_vec, time_stamps, params )
'''Plot A concentration in the reactor'''
import matplotlib.pyplot as plt
fig, ax1 = plt.subplots(1, figsize=(7, 7))
#plt.ylim = (0.0, np.max(u_vec_history[:,0:4]))
#plt.ylim = (0.0, 500)
ax1.plot(time_stamps/tau_0,u_vec_history[:,0],'b-',label='A' )
ax1.set_xlabel(r'Time [$\tau_0$] ($\tau_0=$%4.1f h)'%tau_0,fontsize=16)
ax1.set_ylabel(r'$c_A$ [mol/L]',fontsize=16,color='black')
#ax1.tick_params(axis='y', labelcolor='blue', labelsize=14)
ax1.tick_params(axis='x', labelsize=14)
ax1.legend(loc='best',fontsize=12)
ax1.grid(True)
#ax2 = ax1.twinx()
ax1.plot(time_stamps/tau_0,u_vec_history[:,1],'g-',label='B' )
#ax2.set_ylabel(r'$c_B$ [mol/L]',fontsize=16,color='green')
#ax2.tick_params(axis='y', labelcolor='green', labelsize=14)
ax1.legend(loc='best',fontsize=12)
#ax2.grid(True)
#ax3 = ax1.twinx()
ax1.plot(time_stamps/tau_0,u_vec_history[:,2],'r-.',label='C')
#ax3.set_ylabel(r'$c_C$ [mol/L]',fontsize=16,color='red')
#ax3.tick_params(axis='y', labelcolor='red', labelsize=14)
ax1.legend(loc='best',fontsize=12)
#ax3.grid(True)
#ax3.spines["right"].set_position(("axes", 1.2))
plt.title('Semi-Batch CSTR w/ Heat Exchange Coil',fontsize=20)
fig.tight_layout() # otherwise the right y-label is slightly clipped
plt.show()
print('')
'''Plot reactor temperature'''
import matplotlib.pyplot as plt
fig, ax1 = plt.subplots(1, figsize=(10, 7))
ax1.plot(time_stamps/tau_0,u_vec_history[:,3],'r-',label='$T_c=$ '+str(params['temp_c']) )
ax1.set_xlabel(r'Time [$\tau_0$] ($\tau_0=$%4.1f h)'%tau_0,fontsize=16)
ax1.set_ylabel(r'$T$ [K]',fontsize=16,color='red')
ax1.tick_params(axis='y', labelcolor='red', labelsize=14)
ax1.tick_params(axis='x', labelsize=14)
ax1.legend(loc='upper left',fontsize=12)
ax1.grid(True)
#ax2 = ax1.twinx()
#k_cte = params['k_0']*np.exp(-params['energy_a_over_r']/u_vec_history[:,2])
#ax2.plot(time_stamps/tau_0,k_cte,'-',color='navy',label='$T_c=$ '+str(params['temp_c']) )
#ax2.set_ylabel(r'$k$ [min$^{-1}$]',fontsize=16,color='navy')
#ax2.tick_params(axis='y', labelcolor='navy', labelsize=14)
#ax2.legend(loc='best',fontsize=14)
#ax2.grid(True)
ax3 = ax1.twinx()
q_dot = - params['u_a'] * (u_vec_history[:,3] - params['temp_c']) / 3600.0
ax3.plot(time_stamps/tau_0,q_dot,'-',color='green',label='$T_c=$ '+str(params['temp_c']) )
ax3.set_ylabel(r'$\dot{Q}$ [W]',fontsize=16,color='green')
ax3.tick_params(axis='y', labelcolor='green', labelsize=14)
#ax3.legend(loc='best',fontsize=14)
#ax2.grid(True)
ax3.spines["right"].set_position(("axes", 1.2))
plt.title('Semi-Batch CSTR w/ Heat Exchange Coil',fontsize=20)
fig.tight_layout() # otherwise the right y-label is slightly clipped
#plt.xticks(fontsize=14)
#plt.yticks(fontsize=14)
#plt.legend(loc='best',fontsize=12)
#plt.grid(True)
plt.show()
print('')
'''Create interactive plot'''
import plotly.offline as py
py.init_notebook_mode()
import plotly.graph_objs as go
from ipywidgets import interact
layout=go.Layout(title="Semi-Batch CSTR w/ Heat Exchange Coil",
xaxis={'title':'Time [tau]'},
yaxis=dict(side='left',title='HELLO'),
yaxis2=dict(overlaying='y',anchor='x',side='right',title='HELLO')
)
fig = go.FigureWidget(layout=layout)
scatt = fig.add_scatter()
fig.layout.titlefont.size = 22
fig.layout.titlefont.family = 'Rockwell'
fig.layout.xaxis.title = 'Time [tau]'
fig.layout.yaxis.title = 'c_A [mol/L]'
#fig.layout.yaxis2.title = 'T [K]'
'''Vary the coolant temperature and inflow molarity of A'''
@interact(coolant_T=(290.0, 310.0, 1.0), inflow_A=(0.5, 1.5, 0.05),
select=['Temperature', 'Heat_Pwr', 'Molarity_A', 'Molarity_B', 'Molarity_C'])
def update(coolant_T=290.0, inflow_A=1.0, select='Molarity_A'):
with fig.batch_update():
scatt.x=time_stamps/tau_0
params['temp_c'] = coolant_T
params['c_a_in'] = inflow_A
history = run_cstr( f_vec, time_stamps, params )
if select == 'Molarity_A':
scatt.y=history[:,0]
scatt.line.color='blue'
fig.layout.yaxis.title = 'c_A [mol/L]'
elif select == 'Molarity_B':
scatt.y=history[:,1]
scatt.line.color='green'
fig.layout.yaxis.title = 'c_B [mol/L]'
elif select == 'Molarity_C':
scatt.y=history[:,2]
scatt.line.color='green'
fig.layout.yaxis.title = 'c_C [mol/L]'
elif select == 'Heat_Pwr':
q_dot = - params['u_a'] * (history[:,3] - params['temp_c']) / 3600
scatt.y=q_dot
scatt.line.color='red'
fig.layout.yaxis.title = 'Q [W]'
else:
scatt.y=history[:,3]
scatt.line.color='red'
fig.layout.yaxis.title = 'T [K]'
fig