Virial Coefficient for Classical Gases

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

Introduction

The virial expansion for the mechanical equation of state can be expressed in the form

((p V)/(N (Subscript[k, B]) T))  =  
      1  +  (Subscript[B, 2]) (N/V)  +  ⋯

where

\(\(\(B\_2\)\(\ \)\(=\)\(\ \)\(\(-2\)\ π\ \(∫\_0\%∞\((Exp[\(-\(v[r]\/\(\(k\_\(\(B\)\(\ \)\)\) T\)\)\)]\  - \ 
                1)\)\ \(r\^2\) dr\)\)\(\ 		\)\)\)

is the second virial coefficient due to an intermolecular potential.  [By convention, the first virial coefficient is unity.]   The second virial coefficient has dimensions of volume and is expected to be comparable to the volume, ((Subscript[r, 0])^3), over which the interaction is significant.  Its temperature dependence is most clearly expressed in terms of a reduced temperature variable, τ  =  (Subscript[k, B]) T/(Subscript[v, 0]) where (Subscript[v, 0]) is a suitable measure of the potential depth.

In this notebook, we study second virial coefficient for simple classical gases using a realistic intermolecular  potential.  The intermolecular interaction features a short-range repulsion and a long-range attraction.  The short-range repulsion arises when the electron clouds begin to overlap and there are no shared molecular orbitals which would lead to binding.  Hence, we consider gases to be simple when their constituents do not readily bind.  The long-range attraction arises from the van der Waals force due to mutually induced electric dipole moments; that mechanism leads to an (r^(-6)) potential at large distances.  Therefore, the Lennard-Jones potential,

v(r)  =  
      (Subscript[v, 0]) ( (((((Subscript[r, 0])/r)))^12)  -  
            2  (((((Subscript[r, 0])/r)))^6)  )

provides a realistic empirical model of the intermolecular potential for simple gases.  This potential possesses a minimum at r = (Subscript[r, 0]) with depth -(Subscript[v, 0]), the expected (r^(-6)) asymptotic tail, and a steep small-r repulsion.

To illustrate the roles of various aspects of the potential, we first consider the virial coefficients produced by hard-core and square-well potentials.  We also consider a piecewise constant potential consisting of a hard core followed by a square well, with parameters chosen to approximate the Lennard-Jones potential.  This simple potential is shown to lead to the van der Waals equation of state using suitable approximations and to produce a virial coefficient whose temperature dependence is quite similar to that of the Lennard-Jones potential.

