Figure 6.38:

Coolant temperature at reactor outlet versus temperature at reactor inlet; three steady-state solutions.

Code for Figure 6.38

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_shooting.m - Ammonia tubular reactor: shooting method for Taend
import numpy as np
from scipy.integrate import solve_ivp
from misc import save_ascii, ode15s

P     = 300.       # atm
R     = 82.05e-6   # (m^3 atm)/(mol K)
v0    = 0.05713    # m/sec
A     = 1.         # m^2
UA    = 0.5        # 1/m
delG0 = -4250.     # cal/mol ammonia
delH0 = -12000.    # cal/mol ammonia
Tstd  = 298.       # K
Rg    = 1.987      # cal/mol K
K0    = np.exp(-delG0/(Rg*Tstd))
deltaH = -2.342    # s K/mol
Ta0   = 323.       # K
len_  = np.linspace(0., 12., 100)

k0 = 7.794e11   # reverse rate constant
E  = 20000.     # K

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)

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]

nTs = 100
Tvec   = np.linspace(273., 1000., nTs)
Taend  = np.zeros(nTs)
Trend  = np.zeros(nTs)
xnhvec = np.zeros(nTs)

eps_sq = np.sqrt(np.finfo(float).eps)
for i, T in enumerate(Tvec):
    x0_i = [Nnh0, T, T]
    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, :]
    Taend[i]  = sol.y[2, -1]
    Trend[i]  = sol.y[1, -1]
    xnhvec[i] = sol.y[0, -1] / NT[-1]

table = np.column_stack([Tvec, Ta0*np.ones(nTs), Taend, Trend, xnhvec])
save_ascii('ammonia_shooting.dat', table)