Van der Waals

from Statistical Physics using Mathematica
© James J. Kelly, 1996-2001

Introduction

The van der Waals equation of state

(p + (a/(V^2))) (V - b) == 
      R T

has a surprisingly rich structure, as revealed by the indicator diagram below.  For low density (large molar volume V), the isotherms are very similar to the rectangular hyperbolae of the ideal gas, for which pressure is inversely proportional to volume.  However,  the steep slope in pressure versus volume for high densities is more typical of the nearly incompressible behavior of liquids.  For low temperatures and intermediate densities, the cubic van der Waals equation exhibits a rather interesting nonmonotonic backbending behavior.  The isotherm (Subscript[T, c]) which possesses a simple inflection point for which

\(\((∂p\/∂V)\)\_T = \(\((∂\^2 p\/∂V\^2)\)\_T = \(0\ \ \ ⟹\ \ \ κ\_T = ∞\)\)\)

must be identified with the critical point where the distinction between liquid and vapor subdivisions of the fluid phase loses its meaning.  The infinite compressibility suggests that relatively large density fluctuations can be expected near the critical point.  For temperatures below the critical temperature, the region in which pressure increases with increasing volume is unphysical, corresponding to an unstable situation.  It is in this region that we expect to find a transition between liquid and vapor phases.

Although the van der Waals equation of state may not provide the most accurate description of fluids, it does provide a good qualitative model of the liquid-vapor transition.  We examine some of the properties of the van der Waals equation, with particular attention to the vicinity of the critical point.  The physical origin of the equation itself is considered in more detail within the notebook virial.nb, in which the parameters of the van der Waals equation are related to the properties of the intermolecular interaction by means of the virial expansion.  The parameter a is related to the long-range attraction whereas b is related to the short range repulsion.

Glossary:
	p	= pressure
	V	= molar volume (volume per mole)
	T	= temperature
	R	= molar gas constant
	(Subscript[p, c])	= critical pressure
	(Subscript[V, c])	= critical molar volume
	(Subscript[T, c])	= critical temperature
	a	= parameter representing attractive part of intermolecular interaction
	b	= parameter representing repulsive interaction (excluded volume)
	(Subscript[V, g])	= molar volume for gas phase
	(Subscript[V, ℓ])	= molar volume for liquid phase

Initialization

Defaults, packages, and symbols

ClearAll["Global`*"];
Off[General::spell, General::spell1];

$DefaultFont = {"Times", 12};
$TextStyle = {FontFamily -> "Times", FontSize -> 12, 
          FontSlant -> "Italic"};

Needs["Utilities`Notation`"];

Symbolize[(Subscript[T, x])_];Symbolize[(Subscript[V, x])_];Symbolize[(Subscript[p, x])_]

Symbolize[\(p\&~\)]

Properties of the Van der Waals equation of state

Definition

vdWaalsEq = (p + (a/(V^2))) (V - b) == R T;

\(p\&~ = 
      p /. \(Solve[vdWaalsEq, 
              p]\)[[1]]\  // 
        Simplify\)

