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.
[ ]:
from pyglenn import ThermochemicalCalculator, R
print("Universal gas constant R =", R, "J/(mol.K)")
[ ]:
%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.
[ ]:
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 * (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.
[ ]:
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())
[ ]:
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")
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.
[ ]:
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")
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.
[ ]:
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")
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.
[ ]:
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.