Figure 8.22:

Total concentration of A in the reactor effluent versus particle size.

Code for Figure 8.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
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
# Converted from suspension.m - suspension reactor with particle size variation
import numpy as np
from scipy.optimize import fsolve
from misc import colloc, octave_save

k1      = 1e3      # cm^3/mol min
kma     = 0.01/60  # cm/min
kmb     = 0.01/60  # cm/min
alpha   = 1.
caf     = 1e-3     # mol/cm^3
cbf     = 1e-3     # mol/cm^3
thetabar = 10.     # min

ncolpt = 20
z_col, A_col, B_col, Q_col = colloc(ncolpt-1, 'left')
theta   = z_col / (1. - z_col)
p_rtd   = np.exp(-theta/thetabar) / thetabar
Qz      = Q_col / (1.-z_col)**2

def particles(x, r):
    ca    = x[:ncolpt]
    cb    = x[ncolpt:2*ncolpt]
    cabar = x[2*ncolpt]
    cbbar = x[2*ncolpt+1]

    inta  = Qz @ (ca * p_rtd)
    intb  = Qz @ (cb * p_rtd)
    rx    = k1 * ca * cb
    intrx = thetabar * (Qz @ (rx * p_rtd))
    first = np.concatenate([A_col @ ca, A_col @ cb])
    rhs   = np.zeros(2*ncolpt + 2)
    rhs[:2*ncolpt] = first - np.concatenate([kma*3./r*(cabar-ca) + rx,
                                              kmb*3./r*(cbbar-cb) + rx])
    # initial condition (t=0 BC)
    rhs[0]       = ca[0] - caf
    rhs[ncolpt]  = cb[0]
    # total balances
    rhs[2*ncolpt]   = (caf - inta - intrx
                        - k1*cabar*cbbar*thetabar/alpha - cabar/alpha)
    rhs[2*ncolpt+1] = (cbf - alpha*intb - alpha*intrx
                        - k1*cabar*cbbar*thetabar - cbbar)
    return rhs

r0   = 1.;  rfin = 1e-5;  nrs = 50
rvec = np.logspace(np.log10(r0), np.log10(rfin), nrs)
cabar_v = np.zeros(nrs);  cbbar_v = np.zeros(nrs)
catot   = np.zeros(nrs);  cbtot   = np.zeros(nrs)

x0 = np.concatenate([caf*np.ones(ncolpt), cbf*np.ones(ncolpt), [caf, cbf]])
for i, r in enumerate(rvec):
    x, _, info, _ = fsolve(lambda x: particles(x, r), x0, full_output=True)[:4]
    if info != 1:
        print(f'warning: fsolve failed at r={r}')
    x0 = x.copy()
    ca_s = x[:ncolpt];  cb_s = x[ncolpt:2*ncolpt]
    cabar_v[i] = x[2*ncolpt];  cbbar_v[i] = x[2*ncolpt+1]
    inta  = Qz @ (ca_s * p_rtd)
    intb  = Qz @ (cb_s * p_rtd)
    catot[i] = (inta*alpha + cabar_v[i])/(1.+alpha)
    cbtot[i] = (intb*alpha + cbbar_v[i])/(1.+alpha)

# CSTR solution
def cstr_eqs(x):
    ca, cb = x
    return [caf*alpha/(1.+alpha) - ca - k1*ca*cb*thetabar,
            cbf/(1.+alpha) - cb - k1*ca*cb*thetabar]

y, _, info, _ = fsolve(cstr_eqs, [caf, cbf], full_output=True)[:4]
cacstr = y[0];  cbcstr = y[1]

table1 = np.column_stack([rvec*1e4, 1000.*catot,
                           1000.*cacstr*np.ones(nrs),
                           1000.*caf*alpha/(1.+alpha)*np.ones(nrs)])

# second set: three specific r values
rvec2 = np.array([1e-4, 1e-3, 1e-2])
table2_cols = [theta]
x0 = np.concatenate([caf*np.ones(ncolpt), cbf*np.ones(ncolpt), [caf, cbf]])
for r in rvec2:
    x, _, info, _ = fsolve(lambda x: particles(x, r), x0, full_output=True)[:4]
    if info != 1:
        print(f'warning: fsolve failed at r2={r}')
    x0 = x.copy()
    ca_s = x[:ncolpt];  cb_s = x[ncolpt:2*ncolpt]
    table2_cols.append(1000.*ca_s)
    table2_cols.append(1000.*cb_s)

table2_cols.append(z_col)
table2_cols.append(p_rtd)
table2 = np.column_stack(table2_cols)

octave_save('suspension.dat',
            ('table1', table1),
            ('table2', table2))