Figure 6.30:

Phase portrait of conversion versus temperature at showing stable and unstable limit cycles.

Code for Figure 6.30

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
# Converted from phase_portrait_2.m - Phase portrait with stable+unstable limit cycles (theta=73.1)
import numpy as np
from scipy.integrate import solve_ivp
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    = 73.1,
    T_set    = 321.53,
    c_set    = 0.48995,
    T_fs     = 298.,
    Kc       = 0.,
)
p['C_ps'] = p['C_p'] * p['rho']

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

def reverserhs(t, x, p):
    vals = rhs(t, x, p)
    return [-vals[0], -vals[1]]

eps_sq = np.sqrt(np.finfo(float).eps)

# Find stable limit cycle: integrate forward to settle, then one cycle
x0 = [(1-0.8)*p['c_Af'], 310.]
sol = solve_ivp(lambda t, x: rhs(t, x, p), [0., 40.*p['theta']], x0,
                method='Radau', t_eval=np.linspace(0., 40.*p['theta'], 3),
                rtol=eps_sq, atol=eps_sq)
x0 = list(sol.y[:, -1])
tout = np.linspace(0., 3.*p['theta'], 200)
sol  = solve_ivp(lambda t, x: rhs(t, x, p), [0., 3.*p['theta']], x0,
                 method='Radau', t_eval=tout, rtol=eps_sq, atol=eps_sq)
x    = sol.y.T
u    = (x[:, 1] - p['T_set']) * p['Kc'] + p['T_fs']
conv = (p['c_Af'] - x[:, 0]) / p['c_Af']
stablim = np.column_stack([tout, x, conv, u])

# Find unstable limit cycle: integrate reverse to settle, then one forward cycle
x0 = [(1-0.7)*p['c_Af'], 305.]
sol = solve_ivp(lambda t, x: reverserhs(t, x, p), [0., 40.*p['theta']], x0,
                method='Radau', t_eval=np.linspace(0., 40.*p['theta'], 3),
                rtol=eps_sq, atol=eps_sq)
x0 = list(sol.y[:, -1])
# use rhs (not reverserhs) for the one-cycle trace (matching Octave)
sol  = solve_ivp(lambda t, x: rhs(t, x, p), [0., 3.*p['theta']], x0,
                 method='Radau', t_eval=tout, rtol=eps_sq, atol=eps_sq)
x    = sol.y.T
u    = (x[:, 1] - p['T_set']) * p['Kc'] + p['T_fs']
conv = (p['c_Af'] - x[:, 0]) / p['c_Af']
unstablim = np.column_stack([tout, x, conv, u])

table = np.column_stack([stablim, unstablim])
save_ascii('phase_portrait_2.dat', table)
print(f'max T stablim   = {np.max(stablim[:, 2]):.6g}')
print(f'max conv stablim= {np.max(stablim[:, 3]):.6g}')