Figure 4.14:

Reaching steady state in a CSTR.

Figure 4.14

Code for Figure 4.14

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 strcstr.m
import numpy as np
from scipy.integrate import solve_ivp
from misc import save_ascii

class P:
    pass

p = P()
p.caf = 2.0
p.cao = 0.0
p.k = 0.1
p.theta = 100.0
tf = 120.0

c_ss = p.caf / (1 + p.k*p.theta)

def mat_bal(t, x, p):
    ca = x[0]
    return [(1/p.theta)*(p.caf - ca) - p.k*ca]

t = np.arange(0, tf+0.1, 0.1)
sqeps = np.sqrt(np.finfo(float).eps)
tff = int(10*tf) + 1

# Case 1: cao = 0
sol1 = solve_ivp(lambda t, x: mat_bal(t, x, p), [0, tf], [p.cao],
                 method='Radau', t_eval=t, rtol=sqeps, atol=sqeps)
solution = sol1.y.T

# Case 2: cao = 2
p.cao = 2.0
sol2 = solve_ivp(lambda t, x: mat_bal(t, x, p), [0, tf], [p.cao],
                 method='Radau', t_eval=t, rtol=sqeps, atol=sqeps)
solution1 = sol2.y.T

answer = np.column_stack([t, solution[:, 0],
                           c_ss*np.ones(len(t)),
                           solution1[:, 0]])
save_ascii('strcstr.dat', answer)