Figure 6.7:

Steady-state conversion versus residence time for \Delta H_R = -3 x10^5~kJ/kmol; ignition and extinction points.

Code for Figure 6.7

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
# Converted from s_shape.m - S-shaped multiplicity curve: theta vs conv for single DeltaH
import numpy as np
from scipy.optimize import fsolve
from misc import save_ascii

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.,
    DeltaH_R = -3e5,
)
p['C_ps'] = p['rho'] * p['C_p']

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

x0     = [1., p['T_f']]
nc_As  = 200
c_Avect = np.linspace(0.995*p['c_Af'], 0.005*p['c_Af'], nc_As)
# first row: theta=0, T=Tf, conv=0, info=0
tmp_table = [[0., p['T_f'], 0., 0.]]

for c_A in c_Avect:
    X, info_d, ier, msg = 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']
    tmp_table.append([theta, T, conv, float(ier)])
    x0 = X

table = np.array(tmp_table)
save_ascii('s_shape.dat', table)