Figure 7.10:

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

Figure 7.10

Code for Figure 7.10

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
112
113
114
115
116
117
118
119
120
# Converted from etaordern_small.m - eta vs phi for small 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 = [0., 0.25, 0.5, 0.75, 1.]
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]
    c = max(x[0], 0.)
    if t == 0:
        xdot[1] = p['k'] * c**order
    else:
        xdot[1] = -p['q']/t * x[1] + p['k'] * c**order
    return xdot

def ivp_shifted(t, x, order, p, r0):
    """ODE in shifted coords: t=0 corresponds to r=r0 (dead zone boundary)."""
    xdot = np.zeros(2)
    xdot[0] = x[1]
    c = max(x[0], 0.)
    r = t + r0
    xdot[1] = -p['q']/r * x[1] + p['k'] * c**order
    return xdot

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

table_list = []
nord = len(ordervec)

for j in range(nord):
    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: shifted coordinates, tau=0 at dead zone boundary (r=start=1).
        # Use asymptotic dead zone IC: c ~ A*tau^m, dc/dtau ~ m*A*tau^(m-1)
        # where A is from dominant balance of full ODE (incl. curvature term q/r):
        #   A^(1-n) = k / (m*(m-1+q)),  m = 2/(1-n).
        # Start at tau=eps_start chosen so phi_start slightly exceeds plot range.
        if order == 0:
            # Zeroth order: exact solution available, start from dead zone boundary
            # directly; source term k*c^0=k nonzero even at c=0, so [0,0] works.
            eps_start = 0.
            c_eps = 0.
            dc_eps = 0.
            t2max = 40.
            tsteps2 = np.concatenate([[eps_start], np.logspace(-2, np.log10(t2max), 300)])
        else:
            m = 2.0 / (1.0 - order)
            A_coef = (p['k'] / (m * (m - 1 + p['q'])))**(1.0 / (1.0 - order))
            # phi_crit ~ sqrt((n+1)/2 * A^(n-1) * k) / (q+1) = sqrt(m*(m-1+q)*(n+1)/2) / (q+1)
            phi_crit_est = np.sqrt((order+1)/2 * m * (m - 1 + p['q'])) / (p['q'] + 1)
            # Choose eps_start so phi_start = phi_crit*(1 + 1/eps_start) > 12
            phi_target = 12.0
            eps_start = phi_crit_est / (phi_target - phi_crit_est)
            c_eps = A_coef * eps_start**m
            dc_eps = m * A_coef * eps_start**(m - 1)
            t2max = 200.
            tsteps2 = np.concatenate([[eps_start],
                                      np.logspace(np.log10(eps_start * 1.05),
                                                  np.log10(t2max), 400)])
        x0_s2 = np.array([c_eps, dc_eps])
        sol2 = solve_ivp(lambda t, x: ivp_shifted(t, x, order, p, start),
                         [tsteps2[0], tsteps2[-1]], x0_s2,
                         method='Radau', t_eval=tsteps2,
                         rtol=eps_sq, atol=1e-25)
        ts2 = sol2.t; xs2 = np.asarray(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:
            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.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, eta])
    table_list.append(tab)

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