Part II: Thermodynamics & Equations of State

Mixtures & Solutions

A rigorous treatment of ideal and non-ideal mixtures: Raoult's law, colligative properties, Henry's law, regular solution theory, and their applications to distillation, osmotic processes, and phase equilibria

1. Introduction — Ideal and Non-Ideal Mixtures

A mixture is a system containing two or more chemical species. When we dissolve a solute in a solvent, or combine two miscible liquids, the thermodynamic properties of the resulting solution depend on the nature and strength of the intermolecular interactions between the components. The central quantity governing mixture behavior is the chemical potential $\mu_i$, which determines how the Gibbs energy changes when species $i$ is added to the system.

Ideal Mixtures

In an ideal mixture, the intermolecular interactions between unlike molecules (A–B) are identical to those between like molecules (A–A and B–B). The enthalpy of mixing is exactly zero: $\Delta H_{\text{mix}} = 0$. The only driving force for mixing is the entropy increase associated with the increased number of accessible microstates.

Examples: benzene + toluene, hexane + heptane — structurally similar molecules with nearly identical intermolecular forces.

Non-Ideal Mixtures

Real solutions exhibit deviations from ideal behavior because A–B interactions differ from A–A and B–B interactions. Positive deviations ($\Delta H_{\text{mix}} > 0$) occur when A–B interactions are weaker; negative deviations ($\Delta H_{\text{mix}} < 0$) when A–B interactions are stronger.

Positive: ethanol + hexane. Negative: chloroform + acetone (hydrogen bonding between unlike molecules).

The thermodynamic framework for mixtures relies on partial molar quantities. For any extensive property $Y$ (volume, entropy, Gibbs energy), the partial molar quantity of component $i$ is:

$\bar{Y}_i = \left(\frac{\partial Y}{\partial n_i}\right)_{T,P,n_{j \neq i}}$

The partial molar Gibbs energy is precisely the chemical potential: $\bar{G}_i = \mu_i$.

For an ideal gas mixture, the chemical potential of component $i$ is:

$\mu_i^{\text{ig}}(T,P) = \mu_i^{\circ}(T) + RT \ln\left(\frac{P_i}{P^{\circ}}\right) = \mu_i^{\circ}(T) + RT \ln\left(\frac{y_i P}{P^{\circ}}\right)$

where $y_i$ is the mole fraction in the gas phase, $P_i = y_i P$ is the partial pressure, and $\mu_i^{\circ}(T)$ is the standard chemical potential at pressure $P^{\circ}$.

2. Derivation: Raoult's Law and Ideal Solutions

Raoult's law follows directly from requiring thermodynamic equilibrium between liquid and vapor phases. At equilibrium, the chemical potential of each component must be equal in both phases:

$\mu_i^{\text{liquid}}(T,P) = \mu_i^{\text{vapor}}(T,P)$

For the vapor phase (treated as ideal gas), the chemical potential is:

$\mu_i^{\text{vapor}} = \mu_i^{*,\text{vapor}}(T) + RT \ln\left(\frac{P_i}{P_i^*}\right)$

where $P_i^*$ is the vapor pressure of pure component $i$ at temperature $T$. For the liquid phase of an ideal solution, we define the chemical potential using the convention that intermolecular interactions are identical for all species:

$\mu_i^{\text{liquid}} = \mu_i^{*,\text{liquid}}(T,P) + RT \ln x_i$

where $x_i$ is the mole fraction in the liquid phase and $\mu_i^{*,\text{liquid}}$ is the chemical potential of pure liquid $i$. For the pure component ($x_i = 1$), liquid-vapor equilibrium gives $\mu_i^{*,\text{liquid}} = \mu_i^{*,\text{vapor}}$. Substituting the equilibrium condition:

$\mu_i^{*,\text{liquid}} + RT \ln x_i = \mu_i^{*,\text{vapor}} + RT \ln\left(\frac{P_i}{P_i^*}\right)$

$RT \ln x_i = RT \ln\left(\frac{P_i}{P_i^*}\right)$

$\boxed{P_i = x_i \, P_i^*}$

Raoult's Law: The partial vapor pressure of each component equals its mole fraction times the vapor pressure of the pure component.

Thermodynamics of Ideal Mixing

For $n$ moles of an ideal solution, the Gibbs energy of mixing is obtained by comparing the Gibbs energy of the mixture with that of the separated pure components:

$\Delta G_{\text{mix}} = G_{\text{mixture}} - \sum_i n_i \mu_i^* = \sum_i n_i (\mu_i - \mu_i^*)$

