Figure 8.33:

Conversion of reactant for single, ideal CSTR, and as a function of internal flowrate in a 2-CSTR mixing model.

Code for Figure 8.33

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
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
# Converted from selectivity.m - selectivity in two-CSTR network
import numpy as np
from scipy.integrate import solve_ivp
from scipy.optimize import fsolve
from misc import octave_save

eps_sq = np.sqrt(np.finfo(float).eps)

k1  = 1.;  k2  = 2.
theta1 = 1.;  theta2 = 2.
caf = 1.;  cbf = 1.
rhotest = 1.;  n = 2.;  alpha = 0.1
theta   = theta1 + theta2

def tworeac(x, rho):
    ca1, cb1, ca2, cb2 = x
    r11 = k1*ca1*cb1;  r21 = k2*ca1**n
    r12 = k1*ca2*cb2;  r22 = k2*ca2**n
    R   = np.array([-(r11+r21)*theta1, -r11*theta1,
                    -(r12+r22)*theta2, -r12*theta2])
    flow = np.array([alpha*caf - (alpha+rho)*ca1 + rho*ca2,
                     -(alpha+rho)*cb1 + rho*cb2,
                     -(1.+alpha)*ca2 + (alpha+rho)*ca1 - rho*ca2,
                     cbf + (alpha+rho)*cb1 - rho*cb2 - (1.+alpha)*cb2])
    return flow + R

def onereac(x):
    ca, cb = x
    r1 = k1*ca*cb;  r2 = k2*ca**n
    R  = np.array([-(r1+r2)*theta, -r1*theta])
    flow = np.array([alpha*caf - (1.+alpha)*ca,
                     cbf - (1.+alpha)*cb])
    return flow + R

npts   = 100
rho0   = 0.;  rhofin = 1.
rhovec = np.linspace(rhofin, rho0, npts)

y0 = np.array([0., cbf-caf])
y, _, info, _ = fsolve(onereac, y0, full_output=True)[:4]
conv1  = (alpha*caf - (1.+alpha)*y[0]) / (alpha*caf) * np.ones(npts)
yield1 = ((cbf - (1.+alpha)*y[1]) /
          (alpha*caf - (1.+alpha)*y[0])) * np.ones(npts)

conv2  = np.zeros(npts)
yield2 = np.zeros(npts)
x0 = np.array([0., cbf-caf, 0., cbf-caf])
for i, rho in enumerate(rhovec):
    x, _, info, _ = fsolve(lambda x: tworeac(x, rho), x0, full_output=True)[:4]
    if info != 1:
        print(f'warning: fsolve failed, i={i}')
    x0 = x.copy()
    conv2[i]  = (alpha*caf - (1.+alpha)*x[2]) / (alpha*caf)
    denom     = alpha*caf - (1.+alpha)*x[2]
    yield2[i] = (cbf - (1.+alpha)*x[3]) / denom if abs(denom) > 0 else 0.

table1 = np.column_stack([rhovec, yield1, yield2, conv1, conv2])

# step test results
def onestep(t, x):
    return [(alpha - (1.+alpha)*x[0]) / (theta1+theta2)]

nts   = 100
tfin  = 10.*theta
times = np.linspace(0., tfin, nts)
sol   = solve_ivp(onestep, [0., tfin], [0.], method='Radau',
                  t_eval=times, rtol=eps_sq, atol=eps_sq)
y_sol = sol.y[0]
impy  = np.array(onestep(0., [y_sol]))  # not used in table

rho = 0.
a   = np.array([alpha/theta1, 0.])
B   = np.array([[-(alpha+rho)/theta1, rho/theta1],
                [(alpha+rho)/theta2, -(1.+alpha+rho)/theta2]])
def twostep1(t, x):
    return a + B @ x
sol1 = solve_ivp(twostep1, [0., tfin], [0., 0.],
                 method='Radau', t_eval=times, rtol=eps_sq, atol=eps_sq)
x1   = sol1.y.T

rho = 1.
a   = np.array([alpha/theta1, 0.])
B   = np.array([[-(alpha+rho)/theta1, rho/theta1],
                [(alpha+rho)/theta2, -(1.+alpha+rho)/theta2]])
def twostep2(t, x):
    return a + B @ x
sol2 = solve_ivp(twostep2, [0., tfin], [0., 0.],
                 method='Radau', t_eval=times, rtol=eps_sq, atol=eps_sq)
x2   = sol2.y.T

# impulsive response approximation
# Octave: diag(a)*ones(2,nts) + B*x' gives (2 x nts); transpose to (nts x 2)
rho = 0.
a0  = np.array([alpha/theta1, 0.])
B0  = np.array([[-(alpha+0.)/theta1, 0./theta1],
                [(alpha+0.)/theta2, -(1.+alpha+0.)/theta2]])
impx1 = (a0[:, None] + B0 @ x1.T).T   # shape (nts, 2)
rho = 1.
a1  = np.array([alpha/theta1, 0.])
B1  = np.array([[-(alpha+1.)/theta1, 1./theta1],
                [(alpha+1.)/theta2, -(1.+alpha+1.)/theta2]])
impx2 = (a1[:, None] + B1 @ x2.T).T   # shape (nts, 2)

# Octave: [times y x1(:,2) x2(:,2) impy impx1(:,2) impx2(:,2)]
impy_vec = np.array([(alpha - (1.+alpha)*yv)/(theta1+theta2) for yv in y_sol])
table2   = np.column_stack([times, y_sol, x1[:,1], x2[:,1],
                             impy_vec, impx1[:,1], impx2[:,1]])

octave_save('selectivity.dat',
            ('table1', table1),
            ('table2', table2))