Figure 6.35:

Molar flow of o-xylene versus reactor length for different feed temperatures.

Code for Figure 6.35

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
# Converted from xyleneT.m - o-xylene oxidation PFR temperature profiles
import numpy as np
from scipy.integrate import solve_ivp
from misc import save_ascii

R_tube = 1.25e-2      # tube radius, m
Ac     = np.pi*R_tube**2
Qrho   = 0.0026371    # mass flow, kg/sec
Pf     = 1.01e2       # feed pressure, kN/m^2
Mwf    = 0.98*(0.79*28+0.21*32)+0.02*106.17  # mol wt feed, kg/kmol
Tf     = 625.         # feed temperature, K
Rg     = 8.314        # kJ/(K kmol)
Nf     = Qrho / Mwf   # kmol/sec
E      = 13636.       # K
Tm     = 625.         # K
km     = 2.0822       # 1/sec
Ta     = 625.         # K
Cp     = 0.992        # kJ/kg K
delH   = -1.284e6     # kJ/kmol
U      = 0.373        # kJ/(m^2 sec K)
beta   = delH * Ac / (Qrho*Cp)
Gamma  = 2*np.pi*R_tube*U / (Qrho*Cp)
l      = 1.5          # m

p = dict(Nf=Nf, Pf=Pf, Rg=Rg, E=E, Tm=Tm, km=km, Ta=Ta, beta=beta, Gamma=Gamma, Ac=Ac)

def pfr(t, x, p):
    Nx = x[0]; T = x[1]
    Q  = p['Nf'] / (p['Pf']/(p['Rg']*T))
    cx = Nx / Q
    k  = p['km'] * np.exp(-p['E']*(1./T - 1./p['Tm']))
    rate = k * cx
    return [-rate*p['Ac'],
            -p['beta']*rate + p['Gamma']*(p['Ta']-T)]

npts   = 200
z      = np.linspace(0., l, npts)
yxfeed = 0.019
Nxf    = yxfeed * Nf
Tfeed  = [615., 620., 625., 630.]
nfeed  = len(Tfeed)
yx_all = np.zeros((npts, nfeed))
T_all  = np.zeros((npts, nfeed))

eps_sq = np.sqrt(np.finfo(float).eps)
for i, Tf_i in enumerate(Tfeed):
    x0_i = [Nxf, Tf_i]
    sol  = solve_ivp(lambda t, x: pfr(t, x, p), [0., l], x0_i,
                     method='Radau', t_eval=z, rtol=eps_sq, atol=eps_sq)
    yx_all[:, i] = sol.y[0, :]
    T_all[:, i]  = sol.y[1, :]

table = np.column_stack([z, yx_all, T_all])
save_ascii('xyleneT.dat', table)