$\Delta G_{\text{mix}} = \sum_i n_i RT \ln x_i = nRT \sum_i x_i \ln x_i$

Since $0 < x_i < 1$, each $\ln x_i < 0$, so $\Delta G_{\text{mix}} < 0$: mixing is always spontaneous for ideal solutions. The entropy, enthalpy, and volume of mixing follow:

$\Delta S_{\text{mix}} = -\left(\frac{\partial \Delta G_{\text{mix}}}{\partial T}\right)_P = -nR \sum_i x_i \ln x_i > 0$

$\Delta H_{\text{mix}} = \Delta G_{\text{mix}} + T\Delta S_{\text{mix}} = nRT\sum_i x_i \ln x_i + T\left(-nR\sum_i x_i \ln x_i\right) = 0$

$\Delta V_{\text{mix}} = \left(\frac{\partial \Delta G_{\text{mix}}}{\partial P}\right)_T = 0$

Key Result: Ideal Solution Properties

$\Delta G_{\text{mix}} = nRT \sum_i x_i \ln x_i < 0$ — mixing always spontaneous
$\Delta S_{\text{mix}} = -nR \sum_i x_i \ln x_i > 0$ — entropy always increases
$\Delta H_{\text{mix}} = 0$ — no heat released or absorbed
$\Delta V_{\text{mix}} = 0$ — no volume change on mixing

3. Derivation: Colligative Properties

Colligative properties depend only on the number of solute particles, not their identity. They arise because the solute lowers the chemical potential of the solvent in the liquid phase. We consider a dilute solution with solvent (1) and non-volatile solute (2).

The chemical potential of the solvent in the solution is:

$\mu_1^{\text{soln}}(T,P) = \mu_1^*(T,P) + RT \ln x_1 \approx \mu_1^*(T,P) - RT\,x_2$

Using $\ln x_1 = \ln(1 - x_2) \approx -x_2$ for dilute solutions ($x_2 \ll 1$).

Boiling Point Elevation

At the boiling point, the liquid solvent and vapor are in equilibrium:$\mu_1^{\text{soln}}(T_b) = \mu_1^{\text{vapor}}(T_b)$. For the pure solvent, this equilibrium occurs at $T_b^*$ where $\mu_1^*(T_b^*) = \mu_1^{\text{vapor}}(T_b^*)$. The solute lowers $\mu_1^{\text{soln}}$, so a higher temperature is needed to reach equilibrium with the vapor. Starting from:

$\mu_1^*(T_b) + RT_b \ln x_1 = \mu_1^{\text{vapor}}(T_b)$

Using the Gibbs-Helmholtz equation $\frac{\partial(\mu/T)}{\partial T}\bigg|_P = -\frac{\bar{H}}{T^2}$ and integrating from $T_b^*$ to $T_b = T_b^* + \Delta T_b$:

$\ln x_1 = \frac{\Delta H_{\text{vap}}}{R}\left(\frac{1}{T_b^*} - \frac{1}{T_b}\right) \approx \frac{\Delta H_{\text{vap}} \, \Delta T_b}{R \, (T_b^*)^2}$

Substituting $\ln x_1 \approx -x_2$ and converting mole fraction to molality $m = \frac{x_2}{(1 - x_2) M_1 / 1000} \approx \frac{1000 \, x_2}{M_1}$:

$\boxed{\Delta T_b = K_b \, m \quad \text{where} \quad K_b = \frac{R \, (T_b^*)^2 \, M_1}{1000 \, \Delta H_{\text{vap}}}}$

$M_1$ is the molar mass of the solvent (g/mol), $m$ is the molality (mol/kg). For water: $K_b = 0.512$ K·kg/mol.

Freezing Point Depression

By an analogous derivation, equating the chemical potential of the solvent in solution with that of the pure solid solvent at the freezing point:

$\ln x_1 = -\frac{\Delta H_{\text{fus}}}{R}\left(\frac{1}{T_f^*} - \frac{1}{T_f}\right) \approx -\frac{\Delta H_{\text{fus}} \, \Delta T_f}{R \, (T_f^*)^2}$

$\boxed{\Delta T_f = K_f \, m \quad \text{where} \quad K_f = \frac{R \, (T_f^*)^2 \, M_1}{1000 \, \Delta H_{\text{fus}}}}$

For water: $K_f = 1.86$ K·kg/mol. Note $K_f > K_b$ because $\Delta H_{\text{fus}} < \Delta H_{\text{vap}}$.

Osmotic Pressure

