# Van der Waals equation of state: Maxwell construction

Reference: Daniel V. Schroeder, An Introduction to Thermal Physics, (Addison-Wesley, 2000) – Problems 5.52 – 5.53.

To get a better feel for the van der Waals equation of state, we can look at how the Gibbs free energy varies with the volume and pressure. Schroeder derives an expression for ${G}$ in his equation 5.56, which I will just quote:

$\displaystyle G=-NkT\ln\left(V-Nb\right)+\frac{\left(NkT\right)\left(Nb\right)}{V-Nb}-\frac{2aN^{2}}{V}+c\left(T\right) \ \ \ \ \ (1)$

where ${c\left(T\right)}$ is a function of temperature only, and can be ignored in what follows since we’ll be specifying a particular temperature in each example, so ${c\left(T\right)}$ is effectively a constant in each case.

This equation is easier to work with if we rewrite it in terms of the reduced variables, defined earlier as ${v\equiv V/V_{c}}$, ${t\equiv T/T_{c}}$ and ${p\equiv P/P_{c}}$, with:

 $\displaystyle V_{c}$ $\displaystyle =$ $\displaystyle 3Nb\ \ \ \ \ (2)$ $\displaystyle T_{c}$ $\displaystyle =$ $\displaystyle \frac{8}{27}\frac{a}{bk}\ \ \ \ \ (3)$ $\displaystyle P_{c}$ $\displaystyle =$ $\displaystyle \frac{a}{27b^{2}} \ \ \ \ \ (4)$

If we divide both sides of 1 by ${NkT_{c}}$ and simplify using the above equations, we get

$\displaystyle G_{c}\equiv\frac{G}{NkT_{c}}=-t\ln\left(3v-1\right)+\frac{t}{3v-1}-\frac{9}{4v}+c^{\prime}\left(t\right) \ \ \ \ \ (5)$

where ${c^{\prime}\left(t\right)=c\left(T\right)/NkT_{c}-t\ln\left(Nb\right)}$ is another function that depends only on temperature, so can be ignored in what follows.

To plot ${G_{c}}$ as a function of reduced pressure, we can use the equation derived in the last post:

$\displaystyle p=\frac{8tv^{2}-9v+3}{\left(3v-1\right)v^{2}} \ \ \ \ \ (6)$

We can specify a temperature ${t}$ and then vary ${v}$, using it to calculate ${p}$ and the corresponding ${G_{c}}$ from the last two equations, and thus generate a plot of ${G_{c}}$ versus ${p}$. For ${t=0.95}$ the isotherm looks like this:

If we plot the Gibbs energy, we get:

As we reduce the volume, we move from right to left on the ${pv}$ plot, and the corresponding path on the Gibbs plot is the line coming up from the bottom. At point ${C}$, we can continue to decrease ${v}$, which causes ${G}$ to carry on along the same path above point ${C}$, until it reaches a sharp corner, at which point it decreases until it reaches point ${B}$, which is directly above (that is, it has the same pressure) as point ${C}$. However, note that ${B}$ does not have the minimum Gibbs energy for that pressure, since point ${C}$ is below it. Thus ${B}$ is unstable and does not correspond to a physical state of the substance. Back on the ${pv}$ plot, as we decrease ${v}$ below point ${B}$, the Gibbs path decreases until it reaches the lower sharp corner. This corner corresponds to the minimum between ${A}$ and ${B}$ on the ${pv}$ plot. Again, note that this isn’t the minimum Gibbs energy for that pressure, so this is also unstable. In fact the whole section of the ${pv}$ plot where the pressure decreases with decreasing volume (between the minimum between ${A}$ and ${B}$ and the maximum between ${B}$ and ${C}$) lies on a Gibbs path that is not the minimum energy, so this whole range of volumes is unstable and is not found in the real world.

As we continue to decrease ${v}$ from the minimum point down to point ${A}$, the Gibbs energy returns from the lower sharp corner until it reaches point ${A}$, where the Gibbs energy is equal to the value it had at point ${C}$. Thus points ${A}$ and ${C}$ have the same Gibbs energy at the same pressure, but point ${A}$ has a much lower volume than point ${C}$, so it corresponds to a condensed state (a liquid) while ${C}$ corresponds to a gas.

