Figure 4.20:

Benzene conversion versus reactor volume.

Code for Figure 4.20

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
# Converted from benz_pyrol.m
import numpy as np
from scipy.integrate import solve_ivp
from misc import octave_save

class P:
    pass

p = P()
p.NBf = 60e3
p.R   = 0.08205
p.T   = 1033.0
p.P   = 1.0
p.k1  = 7e5
p.k2  = 4e5
p.K1  = 0.31
p.K2  = 0.48

def benzene_pyrol_rhs(volume, x, p):
    ext1, ext2 = x
    NB = p.NBf - 2*ext1 - ext2
    ND = ext1 - ext2
    NH = ext1 + ext2
    NT = ext2
    Q = p.NBf * p.R * p.T / p.P
    cB = NB/Q; cD = ND/Q; cT = NT/Q; cH = NH/Q
    return [p.k1*(cB*cB - cD*cH/p.K1),
            p.k2*(cB*cD - cT*cH/p.K2)]

x0 = [0.0, 0.0]
vout = np.linspace(0, 1600, 50)

sqeps = np.sqrt(np.finfo(float).eps)
sol = solve_ivp(lambda v, x: benzene_pyrol_rhs(v, x, p),
                [0, 1600], x0, method='Radau',
                t_eval=vout, rtol=sqeps, atol=sqeps)

x = sol.y.T
conv = (2*x[:, 0] + x[:, 1]) / p.NBf
yB   = (p.NBf - 2*x[:, 0] - x[:, 1]) / p.NBf
yD   = (x[:, 0] - x[:, 1]) / p.NBf
yT   = x[:, 1] / p.NBf
yH   = (x[:, 0] + x[:, 1]) / p.NBf

data = np.column_stack([vout, conv, yB, yD, yT, yH])
aux_data = np.array([[403.25, 0.50, 403.25, 0.50],
                     [403.25, 0.0,  0.0,    0.50]])
octave_save('benz_pyrol.dat', ('data', data), ('aux_data', aux_data))