Consider a solution separated from pure solvent by a semipermeable membrane (permeable only to the solvent). At equilibrium, the chemical potential of the solvent must be equal on both sides. The pressure on the solution side exceeds that on the pure solvent side by the osmotic pressure $\pi$:

$\mu_1^*(T, P + \pi) = \mu_1^*(T, P) + RT \ln x_1$

Using $\left(\frac{\partial \mu}{\partial P}\right)_T = \bar{V}_1 \approx V_1^*$ (molar volume of pure solvent):

$V_1^* \, \pi = -RT \ln x_1 \approx RT \, x_2$

For a dilute solution, $x_2 \approx n_2/n_1$ and $n_1 V_1^* \approx V$ (total volume):

$\boxed{\pi = \frac{n_2}{V} RT = M \, RT}$

van't Hoff equation: $M = n_2/V$ is the molar concentration of solute. This has the same form as the ideal gas law!

4. Derivation: Henry's Law and Dilute Solutions

Raoult's law describes the behavior of the solvent in a dilute solution (as $x_1 \to 1$). The solute, however, experiences a very different molecular environment — it is surrounded almost entirely by solvent molecules. Henry's law describes the solute behavior in this limit.

Consider component 2 (solute) at very low concentration ($x_2 \to 0$). Each solute molecule is surrounded by solvent molecules, so the effective interaction it experiences is the A–B interaction, not the B–B interaction. The chemical potential of the solute in the liquid is:

$\mu_2^{\text{soln}}(T,P) = \mu_2^{\dagger}(T,P) + RT \ln x_2$

where $\mu_2^{\dagger}$ is a reference chemical potential that accounts for the solute-solvent interactions (different from the pure component $\mu_2^*$).

Equating with the vapor phase chemical potential at equilibrium:

$\mu_2^{\dagger} + RT \ln x_2 = \mu_2^{\text{vapor},\circ} + RT \ln\left(\frac{P_2}{P^{\circ}}\right)$

$P_2 = x_2 \, \exp\left(\frac{\mu_2^{\dagger} - \mu_2^{\text{vapor},\circ}}{RT}\right) \cdot P^{\circ}$

$\boxed{P_2 = K_H \, x_2}$

Henry's law: The Henry's law constant $K_H = P^{\circ} \exp\left(\frac{\mu_2^{\dagger} - \mu_2^{\text{vapor},\circ}}{RT}\right)$ encodes the strength of the solute-solvent interaction.

The key distinction: Raoult's law uses $P_2^*$ (vapor pressure of pure solute), while Henry's law uses $K_H$ (which depends on the solute-solvent interaction). If A–B interactions are weaker than B–B interactions, $K_H > P_2^*$ (positive deviation). If A–B interactions are stronger, $K_H < P_2^*$ (negative deviation).

Activity Coefficients

For non-ideal solutions, we introduce the activity $a_i$ and activity coefficient $\gamma_i$ to account for deviations from ideal behavior:

$\mu_i = \mu_i^* + RT \ln a_i \quad \text{where} \quad a_i = \gamma_i \, x_i$

$\gamma_i = \frac{a_i}{x_i}$

For an ideal solution, $\gamma_i = 1$ for all components. Raoult's law convention:$\gamma_i \to 1$ as $x_i \to 1$. Henry's law convention:$\gamma_i^H \to 1$ as $x_i \to 0$.

The modified Raoult's law incorporating activity coefficients is:

$P_i = \gamma_i \, x_i \, P_i^*$

The Gibbs-Duhem equation constrains activity coefficients. At constant $T$ and $P$, for a binary mixture:

$x_1 \, d\ln\gamma_1 + x_2 \, d\ln\gamma_2 = 0$

5. Derivation: Regular Solution Theory

Regular solution theory, developed by Joel Hildebrand, provides the simplest model for non-ideal mixtures. It assumes that the excess entropy of mixing is zero ($\Delta S^E = 0$), so all non-ideality comes from a non-zero enthalpy of mixing. The molecules are assumed to mix randomly (as in an ideal solution), but with different interaction energies.

Consider a binary mixture of $N$ total molecules with $N_1 = x_1 N$ of type 1 and $N_2 = x_2 N$ of type 2, each having coordination number $z$. Let $\epsilon_{11}$, $\epsilon_{22}$, and $\epsilon_{12}$ be the pairwise interaction energies. In a random mixture, the number of each type of pair is:

$N_{11} = \frac{z N}{2} x_1^2, \quad N_{22} = \frac{z N}{2} x_2^2, \quad N_{12} = z N \, x_1 x_2$

