Figure 4.18:

Semi-batch reactor polymer content for different monomer addition policies.

Code for Figure 4.18

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
# Converted from filling.m
import numpy as np
from scipy.optimize import fsolve
from misc import save_ascii

k    = 0.1   # min^-1
rhos = 0.9e3 # kg/m^3
rhom = 0.8e3 # kg/m^3
rhop = 1.1e3 # kg/m^3
VR   = 20.0  # m^3
Qf0  = 1.0   # m^3/min
Mm   = 100.0 # kg/kmol
CM0  = 0.0   # kmol/m^3
V0   = 10.0  # m^3

CMf = rhom/Mm           # kmol/m^3
a   = Qf0*CMf           # kmol/min
deltaV = (1/rhop - 1/rhom)*Mm  # m^3/kmol

class P:
    pass
p = P()
p.k = k; p.VR = VR; p.Qf0 = Qf0; p.V0 = V0
p.CMf = CMf; p.deltaV = deltaV

def stopping(time, p):
    return ((p.VR - p.V0) - (p.Qf0 + p.deltaV*p.Qf0*p.CMf)*time
            + p.deltaV*p.Qf0*p.CMf/p.k*(1 - np.exp(-p.k*time)))

x0 = (VR - V0)/Qf0
tstop = fsolve(lambda t: stopping(t[0], p), [x0])[0]

tfin = 50
npts = tfin + 1
t1 = np.linspace(0, tstop, npts)
V1 = V0 + (Qf0 + deltaV*Qf0*CMf)*t1 - deltaV*Qf0*CMf/k*(1-np.exp(-k*t1))
M1 = Qf0*CMf/k*(1-np.exp(-k*t1))
P1 = Qf0*CMf*(t1 - 1/k*(1-np.exp(-k*t1)))*Mm
t2 = np.linspace(tstop, tfin, npts)
Qf1 = Qf0*np.ones(npts)

V20 = V1[-1]; M20 = M1[-1]; P20 = P1[-1]
M2 = M20*np.exp(-k*(t2-tstop))
P2 = P20 + M20*(1-np.exp(-k*(t2-tstop)))*Mm
V2 = V20 + deltaV*M20*(1-np.exp(-k*(t2-tstop)))
Qf2 = np.zeros(npts)

V3 = VR*np.ones(npts)
M3 = M20*np.exp(-k*(deltaV*CMf+1)*(t2-tstop))
Qf3 = -deltaV*k*M3
P3 = P20 + M20*k/(k*(deltaV*CMf+1))*(1-np.exp(-k*(deltaV*CMf+1)*(t2-tstop)))*Mm

table = np.column_stack([
    np.concatenate([t1, t2]),
    np.concatenate([V1, V2]),
    np.concatenate([M1, M2])*Mm,
    np.concatenate([P1, P2]),
    np.concatenate([Qf1, Qf2]),
    np.concatenate([V1, V3]),
    np.concatenate([M1, M3])*Mm,
    np.concatenate([P1, P3]),
    np.concatenate([Qf1, Qf3])
])
save_ascii('filling.dat', table)