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.

2 thoughts on “Van der Waals equation of state: Maxwell construction

  1. Pingback: Helmholtz free energy of a van der Waals fluid | Physics pages

  2. Pingback: Van der Waals fluid at the critical point | Physics pages

Leave a Reply

Your email address will not be published. Required fields are marked *