Figure 8.31:

RTD for the optimal reactor configuration.

Code for Figure 8.31

Text of the GNU GPL.

main.py

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
# Converted from leven_rtd.m - RTD for optimal Levenspiel reactor
import numpy as np
from scipy.integrate import solve_ivp
from scipy.optimize import fsolve
from misc import save_ascii

eps_sq = np.sqrt(np.finfo(float).eps)

a  = 5.;   b  = 0.05
c0 = 5.;   conv = 0.95
d  = np.sqrt((1.-2.*b)**2 - 4.*b*(b+1.))
w1 = (1.-2.*b + d)/(2.*b)
w2 = (1.-2.*b - d)/(2.*b)
x1 = np.sqrt(w1/a);   x2 = np.sqrt(w2/a)

def rxrate(c, a, b):
    return c/(1.+a*c*c) + b*c

fixrate = rxrate(x2, a, b)
c1, _, info, _ = fsolve(lambda c: fixrate - rxrate(c, a, b), c0, full_output=True)[:4]
c1 = float(c1[0])

def pfr_dtdc(c, t):
    return [-1./rxrate(c, a, b)]

ncs  = 10
sol  = solve_ivp(pfr_dtdc, [c0, c1], [0.], method='Radau',
                 t_eval=np.linspace(c0, c1, ncs), rtol=eps_sq, atol=eps_sq)
theta1 = sol.y[0,-1]

c2     = x2
theta2 = (c1-c2)/rxrate(c2, a, b)

c3 = (1.-conv)*c0
sol2 = solve_ivp(pfr_dtdc, [c2, c3], [0.], method='Radau',
                 t_eval=np.linspace(c2, c3, ncs), rtol=eps_sq, atol=eps_sq)
theta3 = sol2.y[0,-1]

# RTD of PFR(theta1) -> CSTR(theta2) -> PFR(theta3)
npts    = 100
thetamin = theta1 + theta3
thetamax = thetamin + 5.*theta2
theta    = np.linspace(thetamin, thetamax, npts)
p        = (1./theta2)*np.exp(-(theta - thetamin)/theta2)
table    = np.vstack([[0., 0.], [thetamin, 0.], np.column_stack([theta, p])])
save_ascii('leven_rtd.dat', table)