Figure 8.12:

Start-up of the tubular reactor; c_A(t,z) versus z for various times, 0\leq t\leq 2.5~min, \Delta t=0.25~min.

Code for Figure 8.12

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
# Converted from dispersedpfrtran.m - transient dispersed PFR (ode15i DAE)
# Uses SUNDIALS IDA via scikits.odes: boundary rows are algebraic constraints.
import numpy as np
from scikits.odes import dae
from misc import colloc, save_ascii

k     = 1./2.
v     = 1./2.
D     = 0.01
caf   = 1.
order = 2
length = 1.
ncolpt = 50

z_col, A_col, B_col, Q_col = colloc(ncolpt-2, 'left', 'right')

def dae_residual(t, y, ydot, res):
    first  = A_col @ y
    Ra     = -2.*k * (y**order)
    second = B_col @ y
    # interior: res = ydot - (-v*first + D*second + Ra)
    res[:] = ydot - (-v*first + D*second + Ra)
    # Danckwerts BCs (algebraic — overwrite residual rows)
    res[0]          = v*y[0] - D*first[0] - v*caf
    res[ncolpt - 1] = first[ncolpt - 1]

tsteps = np.linspace(0., 2.5, 11)
y0     = np.zeros(ncolpt)
ydot0  = np.zeros(ncolpt)

algvar = np.ones(ncolpt)
algvar[0]          = 0.   # algebraic
algvar[ncolpt - 1] = 0.   # algebraic

solver = dae('ida', dae_residual,
             algebraic_vars_idx=np.where(algvar == 0)[0],
             compute_initcond='yp0',
             old_api=False,
             rtol=1e-6, atol=1e-6)

solution = solver.solve(tsteps, y0, ydot0)

sol_y = solution.values.y   # shape (11, ncolpt)
table = np.column_stack([z_col, sol_y.T])
save_ascii('dispersedpfrtran.dat', table)