\(\((\(p\&~\ V\)\/\(R\ T\) /. V -> 1\/x\  // 
          Series[#1, {x, 0, 3}] &)\) /. x -> 1\/V\)

Critical parameters

The critical parameters are obtained by locating the inflection point in the critical isotherm.

\(eq1 = \(T\_c == 
            T /. \(Solve[∂\_V p\&~ == 0, 
              T]\)[[1]]\) /. 
        V -> V\_c\)

\(eq2 = \(T\_c == 
            T /. \(Solve[∂\_{V, 2}p\&~ == 0, 
              T]\)[[1]]\) /. 
        V -> V\_c\)

sol1 = Solve[{eq1, eq2}, {(Subscript[V, c]), 
            (Subscript[T, c])}][[3]]

sol2 = Solve[
          sol1 /. {Rule -> Equal}, {a, 
            b}][[1]]

\(p\_c = \(p\&~ /. {T -> T\_c, V -> V\_c}\) /. sol1\  // 
        Simplify\)

criticalParameters = {(Subscript[p, c]), (Subscript[V, c]), (Subscript[T, c])} /. sol1

With only two parameters, the critical compressibility factor for the van der Waals fluid is a constant, independent of fluid type.

criticalCompressibilityFactor = (((Subscript[p, c]) (Subscript[V, c]))/(R (Subscript[T, c]))) /. 
        sol1

Example: nitrogen

The following parameters for nitrogen are obtained by fitting data near 273 K and are expressed in SI units.

nitrogenParameters = {a -> 0.136, b -> 38.5 (10^(-6)), 
          R -> 8.314};

criticalParameters /. nitrogenParameters

These values may be compared with 
	(Subscript[p, c]) = 3.39 x (10^6) Pa 
	(Subscript[V, c]) = 9.01 (x10^(-5)) (m^3)/mole  
	(Subscript[T, c]) = 126.2  kelvin,
which corresponds to a critical compressibility factor of 0.291.  The van der Waals model misses the critical density, and hence the compressibility factor, by about 25% but does reproduce the critical temperature and pressure fairly accurately.  Obviously, it is not possible to reproduce all three critical parameters for an arbitrary fluid using a model with only two adjustable parameters.

Dimensionless form

It is useful to reduce the van der Waals equation to a universal dimensionless form, by expressing (p,  V,  T)  in units of ((Subscript[p, c]),  (Subscript[V, c]),  (Subscript[T, c])) and eliminating the constants (a,  b).

eq3 = 
      Thread[vdWaalsEq /. {V -> V (Subscript[V, c]), T -> T (Subscript[T, c]), 
                  p -> p (Subscript[p, c])} /. sol2 /. {(Subscript[V, c]) -> 1, 
              (Subscript[T, c]) -> 1, (Subscript[p, c]) -> 1}, Equal]  // 
        Simplify

reducedP = 
      Simplify[p /. Solve[eq3, 
              p][[1]]]

Therefore, we have obtained a universal equation of state for the van der Waals fluid expunged of any specific properties (a and b) of any particular substance.  This is an example of the principle of corresponding states which suggests that: All substances with the same reduced variables occupy corresponding states with the same behavior.

As a simple example of this principle, consider the compressibility factor defined by the ratio ((p V)/(R T)), which has the value unity for ideal gases.  As shown above, the compressibility factor for a van der Waals fluid assumes the universal value (3/8) at the critical point.  Because the van der Waals parameters can be fitted to data obtained far from the critical point, we would be able to predict both the local and compressibility factor for the critical point if real gases satisfied the van der Waals equation.  However, most real gases have critical compressibility factors between about 0.25 and 0.3, considerably below the prediction of the van der Waals equation.  Nevertheless, the narrow observed range supports the principle of corresponding states as a useful approximation despite the numerical inaccuracy of the critical values predicted by the van der Waals equation.  There are other empirical equations of state, such as the Dieterici equation, which obtain more accurate predictions for the critical compressibility factor using only two parameters and there are more complicated multiparameter models which describe both liquid and vapor phases more accurately.  Nevertheless, the van der Waals equation remains useful because of its historical importance and conceptual simplicity.

Virial expansion

(Series[#1, {x, 0, 3}] &)[((reducedP V)/T) /. V -> (1/x)] /. 
      x -> (1/V)

Display isotherms near critical point

Clear[plotPvsV, arrayPvsV]

plotPvsV[x_] := 
        Plot[Evaluate[reducedP /. T -> x], {V, 0.4, 
            4.0}, PlotRange -> {{0.6, 2.4}, {0.4, 2.5}}, 
          DisplayFunction -> Identity];

arrayPvsV = Table[plotPvsV[x], {x, 0.9, 1.1, 0.05}];

\(\(Show[arrayPvsV, PlotRange -> {{0.4, 4}, {0, 2}}, 
        PlotLabel -> "Isotherms near Critical Point", 
        Frame -> True, 
        FrameLabel -> {\*""V/(Subscript[V, c])"", \*""P/(Subscript[P, c])""}, 
        DisplayFunction -> $DisplayFunction];\)\)

Display equilibrium surface near critical point

\(\(Plot3D[Evaluate[reducedP], {V, 0.5, 2.0}, {T, 0.75, 1.25}, 
        AxesLabel -> {\*""V/(Subscript[V, c])"", \*""T/(Subscript[T, c])"", \*""p/(Subscript[p, c])""}, PlotPoints -> 25, 
        PlotRange -> {{0.5, 2.0}, {0.9, 1.1}, {0, 5}}];\)\)

Liquid-Vapor Transition

Maxwell construction for liquid-vapor coexistence

How should the backbending region be interpreted?  In this region, the molar volume appears to be a triple-valued function of pressure.  The solution with highest density (smallest V) can be associated with a liquid phase, whereas the low-density solution naturally represents the vapor phase.  However, the intermediate solution would appear, at first glance, to have no physical significance because it lies on an unphysical curve with negative compressibility that cannot represent an equilibrium state.  Nevertheless, some insight into the liquid-vapor phase transition can be obtained by considering the unphysical portion of the isotherm to also be meaningful.

Each phase can be represented by a Gibbs free enthalpy surface, G(T, p), and the stable phase is the one that lies lowest.  Coexistence between two phases is described by the curve formed by the intersection of the two Gibbs surfaces.  Under carefully controlled conditions, it is often possible to adjust the parameters of a system so that its state remains on a particular Gibbs surface beyond the intersection with the surface belonging to another phase.  Hence, we can sometimes maintain a particular phase for a relatively long time in a metastable state even though a different phase has a lower free enthalpy and thereby contains the true equilibrium state under those conditions.  However, mechanical or other disturbances of sufficient amplitude will generally cause an abrupt, sometimes explosive, transition to the stable phase.  For example, it is often possible to cool a gas without condensation occurring even at pressures as much as four times the equilibrium vapor pressure, a phenomenon which makes possible the operation of the Wilson cloud chamber.  Similarly, small droplets of water can be cooled to as low as -(40^∘) C without freezing if no nucleation centers (such as dust particles) are available.

The true isotherm for liquid-vapor coexistence should be horizontal (isobaric).  For a given temperature and pressure, the fluid can exist anywhere within a range of molar volumes between (Subscript[V, ℓ]) for a saturated (pure) liquid and (Subscript[V, g]) for a saturated (pure) vapor, with intermediate volumes corresponding to mixed phases in which the liquid and vapor phases coexist.  These phases are separated gravitationally in the lab, but the free enthalpy is the same for all liquid-vapor mixtures coexisting at the same temperature and pressure.  Therefore, the Maxwell equal-areas construction replaces the van der Waals isotherm in the coexistence region by a horizontal line for which the vapor pressure is chosen ensure that the free enthalpies of the saturated liquid and vapor phases, as obtained from the van der Waals equation of state, are equal.

Let G represent the molar Gibbs free enthalpy, such that dG  =  -SdT  +  Vdp.  Under isothermal conditions, dT = 0  ⟹ dG  =  Vdp.  Therefore, the requirement that the free enthalpy for saturated liquid and vapor phases must be equal can be satisfied by choosing the vapor pressure so that it divides the van der Waals curve into two lobes of equal area in the coexistence region.

area1 = ((Subscript[V, g]) - 
              (Subscript[V, ℓ])) reducedP /. V -> 
          (Subscript[V, ℓ])

area2 = ((∫_((Subscript[V, ℓ]))^((Subscript[V, g])))reducedP \[DifferentialD]V  // FullSimplify) /. 
        Log[a_]  -  Log[b_] -> Log[(a/b)]

pReduced[TT_, VV_] := 
      reducedP /. {T -> TT, V -> VV}

Since the Maxwell construction is fairly delicate, we must determine approximate roots before invoking the FindRoot function.  This is done by finding the extreme values of the pressure in the physical region and then finding the volumes for the pressure that is halfway between these extremes.  FindRoot can then be used to improve the approximate solutions for (Subscript[V, ℓ]) and (Subscript[V, g]), the liquid and gas volumes.   Nevertheless, it remains difficult to apply this method far from the critical point; a more refined method would be needed to extend the algorithm.

\(maxwellConstruction[TT_] := 
      Module[{vlimits, midP, V1, V2, volumes}, 
        vlimits = Solve[∂\_V pReduced[TT, V] == 0, V]; 
        midP = 0.5\ pReduced[TT, 
                V /. vlimits[[2]]] + 
            0.5\ pReduced[TT, 
                V /. vlimits[[3]]]; 
        V1 = V /. \(Solve[pReduced[TT, V] == midP, 
                V]\)[[1]]; 
        V2 = V /. \(Solve[pReduced[TT, V] == midP, 
                V]\)[[3]]; 
	
        volumes = 
          FindRoot[
            Evaluate[{area1 == area2 /. T -> TT, 
                pReduced[TT, V\_ℓ] == 
                  pReduced[TT, V\_g]}], {V\_ℓ, V1, 0, 1}, {V\_g, V2, 
              1, 1000}]; 
        Join[{p -> pReduced[TT, V\_ℓ /. volumes]}, 
          volumes]]\)

maxwellConstruction[0.95]

\(maxwellPlot[T_, Vmin_, Vmax_, Pmin_, Pmax_] := 
	
      Module[{values, pp, V1, V2}, 
		values = maxwellConstruction[T]; 
        pp = p /. values; V1 = V\_ℓ /. values; 
        V2 = V\_g /. values; 
		Show[Evaluate[plotPvsV[T]], 
          Graphics[Line[{{V1, pp}, {V2, pp}}]], 
          PlotRange -> {{Vmin, Vmax}, {Pmin, Pmax}}, 
          PlotLabel -> "Maxwell Construction", Frame -> True, 
          FrameLabel -> {\*""V/(Subscript[V, c])"", \*""P/(Subscript[P, c])""}, DisplayFunction -> $DisplayFunction];]\)

maxwellPlot[0.90, 0.5, 2.5, 0.4, 0.8]

Develop an interpolating function for the coexistence boundary

The boundary curve which defines the coexistence region is constructed by performing the Maxwell construction for a range of temperatures and constructing an interpolation function that represents the endpoints of the coexistence isobars.

points = 
        Table[{{(Subscript[V, ℓ]), p}, {(Subscript[V, g]), p}} /. maxwellConstruction[T], {T, 
            0.82, 0.99, 0.01}];

coexistencePoints = 
        Sort[Join[{{1.0, 1.0}}, Partition[Flatten[points], 2]]];

boundary = Interpolation[coexistencePoints];

Display the coexistence region

arrayPvsV = Table[plotPvsV[x], {x, 0.8, 1.2, 0.05}];
lines = 
        Table[Line[{{(Subscript[V, ℓ]), p}, {(Subscript[V, g]), p}} /. 
              maxwellConstruction[T]], {T, 0.85, 0.95, 0.05}];

boundaryPlot = 
        Plot[boundary[V], {V, 0.54, 3.7}, 
          PlotStyle -> {{Thickness[0.006], RGBColor[1, 0, 0]}}, 
          DisplayFunction -> Identity];

\(\(Show[arrayPvsV, boundaryPlot, Graphics[{RGBColor[1, 0, 0], lines}], 
        PlotRange -> {{0.4, 4}, {0, 2}}, 
        PlotLabel -> "Liquid-Vapor Transition", Frame -> True, 
        FrameLabel -> {\*""V/(Subscript[V, c])"", \*""P/(Subscript[P, c])""}, 
        DisplayFunction -> $DisplayFunction];\)\)

Within the liquid-vapor coexistence region, the van der Waals isotherms are replaced by isobaric lines between the saturated liquid and vapor volumes.

Problems

Dieterici equation

The Dieterici equation of state,
					\(TraditionalForm\`p\  = \ \(R\ T\)\/\(V - b\)\ Exp[\(-a\)\/\(R\ T\ V\)]\),
provides a useful alternative to the van der Waals equation for nonideal classical gases.  Note that here V is the molar volume, R is the molar gas constant, and a and b are constants which are properties of the particular gas.
a) Derive the expressions for the critical pressure, molar volume, and temperature ((Subscript[p, c]),  (Subscript[V, c]),  (Subscript[T, c])) in terms of a and b.  Evaluate the critical compressibility factor (((Subscript[p, c]) (Subscript[V, c]))/(R (Subscript[T, c]))) and show that it is in better agreement with the data for most gases than the corresponding result for the van der Waals equation.
b) Express the Dieterici equation in reduced form.
c) Deduce the virial coefficients for the Dieterici equation.  Compare the second virial coefficient with the van der Waals result.

Heat of vaporization

The heat of vaporization for a van der Waals fluid can be obtained by integrating the heat transfer along an isotherm using

\(Δ H\  = \ 
      T \(∫\_\(V\_ℓ\)\%\(V\_g\)\(\((∂S\/∂V)\)\_T\) dV\)\)

where (Subscript[V, ℓ]) and (Subscript[V, g]) are the molar volumes for the saturated liquid and vapor phases, as determined by the Maxwell construction.
a) Show that

\(Δ H\  = \ 
      T \(∫\_\(V\_ℓ\)\%\(V\_g\)\(\((∂p\/∂T)\)\_V\) dV\)\)

and evaluate using the van der Waals equation of state.
b) Compute ΔH for water at atmospheric pressure given that a = 0.565 Pa (m^6) (mole^(-2)) and b = 3.1×(10^(-5)) (m^3) (mole^(-1)).  Compare with experimental data.

Thermodynamics of van der Waals gas

a) Evaluate (Subscript[C, p])  -  (Subscript[C, V]) for the van der Waals equation of state.

b) Determine the relationships between T and V and between p and V that apply to isentropic transformations of a gas that obeys the van der Waals equation of state assuming that (Subscript[C, V]) is constant.