08 - Chemical Equilibrium & Gibbs Free Energy
Whether a reaction proceeds — and how far — is governed by the Gibbs free energy. From the standardized enthalpy and absolute entropy pyglenn provides, we can build the standard molar Gibbs energy
the standard reaction Gibbs energy \(\Delta G^\circ(T)\), and the equilibrium constant
We apply this to the industrially important water-gas shift reaction and to high-temperature dissociation, and check consistency with the van’t Hoff equation.
[1]:
from pyglenn import ThermochemicalCalculator, R
print("Universal gas constant R =", R, "J/(mol.K)")
Universal gas constant R = 8.31446261815324 J/(mol.K)
[2]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
plt.rcParams["figure.figsize"] = (8, 4.5)
plt.rcParams["axes.grid"] = True
pd.set_option("display.float_format", lambda v: f"{v:,.3f}")
1. Gibbs energy and reaction helpers
gibbs returns \(G^\circ(T)\) for one species; reaction_props returns \(\Delta H^\circ\), \(\Delta S^\circ\), \(\Delta G^\circ\) and \(K\) for a reaction given as {species: coefficient} dicts.
[3]:
def gibbs(calc, name, T):
p = calc.calculate_properties(calc.get_available_species(name, exact_match=True)[0].id, T)
return p.h_relative - T * p.s # J/mol
def reaction_props(calc, reactants, products, T):
def acc(mix, key, scale=1.0):
tot = 0.0
for name, nu in mix.items():
p = calc.calculate_properties(calc.get_available_species(name, exact_match=True)[0].id, T)
tot += nu * (getattr(p, key) if key != "g" else (p.h_relative - T * p.s))
return tot
dH = acc(products, "h_relative") - acc(reactants, "h_relative")
dS = acc(products, "s") - acc(reactants, "s")
dG = acc(products, "g") - acc(reactants, "g")
K = np.exp(-dG / (R * T))
return dH, dS, dG, K
2. The water-gas shift reaction
Mildly exothermic, so by Le Chatelier’s principle it is favoured at low temperature and \(K\) falls as temperature rises, crossing \(K=1\) near ~1100 K.
[4]:
reac = {"CO": 1, "H2O": 1}
prod = {"CO2": 1, "H2": 1}
rows = []
with ThermochemicalCalculator() as calc:
for T in [500, 700, 900, 1100, 1300]:
dH, dS, dG, K = reaction_props(calc, reac, prod, T)
rows.append({"T [K]": T, "dH [kJ/mol]": dH/1000, "dS [J/mol/K]": dS,
"dG [kJ/mol]": dG/1000, "K": K})
print(pd.DataFrame(rows).set_index("T [K]").to_string())
dH [kJ/mol] dS [J/mol/K] dG [kJ/mol] K
T [K]
500 -39.815 -38.728 -20.451 136.910
700 -37.858 -35.449 -13.043 9.403
900 -35.779 -32.837 -6.226 2.298
1100 -33.773 -30.822 0.132 0.986
1300 -31.910 -29.264 6.134 0.567
[5]:
Tgrid = np.linspace(400, 1500, 80)
with ThermochemicalCalculator() as calc:
K = np.array([reaction_props(calc, reac, prod, T)[3] for T in Tgrid])
fig, ax = plt.subplots()
ax.semilogy(Tgrid, K)
ax.axhline(1.0, color="0.5", ls="--")
ax.set_xlabel("Temperature [K]")
ax.set_ylabel("equilibrium constant $K$")
ax.set_title("Water-gas shift: $K$ decreases with temperature")
plt.show()
# temperature where K = 1 (dG = 0), by sign change
sign = np.sign(np.log(K))
idx = np.where(np.diff(sign) != 0)[0]
if len(idx):
T_cross = np.interp(0.0, [np.log(K[idx[0]]), np.log(K[idx[0]+1])],
[Tgrid[idx[0]], Tgrid[idx[0]+1]])
print(f"K = 1 (dG = 0) near T = {T_cross:.0f} K")
K = 1 (dG = 0) near T = 1096 K
3. The van’t Hoff check
The van’t Hoff equation relates the slope of \(\ln K\) against \(1/T\) to the reaction enthalpy:
Plotting \(\ln K\) vs. \(1/T\) gives a nearly straight line; its local slope recovers \(\Delta H^\circ\) computed directly — a nice internal-consistency check.
[6]:
invT = 1.0 / Tgrid
lnK = np.log(K)
# local slope near the middle of the range
i = len(Tgrid) // 2
slope = (lnK[i+1] - lnK[i-1]) / (invT[i+1] - invT[i-1])
dH_vanthoff = -slope * R
with ThermochemicalCalculator() as calc:
dH_direct = reaction_props(calc, reac, prod, Tgrid[i])[0]
# Linear regression to visualise the van't Hoff slope
fit = np.polyfit(invT, lnK, 1)
lnK_fit = np.polyval(fit, invT)
dH_fit = -fit[0] * R # slope = -dH/R
fig, ax = plt.subplots()
ax.plot(invT, lnK, label="data")
ax.plot(invT, lnK_fit, "--", color="C1",
label=f"linear fit (slope = {fit[0]:.0f} K)")
ax.set_xlabel("1 / T [1/K]")
ax.set_ylabel(r"$\ln K$")
ax.set_title("van't Hoff plot for the water-gas shift")
ax.legend()
plt.show()
print(f"dH from van't Hoff local slope : {dH_vanthoff/1000:7.2f} kJ/mol")
print(f"dH from linear regression : {dH_fit/1000:7.2f} kJ/mol")
print(f"dH computed directly : {dH_direct/1000:7.2f} kJ/mol")
dH from van't Hoff local slope : -35.20 kJ/mol
dH from linear regression : -37.24 kJ/mol
dH computed directly : -35.20 kJ/mol
4. High-temperature dissociation
Diatomic molecules split into atoms at extreme temperatures — a strongly endothermic process (\(\Delta H^\circ \approx +945\) kJ/mol for N₂, \(+498\) kJ/mol for O₂). Their equilibrium constants are astronomically small at moderate temperature and only become appreciable in the thousands of kelvin, which is why the complete-combustion flame temperatures of notebook 06 are over-predictions.
[7]:
dissociations = {
"N2 <-> 2 N": ({"N2": 1}, {"N": 2}),
"O2 <-> 2 O": ({"O2": 1}, {"O": 2}),
}
Tgrid = np.linspace(1000, 6000, 80)
fig, ax = plt.subplots()
with ThermochemicalCalculator() as calc:
for label, (r, p) in dissociations.items():
logK = [np.log10(reaction_props(calc, r, p, T)[3]) for T in Tgrid]
ax.plot(Tgrid, logK, label=label)
dH_N2 = reaction_props(calc, {"N2": 1}, {"N": 2}, 298.15)[0]
ax.axhline(0.0, color="0.6", ls="--")
ax.set_xlabel("Temperature [K]")
ax.set_ylabel(r"$\log_{10} K$")
ax.set_title("Thermal dissociation becomes significant only at very high T")
ax.legend()
plt.show()
print(f"N2 -> 2N bond dissociation enthalpy at 298 K: {dH_N2/1000:.0f} kJ/mol")
N2 -> 2N bond dissociation enthalpy at 298 K: 945 kJ/mol
5. Enthalpy vs. entropy control
\(\Delta G^\circ = \Delta H^\circ - T\Delta S^\circ\) makes the competition explicit: the enthalpy term dominates at low \(T\), the entropy term (weighted by \(T\)) at high \(T\). Dissociation (\(\Delta S^\circ > 0\), more gas moles) is entropy-driven and thus switched on by high temperature, while the exothermic water-gas shift is enthalpy-favoured at low temperature. We visualise the two contributions for N₂ dissociation.
[8]:
Tgrid = np.linspace(1000, 6000, 80)
with ThermochemicalCalculator() as calc:
dH = np.array([reaction_props(calc, {"N2": 1}, {"N": 2}, T)[0] for T in Tgrid]) / 1000
TdS = np.array([T * reaction_props(calc, {"N2": 1}, {"N": 2}, T)[1] for T in Tgrid]) / 1000
dG = np.array([reaction_props(calc, {"N2": 1}, {"N": 2}, T)[2] for T in Tgrid]) / 1000
fig, ax = plt.subplots()
ax.plot(Tgrid, dH, label=r"$\Delta H^\circ$")
ax.plot(Tgrid, TdS, label=r"$T\,\Delta S^\circ$")
ax.plot(Tgrid, dG, "k--", label=r"$\Delta G^\circ$")
ax.axhline(0, color="0.6", lw=0.8)
ax.set_xlabel("Temperature [K]")
ax.set_ylabel("energy [kJ/mol]")
ax.set_title(r"N$_2$ dissociation: $\Delta G^\circ = \Delta H^\circ - T\Delta S^\circ$")
ax.legend()
plt.show()
Summary
\(G^\circ = H^\circ - T S^\circ\) from
pyglenngives \(\Delta G^\circ(T)\) and hence \(K(T)\) for any reaction.The water-gas shift \(K\) falls through 1 near 1100 K; the van’t Hoff slope reproduces the reaction enthalpy.
Dissociation is entropy-driven and only matters at very high temperature, explaining the flame-temperature caveat of notebook 06.
Next: notebook 09 applies temperature-dependent properties to a gas-turbine cycle.