Gibbs’s phase rule on the phase diagram of water

Gibbs’s phase rule, \(F = C - P + 2\) (gibbs_phase_rule()), says how many intensive variables can change while a set of phases coexists. On a one-component phase diagram, single-phase regions are areas (\(F=2\)), two-phase boundaries are curves (\(F=1\)), and three phases meet only at a point (\(F=0\)). The sketch below draws water’s diagram around its triple point (vaporization and sublimation curves from the Clausius-Clapeyron equation, melting line from the Clapeyron slope) and labels each feature with its degrees of freedom.

import matplotlib.pyplot as plt
import numpy as np

from chemistrykit.thermo.systems.phase_equilibria import ClausiusClapeyron, gibbs_phase_rule

T_tp, P_tp = 273.16, 611.657  # water triple point
vaporization = ClausiusClapeyron(delta_h_vap=45050.0, T_ref=T_tp, P_ref=P_tp)
sublimation = ClausiusClapeyron(delta_h_vap=51060.0, T_ref=T_tp, P_ref=P_tp)
# Clapeyron slope of the melting line: dP/dT = dH_fus / (T dV_fus), negative for water.
dH_fus, dV_fus = 6010.0, -1.63e-6
slope = dH_fus / (T_tp * dV_fus)

T_vap = np.linspace(T_tp, 300.0, 100)
T_sub = np.linspace(250.0, T_tp, 100)
P_melt = np.logspace(np.log10(P_tp), 7, 100)
T_melt = T_tp + (P_melt - P_tp) / slope

fig, ax = plt.subplots(figsize=(7, 5))
ax.plot(T_vap, vaporization.pressure(T_vap), label="liquid-vapor (F = 1)")
ax.plot(T_sub, sublimation.pressure(T_sub), label="solid-vapor (F = 1)")
ax.plot(T_melt, P_melt, label="solid-liquid (F = 1)")
ax.plot([T_tp], [P_tp], "ko", label="triple point (F = 0)")
ax.set_yscale("log")
ax.set_xlim(250.0, 300.0)
ax.set_ylim(50.0, 1e7)
for x, y, text in [(258.0, 1e5, "ice"), (288.0, 1e5, "liquid"), (288.0, 200.0, "vapor")]:
    ax.text(x, y, f"{text}\nF = {gibbs_phase_rule(1, 1)}", ha="center")
ax.set_xlabel("T (K)")
ax.set_ylabel("P (Pa)")
ax.set_title("Water near its triple point: F = C - P + 2")
ax.legend(loc="lower right")
fig.tight_layout()

The rule for one and two components, and with a reaction constraint (e.g. \(N_2O_4 \rightleftharpoons 2NO_2\) in the gas phase, where the two species count as only one independent component):

for C in (1, 2):
    print(f"C = {C}: " + ", ".join(f"P = {P}: F = {gibbs_phase_rule(C, P)}" for P in range(1, C + 3)))
print("N2O4/NO2 gas with the reaction at equilibrium: F =", gibbs_phase_rule(n_components=2, n_phases=1, reactions=1))

plt.show()

C = 1: P = 1: F = 2, P = 2: F = 1, P = 3: F = 0
C = 2: P = 1: F = 3, P = 2: F = 2, P = 3: F = 1, P = 4: F = 0
N2O4/NO2 gas with the reaction at equilibrium: F = 2