Figure 6.23:

Eigenvalues of the Jacobian matrix versus reactor conversion in the region of steady-state multiplicity.

Code for Figure 6.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
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
# Converted from st_st_osc_eigs.m - Oscillatory system eigenvalues (s-shaped part)
# Saves rows 1:npts1+npts2 (the s-shaped segment, lower branch)
import numpy as np
from scipy.optimize import fsolve
from misc import save_ascii

p = dict(
    k_m      = 0.004,
    T_m      = 298.,
    E        = 15000.,
    c_Af     = 2.,
    C_p      = 4.,
    rho      = 1000.,
    T_f      = 298.,
    T_a      = 298.,
    DeltaH_R = -2.2e5,
    U        = 340.,
    theta    = 35.,
    T_set    = 321.53,
    T_fs     = 298.,
    Kc       = 0.,
)
p['C_ps'] = p['C_p'] * p['rho']

def st_st_cA(x, c_A, p):
    theta = x[0]; T = x[1]
    k = p['k_m'] * np.exp(-p['E']*(1./T - 1./p['T_m']))
    return [p['c_Af'] - (1+k*theta)*c_A,
            p['U']*theta*(p['T_a']-T) + p['C_ps']*(p['T_f']-T) - k*theta*c_A*p['DeltaH_R']]

def st_st_theta(x, theta, p):
    c_A = x[0]; T = x[1]
    k = p['k_m'] * np.exp(-p['E']*(1./T - 1./p['T_m']))
    return [p['c_Af'] - (1+k*theta)*c_A,
            p['U']*theta*(p['T_a']-T) + p['C_ps']*(p['T_f']-T) - k*theta*c_A*p['DeltaH_R']]

def compute_lamrow(k, c_A, T, theta, p):
    Jac = np.array([[-1./theta - k, -k*c_A*p['E']/(T*T)],
                    [-k*p['DeltaH_R']/p['C_ps'],
                     -p['U']/p['C_ps'] - 1./theta - k*c_A*p['DeltaH_R']/p['C_ps']*p['E']/(T*T)]])
    lam = np.linalg.eigvals(Jac)
    a = lam[0].real; b = lam[0].imag
    c = lam[1].real; d = lam[1].imag
    return [a,b,c,d] if a >= c else [c,d,a,b]

# Lower branch (logspace in conversion)
npts1 = 200
x0 = [1., p['T_f']]
xvect = np.logspace(np.log10(1e-4), np.log10(0.75), npts1)
c_Avect = (1 - xvect) * p['c_Af']
table = []

for c_A in c_Avect:
    X, _, ier, _ = fsolve(lambda x: st_st_cA(x, c_A, p), x0, full_output=True)
    theta = X[0]; T = X[1]
    conv  = (p['c_Af'] - c_A) / p['c_Af']
    k = p['k_m'] * np.exp(-p['E']*(1./T - 1./p['T_m']))
    lamrow = compute_lamrow(k, c_A, T, theta, p)
    table.append([theta, T, conv] + lamrow + [float(ier)])
    x0 = X

conv_c = conv

# Near-extinction patch (npts2=200)
npts2 = 200
xvect2 = np.linspace(conv_c, 0.95, npts2)
c_Avect2 = (1 - xvect2) * p['c_Af']

for c_A in c_Avect2:
    X, _, ier, _ = fsolve(lambda x: st_st_cA(x, c_A, p), x0, full_output=True)
    theta = X[0]; T = X[1]
    conv  = (p['c_Af'] - c_A) / p['c_Af']
    k = p['k_m'] * np.exp(-p['E']*(1./T - 1./p['T_m']))
    lamrow = compute_lamrow(k, c_A, T, theta, p)
    table.append([theta, T, conv] + lamrow + [float(ier)])
    x0 = X

# The file saves only rows 1:npts1+npts2
table = np.array(table)
tmp = table[:npts1 + npts2, :]
save_ascii('st_st_osc_eigs.dat', tmp)