The Maxwell Construction

(14 comments)

The van der Waals equation of state for a gas is: $$ \left(p + \frac{a}{V^2}\right) \left(V - b\right) = RT $$ in terms of the pressure, $p$, molar volume, $V$ and temperature, $T$. $a$ and $b$ are constants which depend on the gas, so it is often useful to recast this equation into reduced form: $$ \left(p_\mathrm{r} + \frac{3}{V_\mathrm{r}^2}\right)\left(V_\mathrm{r}-\frac{1}{3}\right) = \frac{8}{3}T_\mathrm{r}, $$ or equivalently $$ p_\mathrm{r} = \frac{8T_\mathrm{r}}{3V_\mathrm{r}-1} - \frac{3}{V_\mathrm{r}^2}, $$ where the reduced variables $p_\mathrm{r} = p/p_\mathrm{c}$, $V_\mathrm{r} = V/V_\mathrm{c}$ and $T_\mathrm{r} = T / T_\mathrm{c}$ in terms of the critical pressure, volume and temperature: $$ p_\mathrm{c} = \frac{a}{27b^2}, \;\; V_\mathrm{c} = 3b, \;\; k_\mathrm{B}T_\mathrm{c} = \frac{8a}{27b}. $$ Whilst the van der Waals equation does a better job than the ideal gas law of describing the properties of a real gas it suffers from an artefact for $T_\mathrm{r} < 1$ (that is, temperatures $T < T_\mathrm{c}$), as shown in the plot below.

enter image description here

The "loop" in the region between $V_\mathrm{r} \approx 0.6$ and $V_\mathrm{r} \approx 2.5$ is clearly unphysical for it would imply that a gas at a fixed pressure and temperature could have three different equilibrium volumes. In mathematical terms, this region corresponds to values of $p_\mathrm{r}$ where the cubic equation obtained by rearranging the reduced equation of state above, $$ 3p_\mathrm{r}V_\mathrm{r}^3 - (p_\mathrm{r} + 8T_\mathrm{r})V_\mathrm{r}^2 + 9_\mathrm{r} - 3 = 0, $$ has three real roots.

The Maxwell construction deals with this anomaly by recognising that this region corresponds to conditions where the liquid and gas phases coexist: that is, a gas compressed at constant temperature through these volumes will condense with the volume decreasing at constant pressure. Maxwell connected the high and low-volume states with a horizontal line at a fixed pressure such that the areas cut off above and below this line were equal:

enter image description here

The following program produces a pressure-volume plot for $\mathrm{CO_2}$ using the van der Waals equation of state with the Maxwell construction for $T < T_\mathrm{c}$. We first find the maximum and minimum of the "loop" in $p_\mathrm{r}$ and take our initial guess for the "cut-off" pressure, $p_\mathrm{r,0}$, to be the mean. The three volumes corresponding to this pressure: $V_\mathrm{r,min}$, $V_\mathrm{r,0}$ and $V_\mathrm{r,max}$ are obtained by finding the roots of the above cubic equation. We can then consider the exercise of finding a value of $V_\mathrm{r,0}$ which ensures that the area of the two loops to be the same to be equivalent to finding the solution of the equation $$ \int_{V_\mathrm{r,min}}^{V_\mathrm{r,max}} (p_\mathrm{r} - p_\mathrm{r,0})\,\mathrm{d}V_\mathrm{r} = 0. $$ That is, finding the root of the left hand side: this we can do using Newton's method with scipy.optimize.newton.

enter image description here

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad
from scipy.optimize import newton
from scipy.signal import argrelextrema

palette = iter(['#9b59b6', '#4c72b0', '#55a868', '#c44e52', '#dbc256'])

# Critical pressure, volume and temperature
# These values are for the van der Waals equation of state for CO2:
# (p - a/V^2)(V-b) = RT. Units: p is in Pa, Vc in m3/mol and T in K.
pc = 7.404e6
Vc = 1.28e-4
Tc = 304

def vdw(Tr, Vr):
    """Van der Waals equation of state.

    Return the reduced pressure from the reduced temperature and volume.

    """

    pr = 8*Tr/(3*Vr-1) - 3/Vr**2
    return pr


def vdw_maxwell(Tr, Vr):
    """Van der Waals equation of state with Maxwell construction.

    Return the reduced pressure from the reduced temperature and volume,
    applying the Maxwell construction correction to the unphysical region
    if necessary.

    """

    pr = vdw(Tr, Vr)
    if Tr >= 1:
        # No unphysical region above the critical temperature.
        return pr

    if min(pr) < 0:
         raise ValueError('Negative pressure results from van der Waals'
                         ' equation of state with Tr = {} K.'.format(Tr))

    # Initial guess for the position of the Maxwell construction line:
    # the volume corresponding to the mean pressure between the minimum and
    # maximum in reduced pressure, pr.
    iprmin = argrelextrema(pr, np.less)
    iprmax = argrelextrema(pr, np.greater)
    Vr0 = np.mean([Vr[iprmin], Vr[iprmax]])

    def get_Vlims(pr0):
        """Solve the inverted van der Waals equation for reduced volume.

        Return the lowest and highest reduced volumes such that the reduced
        pressure is pr0. It only makes sense to call this function for
        T<Tc, ie below the critical temperature where there are three roots.

        """

        eos = np.poly1d( (3*pr0, -(pr0+8*Tr), 9, -3) )
        roots = eos.r
        roots.sort()
        Vrmin, _, Vrmax = roots
        return Vrmin, Vrmax

    def get_area_difference(Vr0):
        """Return the difference in areas of the van der Waals loops.

        Return the difference between the areas of the loops from Vr0 to Vrmax
        and from Vrmin to Vr0 where the reduced pressure from the van der Waals
        equation is the same at Vrmin, Vr0 and Vrmax. This difference is zero
        when the straight line joining Vrmin and Vrmax at pr0 is the Maxwell
        construction.

        """

        pr0 = vdw(Tr, Vr0)
        Vrmin, Vrmax = get_Vlims(pr0)
        return quad(lambda vr: vdw(Tr, vr) - pr0, Vrmin, Vrmax)[0]

    # Root finding by Newton's method determines Vr0 corresponding to
    # equal loop areas for the Maxwell construction.
    Vr0 = newton(get_area_difference, Vr0)
    pr0 = vdw(Tr, Vr0)
    Vrmin, Vrmax = get_Vlims(pr0)

    # Set the pressure in the Maxwell construction region to constant pr0.
    pr[(Vr >= Vrmin) & (Vr <= Vrmax)] = pr0
    return pr

