Equation of State using Generalized van der Waals Theory

This notebook illustrates the generalized van der Waals theory (gvdW) for the equation of state for interacting particles. Based on the lecture notes, Properties of Molecular Fluids in Equilibrium by Sture Nordholm.



In [1]:

    
from __future__ import division, unicode_literals, print_function
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import quad
from ipywidgets import interact
from math import sqrt, pi
plt.rcParams.update({'font.size': 14, 'figure.figsize': [6.0, 4.0]})

Pair potentials

The particles are here assumed to interact via a Lennard-Jones and a screened Coulomb potential,

$$ \beta w(r) = \frac{\lambda_B z^2}{r} e^{-r/\lambda_D} + 4\beta \varepsilon_{LJ} \left ( \left ( \frac{\sigma}{r}\right )^{12} - \left ( \frac{\sigma}{r}\right )^{6}\right ) $$

where $\lambda_B$ and $\lambda_D$ are the Bjerrum and Debye lengths, respectively. Any potential may in principle be given and must return the energy in units of $k_BT=\beta^{-1}$.



In [2]:

    
# Debye-Huckel
def wDH(r, z, D):
    lB = 7.0 # Bjerrum length, angstroms
    return lB * z**2 * np.exp(-r/D) / r

# Lennard-Jones
def wLJ(r, eps, sigma):
    return 4 * eps * ( (sigma/r)**12 - (sigma/r)**6 )

# Total potential
def w(r, z, D, eps, sigma):
    return wDH(r, z, D) + wLJ(r, eps, sigma)

Interaction parameter

Here we integrate the above pair potential to get the average interaction energy per particle, assuming that the pair correlation function, $g(r)$, can be described by a simple step function, zero when $r<\sigma$, unity otherwise: $$ \hat{a} = -\frac{1}{2} \int_{\sigma}^{\infty} 4\pi w(r) r^2 dr $$

In this Notebook we simply do the integration numerically so that we can use arbitrary pair potentials. From this we calculate the pressure, $p$, versus density, $n$, using,

$$ \beta p_{gvdW} = \frac{1}{v-v_0} - \frac{\hat{a}}{v^2} $$

where $v=1/n$ and $v_0=2\pi\sigma^3/3$ is the particle volume. For reference we'll also plot EOS for an ideal system (van 't Hoff), $\beta p_{ideal}=n$, where $\beta = 1/k_BT$.



In [3]:

    
def ahat(z, D, eps, sigma):
    return -2*pi*quad(lambda r: w(r, z, D, eps, sigma)*r**2, sigma, np.infty, limit=50)[0]

def Pideal(n):
    return n

def Pgvdw(n, z, D, eps, sigma):
    v0 = 2*pi*sigma**3 / 3
    v  = 1 / n
    a  = ahat(z, D, eps, sigma)
    return 1/(v-v0) - a/v**2

Interactive EOS plot



In [4]:

    
def plot_EOS( eps=1.0, sigma=4.0, z=0.0, Cs=0.3 ):
    D = 3.04/sqrt(Cs)
    plt.plot(n, Pideal(n),  'k--', label='ideal', lw=2)
    plt.plot(n, Pgvdw(n, z, D, eps, sigma),  'r-', label='gvdW', lw=2)
    plt.title('Equation of State')
    plt.xlabel(r'Number density, $n$')
    plt.ylabel(r'Pressure, $\beta p$')
    plt.legend(loc=0, frameon=False)
    
n = np.linspace(1e-7, 6e-3, 100)
    
i = interact(plot_EOS,
             eps=(0.0, 10.0, 0.1), sigma=(0, 5, 0.1), z=(0.0, 10, 1.0), Cs=(1e-3, 1.0, 0.1) )
 
i.widget.children[0].description=r'$\beta\varepsilon_{LJ}$'
i.widget.children[1].description=r'$\sigma_{LJ}$'
i.widget.children[2].description=r'$z$'
i.widget.children[3].description=r'$c_s$ (M)'