Liquid State Theory Statistical Mechanics of Liquids and Complex Fluids George Jackson Department of Chemical Engineering Imperial College London South Kensington Campus London SW7 2AZ
[email protected] http://www3.imperial.ac.uk/people/g.jackson
Summary A molecular description of matter using statistical mechanical theories and computer simulation is the key to understanding and predicting the properties of dense fluids and materials. Fluid systems form an integral part of our modern lifestyle from the use of simple solvents in chemical processing to the design of opto-electronic devices with liquid crystalline and polymeric materials. As well as being of inherent scientific interest, the link between the microscopic interactions between small numbers of molecules and the macroscopic properties of the bulk system comprising them is of great industrial relevance. In this part of the course the main theoretical achievements in liquid state theory will be outlined. This with range from the original approach of van der Waals (1873), through the use of correlation functions to describe the properties of fluids (1950s-1980s), to the more sophisticated perturbation theories (1970sdate). The most recent developments in statistical mechanical theories which are used to probe the bulk or interfacial properties of fluids will be reviewed making contact with the description of real systems wherever possible. The focus will be on the structural and thermodynamic properties, with an emphasis on phase equilibria. 2
G. Jackson – Imperial College London
Synopsis of Lectures 1. Equilibrium Structure: 1.1 particle densities and distribution functions 1.2 pair distribution function 2. Thermodynamics: 2.1 using the language of distribution functions 2.2 internal energy 2.3 pressure virial (as opposed to virial expansion) 3. Perturbation Theories: 3.1 van der Waals (1873) 3.2 Zwanzig (1954) 3.3 Carnahan and Starling and the hard-sphere model (1969) 3.4 Wertheim/SAFT (1980-1990s) for associating and non-spherical molecules 3.5 Examples of the description of real complex fluids and mixtures (from water to polystyrene solutions) 3
G. Jackson – Imperial College London
Recap of Statistical Mechanical Concepts Total energy of classical system of N indistinguishable particles: rN r N rN 1 N 2 E (r , p ) = pi + U (r ) ∑ 2m i =1 kinetic + potential energy
(1)
Canonical partition function: QNVT
rN r N rN r N 1 = ...∫ exp(− β E (r , p ))d r d p where β = 1 /(kT ) 3N ∫ N!h rN rN 1 ( ( ) ) = exp − β integrating over momentum (2) U r d r 3N ∫ N!Λ r r Z = NVT3 N configurational integral Z NVT = ∫ exp(− β U (r N ))d r N N!Λ
Helmholtz free energy: r rr r r r N ≡ r1 r2 r3... rN 4
positions
F = −kT ln QNVT
(3)
r r r r r p N ≡ p1 p2 p3... pN momenta
G. Jackson – Imperial College London
Ensemble average of variable in the canonical ensemble: X NVT = ∑ X i Pi i
rN r N rN r N rN r N ...∫ X (r , p )exp(− β E (r , p ))d r d p ∫ classically = rN r N rN r N ∫ ...∫ exp(− β E (r , p ))d r d p rN r N (N ) rN r N rN r N = ∫ ...∫ X (r , p )P (r , p )d r d p (4) (N ) rN r N Pi ∝ P (r , p ) probability of finding a system in a state i
A probability density in the canonical ensemble can thus be defined (cf. equation 2) as rN r N exp(− β E (r , p )) (N ) rN r N P (r , p ) = rN r N rN r N ∫ ...∫ exp(− β E (r , p ))d r d p r r 1 exp(− β E (r N , p N )) (5) = 3N N h QNVT G. Jackson – Imperial College London ! 5
From the definition of the probability distribution function given in equation 5 it follows that rN r N r r r r (N ) rN r N (N ) rN r N P (r , p )d r d p = P (r , p )d r1 d p1... d rN d pN is the probability of finding the system in the phase space rN r N rN r N element d r d p about configuration r p , i.e., particle 1 r r d p in volume element d r1 and momentum element ... 1 r r particle N in volume element d rN and momentum d pN . r r r4 , p4
r r r2 , p2
r r r7 , p7
r r r12 , p12
r r r10 , p10
r r r8 , p8
r r r1 , p1
r r r3 , p3 r r r11 , p11
r r r6 , p6 6
r r r9 , p9
r r rN , pN
r r rn , pn
System configuration rN r N r p (state i)
r r r5 , p5 G. Jackson – Imperial College London
The momentum can again be integrated out into a kinetic contribution to define the probability density solely in terms of the positions (cf. equations 2 and 5) as rN rN exp(− β U (r )) (N ) rN (N ) rN r N P (r ) = ∫ P (r , p )d p = (6) rN rN exp(− β U (r ))d r ∫ r r r r r P ( N ) (r N )d r N = P ( N ) (r N )d r1 ... d rN is now the probability of rN rN finding the system in the volume element d r about r , r r i.e., particle 1 in d r1 ... particle N in d rN , irrespective of momenta. r r r2
r4
r r7
r r12
System configuration
r r10
r r1 r r8
r r9
r r3
r rn
r r11
r r6 7
r rN
r rN
r r5
G. Jackson – Imperial College London
For applications involving identical particles it is convenient (N ) rN to define a generic probability density ρ (r ) such that rN r r (N ) rN (N ) rN ρ (r )d r = ρ (r )d r1 ... d rN is the probability of finding rN rN the system in the volume element d r about r , irrespective of the particle labels. As there are N! possible permutations rN r r ( ( ) ) exp − β U r N N (N ) (N ) (7) ρ (r ) = N ! P (r ) = N ! rN rN ∫ exp(− β U (r ))d r
System configuration r rN regardless of particle labels. 8
G. Jackson – Imperial College London
1. Equilibrium Structure 1.1 Particle Densities and Distribution Functions (N ) rN r N (N ) rN The full probability densities P (r , p ) and ρ (r ) are high-order functions ~ O(N). The properties of the system can however be expressed in terms of lower-order or reduced particle densities and distribution functions. (n) n The reduced generic particle density ρ N (r ) is proportional
r
to the probability of observing any set of n particles in the r configuration r n and is defined by integrating out the coordinates of the remaining N − n particles r N −n N! (n) r n (N ) rN ρ N (r ) = P (r )d r ∫ (N − n )! (8) rN r N −n N ! ∫ exp(− β U (r ))d r = (N − n )! ∫ exp(− β U (rr N ))d rr N
9
G. Jackson – Imperial College London
The factor N ! / (N − n )! is included in equation 8 to take into account the number of ways that n particles can be chosen from the total N in the system. r r4
r r2
rn r Configuration of reduced set of n particles.
r r1 r r3
…
r rn
The other N − n particle positions are averaged out.
The reduced generic particle density is normalised such that
∫ρ 10
(n) N
rn rn (r )d r =
N! (N − n )!
(9) G. Jackson – Imperial College London
The single-particle density (n = 1) of a uniform fluid is simply the number density r N (1) r (10) ρ (r ) d r = = ρ V r (1) r ( ) ρ r d r =N . and is normalised as ∫ When the particles are uncorrelated (e.g., an ideal gas) the n-body particle density is the product of the single-particle (n) r n (1) r densities, ρ (r ) = ∏ ρ (ri ) = ρ n , the last equality in the case of a homogeneous fluid. An n-body distribution function can then be defined in terms of the reduced generic particle density to quantify the extent by which the structure of the system deviates from complete randomness (particle correlations): (n) r n (n) r n r ρ (r ) ρ N (r ) (11) g N( n ) (r n ) = n N = n ρ (1) r ( ) ρ r ∏ i 11
i =1
G. Jackson – Imperial College London
1.2 Pair Distribution Function
The pair distribution function (n = 2) of a uniform fluid can be obtained at once from equations 8 and 11 as rN r N −2 ( 2) r 2 ρ N (r ) N (N − 1) ∫ exp(− β U (r ))d r (2) r r g N (r1 , r2 ) = = (12) 2 2 Z NVT ρ ρ When the fluid is isotropic the pair distribution function r r depends only on the relative separation r12 = r2 − r1 . The probability of finding any one of the N − 1 particles at a relative distance from a given particle is proportional to r rN r N −2 d r1 ∫ exp(− β U (r ))d r r 1 ∫ (2) r (2) r r ρ N (r12 ) = ∫ ρ N (r1 , r12 ) d r1 = (N − 1) N Z NVT (13) rN r N −2 exp (− β U (r ))d r ∫ = (N − 1)V Z NVT 12 G. Jackson – Imperial College London
In the ideal (uncorrelated) limit r N −1 N −1 lim ρ N( 2 ) (r12 ) = =ρ N →∞ V V
(14)
A radial pair distribution function can thus be defined as (2) r r ρ (r ) n( r ) g (r12 ) = g N( 2 ) (r12 ) = N 12 = ideal 12 (15) ρ n ( r12 ) a form commonly employed in molecular simulation. The forms of the function usually used in the theory of radiation scattering are given in terms of appropriate δ functions, which follow directly from equation (12): 1 N N r r r r (2) r r g N (r1 , r2 ) = 2 ∑∑ δ (r1 − ri )δ (r2 − rj )
ρ
1
i =1 j =1
1 g ( r12 ) = ρ N 13
(16)
r r r ∑∑ δ (r12 − rj + ri ) N
N
i =1 j =1
G. Jackson – Imperial College London
The Fourier transform of the pair distribution function r r r r S ( k ) = 1 + ρ ∫ g ( r12 ) exp − ik ⋅ r12 d r12
(
)
(17)
is the so-called static structure factor which can be determined experimentally from neutron or x-ray scattering. The inverse transform gives the distribution function: r r r r r 1 (18) g ( r12 ) = S ( k ) − 1 exp ik ⋅ r12 d k 3 ∫ ρ (2π )
[
] (
)
1st shell
Structure of a fluid
Radial pair distribution function of liquid argon from neutron scattering. See Yarnell et al., PRA (1973); Hansen and McDonald (2006). 14
2nd shell 3rd shell
G. Jackson – Imperial College London
Characterising the Structure of a Material
gas
liquid
solid long-range order
short-range order
Radial pair distribution function of N = 500 Lennard-Jones particles at different densities. See Glotzer et al. (2006). 15
G. Jackson – Imperial College London
2. Thermodynamic Properties 2.1 Using the language of the distribution functions developed in Section 1.
Let us consider a uniform fluid for which the potential energy can be expressed as a sum of pair-wise additive contributions (e.g., hard-sphere, square-well, Lennard-Jones models etc.): N N rN U (r ) = ∑∑ u2 ( rij )
(19)
i =1 j >i
For simplicity we consider intermolecular potentials which have radial symmetry. When particles interact through pairwise-additive forces, the thermodynamic properties can be expressed in terms of integrals over the pair distribution function. 16
G. Jackson – Imperial College London
2.2 Internal Energy
The total energy of the system involves the kinetic (ideal) and configurational contributions (20) E = U ideal + U U ideal =
3 NkT 2
kinetic energy
The configurational part which is due to the interparticle interactions can be expressed as an ensemble average in the usual way in terms of the N-body probability distribution function: rN rN rN U (r )exp(− β U (r ))d r ∫ U = U NVT = rN rN ( ( ) ) exp β − U r d r ∫ (21) rN (N ) rN rN = ∫ U (r )P (r )d r 17
G. Jackson – Imperial College London
In terms of the pair-wise intermolecular potentials (cf. equation 19) the configurational energy (21) can be written as ⎛ N N ⎞ rN rN u2 ( rij ) ⎟⎟ exp(− β U (r ))d r ∫ ⎜⎜⎝ ∑∑ i =1 j >i ⎠ (22) U= Z NVT The double sum over i and j gives rise to N ( N − 1) / 2 pair terms which yield the same results on integration. We therefore focus on a given pair (in this case 12). rN rN ( ) ( ( ) ) u r U r d r exp β − 1 2 12 ∫ U = N ( N − 1) Z NVT 2 rN r N −2 ⎤ ⎡ ( ( ) ) exp β U r d r − r r 1 ∫ ⎢ ⎥ ( ) d r1 d r2 (23) = ∫∫ u2 r12 N ( N − 1) 2 Z NVT ⎢⎣ ⎥⎦ 2 (2) r r ρ g N (r1 , r2 ) 18
G. Jackson – Imperial College London
The integral in equation 23 over all N particles has been split into integrations over 1 and 2 and over the other N − 2 particles to introduce the previously defined pair distribution function (cf. equation 12): rN rN ( ) ( ( ) ) u r exp − β U r d r 1 2 12 ∫ U = N ( N − 1) (24) 2 Z NVT r r 1 2 (2) r r = ρ ∫∫ u2 (r12 )g N (r1 , r2 ) d r1 d r2 2 r r 1 N2 (2) r = u (r )g N (r12 ) d r1 d r12 2 ∫∫ 2 12 2V r 1 N2 (2) r = u2 (r12 )g N (r12 ) d r12 ∫ 2 V = 2πNρ ∫ u2 (r12 )g (r12 ) r122 dr12 19
r r r where r12 = r2 − r1 as
r ∫ d r1 = V
with
r d r12 = 4π r122 dr12 G. Jackson – Imperial College London
The final expression for the configurational energy per particle in terms of the radial pair distribution function is obtained from equation 24 as U u = = 2πρ ∫ u2 (r12 )g (r12 ) r122 dr12 (25) N This so-called energy equation can be interpreted easily as follows: volume of shell
u=∫
1 2
u2 (r12 ) ρ g (r12 ) 4π r122 dr12 mean number of particles in shell
n( r12 )
mean energy of central particle with particles in shell sum over shells 20
G. Jackson – Imperial College London
2.3 Pressure
An expression for the pressure of the system in terms of the radial pair distribution function can be obtained from the thermodynamic definition of the pressure (volume derivative of Helmholtz free energy): ⎛ ∂F ⎞ (26) P = −⎜ ⎟ ⎝ ∂V ⎠ NT The free energy is related to the canonical partition function through F = −kT ln QNVT (cf. equation 3) so ⎛ ∂ ln QNVT ⎞ ⎛ ∂ ln Z NVT ⎞ P = kT ⎜ ⎟ ⎟ = kT ⎜ (27) ⎝ ∂V ⎠ NT ⎝ ∂V ⎠ NT The second equality of equation 27 stems from the fact that the only the configurational contribution to the partition function depends on the volume r r Z NVT = ∫ exp(− β U (r N ))d r N
21
G. Jackson – Imperial College London
We now employ the method of Green (1947) and reexpress the configurational integral1 in terms of coordinates scaled by the box dimension L = V 3 rN rN Z NVT = ∫ exp(− β U (r ))d r rN r r = ∫ ...∫ exp(− β U (r ))d r1...d rN rN = ∫ ...∫ exp(− β U (r ))dx1dy1dz1...dx N dy N dz N r = V N ∫ ...∫ exp(− β U (r N ))dx1*dy1*dz1*...dx*N dy *N dz *N r r = V N ∫ exp(− β U (r N ))d r *N (28) where 1 * xi = xi / L = xi / V 3 1 3
y = yi / L = yi / V * i
dx = dxi / V dy = dyi / V r* r * * * d ri = dxi dyi dzi = d ri / V * i
22
* i
1 3
1 3
z = zi / L = zi / V
dz = dzi / V * i
* i
1 3
1 3
(29) G. Jackson – Imperial College London
In the case of pairwise-additive potentials (cf. equation 19) N N rN U (r ) = ∑∑ u2 ( rij ) i =1 j >i
the configuration integral can be differentiated with respect to volume as 1 rN r *N ∂ Z ⎛ NVT ⎞ N −1 ⎜ ⎟ = NV ∫ exp(− β U (r ))d r ⎝ ∂V ⎠ NT 0 (30) N 1 N N ⎛ ⎞ ∂u2 ( rij ) ∂rij rN ⎜ r *N V ⎟ − dr exp(− β U (r )) ∑∑ ∫ ⎜ ⎟ kT 0 ⎝ i =1 j >i ∂rij ∂V ⎠ where (xi − x j )2 + ( yi − y j )2 + (zi − z j )2 rij =
[
=V 1 3
1 3
]
[(x − x ) + (y − y ) + (z − z ) ] * i
rij = V r
* ij
* 2 j
* i
* 2 j
∂rij
rij*
rij
=
=
* i
* 2 j
(31)
∂V 3V 3V Note invariance of scaled distances to changes in volume. 23
2 3
G. Jackson – Imperial College London
The derivative of the configurational integral can then be written as rN r *N ⎛ ∂Z NVT ⎞ N −1 ( ( ) ) NV β U r d r exp = − ⎜ ⎟ ∫0 V ∂ ⎝ ⎠ NT 1
1 r N ⎛⎜ N N rij ∂u2 ( rij ) ⎞⎟ r *N VN dr exp(− β U (r )) ∑∑ − ∫ ⎜ ⎟ kT 0 ⎝ i =1 j >i 3V ∂rij ⎠
rN rN N = ∫ exp(− β U (r ))d r V r N ⎛⎜ N N ∂u2 ( rij ) ⎞⎟ r N 1 dr exp(− β U (r )) ∑∑ rij − ∫ ⎜ ⎟ 3VkT ∂rij ⎠ ⎝ i =1 j >i = ρ Z NVT
r N ⎛⎜ N N ∂u2 ( rij ) ⎞⎟ r N 1 exp(− β U (r )) ∑∑ rij dr − ∫ ⎜ ⎟ (32) 3VkT ∂rij ⎠ ⎝ i =1 j >i
where we have transformed back to the original variables (cf. equation 29). 24
G. Jackson – Imperial College London
We proceed to obtain an expression for the pressure from equation 32: 1 ⎛ ∂ ln Z NVT ⎞ ⎛ ∂ ln Z NVT ⎞ = ⎜ ⎟ ⎜ ⎟ ⎝ ∂V ⎠ NT Z NVT ⎝ ∂V ⎠ NT r N ⎛⎜ N N ∂u2 ( rij ) ⎞⎟ r N 1 1 =ρ− exp(− β U (r )) ∑∑ rij dr ∫ ⎜ ⎟ (33) ∂rij ⎠ 3VkT Z NVT ⎝ i =1 j >i
As for the internal energy the double sum over i and j gives rise to N ( N − 1) / 2 pair terms which yield the same results on integration so we can choose a given pair (say 12): ⎛ ∂ ln Z NVT ⎞ ⎜ ⎟ =ρ ⎝ ∂V ⎠ NT
rN r N −2 ⎤ ⎡ ( ( ) ) exp β − U r d r r r (34) 1 ∂u2 ( r12 ) ∫ ⎢ N ( N − 1) ⎥ d r1 d r2 − r12 ∫∫ 6VkT Z NVT ∂r12 ⎢ ⎥⎦ ⎣
r r ρ g (r1 , r2 ) Note presence of pair distribution function (cf. equation 12) 2
25
(2) N
G. Jackson – Imperial College London
In terms of the pair distribution function: 1 ∂u2 ( r12 ) ( 2 ) r r r r ⎛ ∂ ln Z NVT ⎞ 2 g N (r1 , r2 ) d r1 d r2 ρ ∫∫ r12 ⎟ =ρ− ⎜ 6VkT ∂r12 ⎝ ∂V ⎠ NT r ∂u2 ( r12 ) ( 2 ) r 1 2 =ρ− g N (r12 ) d r12 ρ ∫ r12 ∂r12 6kT ∂u2 ( r12 ) 2π 2 =ρ− g (r12 ) r122 dr12 ρ ∫ r12 ∂r12 3kT The pressure is then given by
(35)
∂u2 ( r12 ) 2π 2 ⎛ ∂ ln Z NVT ⎞ ρ ∫ r12 P = kT ⎜ g (r12 ) r122 dr12 (36) ⎟ = ρkT − ∂r12 3 ⎝ ∂V ⎠ NT This is usually referred to as the pressure equation or the virial equation as it corresponds to an average of the virial function r r ∂u ( r ) r1 ⋅ F12 = r12 2 12 ∂r12 26
G. Jackson – Imperial College London
We have now demonstrated the usefulness of the pair distribution function g (r12 ) in expressing the thermodynamic properties such as the energy and pressure of a many-particle system in terms of onedimensional integrals rather than the N-dimensional integrals inherent in the full ensemble averages (cf. equation 4). Though the expressions appear simpler in form, the difficulty is now shifted to determining the distribution function. The other thermodynamic properties can be obtained from the standard thermodynamic relations. For example, the Helmholtz free energy can be obtained by integrating equation 26.
27
G. Jackson – Imperial College London
3. Perturbation Theories As we have seen in the previous sections the structural and thermodynamic properties of fluids of particles interacting through pair-wise forces can be determined from a knowledge of two-body distribution functions. Unfortunately, an analytical form of these functions is only available for a very limited number of systems (e.g., particles interacting through a hardsphere pair potential). The problem is even more complicated if there are many-body forces or if the intermolecular potential is not spherically symmetrical. Perturbation theories provide a means by which the properties of a system (e.g., the distribution functions, free energy or pressure) can be represented as a perturbation from those of a reference system with known properties. Such approaches therefore represent some of the most versatile, accurate and powerful theories of the liquidG.state. 28 Jackson – Imperial College London
The intermolecular pair potential is often found to separate naturally into a steep short-range repulsion and a smoothly varying long-range attraction. The structure of most simple liquids is largely determined by the molecular packing which is dominated by the repulsive interactions. Structure factor of a Lennard-Jones fluid (curve) compared with that an equivalent hard-sphere system (points). See Verlet (1968); Hansen and McDonald (2006).
The attractive interactions may thus be treated as a uniform background potential that contributes to the configurational energy of the fluid but does not affect its structure. 29
G. Jackson – Imperial College London
3.1 van der Waals (1873)
Little was know of the form of intermolecular potentials at the time that van der Waals developed his equation of state for fluids. As Logan showed in his lectures the assumptions made by van der Waals in developing his theory of liquids involve the use of a mean-field approximation (uniform background attractive potential). Here, the total configurational energy is approximated as a sum of one-body potentials which only depends on the position of a given molecule: N r (37) U ≈ ∑ϕ i (ri ) i =1
The configurational integral can thus be expressed as rN rN Z NVT = ∫ exp(− β U (r ))d r
[
30
r r ≈ ∫ exp(− β ϕ (r )) d r
]
N
≈ [(V − Vexc ) exp(− β u )]
N
(38)
G. Jackson – Imperial College London
In equation 38 Vexc is the average volume exclude to a given particle by the other particles and u is the average potential energy experienced by a particle which is assumed to be constant (at a given density). An expression for the pressure can be obtained directly from the configurational integral (cf. equation 27) NkT ∂u ⎛ ∂ ln Z NVT ⎞ P = kT ⎜ −N ⎟ = ∂V ⎝ ∂V ⎠ NT V − Vexc
(39)
It is clear from the expression for the configurational energy (cf. equation 25) that if one assumes g (r12 ) = 1, i.e., no correlations between particles, u = −ρ a 31
with the constant a = −2π ∫ u2 (r12 ) r122 dr12
(40)
G. Jackson – Imperial College London
The volume derivative of the one-particle potential is thus given by ∂u ρ N = 2a= a ∂V V V so that from equation 39 the pressure is simply NkT P= − ρ 2a V − Vexc
(41)
(42)
All that remains is to provide an expression for the excluded volume. If one assumes that the particles are hard-spheres with a diameter σ then the excluded volume 3 4 v = πσ (shaded region below) for a pair of particles is exc 3 . van der Waals approximated the total excluded volume by N V ≈ vexc = Nb (43) exc σ 2 2 where the constant b = πσ 3 3 G. Jackson – Imperial College London 32
Inserting the approximate excluded volume (equation 43) into equation 42 we obtain the familiar form of the van der Waals equation of state in terms of his two constants a and b: NkT P= − ρ 2a (44) V − Nb This can be considered as a perturbative approach in that a uniform attractive background (perturbation) is added to the hard-sphere (reference) term. We will come back to this point in the following section. The hard-sphere contribution to the pressure (or compressibility factor Z) can be expressed in terms of the packing fraction η (the fraction of the volume occupied by the molecules) as PV πσ 3 1 1 Z= = = where η = ρ (45) NkT 1 − ρ b 1 − 4η 6 33
G. Jackson – Imperial College London
3.2 Zwanzig (1954) – The High Temperature Expansion
The basis of modern perturbation theories of liquids can be traced back to the work of Longuet-Higgins, Barker, Pople, and Zwanzig in the early 1950s. Here we follow the treatment of Zwanzig (1954). Let us assume that we can separate the pair potential into reference (0) and perturbative (1) contributions, i.e., u2 ( rij ) = u2( 0 ) ( rij ) + u2(1) ( rij )
(46)
The configurational energy can then can then also be expressed in terms of the two contributions as N N N N rN rN rN (0) (1) U (r ) = ∑∑ u2 ( rij ) + ∑∑ u2 ( rij ) = U 0 (r ) + U1 (r ) (47) i =1 j >i
34
i =1 j >i
G. Jackson – Imperial College London
The Helmholtz free energy of the system of particles interacting via the full potential and the reference potential can be expressed in terms of corresponding partition functions (cf. equations 2 and 3): F = − kT ln Q
F0 = −kT ln Q0
(48)
The difference in free energy between the full and reference system (i.e., the perturbation to the free energy) is Q Z ΔF = F − F0 = − kT ln = − kT ln (49) Q0 Z0 The free energy difference can be expressed in terms of the ratios of the configurational integrals r r r r Z = ∫ exp(− β U (r N ))d r N and Z 0 = ∫ exp(− β U 0 (r N ))d r N (50) as the kinetic (ideal) contributions cancel. 35
G. Jackson – Imperial College London
The ratio of the configurational integrals is rN rN exp(− β U (r ))d r Z ∫ = rN rN Z 0 ∫ exp(− β U 0 (r ))d r rN rN rN exp(− β U 0 (r ))exp(− β U1 (r ))d r ∫ = rN rN ( ( ) ) β − exp U r d r 0 ∫ rN rN rN (N ) = ∫ P0 (U 0 (r ))exp(− β U1 (r ))d r r = exp(− β U1 (r N )) 0
(51)
This corresponds to the average of the Boltzmann factor of the perturbation contribution to the energy with respect to configurations of the reference system; the probability distribution of the reference P0( N ) appears in the integral (cf. equation 6). 36
G. Jackson – Imperial College London
The perturbation to the free energy is obtained from equation 49 as rN ΔF = − kT ln exp(− β U1 (r )) 0
(52)
This important result due to Zwanzig is the starting point for the popular perturbation theories of Barker and Henderson (1967), and of Weeks, Chandler and Andersen (1971). To make the expression more tractable both the exponential and the logarithm can be expanded in powers of the ∞ xi x perturbative energy. Starting with the exponential e = ∑ i =0 i! 1 1 ⎞ ⎛ 2 3 ΔF = − ln ⎜1 − β U1 + [β U1 ] − [β U1 ] + ... ⎟ β 2! 3! ⎠ ⎝ 1
1 2 2 ⎛ = − ln⎜1 − β U1 0 + β U1 β ⎝ 2 1
37
1 3 3 − β U1 0 6
0
⎞ + ... ⎟ 0 ⎠
(53)
G. Jackson – Imperial College London
1 2 1 3 We now expand the logarithm ln (1 + x ) = x − x + x − ... 2 3 1⎛ 1 3 3 ⎞ 2 2 ΔF = − ⎜ − β U1 0 + β U1 − β U1 + ... ⎟ 0 0 β⎝ 6 ⎠ 1 ⎛ 2 2 + ⎜ − β U1 0 + β U1 2β ⎝ −
1 3β
⎛ 2 2 ⎜ − β U1 0 + β U1 ⎝
1 3 3 − β U1 0 6 1 − β 3 U13 0 6
⎞ + ... ⎟ 0 ⎠
2
(54)
3
⎞ + ... ⎟ + ... 0 ⎠
Collecting the terms in powers of the inverse temperature we obtain the so-called high-temperature expansion:
(
)
1 2 2 ΔF = U 1 0 − β U 1 − U 1 0 0 2 1 2 3 + β U1 − 3 U12 U1 0 + 2 U1 0 0 6
(
38
3 0
)+ ...
(55)
G. Jackson – Imperial College London
When one chooses a hard-core reference system (0) and an attractive perturbation (1), as is common practice, the first-order term is the so-called mean-attractive energy. An explicit expression for the mean-attractive energy can be obtained in terms of the pair distribution function of the repulsive reference system following a similar procedure to that employed for the configurational energy in Section 2.2 U1
0
rN rN rN U1 (r )exp(− β U 0 (r ))d r ∫ = rN rN ( ( ) ) exp β U r d r − 0 ∫ ⎛ N N (1) ⎞ rN rN ⎜ ⎟ ( ( ) ) u ( r ) exp − β U r d r ij ⎟ 0 2 ∫ ⎜⎝ ∑∑ i =1 j >i ⎠ = Z0
39
(56)
G. Jackson – Imperial College London
The double sum over i and j again gives rise to N ( N − 1) / 2 pair terms which yield the same results on integration:
U1
0
1 = N ( N − 1) 2
rN rN ∫ u (r12 )exp(− β U 0 (r ))d r (1) 2
Z0
rN r N −2 ⎤ ⎡ exp(− β U 0 (r ))d r r r 1 ∫ (1) ⎥ d r1 d r2 (57) = ∫∫ u2 (r12 )⎢ N ( N − 1) 2 Z0 ⎢⎣ ⎥⎦ r r ρ g (r1 , r2 ) 2
where
(2) 0
rN r N −2 N ( N − 1) ∫ exp(− β U 0 (r ))d r ( 2) r r g 0 (r1 , r2 ) = Z0 ρ2
(58)
is the pair distribution function of the repulsive reference. 40
G. Jackson – Imperial College London
Following a procedure equivalent to that of Section 2.2: U1
0
r r 1 2 (1) (2) r r = ρ ∫∫ u2 (r12 )g0 (r1 , r2 ) d r1 d r2 2 = 2πNρ ∫ u2(1) (r12 )g 0 (r12 ) r122 dr12
(59)
from which it is clear that the first-order term of the perturbation expansion is the average of the attractive part of the potential over configurations of the repulsive reference system. The van der Waals attractive contribution (equation 40) is recovered when g 0 (r12 ) = 1 confirming the perturbative roots of his approach. The second-order term 1 2 2 − β U1 − U1 0 0 2 represents a fluctuation of the attractive energy and is more difficult to evaluate; it depends on high-order distribution functions.
(
41
)
(60)
G. Jackson – Imperial College London
It follows from the ongoing discussion that a knowledge of the thermodynamic properties and structure of the reference system allows one to obtain the properties of the full system (to first-order at the very least). In this regard, the hard-sphere fluid plays a pivotal role as a tractable reference system for which accurate analytical description of the properties is available.
42
G. Jackson – Imperial College London
3.3 Carnahan and Starling (1969) – Hard Sphere Fluid
As is apparent from the previous section the properties of the hard-sphere fluid are central to the development of accurate perturbation approaches for fluids of particles interacting via more realistic intermolecular interactions. The equation of state proposed by Carnahan and Starling (1969) is a key development. The (second-virial) approximation to the hard-sphere repulsive contribution employed by van der Waals (cf. equation 43) fails rapidly at moderate densities because the pair excluded volumes are treated as additive which is clearly only true at very low density. This means that the van der Waals equation of state can not be used to represent the pressure (and other thermodynamic properties) of the hard-sphere fluid at moderate to high densities. 43
G. Jackson – Imperial College London
The coefficients of the virial expansion (as opposed to the virial or pressure equation) of a hard-sphere system can be obtained analytically and numerically to reasonably high order. The compressibility factor of the hard-sphere fluid up to fifth order in density (packing fraction) is accurately represented by
PV = 1 + 4η + 10η 2 + 18.36η 3 + 28.2η 4 + 39.5η 5 + ... Z= NkT
(61)
When the van der Waals expression (equation 45) is expanded as a binomial series only the leading-order term (second virial coefficient) is reproduced exactly.
44
G. Jackson – Imperial College London
If the virial coefficients are approximated by their nearest integer values
Z ≈ 1 + 4η + 10η 2 + 18η 3 + 28η 4 + 40η 5 + ...
(62)
The series can be reproduced as a Padé approximant of the following form:
1 +η +η2 −η3 Z≈ 3 (1 − η )
(63)
This relation (originally proposed by Carnahan and Starling) together with an accurate integral equation solution of the radial pair distribution function (see Appendix D, Hansen and McDonald, 2006) provide a very accurate representation of the properties of the hardsphere fluid for use as a reference in perturbation theories of molecules with more realistic interactions. 45
G. Jackson – Imperial College London
3.4 Wertheim/SAFT – Associating and Non-spherical Molecules
Most molecules of industrial relevance are non-spherical and interact through complex interactions (polar, hydrogen bonding, association etc.). These types of interactions can not be captured adequately by simple intermolecular potentials such as the square-well or Lennard-Jones. Modern theories of the fluid state have been developed within the perturbative formalism of the previous sections to represent the thermodynamic properties and fluid phase equilibria of more complex molecules. The focus in this section is on the theory of associating molecules of Wertheim (1984-86), as implemented in the SAFT equation of state (Jackson et al. 1988, Chapman et al. 1988-1990). 46
G. Jackson – Imperial College London
Model of molecules • Monomeric segments interacting through simple pairwise interactions, e.g., square-well, Lennard-Jones (short-range repulsion and long-range attraction)
•
Nonspherical molecules made up as chains of the fused segments
•
Directional interactions and association (hydrogen bonding, chem. equil.) mediated through additional bonding sites
e.g., butan-1-ol
47
G. Jackson – Imperial College London
Wertheim’s theory of association Wertheim, J. Stat. Phys. 35, 19 (1984); 35, 35 (1984) Wertheim, J. Stat. Phys. 42, 459 (1986); 42, 477 (1986)
Jackson, Chapman, Gubbins, Mol. Phys. 65, 1 (1988) Chapman, Jackson, Gubbins, Mol. Phys. 65, 1057 (1988)
+ + 48
G. Jackson – Imperial College London
Wertheim’s theory of association – Alternative view Galindo, Burton, Jackson, Visco, Kofke, Mol. Phys. 100, 2241 (2002) Hansen, McDonald Theory of Simple Liquids, 3rd Edition (2006) 11.10 Associating Liquids
+ ρ ρ0 ρD
-
total density of segments density of segments as monomers density of segments in dimers
The corresponding mass action equation (cf. chemical equilibria) is
ρ = ρ0 + 2ρ D
49
(64) G. Jackson – Imperial College London
In the case of an ideal dimerising gas (Hill 1956) we can write 1 ρ D (1) = ρ 0 (1) ∫ ρ 0 ( 2 ) [exp (− u (12 ) / kT ) − 1]d ( 2 ) (65) 2 For an interacting system Wertheim (1984) showed that an the expression becomes 1 ρ D (1) = ρ 0 (1) ∫ ρ 0 ( 2 ) g 00 (12 ) [exp (− u (12 ) / kT ) − 1]d ( 2 ) 2 1 = ρ 0 (1) ∫ ρ 0 ( 2 ) g 00 (12 ) f (12 ) d ( 2 ) (66) 2 using an exact cluster expansion. A thermodynamic perturbation theory (TPT1) can be developed with the approximation (67) g 00 (12 ) ≈ g ref (12 )
ρ D (1) ≈ ρ 0 (1) ∫ ρ 0 ( 2 ) g ref (12 ) f (12 ) d ( 2 ) 50
(68)
G. Jackson – Imperial College London
The mass action equation can then be approximated as
ρ (1) = ρ 0 (1) + 2 ρ D (1) ≈ ρ 0 (1) + ρ 0 (1) ∫ ρ 0 ( 2 ) g ref (12 ) f (12 ) d ( 2 )
(69)
A form of the Helmholtz free energy at the TPT1 level that is consistent with this mass action equation is ⎛ δ F [ρ 0 (1) ] ⎞ ⎟⎟ = 0 ⎯ ⎜⎜ ⎯→ mass action equation ⎝ δρ 0 (1) ⎠ ⎛ ⎞ ρ 0 (1) ∫ ⎜⎜⎝ ρ (1) ln ρ (1) − ρ 0 (1) + ρ (1) ⎟⎟⎠ d (1) (70) 1 − ∫ ρ 0 (1) ∫ ρ 0 ( 2 ) g ref (12 ) f (12 ) d ( 2 ) d (1) 2
F assoc . = kT
which follows from the variation of the free energy with the density of monomers. 51
G. Jackson – Imperial College London
In the case of a homogeneous system ρ (1) = ρ ρ 0 (1) = ρ 0 ( 2 ) = ρ 0 X = ρ0 / ρ and the mass action equation (Jackson et al. 1988) can be expressed as a simple quadratic equation: (71) ρ = ρ 0 + ρ 02 ∫ g ref (12 ) f (12 ) d ( 2 ) ≈ ρ 0 + ρ 02 Δ 12
where
Δ 12 = g ref (σ 12 ) K 12 F12
1 X = 1 = X + X ρ Δ 12 or (72) 1 + X ρ Δ 12 The perturbation to the free energy due to association is then ⎛ ⎞ F assoc . ρ 0 (1) = ∫ ⎜⎜ ρ (1) ln − ρ 0 (1) + ρ (1) ⎟⎟ d (1) kT ρ (1) ⎝ ⎠ 2
+
∫ρ
0
(1) ∫ ρ 0 ( 2 ) g ref (12 ) f (12 ) d ( 2 ) d (1)
X 1 ρ 02 1 ρ0 ρ0 (73) X ln = ln − +1− ρ Δ = − + 12 2 2 ρ 2 2 ρ ρ 52
G. Jackson – Imperial College London
For systems with multiple bonding sites one can assume that bonding at a given site is independent of that at another and approximate the free energy with independent terms as
Xa 1⎞ ⎛ ∑A ⎜⎝ ln X a − 2 + 2 ⎟⎠ 1 = 1 + ∑ ρ X b Δ ab
F assoc . = NkT
(74)
XA
(75)
a
One can obtain an expression for the bonding together of monomeric segments into chains by considering the limit of complete association (Chapman et al. 1988). In the following, however, we develop the chain contribution using an alternative approach to provide further insights into the physical nature of the approximations that are employed. 53
G. Jackson – Imperial College London
The contribution to the free energy due to chain formation can be obtained as outlined by Paricaud et al. (2003).
m
We start by considering the formation of a single chain from a monomer fluid N-m segments 1 chain N segments
54
G. Jackson – Imperial College London
The residual (in excess of the ideal contribution) free energy and chemical potential for these two systems are
F
res N
Z = − kT ln NN V
μ res = F Nres− m ,1 − F Nres
Z N − m ,1 = − kT ln N − m + 1 F V ⎛ m − 1 Z N − m ,1 ⎞ ⎟⎟ = − kT ln ⎜⎜ V ZN ⎠ ⎝ res N − m ,1
(76) (77)
where the configuration integrals are defined in the usual way as rN r r ZN = ∫ exp (− U ( r ) / kT )d r1 ... d rN (78) r N − m ,1 r r r Z N − m ,1 = ∫ exp (− U ( r ) / kT )d r1 d rm + 1 ... d rN r r r = V ∫ exp (− U ( r N − m ) / kT )d rm +1 ... d rN (79) 55
G. Jackson – Imperial College London
The monomer m-body distribution function is of the form (cf. equation 11) Vm g (1, 2 ... m ) = ZN
rN r r ∫ exp (− U ( r ) / kT )d rm +1 ... d rN
(80)
An appropriate m-body cavity function can be defined as r y (1, 2 ... m ) = g (1, 2 ... m ) exp (U ( r m ) / kT Vm = ZN
r N −m ∫ exp (− U ( r
V m Z N − m ,1 = =V ZN V
m −1
) r ) / kT )d r
r m + 1 ... d rN
Z N − m ,1 ZN
(81)
The chemical potential of forming a single chain can therefore be expressed in terms of the cavity function as
μ 56
res
⎛ = − kT ln ⎜⎜ V ⎝
m −1
Z N − m ,1 ⎞ ⎟⎟ = − kT ln y (1, 2 ... m ) ZN ⎠
(82)
G. Jackson – Imperial College London
The residual chemical potential of a hard-sphere chain at rolling contact is
μ res = − kT ln y (1, 2 ... m ) = − kT ln g (1, 2 ... m )
(83)
The free energy of forming Nm such chains can therefore approximated as F chain ≈ N m μ res = − N m kT ln g (1, 2 ... m )
(84)
By employing the linear approximation for m-body function g (1, 2 ... m ) ≈ g (σ 12 ) g (σ 23 )... g (σ m −1, m ) g (σ 12 ) ≈ g (σ 23 )... ≈ g (σ m −1, m ) ≈ g ref (σ )
(85) (86)
an expression for chain free energy (TPT1) can be obtained: F chain = − ln g (1, 2 ... m ) ≈ − ln g ref (σ ) m −1 = − ( m − 1) ln g ref (σ ) N m kT (87) 57
G. Jackson – Imperial College London
Statistical Associating Fluid Theory (SAFT)
The SAFT equation of state can be developed within a standard perturbation approach with the association and chain contributions that have just been developed. Original Theory Chapman, Gubbins, Jackson, Radosz, Ind. Eng. Chem. Res., 29, 1709 (1990) SAFT-VR (potentials with variable range) Gil-Villegas, Galindo, Whitehead, Mills, Jackson, Burgess, J. Chem. Phys., 106, 4168 (1997) Numerous other versions (soft-SAFT, PC-SAFT) Blas and Vega (1997), Gross and Sadowski (2001) 58
G. Jackson – Imperial College London
The SAFT free energy is generally expressed as F F ideal F mono . F chain F assoc . = + + + NkT NkT NkT NkT NkT
(88)
which incorporates the ideal, monomer segment, chain and association contributions. The expressions for the SAFT-VR version of the theory are summarised below. The ideal (kinetic) contribution to the free energy is a linear function of the logarithm of the number density F ideal = ln( ρ v / Ω ) − 1 NkT
(89)
where v / Ω is the corresponding de Broglie volume (constant). 59
G. Jackson – Imperial College London
The monomer segment term is obtained from a standard second-order perturbation theory (cf. Section 3.2): Perturbation theory Barker and Henderson (1975) 2 ⎛ 1 ⎞ ⎛ 1 ⎞ mono . hs F = F +⎜ ⎟ F2 ⎟ F1 + ⎜ ⎝ kT ⎠ ⎝ kT ⎠ Hard sphere reference Carnahan and Starling (1969) Mean value theorem a = 1
F hs 4η − 3η 2 = N s kT (1 − η )2
F1 = − 2π N s kT
∞
∫σ
60
(91)
u ( r ) g hs ( r ) r 2 dr
≈ − ρ s a vdw g hs (σ ;η eff . )
Local compressibility approximation (LCA)
(90)
F2 1 hs ∂ a 1* a2 = ≈ K ρs N s kT 2 ∂ρ s
(92) (93)
G. Jackson – Imperial College London
The contribution to the free energy due to association is expressed as described earlier:
F chain = − ( m − 1) ln g mono . (σ ) NkT
(94)
Finally, the chain contribution is given by
Xa 1⎞ ⎛ − ⎟ ⎜ ln X a − ∑ 2 2⎠ a =1 ⎝ 1
(95)
∑
ρ X b Δ ab b =1
(96)
Δ ab ≈ K ab Fab g mono . (σ )
(97)
F assoc . = NkT Xa =
61
1+
s
s
G. Jackson – Imperial College London
The other thermodynamic properties are obtained from standard thermodynamic relations. For example, the pressure and chemical potential can be expressed as the following partial derivatives of the Helmholtz free energy: ⎛ ∂F ⎞ ⎛ ∂F ⎞ μ =⎜ P = −⎜ ⎟ ⎟ (98) ⎝ ∂ N ⎠ VT ⎝ ∂ V ⎠ NT At equilibria the temperature, pressure and chemical potential of the two phases must be the same (solved numerically). Though we have presented the expressions for pure components these can be generalised to mixtures in a straightforward manner (e.g., see Gil-Villegas et al. 1997). The SAFT equation of state provides an accurate description of the thermodynamic properties and fluid phase equilibria of complex fluid mixtures. Some examples of the versatility of the approach are provided in the final section. 62
G. Jackson – Imperial College London
3.5 Examples of the Description of Real Complex Fluids and Mixtures
63
G. Jackson – Imperial College London
Alkanes Vapour-liquid equilibria
Gil-Villegas et al. J. Chem. Phys. (1997)
64
G. Jackson – Imperial College London
Alkanes Vapour-liquid equilibria
McCabe and Jackson PCCP (1999)
65
G. Jackson – Imperial College London
C5H12+PE vapour-liquid equilibria liquid-liquid equilibria
…
PE 100 kg/mol 66
Paricaud et al. Ind. Eng. Chem. Res. (2004) G. Jackson – Imperial College London
Water H2O Vapour-liquid equilibria
Clark et al. Mol. Phys. (2006)
67
G. Jackson – Imperial College London
Water H2O Degree of association
Clark et al. Mol. Phys. (2006)
68
G. Jackson – Imperial College London
Water H2O Vapour-liquid tension
Blas et al. Mol. Phys. (2001); Gloor et al. J. Chem. Phys. (2004) 69
G. Jackson – Imperial College London
HFC-32 Difluoromethane Vapour-liquid equilibria
Galindo et al. J. Phys. Chem. B (1998)
70
G. Jackson – Imperial College London
H2O+HF Vapour-liquid equilibria
Galindo et al. J. Phys. Chem. B (1997)
71
G. Jackson – Imperial College London
H2O+HF Vapour-liquid equilibria
Galindo et al. J. Phys. Chem. B (1997)
72
G. Jackson – Imperial College London
H2O+HF Vapour-liquid equilibria
Galindo et al. J. Phys. Chem. B (1997)
73
G. Jackson – Imperial College London
H2O+HF Vapour-liquid equilibria
Galindo et al. J. Phys. Chem. B (1997)
74
G. Jackson – Imperial College London
H2O+C4E1 Vapour-liquid equilibria Liquid-liquid equilibria
Garcia Lisbona et al. Molec. Phys. (1997)
75
G. Jackson – Imperial College London
H2O+C4E1 Liquid-liquid equilibria
Garcia Lisbona et al. Molec. Phys. (1997)
76
G. Jackson – Imperial College London
H2O+C10E5 Vapour-liquid equilibria Liquid-liquid equilibria
Garcia Lisbona et al. Molec. Phys. (1997)
77
G. Jackson – Imperial College London
H2O+CiEj Vapour-liquid equilibria Liquid-liquid equilibria
Garcia Lisbona et al. JACS (1998)
78
G. Jackson – Imperial College London
H2O+PEG Liquid-liquid equilibria
…
Clark et al. Macromolecules (2008)
PEG 2.18 to 1,020 kg/mol 79
G. Jackson – Imperial College London