The total interaction energy of the mixture is:

$U_{\text{mix}} = N_{11}\epsilon_{11} + N_{22}\epsilon_{22} + N_{12}\epsilon_{12}$

The energy of the unmixed pure components is:

$U_{\text{pure}} = \frac{zN}{2}\left(x_1 \epsilon_{11} + x_2 \epsilon_{22}\right)$

The enthalpy of mixing (at constant pressure, approximately equal to the energy of mixing for condensed phases) is:

$\Delta H_{\text{mix}} = U_{\text{mix}} - U_{\text{pure}} = \frac{zN}{2}\left(\epsilon_{12} - \frac{\epsilon_{11} + \epsilon_{22}}{2}\right) \cdot 2x_1 x_2$

$\Delta H_{\text{mix}} = n \, w \, x_1 x_2$

where the interaction parameter is $w = z N_A \left(\epsilon_{12} - \frac{\epsilon_{11} + \epsilon_{22}}{2}\right)$.

Since the excess entropy is zero for a regular solution, the Gibbs energy of mixing is:

$\boxed{\Delta G_{\text{mix}} = nRT\left(x_1 \ln x_1 + x_2 \ln x_2\right) + n \, w \, x_1 x_2}$

The first term is the ideal entropy contribution; the second is the excess enthalpy.

Phase Separation and Critical Solution Temperature

When the interaction parameter $w$ is sufficiently positive (A–B repulsion), the system can lower its Gibbs energy by separating into two phases. The stability criterion requires the second derivative of $\Delta G_{\text{mix}}$ to be positive:

$\frac{\partial^2 \Delta G_{\text{mix}}}{\partial x_1^2} = nRT\left(\frac{1}{x_1} + \frac{1}{x_2}\right) - 2nw$

The minimum of this expression occurs at $x_1 = x_2 = 1/2$, where it equals $n(4RT - 2w)$. The system becomes unstable when this is negative:

$\boxed{w > 2RT \quad \Longrightarrow \quad \text{Phase separation occurs}}$

The upper critical solution temperature (UCST) is $T_c = w / 2R$. Below this temperature, the mixture separates into two phases.

The activity coefficients in regular solution theory follow from $\ln \gamma_i = \frac{w}{RT} x_j^2$, which is the Margules one-parameter model. The two-parameter Margules model allows asymmetry:

$\ln \gamma_1 = x_2^2 \left[A_{12} + 2(A_{21} - A_{12}) x_1\right]$

$\ln \gamma_2 = x_1^2 \left[A_{21} + 2(A_{12} - A_{21}) x_2\right]$

6. Applications

Distillation

