Figure 7.20:

Dimensionless concentration versus radius for the nonisothermal spherical pellet: lower (A), unstable middle (B), and upper (C) steady states.

Code for Figure 7.20

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
# Converted from whmss.m - Weisz-Hicks multiple steady-state pellet profiles
import numpy as np
from scipy.integrate import solve_ivp, solve_bvp, quad
from scipy.optimize import brentq
from misc import octave_save

Gamma    = 30.
beta     = 0.4
Phiscale = 0.01
tol      = 1e-6

def integrand(c, beta, Gamma):
    return c * np.exp(Gamma*beta*(1-c) / (1 + beta*(1-c)))

value, _ = quad(lambda c: integrand(c, beta, Gamma), 0., 1., limit=200, epsabs=tol)
intercept = np.sqrt(2*value)
Phi = intercept * Phiscale

eps_sq = np.sqrt(np.finfo(float).eps)
rout = np.linspace(0., 3., 100)

def ep_fn(ca):
    denom = 1 + beta*(1 - ca)
    return np.exp(np.clip(Gamma*beta*(1-ca) / np.where(denom > 1e-10, denom, 1e-10), -500., 500.))

def pelletode(r, y):
    ca, dcadr = y[0], y[1]
    ep = np.vectorize(ep_fn)(ca)
    d2 = np.where(np.abs(r) < 1e-10,
                  Phi**2/3. * ca * ep,
                  -2./r * dcadr + Phi**2 * ca * ep)
    return np.vstack([dcadr, d2])

def bc(ya, yb):
    return np.array([ya[1], yb[0] - 1.])

# Shooting to find initial guess for each branch
def pelletode_ivp(r, x):
    ca, dcadr = x[0], x[1]
    denom = 1 + beta*(1 - ca)
    ep = np.exp(Gamma*beta*(1-ca)/denom) if denom > 1e-10 else 0.
    if r < 1e-10:
        return [dcadr, Phi**2/3. * ca * ep]
    return [dcadr, -2./r * dcadr + Phi**2 * ca * ep]

def shoot(c0_val):
    s = solve_ivp(pelletode_ivp, [0., 3.], [c0_val, 0.],
                  method='Radau', rtol=eps_sq, atol=eps_sq)
    return float(s.y[0, -1]) - 1.

brackets = [(1e-11, 1e-10), (0.5, 0.7), (0.90, 0.99)]

results = []
for i, (lo, hi) in enumerate(brackets):
    c0 = brentq(shoot, lo, hi, xtol=1e-12, rtol=1e-12)
    # get shooting profile on fine mesh as initial guess for solve_bvp
    r_dense = np.linspace(1e-8, 3., 500)
    sol_shoot = solve_ivp(pelletode_ivp, [r_dense[0], r_dense[-1]], [c0, 0.],
                          method='Radau', t_eval=r_dense,
                          rtol=1e-10, atol=1e-12)
    y_guess = sol_shoot.y  # shape (2, N)
    sol = solve_bvp(pelletode, bc, r_dense, y_guess, tol=1e-6, max_nodes=50000)
    if not sol.success:
        print(f'Warning: solve_bvp branch {i} did not converge: {sol.message}')
    conc = sol.sol(rout)[0]
    temp = beta * (1. - conc)
    results.append(np.column_stack([rout, conc, temp]))
    print(f'Branch {i}: c(0)={conc[0]:.6f}, c(R)={conc[-1]:.8f}')

octave_save('whmss.dat',
            ('results0', results[0]),
            ('results1', results[1]),
            ('results2', results[2]))