Figure A.5:

Dimensionless concentration versus dimensionless radial position for different numbers of collocation points.

Code for Figure A.5

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
# Converted from pelletcolloc.m - pellet collocation profiles vs npts
# [makes] pelletcolloc.dat
import numpy as np
from scipy.optimize import fsolve
from misc import save_ascii, colloc

Phi = 10.0; n = 1
nint = 100
xint = np.linspace(0, 3, nint)

# Polynomial interpolation (Lagrange)
def polinterp(xcol, xint_arr):
    ncol = len(xcol); nint = len(xint_arr)
    W = np.zeros((nint, ncol))
    for i in range(ncol):
        d = 1.0
        for j in range(ncol):
            if j != i:
                d *= (xcol[i] - xcol[j])
        for k in range(nint):
            n_val = 1.0
            for j in range(ncol):
                if j != i:
                    n_val *= (xint_arr[k] - xcol[j])
            W[k, i] = n_val / d
    return W

# Analytical solution
can      = np.where(xint != 0, 3.0/xint * np.sinh(Phi*xint)/np.sinh(3*Phi), 3*Phi/np.sinh(3*Phi))
etaan    = 1.0/Phi * (1.0/np.tanh(3*Phi) - 1.0/(3*Phi))

colvec = [5, 10, 30, 50]
ncol   = len(colvec)
csave  = np.zeros((nint, ncol))
cerr   = np.zeros((nint, ncol))

for i, npts in enumerate(colvec):
    R0, A, B, Q = colloc(npts - 2, 'left', 'right')
    R = R0*3.0; A = A/3.0; B = B/9.0; Q = Q*3.0

    W = polinterp(R, xint)

    def pellet(c):
        res = B @ c + 2.0*(A @ c)/R - 2.0/(n+1)*Phi**2*(c**n)
        res[0]  = A[0, :] @ c
        res[-1] = 1.0 - c[-1]
        return res

    c0 = 0.5*np.linspace(0, 1, npts)
    c  = fsolve(pellet, c0, full_output=False)
    cint = W @ c
    csave[:, i] = cint
    cerr[:, i]  = np.abs(cint - can) / can

table = np.column_stack([xint, csave, cerr])
save_ascii('pelletcolloc.dat', table)