Figure 8.30:

Inverse of reaction rate versus concentration; optimal sequence to achieve 95% conversion is PFR--CSTR--PFR.

Code for Figure 8.30

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
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
# Converted from leven.m - Levenspiel optimum reactor design
import numpy as np
from scipy.integrate import solve_ivp
from scipy.optimize import fsolve
from misc import octave_save

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):
    return c/(1.+a*c*c) + b*c

c   = np.linspace(0.01, 10., 300)
y   = 1./rxrate(c)
cy  = np.column_stack([c, y])

# find c1: outlet of first PFR, where r(c1) = r(x2)
fixrate = rxrate(x2)
c1, _, info, _ = fsolve(lambda c: fixrate - rxrate(c), c0, full_output=True)[:4]
c1 = c1[0]

# size first PFR: integrate dc/dt = -1/r(c) from c0 to c1
# using c as independent variable
def pfr_dtdc(c, t):
    return [-1./rxrate(c)]

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

# CSTR: c2 = x2
c2     = x2
theta2 = (c1 - c2) / rxrate(c2)

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

# augment cy with special points
cy_extra = np.array([[c0, 1./rxrate(c0)],
                     [c1, 1./rxrate(c1)],
                     [c2, 1./rxrate(c2)],
                     [c3, 1./rxrate(c3)]])
cy = np.vstack([cy, cy_extra])
cy = cy[np.argsort(cy[:,0])]

# optimal: PFR(theta1) -> CSTR(theta2) -> PFR(theta3)

# segregated reactor
def seg(z, x):
    c_val, ctot = x
    if z < 1.:
        t = z/(1.-z)
        # RTD of PFR-CSTR-PFR in series
        if t <= theta1 + theta3:
            pr = 0.
        else:
            pr = np.exp(-(t-(theta1+theta3))/theta2) / theta2
        dcdz  = -rxrate(c_val) / (1.-z)**2
        dctot = pr*c_val / (1.-z)**2
    else:
        dcdz  = 0.
        dctot = 0.
    return [dcdz, dctot]

npts = 50
zseg = np.linspace(0., 1., npts)
sol_seg = solve_ivp(seg, [0., 1.], [c0, 0.], method='Radau',
                    t_eval=zseg, rtol=eps_sq, atol=eps_sq)
csegfin = sol_seg.y[1,-1]

# maximum mixedness
cinf, _, info_c, _ = fsolve(lambda c: c - c0 + rxrate(c)*theta2,
                              c0, full_output=True)[:4]
cinf = cinf[0]

def mm(z, c):
    t = (1.-z)/z
    f = 0. if t <= theta1 + theta3 else 1./theta2
    return [-(f*(c[0]-c0) + rxrate(c[0])) / (z*z)]

# Start from small epsilon: at z=eps, c=cinf gives dc/dz~0 (cinf is equilibrium),
# avoiding the 1/z^2 singularity that trips up Radau at z=0.
z_start = 1e-4
sol_mm = solve_ivp(mm, [z_start, 1.], [cinf], method='Radau',
                   t_eval=np.linspace(z_start, 1., npts), rtol=eps_sq, atol=eps_sq)
cmmfin = sol_mm.y[0,-1]

# build lines table
lines = np.array([
    [c0, 1./rxrate(c0), c1, 1./rxrate(c1), c2, 1./rxrate(c2), c3, 1./rxrate(c3), c1, 1./rxrate(c1)],
    [c0, 0.,            c1, 0.,            c2, 0.,            c3, 0.,            c2, 1./rxrate(c2)]
])

# region polygons
ind1   = (cy[:,0] >= c1) & (cy[:,0] <= c0)
ind3   = (cy[:,0] >= c3) & (cy[:,0] <= c2)
regpfr1  = np.vstack([cy[ind1], [c0, 0.], [c1, 0.], [c1, 1./rxrate(c1)]])
regpfr3  = np.vstack([cy[ind3], [c2, 0.], [c3, 0.], [c3, 1./rxrate(c3)]])
regcstr2 = np.array([[c2, 1./rxrate(c2)],
                     [c1, 1./rxrate(c1)],
                     [c1, 0.],
                     [c2, 0.],
                     [c2, 1./rxrate(c2)]])

octave_save('leven.dat',
            ('cy',      cy),
            ('lines',   lines),
            ('regpfr1', regpfr1),
            ('regcstr2',regcstr2),
            ('regpfr3', regpfr3))