Figure 8.19:

Dimensionless effluent concentration versus dimensionless rate constant for second-order reaction.

Code for Figure 8.19

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
# Converted from mmseg2cstr.m - maximum mixedness and segregated models for 2 CSTRs
# ode15i for mm model → reformulate as standard ODE using a stiff solver
# The implicit equation is: dcdz*z^2 + 4*(1-z)/(2-z)*(c-1) + K*c*c = 0
# → dcdz = -(4*(1-z)/(2-z)*(c-1) + K*c*c) / z^2
import numpy as np
from scipy.integrate import solve_ivp
from misc import save_ascii

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

Kmin = 1e-2;  Kmax = 1e3;  nks = 40
Kseq = np.logspace(np.log10(Kmin), np.log10(Kmax), nks)

def mm_ode(z, c, K):
    if z == 0. or z < 1e-12:
        return [0.]
    return [-(4.*(1.-z)/(2.-z)*(c[0]-1.) + K*c[0]*c[0]) / (z*z)]

def seg_ode(z, x, K):
    c, ctot = x
    if z < 1.:
        t = z/(1.-z)
        dcdt  = -K*c*c / (1.-z)**2
        dctot = 4.*t*np.exp(-2.*t)/(1.-z)**2 * c
    else:
        dcdt = 0.;  dctot = 0.
    return [dcdt, dctot]

npts = 50
z    = np.linspace(0., 1., npts)
c0_seg = np.array([1., 0.])

cmmfin  = np.zeros(nks)
csegfin = np.zeros(nks)

for i, K in enumerate(Kseq):
    cbar = (-1. + np.sqrt(1.+2.*K)) / K

    sol_mm = solve_ivp(lambda z, c: mm_ode(z, c, K),
                       [0., 1.], [cbar], method='Radau',
                       t_eval=z, rtol=eps_sq, atol=eps_sq)
    cmmfin[i] = sol_mm.y[0,-1]

    sol_seg = solve_ivp(lambda z, x: seg_ode(z, x, K),
                        [0., 1.], c0_seg, method='Radau',
                        t_eval=z, rtol=eps_sq, atol=eps_sq)
    csegfin[i] = sol_seg.y[1,-1]

# 2 ideal CSTRs
cbar1 = (-1. + np.sqrt(1.+2.*Kseq)) / Kseq
cbar2 = (-1. + np.sqrt(1.+2.*Kseq*cbar1)) / Kseq
table = np.column_stack([Kseq, cmmfin, csegfin, cbar2])
save_ascii('mmseg2cstr.dat', table)