Figure 7.9:

Effectiveness factor versus Thiele modulus in a spherical pellet; reaction orders greater than unity.

Code for Figure 7.9

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
# Converted from etaordern_large.m - eta vs phi for large orders (n>=1)
import numpy as np
from scipy.integrate import solve_ivp
from scipy.special import iv as besseli
from misc import octave_save

ordervec = [1., 2., 5., 10.]
p = dict(q=2., thiele=10., k=2.)
c0  = 0.01
der = 1e-8
start = 1.

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

def ivp(t, x, order, p):
    xdot = np.zeros(2)
    xdot[0] = x[1]
    if t == 0:
        xdot[1] = p['k'] * x[0]**order
    else:
        xdot[1] = -p['q']/t * x[1] + p['k'] * x[0]**order
    return xdot

def make_g(order, p):
    def g(t, x):
        val = np.sqrt((order+1)/2 * x[0]**(order-1) * p['k']) * t/(p['q']+1)
        return val - p['thiele']
    g.terminal  = True
    g.direction = 0
    return g

table_list = []
no = len(ordervec)

for j in range(no):
    order = ordervec[j]
    if order < 1:
        # Stage 1
        tsteps = np.concatenate([[0.], np.logspace(-2, 5, 200)])
        x0     = np.array([c0, 0.])
        sol    = solve_ivp(lambda t, x: ivp(t, x, order, p),
                           [tsteps[0], tsteps[-1]], x0,
                           method='Radau', t_eval=tsteps,
                           rtol=eps_sq, atol=1e-18)
        ts = sol.t; xs = sol.y.T
        phi = np.sqrt((order+1)/2 * xs[:,0]**(order-1) * p['k']) * ts/(p['q']+1)
        eta = xs[:,1] / (p['k'] * ts * xs[:,0]**order) * (p['q']+1)
        phi[0] = 0.; eta[0] = 1.
        # Stage 2
        tsteps2 = np.concatenate([np.linspace(start, start+2, 30),
                                   np.logspace(np.log10(start+2), np.log10(40.), 20)])
        x0 = np.array([0., der])
        sol2 = solve_ivp(lambda t, x: ivp(t, x, order, p),
                         [tsteps2[0], tsteps2[-1]], x0,
                         method='Radau', t_eval=tsteps2,
                         rtol=eps_sq, atol=1e-18)
        ts2 = sol2.t; xs2 = sol2.y.T
        phidry = np.sqrt((order+1)/2 * xs2[:,0]**(order-1) * p['k']) * (ts2+start)/(p['q']+1)
        etadry = xs2[:,1] / (p['k'] * (ts2+start) * xs2[:,0]**order) * (p['q']+1)
        tab = np.vstack([np.column_stack([phi, eta]),
                         np.column_stack([phidry[1:][::-1], etadry[1:][::-1]])])
    elif order == 1:
        phi_v = np.logspace(-2, 2, 250)
        if p['q'] == 0:
            eta_v = np.tanh(phi_v) / phi_v
        elif p['q'] == 1:
            eta_v = besseli(1, phi_v) / (phi_v * besseli(0, phi_v))
        else:  # sphere
            eta_v = (1./phi_v) * (1./np.tanh(3*phi_v) - 1./(3*phi_v))
        tab = np.column_stack([phi_v, eta_v])
    else:
        x0_  = np.array([c0, 0.])
        t_out = np.concatenate([[0.], np.logspace(-1, 10, 2000)])
        ev    = make_g(order, p)
        sol   = solve_ivp(lambda t, x: ivp(t, x, order, p),
                          [t_out[0], t_out[-1]], x0_,
                          method='Radau', t_eval=t_out,
                          events=ev, rtol=eps_sq, atol=eps_sq)
        ts = sol.t; xs = sol.y.T
        phi = np.sqrt((order+1)/2 * xs[:,0]**(order-1) * p['k']) * ts/(p['q']+1)
        eta = xs[:,1] / (p['k'] * ts * xs[:,0]**order) * (p['q']+1)
        tab = np.column_stack([phi[1:], eta[1:]])
    table_list.append(tab)

octave_save('etaordern_large.dat',
            *[('table%d' % j, table_list[j]) for j in range(no)])