Figure 3.6:

Gibbs energy contours for the pentane reactions as a function of the two reaction extents.

Code for Figure 3.6

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
113
114
115
116
117
118
119
120
121
# Converted from contour_pentane.m
# Note: Uses fmincon which is not directly available in scipy.
# We use scipy.optimize.minimize with constraints instead.
import numpy as np
from scipy.optimize import fsolve, minimize
from misc import octave_save

class P:
    pass

p = P()
p.deltag1 = -3.72  # kcal/mol
p.deltag2 = -4.49  # kcal/mol
p.T = 400.0        # K
p.P = 2.5          # atm
p.R = 1.987        # cal/mol/K
p.K1 = np.exp(-p.deltag1*1000/(p.R*p.T))
p.K2 = np.exp(-p.deltag2*1000/(p.R*p.T))

def G(w, p):
    x, y = w
    z = 0.5 - x - y
    pp = 0.5 + z

    xlogx = 0.0 if x == 0 else abs(x)*np.log(abs(x/pp))
    ylogy = 0.0 if y == 0 else abs(y)*np.log(abs(y/pp))
    zlogz = 0.0 if z == 0 else abs(z)*np.log(abs(z/pp))

    return -p.Gval - (x*np.log(p.K1) + y*np.log(p.K2)) + pp*np.log(p.P) + \
           2*zlogz + xlogx + ylogy

def Gcontour(var, p, ext0, unk):
    w = ext0.copy()
    w[np.array(unk) - 1] = var
    return G(w, p)

def inflate(z, zlast, dels, p, ext0, unk):
    return np.array([
        Gcontour(z, p, ext0, np.array(unk)),
        np.linalg.norm(z - zlast) - dels**2
    ])

p.Gval = 0.0
w0 = np.array([0.2, 0.2])

# Minimize G subject to constraints
from scipy.optimize import minimize
cons = [
    {'type': 'ineq', 'fun': lambda w: w[0]},
    {'type': 'ineq', 'fun': lambda w: w[1]},
    {'type': 'ineq', 'fun': lambda w: 0.5 - w[0] - w[1]}
]
res = minimize(lambda w: G(w, p), w0, method='SLSQP',
               constraints=cons,
               bounds=[(0, 0.5), (0, 0.5)])
w = res.x
print(f'Optimal w = {w}')

Glevels = [0, -1, -2, -2.5, -2.5, -2.53, -2.55, -2.559]
nlev = len(Glevels)
alpha_vec = [0.05, 0.05, 0.075, 0.05, 0.05, 0.05, 0.04, 0.02]
x0vec = np.array([
    [1e-3, 1e-3, 1e-3, 0.133, 0.133, 0.133, 0.133, 0.133],
    [0.02, 0.15, 0.35, 0.3,   0.3,   0.3,   0.3,   0.3  ]
])
delx = np.array([
    [0.01, 0.01, 0.01, 0.01, -0.01, 0., 0., 0.],
    [0.,   0.,   0.,   0.,    0.,   0., 0., 0.]
])
sfinal = [9, 9, 9, 9, 9, 9, 9, 5]

ztables = []
for i in range(nlev):
    p.Gval = Glevels[i]
    ext0 = x0vec[:, i].copy()
    unk = [2]  # 1-indexed
    x0 = np.array([ext0[unk[0]-1]])

    try:
        x_sol = fsolve(lambda x: [Gcontour(x, p, ext0, unk)], x0, full_output=True)
        xsol = x_sol[0]
        ext = ext0.copy()
        ext[np.array(unk) - 1] = xsol
        zlast = ext.copy()
    except:
        zlast = ext0.copy()
        ext = ext0.copy()

    dels = alpha_vec[i]
    svec = np.arange(dels, sfinal[i]+dels, dels)
    narc = len(svec)
    z0 = ext + delx[:, i]
    unk2 = [1, 2]

    pts = [ext.copy()]
    for j in range(narc):
        try:
            z_sol = fsolve(lambda z: inflate(z, zlast, dels, p, ext0, unk2), z0)
            if (z_sol[0] > 0 and z_sol[1] > 0 and
                    0.5-(z_sol[0]+z_sol[1]) > 0):
                pts.append(z_sol.copy())
                z0 = 2*z_sol - zlast
                zlast = z_sol.copy()
            else:
                break
        except:
            break

    if len(pts) > 0:
        ztables.append(np.array(pts))
    else:
        ztables.append(np.array([[np.nan, np.nan]]))

feasbound = np.array([[0., 0.5], [0.5, 0.]])

# Save: each ztable as a named block, then feasbound
pairs = []
for i, zt in enumerate(ztables):
    pairs.append((f'ztable_{i+1}', zt))
pairs.append(('feasbound', feasbound))
octave_save('contour_pentane.dat', *pairs)