Figure 7.28:

Dimensionless equilibrium constant and Thiele modulus versus reactor volume.

Code for Figure 7.28

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
# Converted from fbhw.m - CO oxidation fixed-bed reactor (Hougen-Watson kinetics)
import numpy as np
from scipy.integrate import solve_ivp
from misc import save_ascii

p = dict(
    T    = 838.,
    Rg   = 82.06,
    P    = 1.,
    rhop = 0.68,
    rhob = 0.60,
    Dco  = 0.0487,
    Rp   = 0.1,
    xco  = 0.95,
)
p['ccof'] = 0.167 * (p['P']/(p['Rg']*p['T']))
p['co2f'] = 0.833 * (p['P']/(p['Rg']*p['T']))
p['cIf']  = 0.
p['Kco']  = 8.099e6 * np.exp(409./p['T'])
p['k']    = 1.3828e19 * np.exp(-13500./p['T']) * p['co2f']
p['Qf']   = 792e3
p['Ncof'] = p['Qf'] * p['ccof']
p['No2f'] = p['Qf'] * p['co2f']
p['NIf']  = p['Qf'] * p['cIf']
p['a']    = p['Rp']/3.

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

def pbr(t, x, p):
    Nco = x[0]
    Nt  = p['No2f'] + (Nco + p['Ncof'])/2. + p['NIf']
    cco = p['P']/(p['Rg']*p['T']) * Nco/Nt
    phi = p['Kco']*cco
    Phi = (phi/(1+phi) * np.sqrt(p['k']*p['a']**2 /
                                  (2*p['Dco']*(phi - np.log(1+phi)))))
    eta = 1. if Phi <= 1. else 1./Phi
    return [-p['rhob']/p['rhop'] * eta * p['k'] * cco/(1 + p['Kco']*cco)]

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

npts   = 100
vfinal = 500.
vsteps = np.linspace(0., vfinal, npts)
x0     = np.array([p['Ncof']])

sol = solve_ivp(lambda t, x: pbr(t, x, p), [0., vfinal], x0,
                method='Radau', t_eval=vsteps,
                events=stop_event,
                rtol=eps_sq, atol=eps_sq)
if len(sol.t_events[0]) > 0:
    vplot = np.append(sol.t, sol.t_events[0][0])
    Nco   = np.append(sol.y[0], sol.y_events[0][0, 0])
else:
    vplot = sol.t
    Nco   = sol.y[0]
No2   = p['No2f'] + (Nco - p['Ncof'])/2.
Nco2  = p['Ncof'] - Nco
Nt    = Nco + No2 + Nco2 + p['NIf']
cco   = p['P']/(p['Rg']*p['T']) * Nco/Nt
co2   = p['P']/(p['Rg']*p['T']) * No2/Nt
cco2  = p['P']/(p['Rg']*p['T']) * Nco2/Nt
phi   = p['Kco']*cco
Phi   = (phi/(1+phi) * np.sqrt(p['k']*p['a']**2 /
                                (2*p['Dco']*(phi - np.log(1+phi)))))
eta   = 1./Phi

table = np.column_stack([vplot, cco, co2, cco2, phi, Phi, eta])
save_ascii('fbhw.dat', table)