Distillation exploits the difference in composition between liquid and vapor phases at equilibrium. For a binary ideal mixture, the vapor is enriched in the more volatile component. The McCabe-Thiele method uses the equilibrium curve (from Raoult's law) and operating lines to determine the number of theoretical plates needed for separation.

Azeotropes present a fundamental limit: at the azeotropic composition, liquid and vapor have the same composition, so simple distillation cannot achieve further separation. Examples: ethanol-water (95.6% ethanol at 78.1°C), HCl-water.

Osmotic Pressure in Biology

Osmotic pressure is critical in biology. Cell membranes act as semipermeable barriers. Red blood cells in a hypotonic solution swell and lyse; in hypertonic solution they crenate. Isotonic saline (0.9% NaCl, $\pi \approx 7.7$ atm) maintains cellular integrity.

Reverse osmosis (applying $P > \pi$) is used for desalination. Seawater has $\pi \approx 27$ atm, requiring pressures of 50–80 atm for efficient desalination.

Azeotropes

A minimum-boiling azeotrope occurs when positive deviations from Raoult's law are large enough that the total pressure exceeds both pure component pressures. The mixture boils at a lower temperature than either pure component.

A maximum-boiling azeotrope arises from strong negative deviations (e.g., HNO$_3$–water at 68.4% HNO$_3$). The excess A–B attraction stabilizes the liquid, requiring more energy to vaporize.

Polymer Solutions (Flory-Huggins)

For polymer solutions, the entropy of mixing is drastically reduced because a polymer chain occupies many lattice sites. The Flory-Huggins theory gives:

$\frac{\Delta G_{\text{mix}}}{NkT} = \frac{\phi_1}{r_1}\ln\phi_1 + \frac{\phi_2}{r_2}\ln\phi_2 + \chi\,\phi_1\phi_2$

where $\phi_i$ are volume fractions, $r_i$ are the number of segments, and $\chi$ is the Flory-Huggins interaction parameter. For a polymer ($r_2 \gg 1$), the critical point shifts to very low polymer concentrations.

7. Historical Context

William Henry (1803)

The English chemist William Henry published his law describing the solubility of gases in liquids, showing that the amount of dissolved gas is proportional to its partial pressure above the liquid. His work provided the first quantitative description of gas-liquid equilibria and remains essential in fields from carbonated beverage production to blood gas physiology.

François-Marie Raoult (1887)

Raoult, a French chemist, systematically measured the vapor pressure lowering of solvents upon addition of non-volatile solutes. His meticulous experiments with hundreds of solute-solvent combinations established the empirical law that bears his name. Raoult's work provided a practical method for determining molecular weights of dissolved substances.

Jacobus Henricus van't Hoff (1887)

Van't Hoff derived the osmotic pressure equation $\pi V = nRT$ by analogy with the ideal gas law. He showed that dissolved molecules exert osmotic pressure just as gas molecules exert gas pressure. This work earned him the first Nobel Prize in Chemistry (1901). His insight unified the thermodynamics of solutions with the kinetic theory of gases.

Joel Hildebrand (1916–1970)

Hildebrand developed regular solution theory over decades of work at UC Berkeley. He introduced the solubility parameter $\delta = \sqrt{\Delta U_{\text{vap}}/V_m}$ to predict miscibility: substances with similar $\delta$ values tend to be miscible. His 1936 book “Solubility of Non-Electrolytes” became the standard reference for solution thermodynamics.

8. Python Simulation — T-x Phase Diagram

The following simulation plots a temperature-composition (T-x) phase diagram for a binary mixture, showing the liquid (bubble) and vapor (dew) curves. It includes both the ideal case (Raoult's law) and a non-ideal case with a minimum-boiling azeotrope using the one-parameter Margules model for activity coefficients.

Python

script.py158 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

# ============================================================
# T-x Phase Diagram for Binary Mixture
# Includes ideal (Raoult) and non-ideal (Margules azeotrope) cases
# ============================================================

# Antoine equation: log10(P/mmHg) = A - B/(C + T/degC)
# Parameters for benzene (1) and toluene (2)
def antoine_pressure(T_C, A, B, C):
    """Vapor pressure in mmHg from Antoine equation."""
    return 10.0**(A - B / (C + T_C))

# Antoine parameters (NIST)
# Benzene
A1, B1, C1 = 6.90565, 1211.033, 220.790
# Toluene
A2, B2, C2 = 6.95464, 1344.800, 219.482

# ---- IDEAL CASE (Raoult's Law) ----
# At a given T, find x1 from: P_total = x1*P1*(T) + (1-x1)*P2*(T) = 760 mmHg
P_total = 760.0  # 1 atm in mmHg

# Boiling points of pure components
# Benzene: ~80.1 C, Toluene: ~110.6 C
T_range = np.linspace(79.0, 112.0, 500)

x1_liquid_ideal = np.zeros_like(T_range)
x1_vapor_ideal = np.zeros_like(T_range)

for i, T in enumerate(T_range):
    P1star = antoine_pressure(T, A1, B1, C1)
    P2star = antoine_pressure(T, A2, B2, C2)
    # From Raoult: P = x1*P1* + (1-x1)*P2* = 760
    x1 = (P_total - P2star) / (P1star - P2star)
    x1_liquid_ideal[i] = x1
    # Vapor composition: y1 = x1*P1*/P
    y1 = x1 * P1star / P_total
    x1_vapor_ideal[i] = y1

# Filter valid range
mask_ideal = (x1_liquid_ideal >= 0) & (x1_liquid_ideal <= 1)
T_ideal = T_range[mask_ideal]
x1_liq_ideal = x1_liquid_ideal[mask_ideal]
x1_vap_ideal = x1_vapor_ideal[mask_ideal]

# ---- NON-IDEAL CASE (Margules model with azeotrope) ----
# Use a hypothetical system with Margules parameter w/(RT) that causes azeotrope
# Modified Raoult's law: P_i = gamma_i * x_i * P_i*
# One-parameter Margules: ln(gamma_1) = A * x2^2, ln(gamma_2) = A * x1^2
# Use ethanol(1)-water(2) like system with synthetic Antoine params

# Synthetic Antoine parameters for a system with azeotrope
# Component 1 (more volatile): bp = 78 C
# Component 2 (less volatile): bp = 100 C
def P1_star(T):
    return antoine_pressure(T, 8.11220, 1592.864, 226.184)  # ethanol-like

def P2_star(T):
    return antoine_pressure(T, 8.07131, 1730.630, 233.426)  # water-like

# Margules parameter (dimensionless A = w/RT, roughly constant)
A_margules = 1.7  # strong positive deviation -> azeotrope

T_range2 = np.linspace(76.0, 101.0, 500)

x1_liquid_azeo = []
x1_vapor_azeo = []
T_azeo_list = []

for T in T_range2:
    P1s = P1_star(T)
    P2s = P2_star(T)
    # Scan x1 from 0 to 1 and find where P_total = 760
    x1_arr = np.linspace(0.001, 0.999, 2000)
    x2_arr = 1.0 - x1_arr
    gamma1 = np.exp(A_margules * x2_arr**2)
    gamma2 = np.exp(A_margules * x1_arr**2)
    P_calc = gamma1 * x1_arr * P1s + gamma2 * x2_arr * P2s

# Find crossings where P_calc = 760
    diff = P_calc - P_total
    sign_changes = np.where(np.diff(np.sign(diff)))[0]

for idx in sign_changes:
        # Linear interpolation
        x1_val = x1_arr[idx] + (x1_arr[idx+1] - x1_arr[idx]) * (-diff[idx]) / (diff[idx+1] - diff[idx])
        x2_val = 1.0 - x1_val
        g1 = np.exp(A_margules * x2_val**2)
        g2 = np.exp(A_margules * x1_val**2)
        y1 = g1 * x1_val * P1s / P_total
        x1_liquid_azeo.append(x1_val)
        x1_vapor_azeo.append(y1)
        T_azeo_list.append(T)

x1_liquid_azeo = np.array(x1_liquid_azeo)
x1_vapor_azeo = np.array(x1_vapor_azeo)
T_azeo_list = np.array(T_azeo_list)

# Sort for plotting
sort_liq = np.argsort(x1_liquid_azeo)
sort_vap = np.argsort(x1_vapor_azeo)

# ---- PLOTTING ----
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))

