Figure 7.26:

Molar flow of A versus reactor volume for second-order, isothermal reaction in a fixed-bed reactor.

Code for Figure 7.26

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
# Converted from fborder2.m - second-order fixed-bed with two eta approximations
import numpy as np
from scipy.integrate import solve_ivp
from misc import octave_save

p = dict(
    k      = 2.25e5,
    Nafin  = 10.,
    T      = 550.,
    rhop   = 0.68,
    rhob   = 0.60,
    Da     = 0.008,
    Rp     = 0.45,
    Rg     = 82.06,
    P      = 4.,
    n      = 2.,
    xa     = 0.75,
    nvs    = 100,
    vfinal = 5e5,
)

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

def pbr(t, x, p, eta_fn):
    Na  = x[0]
    ca  = p['P']/(p['Rg']*p['T']) * Na/(2*p['Nafin'])
    Phi = p['Rp']/3 * np.sqrt((p['n']+1)/2 * p['k'] * ca / p['Da'])
    return [-p['rhob']/p['rhop'] * eta_fn(Phi) * p['k'] * ca**p['n']]

def stop_event(t, x):
    return x[0] - (1 - p['xa'])*p['Nafin']
stop_event.terminal  = True
stop_event.direction = 0

vfinal = p['vfinal']
nvs    = p['nvs']
vsteps = np.linspace(0., vfinal, nvs)
x0     = np.array([p['Nafin']])

# Exact eta
eta_exact = lambda Phi: (1./Phi) * (1./np.tanh(3*Phi) - 1./(3*Phi))
sol = solve_ivp(lambda t, x: pbr(t, x, p, eta_exact), [0., vfinal], x0,
                method='Radau', t_eval=vsteps,
                events=stop_event,
                rtol=eps_sq, atol=eps_sq)
if len(sol.t_events[0]) > 0:
    t_end1 = sol.t_events[0][0]; y_end1 = sol.y_events[0][0]
    t1 = np.append(sol.t, t_end1); y1 = np.vstack([sol.y.T, y_end1])
else:
    t1 = sol.t; y1 = sol.y.T; t_end1 = sol.t[-1]
vout  = t1 / 1000.
VR    = t_end1 / 1000.
x_sol = y1

# Asymptotic eta = 1/Phi
eta_asy = lambda Phi: 1./Phi
sol2 = solve_ivp(lambda t, x: pbr(t, x, p, eta_asy), [0., vfinal], x0,
                 method='Radau', t_eval=vsteps,
                 events=stop_event,
                 rtol=eps_sq, atol=eps_sq)
if len(sol2.t_events[0]) > 0:
    t_end2 = sol2.t_events[0][0]; y_end2 = sol2.y_events[0][0]
    t2 = np.append(sol2.t, t_end2); y2 = np.vstack([sol2.y.T, y_end2])
else:
    t2 = sol2.t; y2 = sol2.y.T; t_end2 = sol2.t[-1]
voutasy = t2 / 1000.
VRasy   = t_end2 / 1000.
xasy    = y2

table    = np.column_stack([vout, x_sol])
tableasy = np.column_stack([voutasy, xasy])
Naout  = (1 - p['xa']) * p['Nafin']
Natop  = (1 - p['xa'] + 0.10) * p['Nafin']
dashedlines = np.array([
    [0,       Naout, VRasy, 0,     VR,    0,     VR,    2, VRasy, 2],
    [1.1*VR,  Naout, VRasy, Natop, VR,    Natop, VR,    3, VRasy, 3],
])

octave_save('fborder2.dat',
            ('table', table), ('tableasy', tableasy), ('dashedlines', dashedlines))