The van der Waals equation of state

Learning Scientific Programming with Python (2nd edition)

P6.4.3: The van der Waals equation of state

Question P6.4.3

The van der Waals equation of state may be written as follows to give the pressure, $p$, of a gas from its molar volume, $V$, and temperature, $T$: $$ p = \frac{RT}{V - b} - \frac{a}{V^2}, $$ where $a$ and $b$ are molecule-specific constants and $R=8.314\;\mathrm{J\,K^{-1}\,mol^{-1}}$ is the gas constant. It can readily be rearranged to yield the temperature for a given pressure and volume, but its form giving the molar volume in terms of pressure and temperature is a cubic equation: $$ pV^3 - (pb + RT)V^2 + aV - ab = 0 $$ Of the three roots to this equation, below the critical point, $(T_\mathrm{c}, p_\mathrm{c})$ all are real: the largest and smallest give the molar volume of the gas phase and liquid phase respectively; above the critical point, where no liquid phase exists, only one root is real and gives the molar volume of the gas (also known in this region as a supercritical fluid). The critical point is given by the condition $(\partial p / \partial V)_T = (\partial^2 p / \partial V^2)_T = 0$ and for a van der Waals gas is given by the formulae $$ T_c = \frac{8a}{27Rb}, \quad \quad p_c = \frac{a}{27b^2} $$

For ammonia the van der Waals constants are $a = 4.225\;\mathrm{L^2\,bar\,mol^{-2}}$ and $b = 0.03707\;\mathrm{L\,mol^{-1}}$.

(a) Find the critical point of ammonia, and then determine the molar volume at room temperature and pressure, $(298\;\mathrm{K}, 1\;\mathrm{atm})$ and at $(500\;\mathrm{K}, 12\;\mathrm{MPa})$.

(b) An isotherm is the set of points $(p, V)$ at a constant temperature satisfying an equation of state. Plot the isotherm ($p$ against $V$) for ammonia at 350 K using the van der Waals equation of state and compare it with the 350 K isotherm for an ideal gas, which has the equation of state $p=RT/V$.

Solution P6.4.3

(a) First write the van der Waals constants in SI units and calculate the critical temperature and pressure, $(T_\mathrm{c}, p_\mathrm{c})$:

In [x]: a = 4.225 / 1.e6 * 1.e5 # m^6.Pa.mol^-2
In [x]: b = 0.03707 / 1.e3      # m^3.mol^-1
In [x]: R = 8.314               # J.K^-1.mol^-1
In [x]: Tc, pc = 8*a/27/R/b, a/27/b**2
In [x]: print('Tc = {:.1f} K, pc = {:.1f} MPa'.format(Tc, pc/1.e6))
Tc = 406.2 K, pc = 11.4 MPa

Next define a function returning a Polynomial object representing the equation of state:

In [x]: import numpy as np
In [x]: Polynomial = np.polynomial.Polynomial
In [x]: def get_poly(a, b, T, p):
   ...:     return Polynomial([-a*b, a, -(p*b + R*T), p])

Room temperature (298 K) and pressure ($1\;\mathrm{atm} = 101325\;\mathrm{Pa}$) are below the critical point, and so we expect three roots, from which we want the smallest and largest:

In [x]: poly = get_poly(a, b, 298., 101325.)
In [x]: roots = poly.roots()
In [x]: roots.sort()
In [x]: Vliq, Vgas = roots[0], roots[-1]
In [x]: print('V(liq) = {:.1f} cm^3.mol^-1\n'
              'V(gas) = {:.1f} dm^3.mol^-1'.format(Vliq*1.e6, Vgas*1.e3))
V(liq) = 54.4 cm^3.mol^-1
V(gas) = 24.3 dm^3.mol^-1

At $(500\;\mathrm{K}, 12\;\mathrm{MPa})$ we are above the critical point and there is only one real root:

In [x]: poly = get_poly(a, b, 500, 12.e6)
In [x]: roots = poly.roots()
In [x]: Vgas = roots[roots.imag==0][0].real
In [x]: print('V(gas) = {:.1f} cm^3.mol^-1'.format(Vgas*1.e6))
V(gas) = 271.5 cm^3.mol^-1

(b) 350 K is below the critical point, so ammonia will condense under sufficient compression. The molar volume will therefore span several orders of magnitude between $10^{-5}$ and $10^{-2}\;\mathrm{m^3}$ and a log scale is appropriate.

In [x]: import matplotlib.pyplot as plt
In [x]: plt.plot(V, R*T/V, label='Ideal gas')
In [x]: plt.plot(V, R*T/(V-b) - a*V**-2, label='van der Waals')
In [x]: plt.ylim(0, 1.e7)
In [x]: plt.xlabel(r"$V\;/\mathrm{m^3\,mol^{-1}}$")
In [x]: plt.ylabel(r"$p\;/\mathrm{Pa}$")
In [x]: plt.xscale('log')
In [x]: plt.legend()
In [x]: plt.show()

350 K isotherms for the ideal gas and van der Waals equations of state for ammonia.