Perturbation and Variational Approaches to Equilibrium Thermodynamics of Gases, Liquids, and Phase I ransitions Perturbation and variational approaches to equation of state and other equilibrium thermodynamic properties of simple liquids and phase transitions, vapor-liquid and liquid-solid, are reviewed and discussed G. Ali Mansoori and Frank B. Canfield and variational approaches to equation Poferturbation state and other equilibrium thermodynamic properties of simple liquids and phase transitions are reviewed and discussed. Perturbation expansion of Zwanzig is introduced in a generalized form for the Helmholtz free energy of an original system with respect to the properties of an arbitrary reference system. Perturbation approaches of Smith and Alder, McQuarrie and Katz, Barker and Henderson, and variational approaches of Kozak and Rice, Mansoori and Canfield (M&C), and Bazarov are introduced and their relations with each other and with the perturbation expansion of Zwanzig are explained. The generalization of the variational approach of M&C for binary fluid mixture is presented, and the related variational parameters are characterized. Methods for prediction of melting transition by perturbation and variational approaches are presented, and it is shown that by the variational approaches of M&C, as well as by those of Bazarov, that the liquid-solid phase transition properties can be calculated in a purely predictive sense. The equation of state for a substance, solid, liquid, or gas, enables one to derive configurational thermodynamic properties such as enthalpy from readily observa12
INDUSTRIAL A N D ENGINEERING CHEMISTRY
ble behavior of the substance, its PVT surface. The practical importance of thermodynamic properties in performing design calculations has stimulated continuing interest arnong engincers in equations of state, but as yet, no wholly satisfactory equation has been developed. After the early work of scientific interest on equations of state for gases at low pressure which ultimately led to the ideal gas law, the most significant development, also in the interest of science, was van der Waals equation of state in 1873 for gases and liquids. Thus, in 1873, a reasonably accurate equation of state for both gases and liquids was available before the practical requirements of industry demanded it. The first half of the present century showed a reversal of this happy situation, and motivation began to be provided by practical as well as scientific considerations. By use of tools provided by statistical mechanics, the virial equation of state, long used for pure gases, was given a firm theoretical foundation, and methods for rigorous treatment of gas mixtures via the virial equation were developed. Most equations of state for simultaneous description of gases and liquids in practical use today show the influence of the virial equation on their development, although a variety of terms is added to ac-
count for the lack of convergence of the virial equation for the liquid state. Probably the most widely used equation of state today, indeed the standard of reference often used when proposing new equations of state, is the Benedict-Webb-Rubin equation. Both the form of this equation and the recipes for estimating mixture behavior from known behavior of its components are based upon virial expansion. While the virial expansion has been valuable in treating gases and gaseous mixtures, its failure to converge in the liquid state forces one to look for a more complete description, either empirically or theoreticaliy. The theoretical problem has been approached in several ways, and these are reviewed briefly below. The main emphases of this review are perturbation and variational approaches to evaluate the partition function, hence thermodynamic properties and equation of state, because these approaches show considerable promise for practical application. As a point of reference to the reader, this review is based largely on references available to the authors as of June 1969. Interpretive Approaches (2)
Interpretive approaches to the equation of state start from a n approximate description of the structure of the system. These approaches are called “lattice” theories because the structure customarily is assumed to resemble the regular lattice structure of crystalline solids. Assumptions regarding structure must be guided both by physical reality and by the ability to calculate the partition function resulting from the assumed structure. For crystalline solids, these requirements are harmonious because solids are known to have regular structures which are disturbed only slightly by thermal motions and such regular, more or less static, structures lead to a partition function which can be evaluated. Liquids, on the other hand, offer a serious challenge to the viability of interpretive approaches because their structure, if it can be called that, is continually changing and can be visualized only on an instantaneous basis. T o account for physical reality in liquids, a static average over the instantaneous structure suggests itself and such a n approach may lead to a satisfactory “lattice” theory for liquids. This is true only if one can solve the mathematical problems in the form of very difficult and complicated combinatorial versions (17), which is still a difficult task. Predictive Approaches (58)
Predictive approaches place emphasis, initially, on the process by which the intermolecular forces determine the structure, in the hope that a correct mathematical description of this process will lead to equations whose solutions describe the actual structure. Theories of this class are called “distribution function” theories because the equations involve distribution functions specifying the probability of finding sets of molecules in particular configurations. T h e three theories, customarily called Born-Green
(BG) (12), hypernetted chain (HNC), and PercusYevick (PY) (52), of which any discussion of the theory of fluids must take account, had very different origins, and each requires different initial assumptions to get tractable results. The first, BG, which was proposed in 1936, requires the direct assumption that the triplet distribution function can be expressed as a product of pair distribution function for the three pairs involved, that is, the superposition approximation. Although the HNC and PY theories bypass direct assumptions regarding the triplet distribution function in their initial formulation, attempts to improve these theories ultimately lead to approximation or calculation of the triplet distribution function. T h e PY theory, while capable of producing an analytic and approximately correct equation of state for hard spheres (62, 65), must be improved before it is equally successful for more realistic systems. Rice and coworkers (67) have replaced the superposition approximation in the BG theory with a more adequate approximation, but extensive calculations have demonstrated that none of the predictive theories are capable of predicting the properties of simple liquids and phase transitions at the present time. T h e triplet distribution function, complicated and unknown, stands in the way of successful predictive approaches to the equation of state and other equilibrium properties of dense gases and liquids. Perturbation and Variational Approaches
Variational and perturbation approaches, the subjects of this review, are neither predictive nor interpretive theories. They are mathematical * means of expanding the configurational partition function of an original system around a relatively simple reference system whose properties are known. The reference system can be: a real substance whose properties are known (55, 69) or can be a hypothetical system, such as hard spheres, as long as relationships are available to describe its behavior. This review is concerned only with the second approach where properties of the reference system can be calculated from a statistical mechanic theory. T h a t is, by defining a particular intermolecular potential function for the reference system, one must be able to calculate other properties of the reference system such as Helmholtz free energy, compressibility, and the radial distribution function. The reference system should resemble the original system as close as possible to speed the convergence of the power series of the partition function of the original system, and could be the result of a predictive or an interpretive theory. T h e equation of state for a gas at very high temperatures is determined largely by the repulsive forces acting between molecules, that is, a real gas at very high temperatures behaves much like an assembly of hard spheres. Thus, it is reasonable to expect that the equation of state a t somewhat lower temperatures could be obtained by treating the attractive forces as perturbations about the repulsive forces. If the perturbational approach is valid, the equation of state for gases at moderate temperatures will be of the following form: VOL. 6 2
NO. 8
AUGUST 1970
13
Pv RT
- = ao(u)
N
a&) + adu> -+ + ... RT (RT)
U
p 2
Perturbation Approach of Zwanzig (68)
T h e thermodynamic properties of a system may be obtained from its partition function. The partition function, Q, for a classical system with N particles and total potential energy, U(rl, . . . , rN) = U, of interaction of all the N ( N - 1)/2 pairs of particles in the system, in equilibrium and a t temperature T , is (22)
and
2 =
s s ... v
e-Budrl
,I
=
h/(2
. . . drN
nrnkT)lf2
p
=
uljo (11) >j-1 where u z j is the potential energy function of interaction between particles i a n d j . For simplicity, let us define
l/kT
Xil
=
i
uij
-
u 2 J.o
- (u j=1
Zwanzig (68),who originally applied this idea, demonstrated that the available data on argon and nitrogen for the range of 0” to 150°C, and u p to a density of six or seven hundred amagat units support this equation of state.
Q = Z/N!P
N
C
=
= -
4v4 ( N N-! 4) !
ssssx
(6)
Zwanzig (68) was able to find the following infinite series for the difference between the Helmholtz free energies of the original and the reference systems :
where and in general w1 =j !
2 n,
a1 =
(U1)O
(8)
(- I)”nn-l.(C n, - I ) !
( Zrne = j)
u
and u 1 = - Uo (10) I n the case where only the two-body interactions are considered, one can write for the potential energy functions the following relationships : 14
I N D U S T R I A L A N D E N G I N E E R I N G CHEMISTRY
By a similar procedure as shown above, it is possible to find the relation for w3, a4,. . . with respect to the correlation functions. T h e approach is straightforward, but the results are lengthy. As is well known (52), the second correlation function, g(2)(r1,r2), is quite complicated and its analytic functional form is not available. For systems of particles with simple potential energy functions, such as hard disks and hard spheres, some satisfactory approximations are available for pair correlation functions. T h e terms w2, ( 3 3 , . . . involve also the third, fourth, fifth, . . . correlation functions for which no satisfactory approximations are available. Their direct evaluation from Monte
Carlo or molecular dynamics calculations would be prohibitively difficult. Thus, Zwanzig truncated Equation 7 after the first term
- Fo
F
=
~1
+ O(l/kT)
(18)
Zwanzig applied this perturbation equation to the high temperature behavior of a gas of molecules that interact according to Lennard- Jones potential, and for the reference system he used the hard-sphere potential
u(r) = 4 uo(r)
-
e[(;)"
=
t
for r 0 for r
I');(
5
>
B*/c3 = B / b C =
1
+
+
(1/T*ce){-2.667
1.333[(l/c6) - 111 (25)
C*/c6 = C/b&' =
0.625
+ (1/T*c6){-1.247 + 1.333[(1/c6) - 111
D*/c'
=
(26)
D/bo%' =
0.287
+ ( l / T * ~ ~ ) { 0 . 0 8 9+9 0.637[(1/c6)
where
bo
-
111 (27)
2/3 rNg3
For high densities, Smith and Alder (60) used the configurational partition function due to free-volume theory of Lennard-Jones and Devonshire (33). For a Lennard- Jones potential with a hard-sphere cutoff, the partition function due to this free-volume theory is
d d
where d was the diameter of the molecules of the hardsphere reference system and d/u = c was a number between zero and one. By considering the definition of compressibility, PV/Nk T, with respect to the Helmholtz free energy according to classical thermodynamics
-PV -
with
we get from Equation 18, by considering Equations 16, 19, and 20
m d y > = [(I
PV
2
k T bp (2 ~p
4 e[
(:y
ao*
- (:>"]g~(r)r"dr}
(22)
-
I n deriving Equation 22, N ( N - 1) is approximated as N 2 [which is sufficiently accurate for N O(lOZ3)], go@)(r)is simply shown by go(r), the system of coordinates is brought to the center of molecule number 1 (which replaces drldrz by drldrlz), the integration over drl is performed (which gives V , the total volume of the system), drlz is replaced by 4 nr1Z2dr12, and the dummy indices 1, 2 are discarded; p = N / V . Equation 22 is the equation of state of a Lennard-Jones 12-6 gas, which is a result of the first-order perturbation theory from the hard-sphere reference system, where (PV/Nk T)o is the hard-sphere equation of state with the diameter of hard spheres equal to d = cu. For the reference hard-sphere system, Zwanzig used the solution of BG equation (26) for the radial distribution function and the equation of state. Smith and Alder (60) by joining Equation 22 to the virial expansion
PV -- - 1 Nk T
+ B ( T ) / u + C(T)/u2 + D(T>/u3+ . . .
and
q
+ p2go,z(r) 3- .
.I
p* = PU3
- (c/ao*)I"
=
where ao* is the reduced nearest-neighbor distance in a face-centered cubic lattice. T h e series in powers of 1 / T * was obtained by the expansion of the exponential in Equation 28. T h e first perturbation term could be calculated analytically, which then leads to the following equation of state :
PV __
q-112
Nk T
12
- (1
T"
12 -
T*
+ 24T
-j; (p*4
- 1.2 p*Z) +
- q-"">[Ll(q)P*4
- 2 Ml(q)P*"]
+
+ Q - ~ / ~ ) J ~ ( Q-) 2(1 P * ~+ q-1'2)Mz(~)p*2] (29)
[(3
where
+ (1/128) l t ( q / 2 ) + (2/729) l t ( q / 3 ) Mt(q) = m d d + (1/16) mi(q/2) + (2/27) mt(q/3) Lt(q)
= lS(9.1
and where (23)
and the expansion of the hard-sphere radial distribution function (46)
+
- Y)41 - 1
= ao/u = (dT/p*)1/3
Iz(P) =
gdr; P) = exp[-uo(r)/kT1 X [1 pgo,i(r)
+Y)/O
1
-I- 4.2 q -I- 3 q2 (1
- 9)'
and
m ( q ) = 11/(1 (24)
were able to find the functional forms of the first few virial coefficients with respect to T* = kT/e and c = d/a,
+ (1/3) q3 -
- qY1
-1
Even though the calculations by Zwanzig were not as successful as expected (probably because of the computational errors) (4,5 ) , Smith and Alder showed that it was possible to evaluate theoretically the thermodynamic VOL. 6 2
NO. 8
AUGUST 1970
15
properties for a potential such as the Lennard- Jones potential at all densities. Smith and Alder's calculations lead to an expansion of the thermodynamic properties in powers of the reciprocal temperature. The convergence of the expansion on the equation of state was such that the first two terms were able to approximate P V / N k T for a reduced temperature, T* = k T / e , greater than 2.0 to within 0.03 unit u p to almost solid densities.
to Equation 22 by replacing d with a and 4 e [ ( v / r ) 1 2( ~ / r ) with ~ ] -e(cr/r)6. For the calculation of (dQ)/ b(l/n), they used the method of Rowlinson (56). Finally for the equation of state, McQuarrie and Katz derived the following relation for the Lennard- Jones n - 6 potential
Perturbation Approach of McQuarrie and Katz (43)
I n the perturbation approach of Zwanzig, which is presented above, there is no contribution by the repulsive part of the potential function of the original system in the calculations. If we had assumed a modified Lennard- Jones potential function for the original system such that
( u (.r.) =
m
{u(r) = 4
-
e[(:)"
6 ):(]
for r
d
7
where
and (30) y = -ln pa
we would have obtained the same results for the equation of state as Equation 22 and the results would be expected to be more nearly correct for Equation 30 than for Equation 19. So actually Equation 22 is the equation of state for a modified Lennard-Jones potential. The repulsive part of a Lennard-Jones potential and also real intermolecular potentials, though steep, is not infinitely steep as implied in Equation 30. I n consequence, numerical results are extremely sensitive to the value selected for the diameter, d, of the hard-sphere reference system. T o solve this problem, McQuarrie and Katz (43), considered a system with a pairwiseadditive intermolecular potential as follows = a[(;)"
-
(y]
- 6)]
-
nY'(n-O)e
(32)
From Equations 31 and 32 it is evident that u,(r) = 4 e[(:)" n -12
-
(5)6]
pa
J
1nx
0
. exp[-ppax]dx
For the hard-sphere radial distribution function, they used the numerical table derived by Throop and Bearman (63) from the PY equation, and for the hardsphere compressibility, (PV/:VkT)0, they used the Pad6 approximant of Ree and Hoover (57), which will be presented later in this review. The equation of state (Equation 36) reliably produces temperatures as low as T* = 3. While the method of McQuarrie and Katz (like the previously mentioned perturbation methods) cannot predict the low temperature behavior, it is independent of any adjustable parameter. Hard-Sphere Reference System
(31)
with a = [n/(n
-
(33)
I n the perturbation approaches that are reviewed above, and also in the perturbation and variational approaches that will be introduced later for calculating the properties of fluids (gases and liquids), a hardsphere fluid is always used for the reference system. This is because only for a hard-sphere fluid are satisfactory relations for the compressibility and the radial distribution function available. One such relation for hard-sphere compressibility is the Pad6 approximant of Ree and Hoover (51)
(%)o
and
=
+
+
+ 0.277264 q2) + 1.301008 a')
4 ~ ( l 0.254028 7 (1 - 2.245972 r7
(37) where
By using this characteristic of u,(r) as given by Equation 33, McQuarrie and Katz rigorously expanded the partition function in powers of 1/ n
I n this relationship, Q ( n = ) is the partition function of Equation 34. The resulting equation of state is similar 16
I N D U S T R I A L A N D E N G I N E E R I N G CHEMISTRY
7 =
apd3/6
(38)
Equation 37 is in very good agreement with the machinecalculated data for hard spheres. Solution of the Percus-Yevick integral equation for the hard-sphere system produces two equations of state depending on whether for pressure the Ornstein-Zernike (47, 62, 65), or virial theorem (62) approximations are used, as follows:
=
1+17+v2 (1 - 17)3 (from Ornstein-Zernike equation)
(%&
’ -k
+
=
(1 -
(39)
(from virial theorem)
v2
712
(40) Thiele (62) has shown that the values of compressibilities from virial theorem equation are slightly less, and from Ornstein-Zernike equation are slightly more than the values calculated by Monte Carlo and molecular dynamic methods. Accordingly it is concluded (40) that the arithmetic average of the above two equations should be better than either of the two
(ao 1
+ 17 +
=
92
(1 -
- 3 173~2
(41)
v)3
While Equation 41 is slightly inferior to Equation 37 in matching machine-calculated data, it has the advantage of being simple. Also there is no distinction between Equations 41 and 37 in graphical representations (40)with scale sizes as used for the graphs in this report. Recently, Carnahan and Starling (74) showed that 1
+ 17 + v 2 - v3
(g)O =
(1 -
(42)
d3
is in a better agreement with the machine-calculated data for hard spheres than even the Pad6 approximant of Ree and Hoover, as shown by Equation 37, while Equation 42 has a much simpler form than the Pad6 approximant. For the calculation of the hard-sphere radial distribution function, RDF, there are several approaches available (26, 62, 65). The solution of PY integral equation for the hard-sphere R D F (62, 65) is probably closer to the exact values than any other approach. While the analytic relation for R D F as derived from PY integral equation is not available, Wertheim (65) has calculated the analytical relation for the Laplace transform, G(S), ofgo(x)x; x = r / d , as follows G(S) =
lm
e-sz.go(x)xcix
(40). For x > 2, Throop and Bearman (63) have calculated go(x) numerically from the above relationships and have presented it in a table as a function of x and pd3. By having the equation of state and the R D F of the hard spheres in hand, we can calculate other equilibrium properties of hard spheres such as Helmholtz free energy, pressure, internal entropy, and internal energy. Examination of Adequacy of Perturbation Expansion of Zwanzig
As was explained before, w2, the coefficient of - (1 /k T ) in Equation 7, which is shown by Equation 17, is a complicated function of the second, third, and fourth distribution functions. These distribution functions by themselves are very complicated and unknown functions. Consequently, Zwanzig and the other investigators, forced to truncate Equation 7 after the first term as shown by Equation 18, found the resulting equation of state to be valid at moderately high temperatures but not at conditions characteristic of the liquid state. For this purpose Barker and Henderson ( 4 ) ,realizing that failure at low temperature could result either from inadequacies in the perturbational treatment of attractive forces or in the treatment of the finite steepness of the repulsive forces, decided to use the perturbation theory of Zwanzig for the square-well potential function. I n this case, the effect of the attractive forces is not complicated by the “softness” of the repulsive part of the potential, which is in fact infinitely steep. Furthermore, Monte Carlo (54)and molecular-dynamics calculations provided Barker and Henderson with machinecalculated data. Then the perturbation theor‘y could be compared with these data without the uncertainty due to the lack of knowledge of the true potential function that is inevitable in applications to real fluids. Also, they were able to find two different approximate, but manageable, relations for w z which enabled them to test the convergence of Equation 7 as follows: macroscopic compressibility (m.c.) approximation
(43)
and local compressibility (1.c.) approximation :
(44)
The other advantage of choosing square-well potential was that no adjustable parameter was needed.
where
and S(S) = (1 - 1 7 ) ’ ~ ~
+ 6 ~ ( -1 7)s2 + 18 7’s
(49)
-
+
12 ~ ( 1 2 17)
(46)
where r) is defined by Equation 38. For 1 < x < 2, go(x) is calculated analytically by Wertheim (65),which is shown in a form free of complex functions in reference
u(r) =
[u
-E
(r
< u)
(u < r < Xu) (r > Xu)
By using the 1.c. approximation, Barker and Henderson obtained for the difference of the Helmholtz free energy VOL. 6 2
NO. 8 A U G U S T 1 9 7 0
17
3
0
0.2
0
0
04
VV
Figure 7. Compressibilities of square-well @id, calculated by Barker and Henderson for dzyerent e / k T ’ s as labeled from Equation 50. Thepoints are Monte Carlo values of Rotenbeig (54)
Figure 2. Compressibilities of square-well Jluid, calculated by Barker and Henderson for dtj’erent e / k T ’ s as labeled from Equation 50. Tfie points are molecular djnamic calculations of Alder and ?VainWright ( 4 )
of the original square-well fluid and reference hardsphere fluid :
( 3 ) as compared with the machine calculated data of Alder. T h e analysis in Table I indicates that the perturbation expansion of Zwanzig is adequate, and in the case of more realistic potentials it is necessary only to choose a suitable diameter for the hard-sphere reference system for every temperature and/or for ever)- density, to speed up the convergence of Equation 7 such that the first few terms of this series are sufficient for the representation of the true value of the series. The perturbation and variational approaches, which will be presented below are basically designed for the solution of this problem.
where
77 = 7rpd3/6
d =
a
To arrive a t this relationship, the hard-sphere compressibility Equation 39 was used. For X < 2 (which was the case of their calculation, X = 1.5), analytic expressions for po(r) were available (65). T h e equation of state can be obtained by straightforward numerical differentiation of Equation 50 by the use of thermodynamic Equation 21. Equation 37 was used for the hardsphere pressure required for the zero-order term. The difference between 1.c. and m.c. approximations for the second term, as shown by Equations 47 and 48, was not sufficiently large to affect their graphical comparisons with machine-calculated data. Their results indicated that the useful convergence of the perturbation expansion of Zwanzig extends to very low temperatures for the square-well potential, so that the truncated expansion with first- and second-order terms gives a good approximation to the correct pressures at all temperatures relevant for liquids and compressed gases. Also, their results suggested that their two approximations to the second-order term were both reasonably accurate (3, 7, 67). Figures 1 and 2 show the equation of state o f a squarewell fluid, calculated by Barker and Henderson, by using Equation 50 and the thermodynamic Equation 21. Also Table I shows the critical properties of the square-well fluid calculated by Barker and Henderson 18
I N D U S T R I A L A N D E N G I N E E R I N G CHEMISTRY
Pertubation Approach of Barker and Henderson
For more realistic potentials than the square-well potential, Barker and Henderson (3, 5) presented the following approach. They considered an arbitrary potential function u ( r ) and defined a modified potential function v(d, u,a , y ;r ) by Equation 51 : u(d,
U,
a, 7 ; r ) = u [ d
+ (r - d ) / a ] ,d +
(r - d ) / a = O,u 3.95. Figure 5 shows the densities of coexisting phases, and Table 11 shows the critical constants for Lennard- Jones 12 :6 potential as calculated by Barker and Henderson. Overall, with d calculated from Equation 57, Barker and Henderson were able to extend the applicability of the perturbation equation of state of Zwanzig, Equation 7, to conditions characteristic of liquids and to reduced densities near the solid phase. Where three-body forces are also considered, Helmholtz free energy takes the following form (6, 7) xu(.)
- - Fo - 2 a p P
+
2.70797 1.68918 p’ - 0.31570 p “ 1 - 0.59056 p ‘ 0.20059 p “
+
(62)
with p‘ = pd3. Equation 60 is used to calculate equilibrium thermodynamic properties of argon (7). The values of compressibility, internal energy, and critical constants calculated based on Equation GO are superior to the case where only two-body forces are considered in comparison with experimental argon data. This is true when for the two-body potential function a more realistic relation than the LennardJones function, such as the one derived by Barker and Pompe €or argon ( 6 ) , is used. When the LennardJones function together with the three-body potential (Equation 61) are used, the calculated values for the thermodynamic properties of the liquid argon deteriorate more than when only the Lennard-Jones potential is used. This could be partially due to the LennardJones potential function not being a true two-body potential function, and it could be considered as an “effective potential function” that also automatically compensates for three-body forces [see p. 4719 of Ref. (5)1. Accordingly in the case of Figure 4, where the results of perturbation calculations are compared with the experimental argon data, it is not implied that argon is exactly a substance with only two-body forces, but that E and u used for argon in Figure 4, which are the same values used by Barker e t al. (6) in Equation 60, are “effective” force constants. Kim et al. (24), following Green (79),in expanding the Helmholtz free energy of a nonclassical system with respect to Plank’s constant, h,
F
=
F(classica1)
+ +...
h2 Np 24 r m k T
g(r)V2u(r)r2dr
(63)
im
go(d, r)u(r)r2dr TABLE II. CRITICAL CONSTANTS FOR 12:6 POTENTIAL
(p2P/6)
ss
(123)w(123)dr~dr3 (60)
where g0(~)(123) is the third correlation function a n d ~ ( 1 2 3 )is the three-body potential function. T h e most important part of the three-body potential function available is probably the triple-dipole dispersion interaction 20
d3
100
T*
Nk T
+ 3 cos 81 cos 8 2 cos 8 3 ) / ( r 1 r z r 3 ) ~ (61)
INDUSTRIAL A N D ENGINEERING CHEMISTRY
TC* PO*
pc * P,VJNkT,
Mansoori and Canfield (Var.
Barker and Henderson (Pert.
APP.1 1.36 0,325 0.165 0.37
AM.) 1.35 0.30 0.14 0.34
Machine Calculations (64) 1.32-1.36 0.32-0.36 0.13-0.17 0.30-0.36
Experimental (Argon) (55) 1.26 0.316 0.117 0.233
0.2
0.4
0.6
0.8
1.0
YC
T*
Figure 4. Compressibilities of Lennard- Jones 12 : 6 Juid calculated by perturbation approach of Barker and Henderson (see text) for digerent kT/e’s as labeled. The points given by 0,0, and are machinecalculated values; the points given by 0 were calculated by using a jive-term virial expansion, and the points given by A, m, and V are experimental values for argon ( 4 )
Figure 5. Densities of coexisting phases for Lennard-Jones 12:6 Juid calculated by perturbation approach of Barker and Henderson. and are machine-calculated values and the T h e points given by points given by 0 and are experimental values for liquid and solid argon, respectively (4)
and using the Helmholtz free energy derived by perturbation approach of Barker and Henderson for F(classical), were able to calculate the thermodynamic properties of nonclassical simple substances, helium, hydrogen, and neon, by perturbation method. Of course, g ( r ) is not known for any potential function other than the hard spheres, consequently they also use go(r) instead of g ( r ) , while in fact go(r) is the first term in the perturbation expansion of g ( r ) , with a hard-sphere reference system. For the potential function, LennardJones 1 2 :6 potential function is used, and for the hardsphere compressibility, the Pad6 approximant of Ree and Hoover (51) is applied. This approach gives good results for gaseous helium and hydrogen and satisfactory results for liquid neon (24). The reason for applying this approach for only high temperature helium and hydrogen is that the force constants available for these substances are only reliable at high temperatures ( 2 3 ) .
and its deterioration a t low temperatures, can probably be explained in light of the following variational approaches. Variational Approach of Mansoori and Canfield
(37, 38, 39) According to Taylor’s theorem (73) the exponential in the right-hand side of Equation 5 can be expanded around w1 and then truncated as follows Q/Qo
exp(-Pwd
+ 2P2! exp(-PPEd((U1 -
-
t ~ ) (64) ~ ) ~
where
u1 < E1 < w1 and Q/Qo
=
exp(-PwJ
+ 2-P2! exp(-Pwd((Ul
- wJ2)o
-
P3 e x p ( - h ) ( ( U 1 - ~ 1 ) -~ ) ~ 3!
Variational Approach of Kozak and Rice (27)
Kozak and Rice extended the perturbation equation of state that was derived by Smith and Alder (60) by also considering the higher order virial coefficients. Then they considered the parameter c = d / u as an equilibrium property of the thermodynamic system. I n other words, they required that at equilibrium the Helmholtz free energy should be minimum with respect to c. Even though they obtained satisfactory results a t high temperatures, they indicated that at low temperatures the equation of state deteriorates in quality. The parameter c has basically mathematical significance with little physical meaning. T h e success of the variational method of Kozak and Rice at high temperatures,
=
P4 - exp(-Pb)((U~ 4!
$d4)0
(65)
where u 1
< t 2 < w1
I n both Equations 64 and 65 the last terms on the righthand sides are always positive
P2
- exp(-PEd((U1
2!
3
-
6d2)0
-
E z ) ~ ) 3~ 0
0
and
P4
- exp(-Ptz)((U,
4!
VOL. 6 2 NO. 8 A U G U S T 1970
21
consequently, we get the following inequalities
. 3 & - [ ( 3P
+
2
2 P2
where Inequality 66 was originally derived by Lukes and Jones (36). By using the definition of Helmholtz free energy with respect to partition function, F = -kT In Q, the following inequalities are derived
F - Fo
o T* bc2
cm1)
-
(84)
where cml is calculated by Equations 82 and 83. Figure 8 is the result of the following relationship for Helmholtz free energy
F*(p*, T * )
Fo*(p*, T*; cm1) -Ir l ( p * ; cmd T"
+
r & * ; cmd T#
(85)
I t is also possible to use the following relationship for free energy
F*(p*, T*) S Fo*(p*, T*; cm1) T h e only variable a t our disposal in the right-hand side of the above relations is c = d / u . I t is evident that if the right-hand side of Inequality 80 is minimized with respect to c, this inequality will be closer to equality. Then it is necessary and sufficient to have
+ r l ( p *T*;
+
Numerical calculations (37, 38) indicated that Equations 85 and 86 diverge a t low densities. This is in accordance with the approximations introduced in Equations 85 and 86 that are expected to be valid only at high densities ( 3 ) . Figure 9 shows the densities of coexisting liquid and vapor phases as calculated based on Equation 84. Table I1 shows the critical constants of 12 : 6 fluid as calculated by present variational approach based on Equation 84 for free energy, perturbation approach of Barker and VOL. 6 2 NO. 8 A U G U S T 1 9 7 0
23
u t j is the potential function for molecules i a n d j in the
original system and ulJo is the potential function for molecules i a n d j in the reference system. I t is interesting to note that the form of the composition dependence in Equation 87 is equivalent to the successful van der Waals conformal solution approximation (29, 32). By considering a binary mixture of hard spheres for the reference system, Inequality 87 can be written as
-_Fo
+ 16 a p p
[T1d1 13e11 x12
u11*(x)gllo(x)x2dx
+
Nk T
(88)
Jrn u~~*(x)gl~~(x)x~dx]
0
2
4
6
8
10
KT/c
Figure 9. Densities of coexisting phases for Lennard-Jones 12:6 j u i d calculated by considering Equation 84 for Helmholtz free energy. The machine-calculated and experimental points are the same as in Figure 5
where d l l and d22 are the diameters of the molecules of the components 1 and 2 in the reference binary mixture of hard spheres, respectively, and diz =
Henderson, and which are compared with the machinecalculated data for Lennard-Jones 1 2 : 6 fluid, and experimental data of argon. As mentioned previously, Rozak and Rice, considering c = d/u an equilibrium property of the thermodynamic system, minimized free energy with respect to c, to get to the value of the free energy a t equilibrium. By considering the fact that parameter c has basically mathematical significance, we may evaluate the variational approach of Kozak and Rice in the light of Inequality 68, which is a mathematical inequality. At high temperatures,'terms of the order of T*-', and higher, are negligible compared to the zeroth term, FO w1, but a t low temperatures the effect of terms of the order of T*-l, and higher become significant. Kozak and Rice got satisfactory results for compressibility a t high temperatures because they were, in effect, using Inequality 68. At lower temperatures, where their inequality does not represent Inequality 68, calculations were not satisfactory.
+
dii
+
d22
2
If we consider Lennard-Jones potential function for the two components of the original system, we will have
and
For the hard-sphere radial distribution functions in the binary mixture, g l l , 9 2 2 , g l z , no analytic relationships are available. Lebowitz (30) has derived Laplace transforms of these functions analytically, from the PY equation for binary hard spheres. Now by a similar procedure as in the case of pure components and by using the method of Frisch et al. (It?), we can reformulate Inequality 88 with respect to the Laplace transforms. T h e process of reformulation is very lengthy (37) and involves numerous algebraic manipulations. I n this publication we report only the final working formulas. Let us define
Variational Approach of Mansoori and Canfield for Binary Fluid Mixture
For a binary fluid mixture, by considering only twobody forces and introducing the radial distribution functions of the two components (49) in the reference system, Inequality 68 will be as follows (37) Now Inequality 88 can be written in the following dimeruionless form : x 22
5
+
~22~(r12)922*(r12)
1
x 1 x ~ u 1 2 ~ ( r 1 2 ) g 1 2 ~ ( drldrz rl~)
(87)
where x1 and x2 are the mole fractions of the two components with x z = 1 - x1, and uijl = uti - u d J owhere 24
I N D U S T R I A L A N D E N G I N E E R I N G CHEMISTRY
and
where
(95)
and where
D*(s)
=
H - Ll*(s)e"
where y i = 7rpxi/6. Lebowitz and Rowlinson ( 3 7 ) show that Equations 98 and 99 are in very good agreement (with P from Equation 98 slightly above and P from Equation 99 slightly below) with Monte Carlo computations of Alder (7) for several values of d11/d2z and xl/xz, and a large range of 4. This indicates that the arithmetic average of Equations 98 and 99 will be a better approximation to the equation of state than either one alone.
5%
- L2*(s)eS/&I+
-
=
S*(s)eS(*+l
+ Bis + H Lz*(s) = 12 7iiAzs2 + Bzs + H L*(s) = 12 qzzAis2
(96)
s*( 5 )
Integration of Equation 100 over density gives the following relationship for the Helmholtz free energy of a binary hard-sphere mixture (37)
and
H = 36 77111722p1(1 -
t =
711
+
Pd2
(+ +
1722
A1 =
W 1
B1 = and
'>1
pi-
z
Y =
x11*?
+
Z =
XZ,
XIP?
and W = x m ( l -
F*(P*, Ti*, xi)
~
+
XZ,
)
~
1
I Fo*(P*, Ti*, xi;
+
cii, ~22) Cl1,
c 2 2 ) / ~ 1 * (102)
If we assume the following relationships for interaction parameters ,712 and €12:
then Solution of Percus-Yevick integral equation for a system consisting of a mixture of hard spheres, as for a pure component system, also produces two different equations of state whether for pressure the OrnsteinZernike or virial theorem relations are used (30),that is,
(101)
and Fid is the free energy of a mixture of ideal gases with mole fractions x1 and x2. Inequality 93 can be written in the following form
rl*(P*, xl; Since
(1 - [)-2
where
A2 = A3 =
p1
,712
=
+
,711
,722.
2
'
we will have c12
+ c22)/(oL + 1)
= (arc11
and €12/Ell
=
7-1/2
To bring the Inequality 102 closer to equality, it is necessary and sufficient to choose c11 and c 2 2 such that
b~,,* arl* 1 -+----=0 bcii bcii Ti*
-a + ~ -" -* -mi*1 bc22
3622
VOL. 6 2
Ti*
NO.
a
(104)
- 0 AUGUST 1970
25
with the conditions that
(105) and
Variational calculation on Inequality 102, with the variational parameters cll and c22, is expected to predict the properties of binary mixtures as w7ell as Inequality 74 does for pure components (by putting U I I = u22 and dll = d22, Inequality 93 becomes identical to Inequality 74). Recent machine-calculated data by Singer (59) and McDonald (41) for Lennard-Jones mixtures, with application of Equations 103 for unlike pair interaction parameters, will permit a test of this proposition. Equations 103, known to be in error for real fluids, will not compromise the comparison adequately. Discussion
I n the perturbation and variational approaches developed earlier, the general idea was to consider a rnodified potential, that is, u ( r ) = uo(r)
+
(106)
UI(T)
where for r for r and