Vr = np.linspace(0.5, 3, 500)

def plot_pV(T):
    Tr = T / Tc
    c = next(palette)
    ax.plot(Vr, vdw(Tr, Vr), lw=2, alpha=0.3, color=c)
    ax.plot(Vr, vdw_maxwell(Tr, Vr), lw=2, color=c, label='{:.2f}'.format(Tr))

fig, ax = plt.subplots()

for T in range(270, 320, 10):
    plot_pV(T)

ax.set_xlim(0.4, 3)
ax.set_xlabel('Reduced volume')
ax.set_ylim(0, 1.6)
ax.set_ylabel('Reduced pressure')
ax.legend(title='Reduced temperature')

plt.show()
Current rating: 4.4

Comments

Comments are pre-moderated. Please be patient and your comment will appear soon.

Chuangde 5 years ago

Hello,it's really a good job. If i want to use another EOS, like Peng-Robison eos or Carnahan-Starling eos, to determine the co-existence density of gas and liquid, I just replace the vdw eos by other eos,isn't right? Thank you!

Link | Reply
Current rating: 5

christian 5 years ago

In principle, yes, but you might have a bit more work to do with more complex equations of state. Let me know how you get on!

Link | Reply
Current rating: 4.8

Francis 3 years, 11 months ago

Great work. Thank you and keep it up.

Link | Reply
Currently unrated

RN 3 years, 6 months ago

Hello.
I really like this work. I would like to extend the plot to lower temperature. Unfortunately, when try to get lower values, the codes is not working. How could we solve this ?
Thank you by advance.

Link | Reply
Currently unrated

christian 3 years, 6 months ago

Thank you – can you be a bit more specific about how it's failing at low T? At very low T the substance is not likely to be well described by a gaseous equation of state, and I wonder if the loops go a bit crazy.

Link | Reply
Currently unrated

Robert Topper 3 years, 2 months ago

Nice project!

If the pressure ever dips below 0, the method fails. This is a limitation of the van der Waals equation.

If the temperature gets too low, the first term in the van der Waals equation (which is always >0) would become smaller than the second term (which is always < 0). This will result in the pressure being (unphysically) negative. The code could monitor for this situation by checking the sign of the pressure at the minimum and making sure that it's > 0 . If <0, an error message should probably be generated.

Link | Reply
Currently unrated

christian 3 years, 2 months ago

Yup - that would be useful. I've added a check in the vdw_maxwell function.

Thanks for getting in touch,
Christian

Link | Reply
Current rating: 5

Robert Topper 3 years, 2 months ago

Glad to help!

Link | Reply
Currently unrated

Nonso 3 years, 3 months ago

Hello Christian,

Thank you for the article!
Basically, helped me understand what a maxwell construction is and why it's actually needed.
I'm taking a Thermodynamics class this semester and I literally didn't understand the explanation from class.

A little thing I wanted to point out though, the pressure dependence in the VDW equation has a plus sign, not a minus. It doesn't affect anything cos I see you used the correct reduced VDW equation but it got me confused at first.

Link | Reply
Currently unrated

christian 3 years, 3 months ago

You're absolutely right: thank you for letting me know about this error. I've fixed it now.
Glad you found the article useful!
Christian

Link | Reply
Current rating: 5

LarP 3 years ago

Thank you very much for your script, I found it very useful. The only thing I can't do is printing the value of the volume that makes the two areas equal, as well as the associated pressure. Could you help me with it?

Thanks again!

Link | Reply
Currently unrated

christian 3 years ago

Glad you found it interesting!

Isn't the area you want just the integral of pr from Vr0 to Vrmax (minus the rectangle under pr0)?

A = quad(lambda vr: vdw(Tr, vr) - pr0, Vr0, Vrmax)[0]

(after the line Vrmin, Vrmax = get_Vlims(pr0) has obtained the converged Maxwell construction).

Link | Reply
Current rating: 5

Abdoulaye 1 year, 5 months ago

Hi,

Are you using the actual derivative or secant method when iterating throughout Newton method?

Link | Reply
Currently unrated

christian 1 year, 5 months ago

Hi Abdoulaye,
As you can see from the code, I do not calculate the derivative explicitly and pass it to `newton`, so this library function will use the secant method.
Best wishes,
Christian

Link | Reply
Currently unrated

New Comment

required

required (not published)

optional

required