# Plot 1: Ideal case (Benzene-Toluene)
ax1.plot(x1_liq_ideal, T_ideal, 'b-', linewidth=2, label='Liquid (bubble curve)')
ax1.plot(x1_vap_ideal, T_ideal, 'r-', linewidth=2, label='Vapor (dew curve)')
ax1.fill_betweenx(T_ideal, x1_liq_ideal, x1_vap_ideal, alpha=0.1, color='purple')
ax1.set_xlabel('Mole fraction of benzene, $x_1$ / $y_1$', fontsize=12)
ax1.set_ylabel('Temperature (\u00b0C)', fontsize=12)
ax1.set_title('Ideal System: Benzene-Toluene\n(Raoult\'s Law)', fontsize=13, fontweight='bold')
ax1.legend(fontsize=10)
ax1.set_xlim(0, 1)
ax1.set_ylim(78, 112)
ax1.grid(True, alpha=0.3)
ax1.text(0.5, 95, 'Liquid + Vapor', fontsize=11, ha='center', style='italic', color='purple')
ax1.text(0.15, 108, 'Liquid', fontsize=11, ha='center', color='blue')
ax1.text(0.85, 82, 'Vapor', fontsize=11, ha='center', color='red')

# Plot 2: Non-ideal case with azeotrope
ax2.plot(x1_liquid_azeo[sort_liq], T_azeo_list[sort_liq], 'b-', linewidth=2, label='Liquid (bubble curve)')
ax2.plot(x1_vapor_azeo[sort_vap], T_azeo_list[sort_vap], 'r-', linewidth=2, label='Vapor (dew curve)')
ax2.set_xlabel('Mole fraction of component 1, $x_1$ / $y_1$', fontsize=12)
ax2.set_ylabel('Temperature (\u00b0C)', fontsize=12)
ax2.set_title('Non-Ideal System: Minimum-Boiling Azeotrope\n(Margules Model, A = 1.7)', fontsize=13, fontweight='bold')
ax2.legend(fontsize=10)
ax2.set_xlim(0, 1)
ax2.grid(True, alpha=0.3)

