Figure 4.22:

Molar flowrate of ethane, ethylene and NO versus reactor volume for ethane pyrolysis example.

Code for Figure 4.22

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

A0 = np.array([1e14, 1e12, 3e14, 3.4e12, 1e12, 1e13])
Ea = np.array([217600., 0., 165300., 28500., 0., 200800.])
R  = 8.3144    # J/gmol-K
R1 = 82.057    # cc-atm/gmol-K
T  = 1050.0    # K
P  = 1.0       # atm
y1f = 0.95; y2f = 0.05
C1f = y1f*P/(R1*T)
C2f = y2f*P/(R1*T)
Qf  = 600.0    # cc/sec
N1f = C1f*Qf; N2f = C2f*Qf

class Params:
    pass
p = Params()
p.P = P; p.T = T; p.R1 = R1
k_vals = A0 * np.exp(-Ea/(R*T))
p.k1  = k_vals[0]; p.k_1 = k_vals[1]
p.k2  = k_vals[2]; p.k3  = k_vals[3]
p.k4  = k_vals[4]; p.k_4 = k_vals[5]

def rxrate(v, N, p):
    C = N/np.sum(N) * p.P/(p.R1*p.T)
    r1 = p.k1*C[0]*C[1] - p.k_1*C[2]*C[3]
    r2 = p.k2*C[2]
    r3 = p.k3*C[0]*C[4]
    r4 = p.k4*C[1]*C[4] - p.k_4*C[3]
    return np.array([
        -r1 - r3,
        -r1 - r4,
         r1 - r2 + r3,
         r1 + r4,
         r2 - r3 - r4,
         r2,
         r3
    ])

N0 = np.array([N1f, N2f, 0., 0., 0., 0., 0.])
v = np.linspace(0, 1500, 200)
sqeps = np.sqrt(np.finfo(float).eps)
sol = solve_ivp(lambda v, x: rxrate(v, x, p), [0, 1500], N0,
                method='Radau', t_eval=v, rtol=sqeps, atol=sqeps)
N = sol.y.T
table = np.column_stack([v, N])
save_ascii('ethane_NO_1.dat', table)