In fact, this phase transformation occurs for only one pressure, as we can see from the following argument. At constant temperature, the thermodynamic identity states that

$\displaystyle \left(\frac{\partial G}{\partial P}\right)_{T}=V \ \ \ \ \ (7)$

As we go around the triangular loop in the Gibbs plot above, the overall change in ${G}$ is zero, so if we integrate it around the loop, we must get zero:

$\displaystyle 0=\oint dG=\oint\left(\frac{\partial G}{\partial P}\right)_{T}dP=\oint V\;dP \ \ \ \ \ (8)$

If we divide this equation by ${V_{c}P_{c}}$ it changes nothing, since both these quantities depend only on ${T}$, which is held constant. Thus

$\displaystyle \oint v\;dp=0 \ \ \ \ \ (9)$

To do this integral directly requires inverting 6 to get ${v}$ as a function of ${p}$. However, this requires solving a cubic equation in ${v}$, which gives 3 different roots. Just for completeness, I did actually do this using Maple, and if we plot the result we get, for ${t=0.95}$:

The red, blue and green curves represent contributions from each of the three roots. As we go round the triangular loop above, we go from point ${A}$ to ${B}$ to ${C}$. Requiring this integral to be zero is equivalent (as Schroeder shows) to the area bounded by the vertical line and curve between ${B}$ and ${C}$ being equal to the area bounded by the vertical line and curve between ${A}$ and ${B}$. This is known as the Maxwell construction (yes, it’s the same Maxwell that devised Maxwell’s equations in electrodynamics).

The problem is then to find the value ${p_{0}}$ of ${p}$ that makes these two areas equal. Doing this requires moving the yellow line in the plot left or right until the areas are equal. Solving the problem analytically is a bit trickier, of course, since we’re dealing with a cubic equation.

It’s a lot easier to solve if we revert to the traditional ${pv}$ plot above. The area bounded by the vertical line and curve between ${B}$ and ${C}$ is then given by

$\displaystyle A_{BC}=\int_{v_{B}}^{v_{C}}p\;dv-p_{0}\left(v_{C}-v_{B}\right) \ \ \ \ \ (10)$

That is, it’s the area under the ${p}$ curve between ${v_{B}}$ and ${v_{C}}$ minus the area of the rectangle under the horizontal line between the same two volumes. By the same reasoning, the area in the ${vp}$ plot bounded by the vertical line and curve between ${A}$ and ${B}$ is

$\displaystyle A_{AB}=p_{0}\left(v_{B}-v_{A}\right)-\int_{v_{A}}^{v_{B}}p\;dv \ \ \ \ \ (11)$

Surprisingly, the integral turns out to be fairly simple, using Maple:

 $\displaystyle \int p\;dv$ $\displaystyle =$ $\displaystyle \int\frac{8tv^{2}-9v+3}{\left(3v-1\right)v^{2}}dv\ \ \ \ \ (12)$ $\displaystyle$ $\displaystyle =$ $\displaystyle \frac{3}{v}+\frac{8}{3}t\ln\left(3v-1\right) \ \ \ \ \ (13)$

The problem now becomes:

1. Specify the temperature.
2. Choose a starting value for the pressure ${p_{0}}$. This can be estimated from the ${pv}$ plot by judging by eye what value of ${p_{0}}$ makes the two areas roughly equal.
3. Solve 6 to find the 3 volumes ${v_{A}}$, ${v_{B}}$ and ${v_{C}}$ for this value of ${p_{0}}$. This requires solving a cubic equation in ${v}$.
4. Calculate the two areas using 10 and 11.
5. Adjust ${p_{0}}$ and repeat from step 3 until we get convergence.

I did this using a Maple procedure, with the best result (for ${t=0.95}$):

$\displaystyle p_{0}=0.8119 \ \ \ \ \ (14)$

As a check, we can compare this with the Gibbs energy plot above and note that the point ${A,C}$ does indeed occur at this pressure.

We can repeat the procedure for ${t=0.8}$ and get the result:

$\displaystyle p_{0}=0.3834 \ \ \ \ \ (15)$

The Gibbs plot for this case is:

Again, we can see that the crossover between the two phases occurs at around ${p_{0}=0.3834}$ so it looks like all is well.