Glossary:
	p	= pressure
	V	= volume
	𝒩	= particle number
	(Subscript[k, B])	= Boltzmann constant
	T	= temperature
	v	= volume per particle (V/N)
	(Subscript[v, 0])	= potential depth parameter
	(Subscript[r, 0])	= potential range parameter 
	(Subscript[r, c])	= hard core radius
	(Subscript[B, 2])	= second virial coefficient
	τ	= reduced temperature ((Subscript[k, B]) T / (Subscript[v, 0]))

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`"];
Needs["Miscellaneous`PhysicalConstants`"];
Needs["Miscellaneous`Units`"];
Needs["Graphics`Graphics`"]

Symbolize[(Subscript[k, B])];Symbolize[(Subscript[B, 2])];Symbolize[(Subscript[v, 0])];Symbolize[(Subscript[r, 0])];Symbolize[(Subscript[r, c])];

Definitions for virial coefficients

mayer[potential_] := 
      Exp[-(potential/((Subscript[k, B]) T))] - 1

B2[T_, 
        potential_] := -2 π (∫_0^∞) (r^2) mayer[potential] \[DifferentialD]r

\(B2Plot[virialCoefficient_, opts___Rule] := 
      LogLinearPlot[virialCoefficient, {τ, 1, 100}, opts, 
        PlotRange -> {All, {\(-5\), 2}}, GridLines -> Automatic, 
        Frame -> True, 
        FrameLabel -> {\*""(Subscript[k, B])T/(Subscript[v, 0])"", \*""(Subscript[B, 2])/((Subscript[r, 0])^3)""}, 
        PlotLabel -> "Virial Coefficient"]\)

Error messages

The numerical integrations sometimes produce error messages at the extremes of their range.  We choose to disable these messages because we know that the integrands are well-behaved and the warnings harmless.

Off[NIntegrate::"inum", General::"unfl", 
      NIntegrate::"slwcon"]

Simple potentials

Hard core

B2hardcore[(Subscript[r, c]) : _] = -2 π (∫_0^((Subscript[r, c])))(r^2) (-1) \[DifferentialD]r

A repulsive hard-core potential produces a positive temperature-independent virial coefficient.

Simple square well

B2square[(Subscript[v, 0]) : _, (Subscript[r, 0]) : _] = -2 π (∫_0^((Subscript[r, 0])))(r^2) mayer[-(Subscript[v, 0])] \[DifferentialD]r

squarePlot = 
        B2Plot[B2square[(Subscript[v, 0]), (Subscript[r, 0])]/((Subscript[r, 0])^3) /. {T -> ((τ (Subscript[v, 0]))/(Subscript[k, B]))}];

An attractive square-well potential produces a negative virial coefficient which is strongest at low temperature but less effective at high temperature.

Lennard-Jones potential

The Lennard-Jones "6-12" potential contains a short-range repulsion and long-range attraction.  The potential minimum is located at r = (Subscript[r, 0]) with depth -(Subscript[v, 0]).  Thus, since the virial coefficient for this model depends upon only two parameters, it should be possible to plot (Subscript[B, 2])[τ]/((Subscript[r, 0])^3) for all simple gases together, where τ  =  (Subscript[k, B]) T/(Subscript[v, 0]) is the temperature in units of the potential depth and where (Subscript[B, 2]) is expressed in units of molecular volume ((Subscript[r, 0])^3).  In fact, except for very light gases where quantum effects are more important, this model provides a remarkably accurate description of the data over a wide range of temperatures.

Numerical evaluation of virial coefficient

lennard[r_, (Subscript[v, 0]) : _, (Subscript[r, 0]) : _] = 
        (Subscript[v, 0]) ((((((Subscript[r, 0])/r)))^12) - 2 (((((Subscript[r, 0])/r)))^6));

B2lennard[τ_] := -2 π NIntegrate[
          (r^2) (mayer[lennard[r, (Subscript[v, 0]), (Subscript[r, 0])]] /. {(Subscript[r, 0]) -> 1, 
                  T -> ((τ (Subscript[v, 0]))/(Subscript[k, B]))}), {r, 
            0, ∞}]

lennardPlot = B2Plot[B2lennard[τ]];

The negative contribution for low temperature evidently originates from the intermediate-ranged attractive part of the potential.  At higher temperatures the positive contribution due to the strong short range repulsion produces a sign change and eventually dominates when the temperature is much larger than the attractive well depth.

Compare with hard core repulsion plus square well attraction

The main features of the Lennard-Jones potential can be approximated using a piecewise constant potential as shown.

Clear[corePlusWell];
corePlusWell[r_ /; r < 0.85] = 1000. ;
corePlusWell[
          r_ /; 0.85 ≤ r ≤ 
              1.5] = -0.75;
corePlusWell[r_ /; r > 1.5] = 0;

\(\(potplot = 
        Plot[{lennard[r, 1, 1], corePlusWell[r]}, {r, 0.1, 5}, 
          PlotRange -> {{0, 2}, {\(-1.5\), 2}}, 
          PlotStyle -> {{}, Dashing[{0.02, 0.02}]}, Frame -> True, 
          FrameLabel -> {\*""r/(Subscript[r, 0])"", \*""ν/(Subscript[ν, 0])""}, PlotLabel -> "Potential", 
          DisplayFunction -> Identity];\)\)

B2CorePlusWell[(Subscript[v, 0]) : _, (Subscript[r, 0]) : _, (Subscript[r, c]) : _] = -2 π ((∫_0^((Subscript[r, c])))(r^2) (-1) \[DifferentialD]r + (∫_((Subscript[r, c]))^((Subscript[r, 0])))(r^2) mayer[-(Subscript[v, 0])] \[DifferentialD]r) // 
          Collect[#1, Exp[((Subscript[v, 0])/((Subscript[k, B]) T))]] & // Simplify

corePlusWellPlot = 
        B2Plot[B2CorePlusWell[(Subscript[v, 0]), (Subscript[r, 0]), (Subscript[r, c])] /. {(Subscript[v, 0]) -> 0.75, 
              (Subscript[r, 0]) -> 1.5, (Subscript[r, c]) -> 0.85, T -> (τ/(Subscript[k, B]))}, 
          PlotStyle -> Dashing[{0.02, 0.02}], 
          DisplayFunction -> Identity];

bplot = 
        Show[{lennardPlot, corePlusWellPlot}, 
          DisplayFunction -> Identity];

Show[GraphicsArray[{potplot, bplot}], 
        PlotLabel -> "Comparison between Lennard-Jones and Toy \
Model", GraphicsSpacing -> {-0.1, 0}, 
        DisplayFunction -> $DisplayFunction];

The approximate potential produces a virial coefficient that has the same qualitative features as the more realistic potential, but is stiffer at high temperatures and more attractive at low temperatures.

Derivation of van der Waals equation

The van der Waals equation of state can be derived as an approximation to the virial expansion for the hard-core + square well potential valid for high temperatures relative to (Subscript[v, 0]).

bHighT = B2CorePlusWell[(Subscript[v, 0]), (Subscript[r, 0]), (Subscript[r, c])]  // 
          Series[#1, {(Subscript[v, 0]), 0, 1}] & // Simplify

vdWaalsParameters = {a -> -(Subscript[k, B]) Coefficient[bHighT, (1/T)], 
          b -> Coefficient[bHighT, (1/T), 0]}  // Simplify

eq1 = 
      (Subscript[k, B]) T == ((Subscript[k, B]) T /. Solve[
                    p == (((N (Subscript[k, B]) T) (1 + ((N (b - (a/((Subscript[k, B]) T))))/V)))/V), 
                    T][[1]] /. {V -> N v}  // Simplify)

\(eq1 /. {x_\/\(b + v\) -> \(x\ \((1 - b\/v)\)\)\/v}\  // 
      Simplify\)

Thus, we have obtained the van der Waals equation in its usual form and can relate its parameters to those of the intermolecular potential.  The parameter b is the volume per molecule excluded by the hard core, which is one-half the volume of a sphere of radius  (Subscript[r, c]) because there are two molecules present for the interaction.  Similarly, the parameter a is the absolute value of the volume integral per molecule of the attractive portion of the potential.  Note that the parameters  (Subscript[r, 0]) and (Subscript[v, 0]) cannot be determined separately since only the volume integral enters at this level of approximation.

Examples

Argon with Lennard-Jones potential

argonValues = {(Subscript[r, 0]) -> ((3.82 Meter)/(10^10)), 
        (Subscript[v, 0]) -> 120 Kelvin BoltzmannConstant}

Convert[(Subscript[v, 0]) /. argonValues, ElectronVolt]

Clear[B2argon];
B2argon[TT_] := 
      Module[{tauArgon}, 
          tauArgon = (((Subscript[k, B]) TT)/(Subscript[v, 0])) /. (Subscript[v, 0]) -> (Subscript[v, 0]) /. 
                argonValues /. {(Subscript[k, B]) -> BoltzmannConstant}; 
          B2lennard[tauArgon] ((((Subscript[r, 0]) /. argonValues))^3)]  // 
        N

B2argon[1000 Kelvin]

B2argon[TKelvin /. TKelvin -> TTT Kelvin]/(Meter^3) /. 
      TTT -> 100

\(\(argonPlot = 
        LogLogListPlot[
          Table[{10\^x, 
              6.022\ 10\^23\ Abs[
                  B2argon[TKelvin /. TKelvin -> 10\^x\ Kelvin]\/Meter\^3]}, {x, 1, 4, 0.01}], PlotJoined -> True, 
          PlotStyle -> AbsoluteThickness[2], PlotPoints -> 50, 
          GridLines -> Automatic, Frame -> True, 
          FrameLabel -> {"T (Kelvin)", \*""(Subscript[B, 2])  ((m^3)/Mole)""}, 
          PlotLabel -> "Virial Coefficient for Argon", 
          DisplayFunction -> Identity];\)\)

Show[argonPlot, PlotRange -> {{1, 4}, {-7, -1}}, 
        DisplayFunction -> $DisplayFunction];

Problems

ν[r]   =  
      (A/(r^m))  -  (B/(r^n))