Figure A.8

Figure A.8:

Magnified view of Figure~{\ref {fig:fb2colloc}}.

Code for Figure A.8

pbrsolve.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 # Converted from pbrsolve.m - packed-bed reactor helper function
import numpy as np
from scipy.integrate import solve_ivp

def pbr(t, x, par):
    Na  = x[0]
    ca  = par['P']/(par['Rg']*par['T']) * Na/(2*par['Nafin'])
    Phi = par['Rp']/3 * np.sqrt((par['n']+1)/2 * par['k'] * ca / par['Da'])
    return [-par['rhob']/par['rhop'] * par['eta'](Phi) * par['k'] * ca**par['n']]

def pbrsolve(par):
    vtotal = par['vfinal'] * np.linspace(0., 1., par['nvs'])
    vsteps = vtotal
    x0     = np.array([par['Nafin']])
    eps_sq = np.sqrt(np.finfo(float).eps)

    def stop_event(t, x):
        return x[0] - (1 - par['xa'])*par['Nafin']
    stop_event.terminal  = True
    stop_event.direction = 0

    sol = solve_ivp(lambda t, x: pbr(t, x, par),
                    [vsteps[0], vsteps[-1]], x0,
                    method='Radau', t_eval=vsteps,
                    events=stop_event,
                    rtol=eps_sq, atol=eps_sq)

    if len(sol.t) == par['nvs']:
        print('hey, did not reach final conversion, increase stopping time')

    return sol.t, sol.y.T

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 # Converted from fb2colloc.m - fixed-bed reactor with pellet collocation
# [depends] ~/ch7/figures/functions/pbrsolve.py
# [makes] fb2colloc.dat
import sys
import numpy as np
from scipy.optimize import fsolve
from scipy.integrate import solve_ivp
sys.path.insert(0, '../../ch7/figures/functions')
from misc import save_ascii, octave_save, colloc
from pbrsolve import pbrsolve

xa    = 0.75
Rg    = 8.314e0   # J/K mol
P_pa  = 4*1.013e5 # N/m^2
T     = 550.0     # K
ctf   = P_pa/(Rg*T)*1e-6  # mol/cm^3
Rp    = 0.45      # cm
a     = Rp/3.0
Nafin = 10.0      # mol/sec
caf   = ctf*0.5
Qf    = Nafin/caf
k1    = 2.25e5    # cm^3/mol/s
Da    = 0.008     # cm^2/s
rhob  = 0.60; rhop = 0.68
epsb  = 1.0 - rhob/rhop

npts  = 25
R0, A, B, Q = colloc(npts - 2, 'left', 'right')
R = R0*Rp; A_m = A/Rp; B_m = B/Rp**2; Q_m = Q*Rp

def pellet(ca_vec, caf_loc):
    r1  = k1 * ca_vec**2
    Ra  = -r1
    with np.errstate(divide='ignore', invalid='ignore'):
        res = B_m @ ca_vec + 2.0*A_m @ ca_vec/R + Ra/Da
    res[0]  = A_m[0, :] @ ca_vec   # symmetry BC overwrites r=0 row
    res[-1] = caf_loc - ca_vec[-1]
    return res

# Initial pellet profile (at bed inlet)
ca0  = np.logspace(np.log10(caf)-2, np.log10(caf), npts)
ca_in = fsolve(lambda x: pellet(x, caf), ca0, full_output=False)

def bed_rhs(V, y):
    Naf     = y[0]
    cpellet = y[1:]
    caf_loc = Naf / Qf
    dcadr   = A_m[-1, :] @ cpellet
    r1p     = Da/a * dcadr
    RA      = -(1.0 - epsb)*r1p
    dNaf_dV = RA
    pelletres = pellet(cpellet, caf_loc)
    dydt    = np.concatenate([[dNaf_dV], pelletres])
    return dydt

def event_conv(V, y):
    return xa - (1.0 - y[0]/Nafin)
event_conv.terminal  = True
event_conv.direction = 0

nvs    = 100; vfinal = 4e5
vsteps = np.linspace(0, vfinal, nvs)
y0     = np.concatenate([[Nafin], ca_in])

sol = solve_ivp(bed_rhs, [0, vfinal], y0, t_eval=vsteps,
                events=event_conv, method='Radau',
                rtol=1e-10, atol=1.5e-8)

if sol.status != 1:
    print('# Warning: did not reach target conversion, increase vfinal')

nout  = len(sol.t)
Naf   = sol.y[0]
vplot = sol.t / 1000.0   # lit
VR    = vplot[-1]
xa_act = (Nafin - Naf[-1]) / Nafin
tableex = np.column_stack([vplot, Naf])

# Pellet profiles at inlet and outlet
ca_profiles = sol.y[1:, :]
table2 = np.column_stack([R, ca_profiles[:, 0], ca_profiles[:, -1]])

Naout  = (1.0 - xa_act)*Nafin
Natop  = (1.0 - xa_act + 0.10)*Nafin

# Approximate solutions using pbrsolve from ch7
par = dict(
    k=k1, Nafin=Nafin, T=T, rhop=rhop, rhob=rhob,
    Da=Da, Rp=Rp, Rg=82.06, P=4.0, n=2,
    xa=xa_act, nvs=nvs, vfinal=vfinal,
    eta=lambda Phi: 1.0/Phi * (1.0/np.tanh(3.0*Phi) - 1.0/(3.0*Phi))
)
vap1, xap1 = pbrsolve(par)

par['eta'] = lambda Phi: 1.0/Phi
vap2, xap2 = pbrsolve(par)

vap1   = vap1 / 1000.0   # cm^3 -> L
VRap1  = vap1[-1]
vap2   = vap2 / 1000.0
VRap2  = vap2[-1]
tableap1 = np.column_stack([vap1, xap1[:, 0]])
tableap2 = np.column_stack([vap2, xap2[:, 0]])

dashedlines = np.array([
    [0.0,      Naout, VR,    0.0,  VRap1, 0.0,  VRap2, 0.0],
    [1.1*VR,   Naout, VR,    Natop, VRap1, Natop, VRap2, Natop]
])

octave_save('fb2colloc.dat',
            ('tableex',    tableex),
            ('tableap1',   tableap1),
            ('tableap2',   tableap2),
            ('dashedlines', dashedlines))