4. Non-ideal flow reactor
Written on July, 2025 by Bolton TranThis post in under construction.
This reactor problem is adapted from Fogler-Elements of Chemical Reaction Engineering 6th Ed, example 18-4. A COMSOL solver for this problem is also available here for comparison. Although Fogler’s and the COMSOL solution are at steady-state, I change this to a non-steady-state problem for the extra challenge.
Reactor settings:
1) This is a tubular laminar flow reactor (LFR). $r$ is radial width and $z$ is axial length of the reactor.
2) The reaction is: $A + B \overset{k_1}{\rightarrow} C$ and is exothermic.
3) The reactor operates at non-isothermal condition with a water cooling jacket.
4) Beside non-ideal flow profile (not plug flow), there is also axial and radial dispersion of mass and heat.
Design equations:
There are three key design partial differential equations (PDEs):
1) Mole balance of $A$
\begin{equation}
\frac{\partial c}{\partial t} = -2\frac{v_0}{A_c}\left[1 - \left(\frac{r}{R}\right)^2 \right] \frac{\partial c}{\partial z}
+ D_d \left[ \frac{\partial^2c}{\partial z^2} + \frac{\partial^2c}{\partial r^2} + \frac{1}{r}\frac{\partial c}{\partial r} \right]
- k c
\end{equation}
2) Energy balance of reactor.
\begin{equation}
\frac{\partial T}{\partial t} = -2\frac{v_0}{A_c}\left[1 - \left(\frac{r}{R}\right)^2 \right] \frac{\partial T}{\partial z}
+ \frac{\kappa}{\rho C_p} \left[ \frac{\partial^2 T}{\partial z^2} + \frac{\partial^2 T}{\partial r^2} + \frac{1}{r}\frac{\partial T}{\partial r} \right]
- \frac{\Delta H_{\text{rxn}}}{\rho C_p} k c
\end{equation}
3) Energy balance of cooling jacket.
\begin{equation}
\frac{\partial T_c}{\partial t} = -v_c\frac{\partial T_c}{\partial z}
+ \frac{2\pi R h_c}{A_c \rho_c C_{pc}} (T(z, r=R) - T_c)
\end{equation}
TO ADD: Boundary conditions
TO ADD: Method of lines (MOL) to discretize all but time dimension; convert PDEs into system of ODEs
TO ADD: Finite difference and finite element approach
PYTHON CODE DUMP: more explaination to come
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
#Constants
kB=8.3144/1000 #kJ/mol.K
h=4.135667696E-15*96.49/3600 #kJ/mol.hr
#Reactor setup
L=4 #length, m
R=1 #radius, m
Ac=np.pi*R**2 #cross-section area, m^2
v0=0.1*3600 #flow rate, m^3/hr
cA0=1000 #mol/m^3
#Dispersion coef
Dd=1E-9*3600 #m^2/hr
#Reaction barrier & energy balance params
EA=160 #kJ/mol
dH=-80 # kJ/mol
rho=1000 # kg/m^3
Cp=4.18 # kJ/kg.K
kappa=0.559/1000 #kW/m.K
#Coolant jacket
vc=0.1/1000*3600/(Ac/5) #1D velocity m/hr
rhoc=1000 #kg/m^3
Cpc=4.18 #kJ/kg.K
hc=10000/1000 #kW/m2.K
#Length discretization
Nz=80 #discretization size
dz=L/(Nz-1) #width
z=np.linspace(0, L, Nz)
#Width discretization
Nr=20 #discretization size
dr=R/(Nr-1) #width
r=np.linspace(0, R/2, Nr) #half-tube
#Set up ODEs:
def rhs(t, yflat):
#--Initialize grids--
cflat=yflat[:Nz*Nr]
Tflat=yflat[Nz*Nr:-Nz]
Tcflat=yflat[-Nz:]
c=cflat.reshape((Nz, Nr))
T=Tflat.reshape((Nz, Nr))
Tc=Tcflat.reshape(Nz)
dcdt=np.zeros((Nz, Nr))
dTdt=np.zeros((Nz, Nr))
dTcdt=np.zeros(Nz)
#--Setup ODEs in grid elements--
for i in range(1, Nz-1):
#Coolant grid
dTcdz=(Tc[i]-Tc[i-1])/dz
dTcdt[i]=-vc*dTcdz + 2*np.pi*R/Ac*hc/(rhoc*Cpc)*(T[i, -1]-Tc[i])
for j in range(1, Nr-1):
#Concentration grid
dcdz=(c[i, j]-c[i-1, j])/dz
d2cdz2=(c[i+1, j] - 2*c[i, j] + c[i-1, j])/dz**2
dcdr=(c[i, j]-c[i, j-1])/dr
d2cdr2=(c[i, j+1] - 2*c[i, j] + c[i, j-1])/dr**2
#Temperature grid
dTdz=(T[i, j]-T[i-1, j])/dz
d2Tdz2=(T[i+1, j] - 2*T[i, j] + T[i-1, j])/dz**2
dTdr=(T[i, j]-T[i, j-1])/dr
d2Tdr2=(T[i, j+1] - 2*T[i, j] + T[i, j-1])/dr**2
#rate constant
k=(kB*T[i, j]/h)*np.exp(-EA/kB/T[i, j])
#Plug-flow
#dcdt[i, j]=-v0/Ac*dcdz - k*c[i, j] + Dd*d2cdz2 + Dd*(d2cdr2 + 1/(dr*j)*dcdr)
#Laminar-flow
dTdt[i, j]=-v0/Ac*2*(1-(dr*j/R)**2)*dTdz + (-dH)*k*c[i, j]/(rho*Cp) + kappa/(rho*Cp)*(d2Tdz2 + d2Tdr2 + 1/(dr*j)*dTdr)
dcdt[i, j]=-v0/Ac*2*(1-(dr*j/R)**2)*dcdz - k*c[i, j] + Dd*(d2cdz2 + d2cdr2 + 1/(dr*j)*dcdr)
#--Boundary conditions--
#-Mass transport for c-
#Boundary condition at r=0 (first j index) dcdr=0
dcdz=(c[1:-1, 0]-c[:-2, 0])/dz
d2cdz2=(c[2:, 0] - 2*c[1:-1, 0] + c[:-2, 0])/dz**2
d2cdr2=(c[1:-1, 2] - 2*c[1:-1, 1] + c[1:-1, 0])/dr**2 #approx as next node j=1
k=(kB*T[1:-1, 0]/h)*np.exp(-EA/kB/T[1:-1, 0])
dcdt[1:-1, 0]=-v0/Ac*2*dcdz - k*c[1:-1, 0] + Dd*(d2cdz2 + d2cdr2)
#Boundary condition at r=R (last j index): dcdr=0
d2cdz2=(c[2:, -1] - 2*c[1:-1, -1] + c[:-2, -1])/dz**2
d2cdr2=(c[1:-1, -1] - 2*c[1:-1, -2] + c[1:-1, -3])/dr**2 #approx as last node j=-2
k=(kB*T[1:-1, -1]/h)*np.exp(-EA/kB/T[1:-1, -1])
dcdt[1:-1, -1]= -k*c[1:-1, -1] + Dd*(d2cdz2 + d2cdr2 + 1/R*0)
#Boundary condition at z=0 (first i index) dcdt=0
dcdt[0,:]=0
#Cheat - Boundary condition at z=L (last i index) dcdz=0
dcdt[-1,:]=dcdt[-2,:]
#-Heat transfer for T-
#Boundary condition at r=0 (first j index) dTdr=0
dTdz=(T[1:-1, 0]-T[:-2, 0])/dz
d2Tdz2=(T[2:, 0] - 2*T[1:-1, 0] + T[:-2, 0])/dz**2
d2Tdr2=(T[1:-1, 2] - 2*T[1:-1, 1] + T[1:-1, 0])/dr**2 #approx as next node j=1
k=(kB*T[1:-1, 0]/h)*np.exp(-EA/kB/T[1:-1, 0])
dTdt[1:-1, 0]=-v0/Ac*2*dTdz + (-dH)*k*c[1:-1, 0]/(rho*Cp) + kappa/(rho*Cp)*(d2Tdz2 + d2Tdr2)
#Boundary conditions at r=R (last j index): dTdr=heat transfer
d2Tdz2=(T[2:, -1] - 2*T[1:-1, -1] + T[:-2, -1])/dz**2
dTdr_wall=-hc/kappa*(T[1:-1, -1] - Tc[1:-1]) #heat transfer with cooling jacket
d2Tdr2=(T[1:-1, -1] - 2*T[1:-1, -2] + T[1:-1, -3])/dr**2 #approx as last node
k=(kB*T[1:-1, -1]/h)*np.exp(-EA/kB/T[1:-1, -1])
dTdt[1:-1, -1]=(-dH)*k*c[1:-1, -1]/(rho*Cp) + kappa/(rho*Cp)*(d2Tdz2 + 1/R*dTdr_wall)
#Boundary condition at z=0 (first i index) dTdt=0
dTdt[0,:]=0
#Cheat condition at z=L (last i index) dTdz=0
dTdt[-1,:]=dTdt[-2,:]
#-Jacket Tc-
#Boundary condition at z=0 (first i index) dTcdt=0
dTcdt[0]=0
#Boundary condition at z=L (last i index) dTcdz=0
dTcdt[-1]=2*np.pi*R/Ac*hc/(rhoc*Cpc)*(T[-1, -1]-Tc[-1])
#--Flatten and return--
dcdtflat=dcdt.flatten()
dTdtflat=dTdt.flatten()
dTcdtflat=dTcdt.flatten()
return np.concatenate([dcdtflat, dTdtflat, dTcdtflat])
#Initialize (initial condition)
c0=np.zeros((Nz, Nr))
c0[0, :]=cA0 #open boundary at z=0
cflat0=c0.flatten()
T0=np.ones((Nz, Nr))*550
Tflat0=T0.flatten()
Tc0=np.ones(Nz)*270
Tcflat0=Tc0.flatten()
yflat0=np.concatenate([cflat0, Tflat0, Tcflat0])
#Time discretization
time_min=1E-3
time_max=1E-0
logtime=np.linspace(np.log10(time_min), np.log10(time_max), 100)
#Run solver
sol=solve_ivp(rhs, t_span=(0, 1), y0=yflat0, t_eval=10**logtime, method='BDF')
cs=sol.y[:Nz*Nr].reshape((Nz, Nr, len(logtime)))
Ts=sol.y[Nz*Nr:-Nz].reshape((Nz, Nr, len(logtime)))
Tcs=sol.y[-Nz:].reshape((Nz, len(logtime)))
#Plot
t_index=-1
#conc
c_plot=cs[:, :, t_index] # shape (Nz, Nr)
c_plot_full=np.concatenate((c_plot[:, ::-1], c_plot), axis=1)
#temp
T_plot=Ts[:, :, t_index]
T_plot_full=np.concatenate((T_plot[:, ::-1], T_plot), axis=1)
#grid
rfull=np.concatenate((-r[::-1], r))
Rgrid, Zgrid=np.meshgrid(rfull, z)
fig, axes=plt.subplots(1,2,figsize=(6, 4))
cmesh=axes[0].pcolormesh(Rgrid, Zgrid, c_plot_full, shading='auto',vmin=0, vmax=1000)
axes[0].set_aspect('equal')
fig.colorbar(cmesh, ax=axes[0])
Tmesh=axes[1].pcolormesh(Rgrid, Zgrid, T_plot_full, shading='auto')
axes[1].set_aspect('equal')
fig.colorbar(Tmesh, ax=axes[1])
plt.show()