Figure 6.40:

Ammonia mole fraction versus reactor length.

Code for Figure 6.40

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
# Converted from ammonia_profile.m - Ammonia tubular reactor profiles (3 steady states)
import numpy as np
from scipy.integrate import solve_ivp
from scipy.optimize import fsolve
from misc import save_ascii

P     = 300.
R     = 82.05e-6
v0    = 0.05713
A     = 1.
UA    = 0.5
delG0 = -4250.
delH0 = -12000.
Tstd  = 298.
Rg    = 1.987
K0    = np.exp(-delG0/(Rg*Tstd))
deltaH = -2.342
Ta0   = 323.
len_  = np.linspace(0., 12., 250)
k0    = 7.794e11
E     = 20000.

Q0   = v0 * A
Nt0  = Q0 * P / (R * Ta0)
xn0  = 0.985*0.25;  xh0 = 0.985*0.75;  xnh0 = 0.015
Nnh0 = xnh0*Nt0;  Nh0 = xh0*Nt0;  Nn0 = xn0*Nt0

p = dict(P=P, R=R, A=A, UA=UA, delH0=delH0, E=E, k0=k0, K0=K0, Tstd=Tstd,
         Rg=Rg, deltaH=deltaH, Ta0=Ta0, Nnh0=Nnh0, Nh0=Nh0, Nn0=Nn0)

eps_sq = np.sqrt(np.finfo(float).eps)

def ivpjbr(t, x, p):
    Nnh = x[0]; T = x[1]; Ta = x[2]
    RT   = p['R']*T
    Nn   = p['Nn0'] - 0.5*(Nnh - p['Nnh0'])
    Nh   = p['Nh0'] - 1.5*(Nnh - p['Nnh0'])
    Nt   = Nn + Nh + Nnh
    pn   = (Nn/Nt)*p['P'];  ph = (Nh/Nt)*p['P'];  pnh = (Nnh/Nt)*p['P']
    qadd = (Ta - T)*p['UA']
    k    = p['k0']*np.exp(-p['E']/T)
    K    = p['K0']*np.exp(-p['delH0']/p['Rg']*(1./T - 1./p['Tstd']))
    rxrate = k/RT * (K**2 * pn*ph**1.5/pnh - pnh/ph**1.5)
    return [p['A']*2*rxrate, qadd - rxrate*p['deltaH'], qadd]

def feed_condition(T0, p):
    x0_i = [Nnh0, T0[0], T0[0]]
    sol  = solve_ivp(lambda t, x: ivpjbr(t, x, p), [len_[0], len_[-1]], x0_i,
                     method='Radau', t_eval=len_, rtol=eps_sq, atol=eps_sq)
    return [p['Ta0'] - sol.y[2, -1]]

T0guesses = [300., 600., 900.]
# table1: [len, x1_cols, xnh1, x2_cols, xnh2, x3_cols, xnh3]
table1 = len_.reshape(-1, 1)

for T0g in T0guesses:
    T0sol = fsolve(lambda T0: feed_condition(T0, p), [T0g])
    x0_i  = [Nnh0, T0sol[0], T0sol[0]]
    sol   = solve_ivp(lambda t, x: ivpjbr(t, x, p), [len_[0], len_[-1]], x0_i,
                      method='Radau', t_eval=len_, rtol=eps_sq, atol=eps_sq)
    NT    = Nn0 + Nh0 + 2*Nnh0 - sol.y[0, :]
    xnh   = sol.y[0, :] / NT
    table1 = np.column_stack([table1, sol.y.T, xnh])

save_ascii('ammonia_profile.dat', table1)