Figure 7.24:

Concentration profiles of products.

Figure 7.24

Code for Figure 7.24

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
# Converted from multi_log.m - multi-component pellet with log transformation
import numpy as np
from scipy.optimize import fsolve
from misc import colloc, octave_save

p = dict(
    Rg  = 8.314,
    P   = 1.013e5,
    T   = 550.,
    Rp  = 0.175,
)
p['cf_tot'] = p['P']/(p['Rg']*p['T'])*1e-6
p['cf_a']   = p['cf_tot']*0.02
p['cf_b']   = p['cf_tot']*0.03
p['cf_c']   = p['cf_tot']*5e-4

p['k_k10']  = 6.802e16*2.6e6*80./100*0.05/100
p['k_k20']  = 1.416e18*2.6e6*80./100*0.05/100
p['k_E1']   = 13108.
p['k_E2']   = 15109.
p['k_Ka0']  = 8.099e6
p['k_Kc0']  = 2.579e8
p['k_Ea']   = -409.
p['k_Ec']   = 191.

p['D_a']  = 0.0487; p['D_b'] = 0.0469; p['D_c'] = 0.0487
p['km_a'] = 0.4*9.76; p['km_b'] = 0.4*10.18; p['km_c'] = 0.4*9.76

npts = 40
R_col, A_col, B_col, Q_col = colloc(npts-2, 'left', 'right')
R_col = R_col * p['Rp']
A_col = A_col / p['Rp']
B_col = B_col / p['Rp']**2
Q_col = Q_col * p['Rp']

Aint = A_col[1:npts-1, :]
Bint = B_col[1:npts-1, :]
Rint = R_col[1:npts-1]

p['npts'] = npts
p['col_R'] = R_col; p['col_A'] = A_col; p['col_B'] = B_col
p['col_Aint'] = Aint; p['col_Bint'] = Bint; p['col_Rint'] = Rint

def pellet(x, p):
    n   = p['npts']
    za  = x[:n]; zb = x[n:2*n]; zc = x[2*n:3*n]
    ca  = np.exp(za); cb = np.exp(zb); cc = np.exp(zc)
    k1  = p['k_k10']*np.exp(-p['k_E1']/p['T'])
    k2  = p['k_k20']*np.exp(-p['k_E2']/p['T'])
    Ka  = p['k_Ka0']*np.exp(-p['k_Ea']/p['T'])
    Kc  = p['k_Kc0']*np.exp(-p['k_Ec']/p['T'])
    den = (1 + Ka*ca + Kc*cc)**2
    r1  = k1*ca*cb/den; r2 = k2*cc*cb/den
    Ra  = -r1; Rb = -0.5*r1 - 4.5*r2; Rc_ = -r2

    retval = np.zeros(3*n)
    # A
    retval[0]   = p['col_A'][0,:]@za
    retval[1:n-1] = (p['col_Bint']@za +
                     (p['col_Aint']@za) * (p['col_Aint']@za + 2./p['col_Rint']) +
                     Ra[1:n-1]/(p['D_a']*ca[1:n-1]))
    dzadr = p['col_A'][-1,:]@za
    retval[n-1] = p['D_a']*dzadr - p['km_a']*(p['cf_a']/ca[-1] - 1)
    # B
    retval[n]   = p['col_A'][0,:]@zb
    retval[n+1:2*n-1] = (p['col_Bint']@zb +
                          (p['col_Aint']@zb) * (p['col_Aint']@zb + 2./p['col_Rint']) +
                          Rb[1:n-1]/(p['D_b']*cb[1:n-1]))
    dzbdr = p['col_A'][-1,:]@zb
    retval[2*n-1] = p['D_b']*dzbdr - p['km_b']*(p['cf_b']/cb[-1] - 1)
    # C
    retval[2*n]   = p['col_A'][0,:]@zc
    retval[2*n+1:3*n-1] = (p['col_Bint']@zc +
                            (p['col_Aint']@zc) * (p['col_Aint']@zc + 2./p['col_Rint']) +
                            Rc_[1:n-1]/(p['D_c']*cc[1:n-1]))
    dzcdr = p['col_A'][-1,:]@zc
    retval[3*n-1] = p['D_c']*dzcdr - p['km_c']*(p['cf_c']/cc[-1] - 1)
    return retval

# Initial guess (linear from small value to bulk)
zain  = np.log(1e-9*p['cf_a']); zaout = np.log(0.75*p['cf_a'])
za0   = zain + R_col/p['Rp']*(zaout - zain)
zbin  = np.log(0.75*p['cf_b']); zbout = np.log(p['cf_b'])
zb0   = zbin + R_col/p['Rp']*(zbout - zbin)
zcin  = np.log(1e-6*p['cf_c']); zcout = np.log(p['cf_c'])
zc0   = zcin + R_col/p['Rp']*(zcout - zcin)
z0    = np.concatenate([za0, zb0, zc0])

eps_sq = np.sqrt(np.finfo(float).eps)
z = fsolve(lambda x: pellet(x, p), z0, full_output=False)

za = z[:npts]; zb = z[npts:2*npts]; zc = z[2*npts:3*npts]
ca = np.exp(za); cb = np.exp(zb); cc = np.exp(zc)
dzadr = A_col[-1,:]@za; dzbdr = A_col[-1,:]@zb; dzcdr = A_col[-1,:]@zc
dcadr = dzadr*ca; dcbdr = dzbdr*cb; dccdr = dzcdr*cc

# Compute product concentrations D and E
p['km_d'] = p['km_a']; p['km_e'] = p['km_b']
p['D_d']  = p['D_a'];  p['D_e']  = p['D_b']
p['cf_d'] = 0.;        p['cf_e'] = 0.

dcddr = (-p['D_a']*dcadr[-1] - 3*p['D_c']*dccdr[-1]) / p['D_d']
dcedr = (-3*p['D_c']*dccdr[-1]) / p['D_e']
cdR   = p['cf_d'] - p['D_d']*dcddr/p['km_d']
ceR   = p['cf_e'] - p['D_e']*dcedr/p['km_e']
concd = (p['D_a']*(ca[-1]-ca) + 3*p['D_c']*(cc[-1]-cc))/p['D_d'] + cdR
ce    = 3*p['D_c']*(cc[-1]-cc)/p['D_e'] + ceR

table   = np.column_stack([R_col, ca, cb, cc, concd, ce, dcadr, dcbdr, dccdr])
table_2 = np.array([[p['Rp'], p['cf_a'], p['cf_b'], p['cf_c'], p['cf_d'], p['cf_e']]])

octave_save('multi_log.dat', ('table', table), ('table_2', table_2))