# Find and mark azeotrope point (where x1 = y1)
diff_xy = np.abs(x1_liquid_azeo - x1_vapor_azeo)
# Only consider points where we have both
if len(diff_xy) > 0:
    min_idx = np.argmin(diff_xy)
    azeo_x = x1_liquid_azeo[min_idx]
    azeo_T = T_azeo_list[min_idx]
    ax2.plot(azeo_x, azeo_T, 'k*', markersize=15, zorder=5)
    ax2.annotate(f'Azeotrope\n$x_1$ = {azeo_x:.2f}, T = {azeo_T:.1f}\u00b0C',
                xy=(azeo_x, azeo_T), xytext=(azeo_x + 0.15, azeo_T + 3),
                fontsize=9, arrowprops=dict(arrowstyle='->', color='black'),
                bbox=dict(boxstyle='round,pad=0.3', facecolor='yellow', alpha=0.8))

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight')
print("Phase diagrams generated successfully.")
print(f"\nIdeal system (Benzene-Toluene):")
print(f"  Benzene bp: {T_ideal[np.argmin(np.abs(x1_liq_ideal - 1.0))]:.1f} C")
print(f"  Toluene bp: {T_ideal[np.argmin(np.abs(x1_liq_ideal - 0.0))]:.1f} C")
if len(diff_xy) > 0:
    print(f"\nNon-ideal system (Margules azeotrope):")
    print(f"  Azeotrope composition: x1 = {azeo_x:.3f}")
    print(f"  Azeotrope temperature: {azeo_T:.1f} C")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

9. Fortran Simulation — Colligative Properties

The following Fortran program computes the boiling point elevation, freezing point depression, and osmotic pressure for various solute concentrations in water. It uses the fundamental equations derived in Section 3 with accurate thermodynamic data for water as the solvent.

Python

script.py161 lines

import subprocess, sys, os, csv, io
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import numpy as np

# --- Embedded Fortran source ---
fortran_source = r"""
program colligative_properties
  implicit none
  integer, parameter :: dp = selected_real_kind(15, 307)

! Physical constants
  real(dp), parameter :: R = 8.314d0        ! J/(mol K)

! Solvent properties (water)
  real(dp), parameter :: M1 = 18.015d0      ! g/mol molar mass of water
  real(dp), parameter :: Tb_star = 373.15d0  ! K boiling point
  real(dp), parameter :: Tf_star = 273.15d0  ! K freezing point
  real(dp), parameter :: dH_vap = 40670.0d0  ! J/mol enthalpy of vaporization
  real(dp), parameter :: dH_fus = 6010.0d0   ! J/mol enthalpy of fusion

! Derived ebullioscopic and cryoscopic constants
  real(dp) :: Kb, Kf

! Variables
  integer, parameter :: n_points = 20
  real(dp) :: m, dTb, dTf, pi_osm, c_molar
  real(dp) :: V_water_L
  integer :: i

! Compute Kb and Kf from first principles
  ! Kb = R * Tb*^2 * M1 / (1000 * dH_vap)
  Kb = R * Tb_star**2 * M1 / (1000.0d0 * dH_vap)
  ! Kf = R * Tf*^2 * M1 / (1000 * dH_fus)
  Kf = R * Tf_star**2 * M1 / (1000.0d0 * dH_fus)

write(*,'(A)') '======================================================'
  write(*,'(A)') '  COLLIGATIVE PROPERTIES OF AQUEOUS SOLUTIONS (Fortran)'
  write(*,'(A)') '======================================================'
  write(*,'(A)')
  write(*,'(A,F10.4,A)') '  Ebullioscopic constant Kb = ', Kb, ' K kg/mol'
  write(*,'(A,F10.4,A)') '  Cryoscopic constant    Kf = ', Kf, ' K kg/mol'
  write(*,'(A)')
  write(*,'(A)') '  Verification:'
  write(*,'(A,F10.4,A)') '    Kb (literature) = 0.5120 K kg/mol'
  write(*,'(A,F10.4,A)') '    Kf (literature) = 1.8600 K kg/mol'
  write(*,'(A)')

write(*,'(A)') 'DATA_START'
  write(*,'(A)') 'molality,dTb_K,dTf_K,pi_atm,Tb_C,Tf_C'

do i = 1, n_points
    ! Molality from 0.05 to 5.0 mol/kg
    m = 0.05d0 + (5.0d0 - 0.05d0) * dble(i - 1) / dble(n_points - 1)

! Boiling point elevation
    dTb = Kb * m

! Freezing point depression
    dTf = Kf * m

! Osmotic pressure: pi = M*R*T where M ~ m * density_water
    ! For dilute solutions, molarity ~ molality * density(water) ~ m * 1.0
    c_molar = m * 1.0d0  ! approximate: molarity ~ molality for dilute aq.
    pi_osm = c_molar * R * 298.15d0 / 101325.0d0 * 1000.0d0  ! convert Pa to atm
    ! More precisely: pi = c * R * T, with c in mol/m^3
    ! c [mol/m^3] = m [mol/kg] * 1000 [kg/m^3]
    ! pi [Pa] = c * R * T = m * 1000 * 8.314 * 298.15
    pi_osm = m * 1000.0d0 * R * 298.15d0 / 101325.0d0  ! in atm

write(*,'(F8.4,A,F10.6,A,F10.6,A,F10.4,A,F10.4,A,F10.4)') &
      m, ',', dTb, ',', dTf, ',', pi_osm, ',', &
      (Tb_star + dTb - 273.15d0), ',', (Tf_star - dTf - 273.15d0)
  end do

write(*,'(A)') 'DATA_END'
  write(*,'(A)')
  write(*,'(A)') '  Physical interpretation:'
  write(*,'(A)') '  - dTb increases linearly with molality (dilute limit)'
  write(*,'(A)') '  - dTf is ~3.6x larger than dTb (dH_fus << dH_vap)'
  write(*,'(A)') '  - Osmotic pressure grows rapidly: 1 mol/kg -> ~24.4 atm'

end program colligative_properties
"""

