Figure 6.5:

Steady-state conversion versus residence time for different values of the heat of reaction.

Code for Figure 6.5

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
# Converted from st_st.m - Steady state multiplicity for multiple DeltaH values
import numpy as np
from scipy.optimize import fsolve
from misc import octave_save

p = dict(
    k_m  = 0.001,
    T_m  = 298.,
    E    = 8000.,
    c_Af = 2.,
    C_p  = 4.,
    rho  = 1000.,
    T_f  = 298.,
    T_a  = 298.,
    U    = 0.,
)
p['C_ps'] = p['rho'] * p['C_p']

def st_st_cA(x, c_A, DeltaH_R, 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*DeltaH_R]

nc_As   = 500
c_Avect = np.linspace(0.995*p['c_Af'], 0.002*p['c_Af'], nc_As)
DelHvec = [-3e5, -2e5, -1e5, -0.5e5, 0., 5e4]
nH      = len(DelHvec)

table_list = []
for j, DeltaH_R in enumerate(DelHvec):
    tmp_table = np.zeros((nc_As, 4))
    x0 = [1., p['T_f']]

    for i, c_A in enumerate(c_Avect):
        X, info_d, ier, msg = fsolve(lambda x: st_st_cA(x, c_A, DeltaH_R, p), x0, full_output=True)
        theta = X[0]; T = X[1]
        conv  = (p['c_Af'] - c_A) / p['c_Af']
        tmp_table[i, :] = [theta, T, conv, float(ier)]
        x0 = X

    table_list.append(tmp_table)

# Save each table separately since octave_save doesn't support cell arrays
vars_to_save = tuple((f'table{j}', t) for j, t in enumerate(table_list))
octave_save('st_st.dat', *vars_to_save)