Figure 4.23:

Molar flowrate of ethane versus reactor volume for inlet temperatures of 1000, 1050 and 1100~K.

Code for Figure 4.23

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 ethane_NO_2.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
R1 = 82.057
Tvec = np.array([1000., 1050., 1100.])
P  = 1.0
y1f = 0.95; y2f = 0.05
Qf = 600.0

class Params:
    pass
p = Params()
p.P = P; p.R1 = R1

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
    ])

v = np.linspace(0, 1500, 200)
sqeps = np.sqrt(np.finfo(float).eps)
table = np.zeros((len(v), len(Tvec)+1))
table[:, 0] = v

for j, T in enumerate(Tvec):
    p.T = T
    N1f = y1f*P/(R1*T)*Qf
    N2f = y2f*P/(R1*T)*Qf
    N0 = np.array([N1f, N2f, 0., 0., 0., 0., 0.])
    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]
    sol = solve_ivp(lambda v, x: rxrate(v, x, p), [0, 1500], N0,
                    method='Radau', t_eval=v, rtol=sqeps, atol=sqeps)
    table[:, j+1] = sol.y[0]  # store ethane

save_ascii('ethane_NO_2.dat', table)