# Write Fortran source
with open('colligative.f90', 'w') as f:
    f.write(fortran_source)

# Compile
result = subprocess.run(['gfortran', '-O2', '-o', 'colligative', 'colligative.f90'],
                       capture_output=True, text=True)
if result.returncode != 0:
    print("Compilation error:", result.stderr)
    sys.exit(1)

# Run
result = subprocess.run(['./colligative'], capture_output=True, text=True)
output = result.stdout
print(output.split('DATA_START')[0])

# Parse CSV data
lines = output.split('DATA_START')[1].split('DATA_END')[0].strip().split('\n')
header = lines[0]
data = []
for line in lines[1:]:
    vals = [float(x.strip()) for x in line.split(',')]
    data.append(vals)

data = np.array(data)
molality = data[:, 0]
dTb = data[:, 1]
dTf = data[:, 2]
pi_atm = data[:, 3]
Tb_C = data[:, 4]
Tf_C = data[:, 5]

# Plot results
fig, axes = plt.subplots(1, 3, figsize=(16, 5))

# Plot 1: Boiling point elevation
axes[0].plot(molality, dTb, 'r-o', markersize=4, linewidth=2, label='$\\Delta T_b$')
axes[0].set_xlabel('Molality (mol/kg)', fontsize=12)
axes[0].set_ylabel('$\\Delta T_b$ (K)', fontsize=12)
axes[0].set_title('Boiling Point Elevation', fontsize=13, fontweight='bold')
axes[0].grid(True, alpha=0.3)
axes[0].legend(fontsize=10)

# Plot 2: Freezing point depression
axes[1].plot(molality, dTf, 'b-s', markersize=4, linewidth=2, label='$\\Delta T_f$')
axes[1].set_xlabel('Molality (mol/kg)', fontsize=12)
axes[1].set_ylabel('$\\Delta T_f$ (K)', fontsize=12)
axes[1].set_title('Freezing Point Depression', fontsize=13, fontweight='bold')
axes[1].grid(True, alpha=0.3)
axes[1].legend(fontsize=10)

# Plot 3: Osmotic pressure
axes[2].plot(molality, pi_atm, 'g-^', markersize=4, linewidth=2, label='$\\pi$')
axes[2].set_xlabel('Molality (mol/kg)', fontsize=12)
axes[2].set_ylabel('Osmotic Pressure (atm)', fontsize=12)
axes[2].set_title('Osmotic Pressure at 25\u00b0C', fontsize=13, fontweight='bold')
axes[2].grid(True, alpha=0.3)
axes[2].legend(fontsize=10)

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight')

# Print summary
print("\nResults plotted successfully.")
print(f"At m = 1.0 mol/kg:  dTb = {np.interp(1.0, molality, dTb):.4f} K, "
      f"dTf = {np.interp(1.0, molality, dTf):.4f} K, "
      f"pi = {np.interp(1.0, molality, pi_atm):.2f} atm")
print(f"At m = 0.1 mol/kg:  dTb = {np.interp(0.1, molality, dTb):.4f} K, "
      f"dTf = {np.interp(0.1, molality, dTf):.4f} K, "
      f"pi = {np.interp(0.1, molality, pi_atm):.2f} atm")

# Cleanup
os.remove('colligative.f90')
os.remove('colligative')

Click Run to execute the Python code

Code will be executed with Python 3 on the server

← Part II Overview Physical Chemistry Home →