- Email: [email protected]
Equations of state from generalized perturbation theory
127
Fluid Phase Equilibria, 67 (1991) 127-150 Elsevier Science Publishers B.V., Amsterdam
Equations of state from generalized perturbation Part II. The Lennard-Jones
theory.
fluid
Gary M. Sowers ’ and Stanley I. Sandler Department of Chemical Engineering, University of Delaware, Newark, DE 19716 (USA) (Received
August 20, 1990; accepted
in final form June 12, 1991)
ABSTRACT Sowers, G.M. and Sandler, S.I., 1991. Equations of state from generalized perturbation theory. Part II. The Lennard-Jones Fluid . Fluid Phase Equilibria, 67: 127-150. We use perturbation theory to derive extensions of our equations of state used previously for the hard-core Lennard-Jones fluid to the more realistic soft-core Lennard-Jones fluid. We find that the use of a “hybrid” Barker-Henderson diameter (see Aim and Nezbeda, 1987) for the hard sphere reference fluid leads to perturbation free energies for the soft-core LennardJones fluid that are similar to the perturbation free energies for a hard-core Lennard-Jones fluid. For the range of temperatures and densities considered here, one of our equations of state is more accurate than a 33-parameter equation (see Ada&i et al., 1988) even though it has only four adjustable parameters. The other equation of state is a straightforward, theoretically based extension of a square-well equation of state (see Lee et al., 1985) to the Lennard-Jones fluid. Both equations of state can be used to correlate the phase diagrams and supercritical PVT data of noble gases and methane, and perform especially well in the correlation of saturated liquid volumes.
INTRODUCTION
The development of a simple, analytic and accurate equation of state for a realistic fluid from statistical mechanics is a difficult problem, even for monatomic fluids interacting with pair potentials as simple as that of Lennard-Jones. Since the thermodynamic and structural properties of a simple fluid can be computed from perturbation and integral equation theories with an acceptable degree of accuracy, some regard the problem of i Present
address:
037%3812/91/$03.50
Dow Chemical
Company,
Plaquemine,
0 1991 Elsevier Science Publishers
LA 70765, USA. B.V. All rights reserved
128
describing these properties from statistical mechanics as solved (Row~nson and Swinton, 1982). However, in engineering applications accurate and simple equations of state with few or no adjustable parameters are required, and perturbation or integral equation theories are of little direct use because they can be mathematically involved and at present are applicable only to relatively simple model fluids. As pointed out by Nezbeda and Aim (1989), a simple, analytic and completely predictive equation of state for the Lennard-Jones fluid that can accurately describe both PVT properties and the phase diagram has not yet been derived from statistical mechanics. There is no theoretically based equation of state that can predict the thermodynamic properties and phase behavior of the Lennard-Jones monatomic fluid more accurately than the empirical 33-parameter equation of state of Nicolas et al. (1979). Although simulation techniques for computing single phase thermodyn~c properties of simple fluids are quite well developed and an improved method of computing the phase diagram of a model fluid from simulation has recently been devised (Panagiotopoulos, 1987), these methods do not provide an equation of state in analytic form. Therefore, from an engineering point of view, the problem of developing an analytic equation of state for the Lennard-Jones fluid from statistical mechanics is still unsolved. Our approach to equation of state development has involved two statistical m~hanical tools; generalized van der Waals theory (Sandler, 1985) and perturbation theory (Hansen and McDonald, 1986). Generalized van der Waals theory provides an expression for the canonical partition function of monatomic or rigid molecular fluids and their mixtures, and can be extended to fluids of flexible chain molecules using the perturbed-hard-chain partition function of Beret et al. (1975). Generalized perturbation theory, based on the h-expansion, was presented in an earlier paper (Sowers and Sandler, 1991). We used the expression for the confi~rational Hehnholtz free energy from generalized perturbation theory as a starting point for our study of the hard-core Lennard-Jones (HCLJ) fluid r
co
-I5y %cLJ(r)= &[( i
_ (;,“I
;:;lz:”
This fluid differs from the soft-core Lermard-Jones potential by having an infinitely steep repulsive part at r = u. It is a useful potential for evaluating perturbation theories because simpler potentials with hard cores, such as the hard-sphere or square-well potential, can be used as reference systems without ambiguity. In contrast, the reference system in perturbation theories for the soft-core Lennard-Jones fluid is softly repulsive, and is typically
129
approximated by a hard-sphere potential with a state-dependent diameter. In our earlier work we presented two equations of state for the HCLJ fluid that fit our computer simulation compressibilities to within 3%. Both of these equations of state were developed from the generalized perturbation theory expression for the configurational Helmholtz free energy PA CoNF = ,AtoNF + PA, N
N
N
(2)
where flAgoNF/N is the configurational free energy of the reference system in the perturbation theory, p = l/kT and /3A,/N is the generalized perturbation free energy (see the third section). The first equation of state results from using a hard-sphere reference system for pAzoNF/N, while the second is novel because it was developed using a square-well reference system. In this paper, we apply generalized perturbation theory to the more realistic soft-core Lennard-Jones fluid, and show how our HCLJ equations of state can be extended to this fluid in a simple and systematic way. In the next section we present a brief review of generalized perturbation theory, while in the third and fourth sections we discuss the existing perturbation theories for the Lennard-Jones fluid that use hard-sphere and square-well reference systems respectively, and how these theories can be used to develop equations of state. The use of a hard-sphere reference system in a perturbation theory of the Lennard-Jones fluid leads to an equation of state very similar to the equation of state for the modified exponential-six fluid proposed by Dodd and Sandler (1989), and to the perturbed-hard-sphere equations of state studied by Aim and Nezbeda (1983, 1984, 1987) and others. In the third section we discuss the choices for a state-dependent hard-sphere diameter in perturbation theory with a hardsphere reference, and the effect of this choice on the equation of state. We present a new perturbed-hard-sphere equation that fits our computer simulation results for the compressibility of the Lennard-Jones fluid to within 2% error. We show, in the fourth section, how the use of perturbation theory leads to a straightforward extension to the Lennard-Jones fluid of the lattice gas model of the square well fluid derived by Lee et al. (1985). This new equation of state is also in good agreement with computer simulation results. In the following section we show that these equations can also be used to correlate the phase diagrams and supercritical PVT data of the noble gases and methane to a high degree of accuracy. THEORY
The generalized perturbation theory expression for the configurational Helmholtz free energy of the Lennard-Jones fluid is derived by first dividing
130
the intermolecular
pair potential into reference and perturbation
U(T) = z+(r) + Xw(r)
parts (3)
where ~a( r) is the reference potential and w(r) = uu( r) - t+(r) is the perturbation potential. When h = 0, the potential u(r) = Q(T), and when X = 1 we recover the Lennard-Jones potential Us from eqn. (3). At fixed N, V and T the change in the configurational Helmholtz free energy of a fluid whose pair potential is u(r) resulting from a change in X is
CONF PA
1 a ( /3AC0NF)
IA=1 =PACoNFI~_O + J0
dX
(4
ah
with
a ( fiACoNF) = ax
(5)
2nNpSjOOw(r)g(r;N,V,T,h)r2dr 0
Defining the generalized perturbation
free energy as
we can write the total configurational
Helmholtz free energy as
PA CoNF = ,8A;oNF + PA, N
N
(7)
N
Equation (7) may appear to be of little value because it requires knowledge of the radial distribution function g( r; N, V, T, X) throughout the entire range 0 < X < 1 to compute the configurational Helmholtz free energy of the fluid of interest. Such information is very difficult to obtain from simulation or from theory. An equation of state for a model fluid can be derived from perturbation theory by making assumptions that allow the terms in the perturbation expansion to be expressed analytically (see, for example, Aim and Nezbeda, 1983; Ponce and Renon, 1976), and limited computer simulation data can then be used to determine the values of the adjustable parameters in the equation of state. We use generalized perturbation theory here in a similar way to extend our equations of state for the HCLJ fluid to the soft-core Lennard-Jones fluid.
APPLICATION
OF GPT USING
A HARD-SPHERE
REFERENCE
SYSTEM
Perturbation theories for the Lennard-Jones fluid using a hard-sphere reference fluid originated with Zwanzig’s (1954) pioneering work. In this section we consider only perturbation theories for the Lennard-Jones fluid
131
in which the reference system is approximated by a hard-sphere fluid whose diameter depends on temperature or on both density and temperature. Two such theories are the perturbation expansions derived by Barker and Henderson (1967) and Weeks et al. (1971) respectively. In Barker-Henderson theory the Hehnholtz free energy is expanded in a double Taylor series in two parameters, a steepness parameter and a well depth parameter, and truncated after the second-order term. An expression for the temperaturedependent Barker-Henderson (BH) diameter of the hard-sphere reference fluid is found by setting one of the first-order terms in the expansion equal to zero, which leads to
dBH(T)= Lo1 - exp[-h_.dr)]dr where uw ( Y) is the Lennard-Jones potential. In Weeks-Chandler-Anderson (WCA) theory the potential is divided into reference and perturbation parts at the potential minimum. A first-order h-expansion of the Helmholtz free energy of the Lennard-Jones fluid is followed by a blip function expansion of the free energy of the soft repulsive reference system about a hard-sphere reference (Anderson et al., 1971). The reference system is approximated by a hard-sphere fluid whose diameter is determined by setting the first-order term in the blip function expansion equal to zero. The WCA diameter is computed by solving the equation
--Puo(r)]JoODyHs(r;dwC*)oxp[
exp[-PuH,(r;dwc,)l}v2dr=0
(9)
for dCcAP~J3, whereyHS
is the background correlation function of the hard-sphere fluid, uO is the pair potential of the soft repulsive reference system and uns( r; dWCA) is the hard-sphere potential with diameter dWCA. Finally, we can define a temperature-dependent hard-sphere diameter by dividing the potential into reference and perturbation parts at the potential minimum as in the WCA theory, and retaining only the zeroth-order term in a Taylor series expansion of r2y,, about r = d (Verlet and Weis, 1972). The resulting temperature-dependent diameter has been referred to as a “hybrid” Barker-Henderson (HBH) diameter (Aim and Nezbeda, 1987), dHBH(T)
=
1
2”601- exp[ -p{ uU(r) 0
+ E}]dr
00)
The free energy of the soft repulsive reference system in each of these perturbation theories for the Lennard-Jones fluid is then CONF
P-&l N
=
~AC$NF
N
(11)
132
where the right-hand side is the configurational Hehnholtz free energy of a hard-sphere fluid whose diameter is computed from the appropriate one of eqns. (S)-(10). If the total configurational Helmholtz free energy /UCoNF/IV is known, eqn. (7) together with eqn. (11) and one of eqns. (8)-(10) can be used to compute the generalized perturbation free energies &4,/N, as follows. Subtracting the relations PA
CONF
N
PZ-1 = o Tdp’
(12)
“2,-l = aP’dd
(13)
/
and CONF
PA0 N
J
we find, by comparison with eqn. (7) that
PA,=
N
“Z-Z,
J p’dp’ o
(14)
where 2 is the compressibility of the Lennard-Jones fluid and 2, is the compressibility of the hard-sphere reference system. Since none of our simulation points are in the two-phase region, we can use the simulation compressibilities in eqn. (14) to compute PAP/N. There is a large body of simulation data available for the Lennard-Jones fluid, but previous studies (Nicolas et al., 1979; Adachi et al., 1988) have not included state points over a wide enough density and temperature range for our purposes. Consequently, we have done additional simulations of the Lennard-Jones fluid in the temperature range 1.35 < kT/c < 6.0 and 0 < pa3 < 0.9. A brief discussion of the simulation method and the results for compressibility and configurational internal energy are given in the Appendix. Although we use only our simulation data to develop the equations of state presented here, we later compare the resulting equations of state with the simulation data of Ada&i et al. (1988). We have used the Carnahan-Starling equation of state z
=
CS
l+77+772-v3 (1 -
d3
(15)
to calculate the hard-sphere compressibility 2,. We use 17 to denote the packing fraction of the hard-sphere fluid, defines as 7 = ( 7r/6)pd3 where d is the hard-sphere diameter. The generalized perturbation free energies /3A,/N calculated from eqn. (14) using eqns. (8), (9) and (10) for the hard-sphere diameter are plotted at a representative temperature kT/e = 1.75 in Fig. 1. The dependence of PAP/N on density is similar at the other five isotherms studied.
133 0.0
a
1
En
la
-1.0
0
4
t
4
-2.0
4
d
0
4 0
-3.0
BHdiameter
cl
4
0
A WCA diameter
0 9
+ HBH diameter
’
-5.0
I 0.2
00
0.4
0.6
0.8
1.0
PO3
Fig. 1. Generalized perturbation free energies computed using a hard-sphere reference fluid with a state-dependent diameter at a representative temperature kT/c = 1.75.
We see from Fig. 1 that /?A,/N is almost, but not quite, a linear function of density regardless of the choice of hard-sphere diameter. However, if the Barker-Henderson diameter is used in the computation of the generalized perturbation free energy, the deviation from linearity is smaller. As a first approximation, we will assume that the generalized perturbation free energies computed using the BH diameter can be fitted with a linear function of density
PAP d
nx
N
=
_n&p*
06)
where p* = pa3 and a is an adjustable parameter determined by fitting to simulation data. When this model for /3A,/N, which we have denoted by &4r,JN, is combined with the Camahan-Starling equation of state for the hard-sphere reference free energy PA,‘ONF/N, the resulting equation of state is z=
1+9+q2-v3 -a&p* (1 -
Id3
07)
which has been referred to as a perturbed-hard-sphere (Aim and Nezbeda, 1987) or Camahan-Starling-van der Waals equation of state. The adjustable parameter a is fitted to the computer simulation compressibilities. In all of our calculations we use the following approximant, accurate to within 0.02%, for the Barker-Henderson diameter: &-i(T) u
1.013/m = 4.142 x 1O-2 + ,I-
in order to avoid repeated numerical integration in eqn. (8).
(18)
134 TABLE 1 Comparison of eqn. (17) and computer simulation cbmpressibilities Data source
a
AZ%
Appendix 1 Ada&i et al. (1988)
6.693 6.482
4.59 9.77
We have determined the parameter a using two different sets of computer simulation data for the compressibility of the Lennard-Jones fluid. The first data set is our simulation data at six supercritical isotherms tabulated in Appendix 1. The second data set is taken from Adachi et al. (1988). State points that may lie in the metastable or two-phase region were removed from Ada&i’s data set by deleting all points below the critical temperature and pressure. The values of the parameter a in eqn. (17) were fitted to each of these data sets and are shown in Table 1. The error between the equation of state and the simulation compressibilities is defined as (19) where ZSrM is the simulation compressibility, ZEoS is the equation of state compressibility and N is the number of state points. We also tested the Nicolas equation of state for the Lennard-Jones fluid using the parameters taken from Ada&i et al. (1988). Although the functional form of the equation was proposed by Nicolas et al. (1979), we will refer to the equation of state with Adachi’s parameter values as the Ada&i equation to avoid confusion. The Ada&i equation produces errors of 3.9% and 3.8% when compared to the subset of Ada&i’s data used here and our simulation data, respectively. If the 33 parameters in the Ada&i equation of state are fitted to the subset of Ada&i’s data used here, the error between the equation of state and the simulation data decreases to 2.26%. The error between eqn. (17) and the simulation data is considerably larger than the error from the Ada&i equation, although eqn. (17) is much simpler. One can improve the performance of this equation of state, without sacrificing too much simplicity, in the,following way. The assumption made in eqn. (16) is that PAP/N is a linear function of density. We improve on this assumption by calculating the correction term j?ACoRR/N necessary to account for the non-linearity in PAP/N
PA‘ORR = MCoNF _ PArF N
N
where Zcs is the compressibility
_+
= /P
,=a
= -
0
of a hard-sphere
dp,
+
+
p
fluid computed
p*
(20)
from the
135
c x
T’=iOO T’=600
-0.15 0.0 0.1
0.2 0.3
Fig. 2. Correction
0.4 0.5
0.6
0.7
PO3 term j3ACoRR/N
0.6 0.9
from eqn. (20) using the BH diameter
in Z,,.
Car&ran-Starling equation of state. This correction term is plotted against density at six isotherms in Fig. 2. We find that it is a complicated function of density and is difficult to fit with a simple function. However, if we use the HBH diameter in Zc, in eqn. (20), the correction term has the form shown in Fig. 3 (the curves have similar shapes if the WCA diameter is used). This correction term can be fitted with the following function of density, which we have used for the HCLJ fluid (Sowers and Sandler, 1991) :
03
8
0.2
2 01
0.0 0.0
0.2
0.4
0.6
0.6
I.0
Fig. 3. Correction term BACoRR /N from eqn. (20) using the HBH diameter lines are the correlations of eqn. (21).
in Z,,.
The solid
136
where pt = 1.0 is the highest reduced density we have considered, and bi, b2 and c are adjustable parameters. The correlations of /UCoRR/N from this model are the solid lines in Fig. 3. The resulting model for the generalized perturbation theory free energies using a hard-sphere reference fluid with a temperature-dependent diameter is CORR
PAPmcd
-=
PAN
--a&p*+
N
(22)
The equation of state for the Lennard-Jones energy, then becomes ==
fluid, obtained from the free
l+7)+Y2-v3 -u&p* (1 -
d3
-cjb,i;Ii,+b2[3(&)2+
G])
x (2p*2 - P*,P*) exp[ -cp*(p* where 17= ( a/6)pd&(
4-m,(T) = (I
- P*,)]
T). For d,,,(T)
(23)
we have used an approximant
21’63m
(24)
0.11 + QE
which is accurate to within 0.05%, to avoid the numerical integration in eqn. (10). The four adjustable parameters in the equation of state (a, b,, b, and c), and the errors between the equation of state and the simulation data are shown in Table 2. The error between eqn. (23) and the simulation compressibilities is almost within the simulation error. It is lower than the error produced by the Ada&i equation of state for the data sets studied here, even though eqn. (23) has only four adjustable parameters, compared with the 33 adjustable parameters in the Ada&i equation. We should point out, however, that the range of density and temperature over which the Ada&i equation is applicable is somewhat larger than the range considered here. The Ada&i equation covers the temperature range 0.5 < kT/c < 6.0 and the density range
TABLE 2 Comparison of eqn. (23) and computer simulation compressibilities Data source
a
b,
Appendix 1 Ada&i et al. (1988)
8.042 7.781
0.731 0.6491
bz
-2.777x10-’ 1.59ox1o-2
C
AZ%
2.000 1.727
1.52 2.16
137 60
N
5.0
0
T’=lOO
4.0
x v X
T-=300 T’=4.00 T’= 6.00
30
0.4
0.2
0.0
0.6
0.6
1.0
Pa3
Fig. 4. Correlations of computer simulation compressibilities
from eqn. (23).
0 < pa3 < 1.2, whereas the parameters in our equations of state were fitted to the region 1.35 < kT/c < 6.0 and 0 < pa3 < 0.9. The correlations of our simulation compressibilities from eqn. (23) are shown in Fig. 4. APPLICATION
OF GPT USING A SQUARE-WELL
REFERENCE
SYSTEM
Both the BH and WCA theories use the hard-sphere fluid as a reference system in the perturbation expansion. Alternatively, since the square-well potential can be used to qualitatively describe the thermodynamic properties of more realistic fluids, one would expect that the perturbation free energy would be smaller when using the square-well potential as a reference system for the Lennard-Jones fluid. A square-well reference perturbation theory of Lennard-Jones fluids was developed by de Lonngi and de1 Rio (1983). They defined a modified potential function v(r;dsw,Rsw,R,,X1,X2)
=
r < d,
UU (z) -e
d,
i %J(W)
RI < r
(25)
where z
=
d,,
+
r-&V
Xl r-R,, x
w=R,+
(26)
2
4 = dsw R, = R,,
+
h(Rm -
dsw)
+ A,( R, - r)
138
and R, = 2’6. The smoothness parameters X, and A, vary the steepness of the repulsive (r < R,) and attractive (r > R,) regions of u(r), respectively. When X, = X, = 0, u(r) is the square-well potential with parameters d R,, and E, and when X, = X, = 1, v( r ) is equal to the Lennard-Jones pL:ntial uLJ(r) with potential parameters E and u. When the configurational free energy of the Lermard-Jones fluid is expanded in a double Taylor series about the square-well reference fluid, corresponding to the point A, = h, = 0, we obtain
PACoNF = bAzoNF N
N
+
il [a(%?‘]+,, (27)
where /3A:ONF/N is the configurational free energy of the square-well reference system. De Lonngi and de1 Rio provide expressions for both firstand second-order terms in the expansion that depend on the radial distribution functions of the square-well reference fluid and their derivatives. The first-order terms A, and A, and the mixed second-order term A,, were set to zero by choosing an effective square-well diameter dsw( T) and an effective square-well range R,,(T) according to
and R,,(T)
= 2%
+ [l -
ex~(~41-‘~~~~~ - exp[-h_h)]d~
(29)
Equation (27) can be written in the form of eqn. (7), provided we make the definition
(30)
1Al=0
PAP -= N
For the compressibility of the reference system Z,, we use the square-well equation of state derived by Lee et al. (1985) ==
1+v+772-773
(l-d3
-
Z, 7&/6
[exp(GkT) - 1117 + [exp( c/2kT)
- l] 9
where n = (r/6)pd&( T). As in the previous section, we approximate the integrals in eqns. (28) and (29) with the following simple analytic expres-
139
-0.2
?=-0.4
2
-0
6
-0.6
0.0
0 +
T’=135 T’=175
0
T’=zoo
k Q
T’=3.00 T’=400
0.2
04
0.6
1.0
0.6
Fig. 5. Generalized perturbation free energies computed using a square-well reference fluid with a state-dependent diameter. The solid lines are the correlations of eqn. (34).
sions, which are accurate to within 0.05%:
(32) %v(T) a
=&w(T)
&w(T)
=
2116+
a
0.347kT/c 0.26 + kT/c
(33)
As a first approximation we set Z,, the maximum coordination number, equal to the constant 18, corresponding to h,, = 1.5. It is important to note that this approximation removes the dependence on X,, from the reference fluid equation of state. The compressibility of the square-well reference system is then given by eqns. (31) and (32). The generalized perturbation free energies PA,/N computed using the square-well reference system are shown in Fig. 5. They can be fitted to an acceptable degree of accuracy by the following quadratic function of density, used for the HCLJ fluid: *=a(&)exp(-E)p*2+[bl(&)+b2(&)l]P*
(34)
where a, b, and b, are adjustable parameters. The correlations of &4,/N from this model are the solid lines in Fig. 5. The equation of state resulting from this model for /3A,/N together with the LLS equation of state for the
140 TABLE 3 Comparison of eqn. (35) and computer simulation compressibilities Data source
a
b,
b,
AZ%
Appendix 1 Ada&i et al. (1988)
- 9.459 - 8.374
- 1.783 - 0.7127
4.581 2.947
3.03 4.18
reference system is Z=
I+17+17*-V3
(1- ?I3 +2n(&)
_
2, [exp(e/2kT) &/6 + [exp(c/XT)
exp( -g)P**+
[4(G)
1177 - I] 17 +4(&)*]p*
(35)
The adjustable parameters in this equation are determined by fitting to simulation compressibilities and are shown in Table 3. The absolute average percent deviations in compressibility, also given in Table 3, are higher than eqn. (23), but still quite reasonable. The correlations of the simulation compressibilities are shown in Fig. 6. An approximation made in eqn. (35) is that the maximum coordination number 2, is a constant. We also examined whether it was possible to improve the performance of eqn. (35) by incorporating &,-dependence into 2,. In the LLS lattice model, 2, is a discontinuous function of the width of the square-well that is calculated from geometry. Since A,, is a function of temperature in the square-well reference perturbation theory, it is desirable to represent 2, as a continuous function of Xs,. One way of doing this is
6.0
5.0
4.0
0 + 0 x v *
I
l-=1.35 T-=1.75 T'=Z.OO T-=300 T-=400 T-=600
N
0.0
0.2
0.4
0.6
0.8
1.0
PO3 Fig. 6. Correlations of computer simulation compressibilities from eqn. (35).
141 to determine 2, by fitting the LLS second virial coefficient to the exact result for the square-well fluid, from which we obtain
&) - l]/[ev( &i=)- 11 -cm =~(&,(T)3- l)[exp(
(36)
(Note that this expression for 2, is undefined at infinite temperature. However, when used in eqn. (35), 2, is multiplied by 9 a d&(T), and this product goes to zero at infinite temperature.) When this model for Z, is used in eqn. (35) and the adjustable parameters are fitted to the simulation data, the performance of the equation of state is no better than when Z, = 18. Therefore we have used Z,,, = 18 in all of our calculations. The advantage of eqn. (35) over eqn. (23) is that it can be written as a polynomial in density. This reduces the computational effort when the equation of state must be solved for the densities of two coexisting phases during a phase equilibrium calculation. Although the deviation between eqn. (35) and the simulation compressibilities is approximately twice that of eqn. (23), the use of the square-well reference system has allowed us to use a very simple model for /3A,/N, and yet develop an accurate equation of state. Other equations of state for the square-well fluid have been proposed (Aim and Nezbeda, 1983; Ponce and Renon, 1976; Lee and Chao, 1988), but none of them offer as good a compromise between simplicity and accuracy as eqn. (31). We used one of these equations of state (Lee and Chao, 1988) for the free energy of the square-well reference fluid, and computed the perturbation free energies &I,/N. The results were similar to the perturbation free energies from the LLS equation, and were fitted with eqn. (34). The resulting
Fig. 7. Phase diagram of the Lennard-Jones fluid. The symbols are the simulation data of Panagiotopoulos (1987). The solid line is the prediction of eqn. (23). The dotted line is the prediction of eqn. (35).
142 TABLE 4 Critical properties of the Lennard-Jones
fluid computed from eqns. (23) and (35)
Equation
K/c
PC03
ZC
Eqn. (23) Eqn. (35) Ada&i et al. (1988) Nicolas et al. (1979)
1.372 1.384 1.273 1.35
0.313 0.308 0.284 0.35
0.363 0.387 0.328 0.300
equation of state is no more accurate than eqn. (35), and slightly more complicated. It also cannot be written as a polynomial in density. Of special interest is the method used in deriving this equation, which is based on making systematic “corrections” to an equation of state for a simpler fluid required to fit a more realistic fluid. In this way we were able to extend the LLS square-well equation of state to the Lennard-Jones fluid by using generalized perturbation theory and incorporating temperature-dependent square-well parameters from the de Lonngi and de1 Rio theory. We have also computed the phase diagram of the Lennard-Jones fluid using eqns. (23) and (35). The values of the adjustable parameters used in both equations were those in the first lines of Table 2 and 3 respectively. A comparison of the equation-of-state phase diagrams and the simulation data of Panagiotopoulos (1987) is shown in Fig. 7. The critical properties computed from both equations are given in Table 4 and compared with the critical properties from Ada&i et al. (1988). From this table and Fig. 7 we can see that both equations of state give nearly identical critical t.emperatures and densities, and both overpredict the critical pressure. Equation (23) gives a slightly more accurate prediction of the phase diagram. APPLICATIONS
TO REAL FLUIDS OF SPHERICAL
MOLECULES
In this section we examine the utility of equations of state developed for the Lennard-Jones fluid in describing real fluids by using eqn. (23) and a modification of eqn. (35) to correlate the phase diagrams and supercritical volumes of methane, argon, krypton and xenon. For these four fluids, we considered the temperature range along the saturation line of 0.5 G T/T, G 0.95, and supercritical state points from T/T, = 1.1 to T/T, = 5.4 for methane, T/T, = 8.6 for argon, T/T, = 5.7 for krypton and T/T, = 4.14 for xenon. The sources of the experimental data were the IUPAC tables (Angus and Armstrong, 1971; Angus et al., 1978) for argon and methane respectively, and Rabinovich et al. (I988) for krypton and xenon. We first correlated the phase diagram and Pn behavior of real fluids using eqns. (23) and (35) with the values of the adjustable parameters used
143 TABLE 5 Equation of state parameters used to fit the phase diagrams of argon, krypton, xenon and methane Equation (23)
Equation (37)
a = 6.700 b, = 1.361 b,= -4.795x10-’ c=2.000
a, = - 1.325 a2 = - 3.191 b, = 6.982 b, = 1.920
X10-’
in
the last section. The values of c and 0 were adjusted to fit the experimental vapor pressures, saturated volumes, and supercritical volumes simultaneously. We found that the fit of supercritical PVT data was acceptable for both equations, but the phase diagrams were not fitted very accurately, especially in the near-critical region. Since the adjustable parameters in Tables 2 and 3 had been fitted to Lennard-Jones compressibilities in the supercritical fluid region, and no information about the phase diagram of the Lennard-Jones fluid was used to determine the parameters, this is not surprising. Also, both equations of state overpredict the critical pressure of the Lennard-Jones fluid with the parameters in Tables 2 and 3, and can be expected to perform badly near the critical point of real fluids. In order to get a more accurate fit of the phase diagram, it was necessary to adjust the parameters in both equations of state. New parameters a, b,, and b, in eqn. (23) were found by fitting the phase diagrams of the noble gases and methane. The value of the c parameter was taken from Table 2 and was not fitted to the phase diagram. The parameters given in Table 5 were used for all four fluids. Once these parameters are fitted to the phase diagram, they were used without change to predict supercritical volumes. Two modifications of eqn. (35) are necessary in order to get an accurate fit of the real fluid phase diagrams. First, the value of Z,,, was set equal to 36, twice the value used to fit Lennard-Jones data. This value of Z,,, is consistent with the value used in our previous work with the square-well fluid and with Kim et al. (1986) for more complex fluids. Second, when Z, is set equal to 36, it is also necessary to modify the temperature dependence of the p*2 term in eqn. (35) in order to obtain an accurate correlation of the phase diagram. The modified equation of state is z=
1+II+172-v3
(l-d3
-
Z, bpWk0 &/6
2(&)2]P*2+
- 11~
+ [ exp( c/2kT) [b2&
- l] 7
+b2(&)‘]p*
(37)
144 TABLE 6 Potential parameters used in eqns. (23) and (37) for methane, argon, krypton and xenon Gas
Equation
r/k(K)
a(A)
Gas
Equation
r/k(K)
a(A)
Methane
Eqn. (23) Eqn. (37) From B,
159.5 164.1 148.2
3.40 3.41 3.82
Argon
Eqn. (23) Eqn. (37) From B2
125.3 128.8 119.8
3.10 3.11 3.40
Krypton
Eqn. (23) Eqn. (37) From B,
174.6 179.5 171.0
3.31 3.32 3.60
Xenon
Eqn. (23) Eqn. (37) From B,
242.3 248.9 221.0
3.61 3.61 4.10
where 2, = 36 and a,, a2, b, and b, are adjustable parameters fitted to the phase diagram. This equation of state differs from eqn. (35) only in the value of Z, and the coefficient of the Pan term. The parameters a,, a,, b, and b, in eqn. (37) were found by fitting the phase diagrams of the noble gases and methane. Again, the parameters given in Table 5 were used for all four fluids, and they are used without change to predict supercritical volumes. The values of the Lennard-Jones potential parameters c/k and u fitted to the phase diagram are shown in Table 6. The Lennard-Jones potential parameters determined from second virial coefficients (Hirschfelder et al., 1954) are included for comparison. The values of the diameter required to fit the phase diagrams of all four fluids are substantially smaller than the values determined from second virial coefficients, while the energy is somewhat larger. In Table 7 we compare eqns. (23) and (37) with two popular engineering equations of state. Also included in the table are results from the original LLS square-well equation of state with e/k and u adjusted to fit the phase diagram, and the maximum coordination number Z,,, set equal to the constant 36. The percentage errors for the vapor pressure, saturated liquid and vapor volumes and supercritical volumes presented in Table 7 are defined for any property X as
(38) where N is the number of experimental data points. We see from Table 7 that both eqns. (23) and (37) provide quite reasonable correlations of phase behavior. Both compare favorably with the empirical Peng-Robinson (1976) and Teja-Pate1 (Pate1 and Teja, 1982) equations in correlations of the vapor pressure and saturated vapor volume, and are better for saturated liquid volumes. In particular, we point out the extremely low errors in saturated liquid volume from eqn. (23).
145 TABLE 7 Comparison of equation of state correlations of the phase diagrams and supercritical volumes of methane, argon, krypton and xenon
v VAP x
VSC%
Gas
Equation
P%
VL'Q%
Argon
Eqn. (23) Eqn. (37) LLS Peng-Robinson Teja-Pate1
1.05 1.47 3.54 0.31 0.75
0.19 1.10 3.57 9.08 3.01
1.81 0.73 ,2.41 1.26 1.20
3.39 1.44 5.37 2.79 0.71
Krypton
Eqn. (23) Eqn. (37) LLS Peng-Robinson Teja-Pate1
1.04 1.19 3.30 0.55 1.13
0.20 0.97 3.44 8.07 3.12
1.64 0.75 2.51 1.81 1.53
3.83 2.31 5.09 3.44 0.84
Xenon
Eqn. (23) Eqn. (37) LLS Peng-Robinson Teja-Pate1
0.98 1.27 3.20 1.53 2.15
0.17 0.92 3.24 7.16 3.47
1.80 0.74 1.99 2.39 2.18
3.98 2.80 5.07 3.43 1.04
Methane
Eqn. (23)
1.38 1.32 4.18 1.54 1.18
0.33 0.91 4.00 8.07 3.54
2.71 2.01 3.17 2.08 1.38
4.14 2.48 4.39 3.51 1.19
2 (37) Peng-Robinson Teja-Pate1
Our purpose in this section was not to present yet another equation of state for atomic fluids, but instead to demonstrate that the equations of state derived from generalized perturbation theory for the Lennard-Jones fluid could be used to correlate phase behavior of real fluids of spherical molecules. Although the adjustable parameters in eqn. (23) and the functional form of part of eqn. (35) had to be changed in order to fit the phase diagrams of the noble gases and methane, one should not conclude from this that the equations of state presented here are poor representations of the Lennard-Jones fluid. Evidence has existed for many years (Barker and Henderson, 1976) that the Lennard-Jones potential is a good effective potential for fluids like argon, but very little quantitative information is available concerning the effect of inaccuracy in the intermolecular potential function on the phase diagram. Also, some researchers have argued (Aim and Nezbeda, 1983) that three-body contributions to the intermolecular potential energy may account for as much as 50% of the total pressure. Clearly, intermolecular potentials more accurate than the Lennard-Jones are
146
necessary to describe the phase diagram of noble gases with a high degree of accuracy.
CONCLUSIONS
We have successfully applied generalized perturbation theory to the development of two new equations of state for the Lennard-Jones fluid. One of these equations (eqn. (23)) is more accurate than the 33-parameter Ada&i equation for the data sets considered here, even though it has only four adjustable parameters. We have shown that the choice of the hard-sphere diameter in a perturbed-hard-sphere equation of state has a small but significant effect on the density dependence of the generalized perturbation free energies. The use of the BH diameter produces generalized perturbation free energies that are almost linear functions of density, and leads to a simple van der Waals-like equation of state. However, when the correction term PA ‘ORR/N representing the deviation PAP/N from linearity is computed, the result is a complex function of density, and is difficult to fit. In contrast, when the HBH or WCA diameter is used in the calculation of PA ‘ORR/, the correction term is simpler and can be fitted with the same function used for the HCLJ fluid. Therefore, the use of either of these two diameters in the calculation of the free energy of the reference system in perturbation theory of the Lennard-Jones fluid leads to behavior of /3A,/N that is quantitatively similar to that of the HCLJ fluid. This suggests that these diameters are more realistic approximations of the “hard-core diameter” for the Lennard-Jones fluid, a conclusion that has been reached earlier by Aim and Nezbeda (1987) in a somewhat different way. The principal new results are eqns. (23) and (35), the latter of which is a straightforward, theoretically based extension of the LLS square-well equation of state to the Lennard-Jones fluid. The importance of these new equations is that their functional forms and their temperature-dependent parameters were found by considering the behavior of PAP/N in generalized perturbation theory. Simulation compressibilities were then used to determine the values of the adjustable parameters in the equations of state. Both equations are able to correlate the phase diagram and supercritical volumes of real atomic fluids, provided the adjustable parameters are fitted to the phase diagrams. While neither of these equations of state is completely predictive for Lennard-Jones or real atomic fluids, both are theoretically based, analytic, and have only a few adjustable parameters. At present, there are no statistical mechanical theories that allow the derivation of a predictive equation of state which is analytic, reasonably simple, and accurate enough for engineering purposes, even for the simple
147
fluids considered in this paper. Perturbation and integral equation theories do give accurate predictions of the PVT properties of the Lennard-Jones fluid, but require complicated numerical solutions. They also fail to give an accurate prediction of the phase diagram of either Lennard-Jones or real atomic fluids near the critical point. For non-spherical or short-chain molecular fluids such as alkanes, the theories are more complicated and less accurate. The generalized perturbation theory used here provides a common framework for developing equations of state for all of these classes of fluids.
ACKNOWLEDGMENTS
This work was supported by grant number FG-85ER13436 from the United States Department of Energy to the University of Delaware. We also thank the Pittsburgh Supercomputer Center, the University of Delaware Computing Center and the National Science Foundation for generous grants of computing time and equipment.
LIST OF SYMBOLS ACONF ACONF
A0
B,P d ECONF
&W
U UO W
YHS =0
z cs =nl
configurational Hehnholtz free energy configurational Hehnholtz free energy of the reference system generalized perturbation free energy second virial coefficient hard-core diameter configurational internal energy radial distribution function range parameter in the square-well potential intermolecular potential intermolecular potential of the reference fluid perturbation potential background correlation function of the hard-sphere fluid compressibility of the reference system compressibility of a hard-sphere fluid computed from the Carnahan-Starling equation maximum coordination number
Greek letters
P 6
l/kT, where k is the Boltzmann constant characteristic energy parameter in an intermolecular potential
148
: P P*
CT
packing fraction perturbation parameter number density reduced number density characteristic length parameter in an intermolecular
potential
REFERENCES
Ada&i, Y., Fijhara, I., Takamiya, M. and Nakanishi, K., 1988. Generalized equation of state for Lennard-Jones fluids. I. Pure fluids and simple mixture. Fluid Phase Equilibria, 39: l-38. Aim, K. and Nezbeda, I., 1983. Perturbed hard-sphere equations of state of real liquids. I. Examination of a simple equation of the second order. Fluid Phase Equilibria, 12: 235-251. Aim, K. and Nezbeda, I., 1984. Perturbed hard-sphere equations of state of real fluids. II. Effective hard-sphere diameters and residual properties. Fluid Phase Equilibtia, 17: 1-18. Aim, K. and Nezbeda, I., 1987. Perturbed hard-sphere equations of state of real fluids. III. Residual parameter up of non-polar liquids. Fluid Phase Equilibria, 34: 171-188. Allen, M.P. and Tildesley, D.J., 1987. Computer simulation of liquids. Clarendon Press, Oxford. Anderson, H.C., Weeks, J.D. and Chandler, D., 1971. Relationship between the hard-sphere fluid and fluids with realistic repulsive forces. Phys. Rev. A, 4: 1597-1607. Angus, S. and Armstrong, B. (Eds.), 1971. International thermodynamic tables of the fluid state, argon. Butterworths, London. Angus, S., Armstrong, B. and de Reuck, K.M. (Eds.), 1978. International thermodynamic tables of the fluid state-5 Methane. Pergamon Press, Oxford. Barker, J.A. and Henderson, D., 1967. Perturbation theory and equation of state for fluids. II. A successful theory of liquids. J. Chem. Phys., 47: 4714-4721. Barker, J.A. and Henderson, D., 1976. What is “liquid”? Understanding the states of matter. Rev. Mod. Phys., 48: 587-671. Beret, S. and Prausnitz. J.M., 1975. Perturbed hard-chain theory: an equation of state for fluids containing small or large molecules. AICHE J., 21: 1123-1132. de Lonngi, D.A. and de1 Rio, F., 1983. Square well perturbation theory of fluids. Mol. Phys., 48: 293-313. Dodd, L.R. and Sandler, S.I., 1989. Monte Carlo study of the Buckingham exponential-six fluid. Mol. Simul., 2: 15-32. Hansen, J.P. and McDonald, I.R., 1986. Theory of simple liquids. Academic Press, London. Hirschfelder, J.O., Curtiss, C.F. and Bird, R.B., 1954. Molecular theory of gases and liquids. John Wiley and Sons, Inc., New York. Kim, C.H., Vimalchand, P., Donohue, M.D. and Sandler, S.I., 1986. Local composition model for chainlike molecules: A new simplified version of the perturbed hard-chain theory. AICHE J., 32: 1726-1734. Lee, R.J. and Chao, K.C., 1988. Equation of state for square-well fluids. Mol. Phys., 65: 1253-1256. Lee, K.-H., Lombardo, M. and Sandler, S.I., 1985. The generalized van der Waals partition function. II. Application to the square-well fluid. Fluid Phase Equilibria, 21: 177-196.
149 Nezbeda, I. and Aim, K., 1989. On the way from theoretical calculations to practical equations of state for real fluids. Fluid Phase Equilibria, 52: 39-46. Nicolas, J.J., Gubbins, K.E., Streett, W.B. and Tildesley, D.J., 1979. Equation of state for the Lennard-Jones fluid. Mol. Phys., 37: 1429-1454. Panagiotopoulos, A.Z., 1987. Direct determination of phase coexistence properties of fluids by Monte Carlo simulation in a new ensemble. Mol. Phys., 61: 813-826. Patel, N. and Teja, A., 1982. A new cubic equation of state for fluids and fluid mixtures. Chem. Eng. Sci., 37: 463-473. Peng, D. and Robinson, D.B., 1976. A new two-constant equation of state. Ind. Eng. Chem. Fundam., 15: 59-64. Ponce, L. and Renon, H., 1976. Analytical equation for the Helmholtz free energy of a pure fluid, using the perturbation theory and a square well potential. J. Chem. Phys., 64: 638-640. Rabinovich, V.A., Vasserman, A.A., Nedostup, V.I. and Veksler, L.S., 1988. Thermophysical properties of neon, argon, krypton and xenon. Hemisphere Publishing, Berlin. Rowlinson, J.S. and Swinton, F.L., 1982. Liquids and liquid mixtures. Butterworths, London. Sandler, S.I., 1985. The general&d van der Waals partition function. I. Basic theory. Fluid Phase Equilibria, 19: 233-257. Sowers, G.M. and Sandler, S.I., 1991. Equations of state from generalized perturbation theory. I. The hard-core Lennard-Jones fluid. Fluid Phase Equilibria, 63: l-25. Verlet, L. and Weis, J., 1972. Equilibrium theory of simple liquids. Phys. Rev. A, 5: 939-952. Weeks, J.D., Chandler, D. and Anderson, H.C., 1971. Role of repulsive forces in determining the equilibrium structure of simple liquids. J. Chem. Phys., 54: 5237-5247. Zwanzig, R.W., 1954. High-temperature equation of state by a perturbation method. I. Nonpolar gases. J. Chem. Phys., 22: 1420-1426.
APPENDIX
We have used the Metropolis Monte Carlo algorithm (Allen and Tildesley, 1987) to simulate a system of 112 atoms interacting with the LermardJones potential. The simulation runs were started from a random initial configuration, averages were recorded after discarding 100 000 configurations, and the lengths of the simulation runs were 3 X lo6 configurations for pa3 G 0.3 and 1.5 X lo6 configurations for pa3 > 0.3. The compressibilities were computed from the relation N-l
N
C
C
r=l
y=r+1
rij
d”(rij)
3NkT + Z,nc
dr ‘J
where the angle brackets denote an ensemble average, and rij is the distance between atom i and atom j. The cut-off distance in the simulation is half of the box length, and long-range corrections Z,, are computed in the usual way. The compressibilities at the isotherms simulated are shown in Table Al. The compressibilities at a reduced density pa3 = 0.05 at all temperatures
150 TABLE Al Compressibilities of the LJ fluid kT/c
PO3
0.05 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
1.35
1.75
2.0
3.0
4.0
6.0
0.846 0.719 0.508 0.367 0.273 0.295 0.543 1.200 2.391 4.408
0.909 0.841 0.733 0.680 0.709 0.860 1.210 1.886 3.042 4.928
0.934 0.889 0.819 0.813 0.880 1.061 1.456 2.156 3.291 4.991
0.988 0.992 1.020 1.108 1.281 1.560 1.982 2.687 3.663 5.115
1.012 1.036 1.119 1.248 1.449 1.762 2.196 2.850 3.782 5.024
1.034 1.082 1.204 1.360 1.601 1.923 2.364 2.949 3.731 4.792
were computed from a virial expansion truncated after the first order term in density. The configurational internal energies were computed from N-l ECONF=
c r=l
N c
U(ri,)+ELRc
j=i+l
i and are shown in Table A2. The error in the simulation data is approximately 1% for the compressibility and 0.2% for the configurational internal energy.
TABLE A2 Configurational internal energies of the LJ fluid PO3
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
kT/r 1.35
1.75
2.0
3.0
4.0
6.0
-
-
-
-
-
- 0.475 - 0.934 - 1.384 - 1.781 - 2.126 - 2.377 - 2.508 -2.444 - 2.086
0.771 1.468 2.114 2.735 3.376 4.024 4.662 5.243 5.687
0.696 1.356 1.979 2.609 3.222 3.846 4.436 4.941 5.264
0.660 1.295 1.910 2.528 3.137 3.738 4.295 4.749 5.047
0.592 1.172 1.740 2.302 2.850 3.366 3.790 4.102 4.556
0.546 1.079 1.605 2.115 2.596 3.025 3.342 3.492 3.443
Fluid Phase Equilibria, 67 (1991) 127-150 Elsevier Science Publishers B.V., Amsterdam
Equations of state from generalized perturbation Part II. The Lennard-Jones
theory.
fluid
Gary M. Sowers ’ and Stanley I. Sandler Department of Chemical Engineering, University of Delaware, Newark, DE 19716 (USA) (Received
August 20, 1990; accepted
in final form June 12, 1991)
ABSTRACT Sowers, G.M. and Sandler, S.I., 1991. Equations of state from generalized perturbation theory. Part II. The Lennard-Jones Fluid . Fluid Phase Equilibria, 67: 127-150. We use perturbation theory to derive extensions of our equations of state used previously for the hard-core Lennard-Jones fluid to the more realistic soft-core Lennard-Jones fluid. We find that the use of a “hybrid” Barker-Henderson diameter (see Aim and Nezbeda, 1987) for the hard sphere reference fluid leads to perturbation free energies for the soft-core LennardJones fluid that are similar to the perturbation free energies for a hard-core Lennard-Jones fluid. For the range of temperatures and densities considered here, one of our equations of state is more accurate than a 33-parameter equation (see Ada&i et al., 1988) even though it has only four adjustable parameters. The other equation of state is a straightforward, theoretically based extension of a square-well equation of state (see Lee et al., 1985) to the Lennard-Jones fluid. Both equations of state can be used to correlate the phase diagrams and supercritical PVT data of noble gases and methane, and perform especially well in the correlation of saturated liquid volumes.
INTRODUCTION
The development of a simple, analytic and accurate equation of state for a realistic fluid from statistical mechanics is a difficult problem, even for monatomic fluids interacting with pair potentials as simple as that of Lennard-Jones. Since the thermodynamic and structural properties of a simple fluid can be computed from perturbation and integral equation theories with an acceptable degree of accuracy, some regard the problem of i Present
address:
037%3812/91/$03.50
Dow Chemical
Company,
Plaquemine,
0 1991 Elsevier Science Publishers
LA 70765, USA. B.V. All rights reserved
128
describing these properties from statistical mechanics as solved (Row~nson and Swinton, 1982). However, in engineering applications accurate and simple equations of state with few or no adjustable parameters are required, and perturbation or integral equation theories are of little direct use because they can be mathematically involved and at present are applicable only to relatively simple model fluids. As pointed out by Nezbeda and Aim (1989), a simple, analytic and completely predictive equation of state for the Lennard-Jones fluid that can accurately describe both PVT properties and the phase diagram has not yet been derived from statistical mechanics. There is no theoretically based equation of state that can predict the thermodynamic properties and phase behavior of the Lennard-Jones monatomic fluid more accurately than the empirical 33-parameter equation of state of Nicolas et al. (1979). Although simulation techniques for computing single phase thermodyn~c properties of simple fluids are quite well developed and an improved method of computing the phase diagram of a model fluid from simulation has recently been devised (Panagiotopoulos, 1987), these methods do not provide an equation of state in analytic form. Therefore, from an engineering point of view, the problem of developing an analytic equation of state for the Lennard-Jones fluid from statistical mechanics is still unsolved. Our approach to equation of state development has involved two statistical m~hanical tools; generalized van der Waals theory (Sandler, 1985) and perturbation theory (Hansen and McDonald, 1986). Generalized van der Waals theory provides an expression for the canonical partition function of monatomic or rigid molecular fluids and their mixtures, and can be extended to fluids of flexible chain molecules using the perturbed-hard-chain partition function of Beret et al. (1975). Generalized perturbation theory, based on the h-expansion, was presented in an earlier paper (Sowers and Sandler, 1991). We used the expression for the confi~rational Hehnholtz free energy from generalized perturbation theory as a starting point for our study of the hard-core Lennard-Jones (HCLJ) fluid r
co
-I5y %cLJ(r)= &[( i
_ (;,“I
;:;lz:”
This fluid differs from the soft-core Lermard-Jones potential by having an infinitely steep repulsive part at r = u. It is a useful potential for evaluating perturbation theories because simpler potentials with hard cores, such as the hard-sphere or square-well potential, can be used as reference systems without ambiguity. In contrast, the reference system in perturbation theories for the soft-core Lennard-Jones fluid is softly repulsive, and is typically
129
approximated by a hard-sphere potential with a state-dependent diameter. In our earlier work we presented two equations of state for the HCLJ fluid that fit our computer simulation compressibilities to within 3%. Both of these equations of state were developed from the generalized perturbation theory expression for the configurational Helmholtz free energy PA CoNF = ,AtoNF + PA, N
N
N
(2)
where flAgoNF/N is the configurational free energy of the reference system in the perturbation theory, p = l/kT and /3A,/N is the generalized perturbation free energy (see the third section). The first equation of state results from using a hard-sphere reference system for pAzoNF/N, while the second is novel because it was developed using a square-well reference system. In this paper, we apply generalized perturbation theory to the more realistic soft-core Lennard-Jones fluid, and show how our HCLJ equations of state can be extended to this fluid in a simple and systematic way. In the next section we present a brief review of generalized perturbation theory, while in the third and fourth sections we discuss the existing perturbation theories for the Lennard-Jones fluid that use hard-sphere and square-well reference systems respectively, and how these theories can be used to develop equations of state. The use of a hard-sphere reference system in a perturbation theory of the Lennard-Jones fluid leads to an equation of state very similar to the equation of state for the modified exponential-six fluid proposed by Dodd and Sandler (1989), and to the perturbed-hard-sphere equations of state studied by Aim and Nezbeda (1983, 1984, 1987) and others. In the third section we discuss the choices for a state-dependent hard-sphere diameter in perturbation theory with a hardsphere reference, and the effect of this choice on the equation of state. We present a new perturbed-hard-sphere equation that fits our computer simulation results for the compressibility of the Lennard-Jones fluid to within 2% error. We show, in the fourth section, how the use of perturbation theory leads to a straightforward extension to the Lennard-Jones fluid of the lattice gas model of the square well fluid derived by Lee et al. (1985). This new equation of state is also in good agreement with computer simulation results. In the following section we show that these equations can also be used to correlate the phase diagrams and supercritical PVT data of the noble gases and methane to a high degree of accuracy. THEORY
The generalized perturbation theory expression for the configurational Helmholtz free energy of the Lennard-Jones fluid is derived by first dividing
130
the intermolecular
pair potential into reference and perturbation
U(T) = z+(r) + Xw(r)
parts (3)
where ~a( r) is the reference potential and w(r) = uu( r) - t+(r) is the perturbation potential. When h = 0, the potential u(r) = Q(T), and when X = 1 we recover the Lennard-Jones potential Us from eqn. (3). At fixed N, V and T the change in the configurational Helmholtz free energy of a fluid whose pair potential is u(r) resulting from a change in X is
CONF PA
1 a ( /3AC0NF)
IA=1 =PACoNFI~_O + J0
dX
(4
ah
with
a ( fiACoNF) = ax
(5)
2nNpSjOOw(r)g(r;N,V,T,h)r2dr 0
Defining the generalized perturbation
free energy as
we can write the total configurational
Helmholtz free energy as
PA CoNF = ,8A;oNF + PA, N
N
(7)
N
Equation (7) may appear to be of little value because it requires knowledge of the radial distribution function g( r; N, V, T, X) throughout the entire range 0 < X < 1 to compute the configurational Helmholtz free energy of the fluid of interest. Such information is very difficult to obtain from simulation or from theory. An equation of state for a model fluid can be derived from perturbation theory by making assumptions that allow the terms in the perturbation expansion to be expressed analytically (see, for example, Aim and Nezbeda, 1983; Ponce and Renon, 1976), and limited computer simulation data can then be used to determine the values of the adjustable parameters in the equation of state. We use generalized perturbation theory here in a similar way to extend our equations of state for the HCLJ fluid to the soft-core Lennard-Jones fluid.
APPLICATION
OF GPT USING
A HARD-SPHERE
REFERENCE
SYSTEM
Perturbation theories for the Lennard-Jones fluid using a hard-sphere reference fluid originated with Zwanzig’s (1954) pioneering work. In this section we consider only perturbation theories for the Lennard-Jones fluid
131
in which the reference system is approximated by a hard-sphere fluid whose diameter depends on temperature or on both density and temperature. Two such theories are the perturbation expansions derived by Barker and Henderson (1967) and Weeks et al. (1971) respectively. In Barker-Henderson theory the Hehnholtz free energy is expanded in a double Taylor series in two parameters, a steepness parameter and a well depth parameter, and truncated after the second-order term. An expression for the temperaturedependent Barker-Henderson (BH) diameter of the hard-sphere reference fluid is found by setting one of the first-order terms in the expansion equal to zero, which leads to
dBH(T)= Lo1 - exp[-h_.dr)]dr where uw ( Y) is the Lennard-Jones potential. In Weeks-Chandler-Anderson (WCA) theory the potential is divided into reference and perturbation parts at the potential minimum. A first-order h-expansion of the Helmholtz free energy of the Lennard-Jones fluid is followed by a blip function expansion of the free energy of the soft repulsive reference system about a hard-sphere reference (Anderson et al., 1971). The reference system is approximated by a hard-sphere fluid whose diameter is determined by setting the first-order term in the blip function expansion equal to zero. The WCA diameter is computed by solving the equation
--Puo(r)]JoODyHs(r;dwC*)oxp[
exp[-PuH,(r;dwc,)l}v2dr=0
(9)
for dCcAP~J3, whereyHS
is the background correlation function of the hard-sphere fluid, uO is the pair potential of the soft repulsive reference system and uns( r; dWCA) is the hard-sphere potential with diameter dWCA. Finally, we can define a temperature-dependent hard-sphere diameter by dividing the potential into reference and perturbation parts at the potential minimum as in the WCA theory, and retaining only the zeroth-order term in a Taylor series expansion of r2y,, about r = d (Verlet and Weis, 1972). The resulting temperature-dependent diameter has been referred to as a “hybrid” Barker-Henderson (HBH) diameter (Aim and Nezbeda, 1987), dHBH(T)
=
1
2”601- exp[ -p{ uU(r) 0
+ E}]dr
00)
The free energy of the soft repulsive reference system in each of these perturbation theories for the Lennard-Jones fluid is then CONF
P-&l N
=
~AC$NF
N
(11)
132
where the right-hand side is the configurational Hehnholtz free energy of a hard-sphere fluid whose diameter is computed from the appropriate one of eqns. (S)-(10). If the total configurational Helmholtz free energy /UCoNF/IV is known, eqn. (7) together with eqn. (11) and one of eqns. (8)-(10) can be used to compute the generalized perturbation free energies &4,/N, as follows. Subtracting the relations PA
CONF
N
PZ-1 = o Tdp’
(12)
“2,-l = aP’dd
(13)
/
and CONF
PA0 N
J
we find, by comparison with eqn. (7) that
PA,=
N
“Z-Z,
J p’dp’ o
(14)
where 2 is the compressibility of the Lennard-Jones fluid and 2, is the compressibility of the hard-sphere reference system. Since none of our simulation points are in the two-phase region, we can use the simulation compressibilities in eqn. (14) to compute PAP/N. There is a large body of simulation data available for the Lennard-Jones fluid, but previous studies (Nicolas et al., 1979; Adachi et al., 1988) have not included state points over a wide enough density and temperature range for our purposes. Consequently, we have done additional simulations of the Lennard-Jones fluid in the temperature range 1.35 < kT/c < 6.0 and 0 < pa3 < 0.9. A brief discussion of the simulation method and the results for compressibility and configurational internal energy are given in the Appendix. Although we use only our simulation data to develop the equations of state presented here, we later compare the resulting equations of state with the simulation data of Ada&i et al. (1988). We have used the Carnahan-Starling equation of state z
=
CS
l+77+772-v3 (1 -
d3
(15)
to calculate the hard-sphere compressibility 2,. We use 17 to denote the packing fraction of the hard-sphere fluid, defines as 7 = ( 7r/6)pd3 where d is the hard-sphere diameter. The generalized perturbation free energies /3A,/N calculated from eqn. (14) using eqns. (8), (9) and (10) for the hard-sphere diameter are plotted at a representative temperature kT/e = 1.75 in Fig. 1. The dependence of PAP/N on density is similar at the other five isotherms studied.
133 0.0
a
1
En
la
-1.0
0
4
t
4
-2.0
4
d
0
4 0
-3.0
BHdiameter
cl
4
0
A WCA diameter
0 9
+ HBH diameter
’
-5.0
I 0.2
00
0.4
0.6
0.8
1.0
PO3
Fig. 1. Generalized perturbation free energies computed using a hard-sphere reference fluid with a state-dependent diameter at a representative temperature kT/c = 1.75.
We see from Fig. 1 that /?A,/N is almost, but not quite, a linear function of density regardless of the choice of hard-sphere diameter. However, if the Barker-Henderson diameter is used in the computation of the generalized perturbation free energy, the deviation from linearity is smaller. As a first approximation, we will assume that the generalized perturbation free energies computed using the BH diameter can be fitted with a linear function of density
PAP d
nx
N
=
_n&p*
06)
where p* = pa3 and a is an adjustable parameter determined by fitting to simulation data. When this model for /3A,/N, which we have denoted by &4r,JN, is combined with the Camahan-Starling equation of state for the hard-sphere reference free energy PA,‘ONF/N, the resulting equation of state is z=
1+9+q2-v3 -a&p* (1 -
Id3
07)
which has been referred to as a perturbed-hard-sphere (Aim and Nezbeda, 1987) or Camahan-Starling-van der Waals equation of state. The adjustable parameter a is fitted to the computer simulation compressibilities. In all of our calculations we use the following approximant, accurate to within 0.02%, for the Barker-Henderson diameter: &-i(T) u
1.013/m = 4.142 x 1O-2 + ,I-
in order to avoid repeated numerical integration in eqn. (8).
(18)
134 TABLE 1 Comparison of eqn. (17) and computer simulation cbmpressibilities Data source
a
AZ%
Appendix 1 Ada&i et al. (1988)
6.693 6.482
4.59 9.77
We have determined the parameter a using two different sets of computer simulation data for the compressibility of the Lennard-Jones fluid. The first data set is our simulation data at six supercritical isotherms tabulated in Appendix 1. The second data set is taken from Adachi et al. (1988). State points that may lie in the metastable or two-phase region were removed from Ada&i’s data set by deleting all points below the critical temperature and pressure. The values of the parameter a in eqn. (17) were fitted to each of these data sets and are shown in Table 1. The error between the equation of state and the simulation compressibilities is defined as (19) where ZSrM is the simulation compressibility, ZEoS is the equation of state compressibility and N is the number of state points. We also tested the Nicolas equation of state for the Lennard-Jones fluid using the parameters taken from Ada&i et al. (1988). Although the functional form of the equation was proposed by Nicolas et al. (1979), we will refer to the equation of state with Adachi’s parameter values as the Ada&i equation to avoid confusion. The Ada&i equation produces errors of 3.9% and 3.8% when compared to the subset of Ada&i’s data used here and our simulation data, respectively. If the 33 parameters in the Ada&i equation of state are fitted to the subset of Ada&i’s data used here, the error between the equation of state and the simulation data decreases to 2.26%. The error between eqn. (17) and the simulation data is considerably larger than the error from the Ada&i equation, although eqn. (17) is much simpler. One can improve the performance of this equation of state, without sacrificing too much simplicity, in the,following way. The assumption made in eqn. (16) is that PAP/N is a linear function of density. We improve on this assumption by calculating the correction term j?ACoRR/N necessary to account for the non-linearity in PAP/N
PA‘ORR = MCoNF _ PArF N
N
where Zcs is the compressibility
_+
= /P
,=a
= -
0
of a hard-sphere
dp,
+
+
p
fluid computed
p*
(20)
from the
135
c x
T’=iOO T’=600
-0.15 0.0 0.1
0.2 0.3
Fig. 2. Correction
0.4 0.5
0.6
0.7
PO3 term j3ACoRR/N
0.6 0.9
from eqn. (20) using the BH diameter
in Z,,.
Car&ran-Starling equation of state. This correction term is plotted against density at six isotherms in Fig. 2. We find that it is a complicated function of density and is difficult to fit with a simple function. However, if we use the HBH diameter in Zc, in eqn. (20), the correction term has the form shown in Fig. 3 (the curves have similar shapes if the WCA diameter is used). This correction term can be fitted with the following function of density, which we have used for the HCLJ fluid (Sowers and Sandler, 1991) :
03
8
0.2
2 01
0.0 0.0
0.2
0.4
0.6
0.6
I.0
Fig. 3. Correction term BACoRR /N from eqn. (20) using the HBH diameter lines are the correlations of eqn. (21).
in Z,,.
The solid
136
where pt = 1.0 is the highest reduced density we have considered, and bi, b2 and c are adjustable parameters. The correlations of /UCoRR/N from this model are the solid lines in Fig. 3. The resulting model for the generalized perturbation theory free energies using a hard-sphere reference fluid with a temperature-dependent diameter is CORR
PAPmcd
-=
PAN
--a&p*+
N
(22)
The equation of state for the Lennard-Jones energy, then becomes ==
fluid, obtained from the free
l+7)+Y2-v3 -u&p* (1 -
d3
-cjb,i;Ii,+b2[3(&)2+
G])
x (2p*2 - P*,P*) exp[ -cp*(p* where 17= ( a/6)pd&(
4-m,(T) = (I
- P*,)]
T). For d,,,(T)
(23)
we have used an approximant
21’63m
(24)
0.11 + QE
which is accurate to within 0.05%, to avoid the numerical integration in eqn. (10). The four adjustable parameters in the equation of state (a, b,, b, and c), and the errors between the equation of state and the simulation data are shown in Table 2. The error between eqn. (23) and the simulation compressibilities is almost within the simulation error. It is lower than the error produced by the Ada&i equation of state for the data sets studied here, even though eqn. (23) has only four adjustable parameters, compared with the 33 adjustable parameters in the Ada&i equation. We should point out, however, that the range of density and temperature over which the Ada&i equation is applicable is somewhat larger than the range considered here. The Ada&i equation covers the temperature range 0.5 < kT/c < 6.0 and the density range
TABLE 2 Comparison of eqn. (23) and computer simulation compressibilities Data source
a
b,
Appendix 1 Ada&i et al. (1988)
8.042 7.781
0.731 0.6491
bz
-2.777x10-’ 1.59ox1o-2
C
AZ%
2.000 1.727
1.52 2.16
137 60
N
5.0
0
T’=lOO
4.0
x v X
T-=300 T’=4.00 T’= 6.00
30
0.4
0.2
0.0
0.6
0.6
1.0
Pa3
Fig. 4. Correlations of computer simulation compressibilities
from eqn. (23).
0 < pa3 < 1.2, whereas the parameters in our equations of state were fitted to the region 1.35 < kT/c < 6.0 and 0 < pa3 < 0.9. The correlations of our simulation compressibilities from eqn. (23) are shown in Fig. 4. APPLICATION
OF GPT USING A SQUARE-WELL
REFERENCE
SYSTEM
Both the BH and WCA theories use the hard-sphere fluid as a reference system in the perturbation expansion. Alternatively, since the square-well potential can be used to qualitatively describe the thermodynamic properties of more realistic fluids, one would expect that the perturbation free energy would be smaller when using the square-well potential as a reference system for the Lennard-Jones fluid. A square-well reference perturbation theory of Lennard-Jones fluids was developed by de Lonngi and de1 Rio (1983). They defined a modified potential function v(r;dsw,Rsw,R,,X1,X2)
=
r < d,
UU (z) -e
d,
i %J(W)
RI < r
(25)
where z
=
d,,
+
r-&V
Xl r-R,, x
w=R,+
(26)
2
4 = dsw R, = R,,
+
h(Rm -
dsw)
+ A,( R, - r)
138
and R, = 2’6. The smoothness parameters X, and A, vary the steepness of the repulsive (r < R,) and attractive (r > R,) regions of u(r), respectively. When X, = X, = 0, u(r) is the square-well potential with parameters d R,, and E, and when X, = X, = 1, v( r ) is equal to the Lennard-Jones pL:ntial uLJ(r) with potential parameters E and u. When the configurational free energy of the Lermard-Jones fluid is expanded in a double Taylor series about the square-well reference fluid, corresponding to the point A, = h, = 0, we obtain
PACoNF = bAzoNF N
N
+
il [a(%?‘]+,, (27)
where /3A:ONF/N is the configurational free energy of the square-well reference system. De Lonngi and de1 Rio provide expressions for both firstand second-order terms in the expansion that depend on the radial distribution functions of the square-well reference fluid and their derivatives. The first-order terms A, and A, and the mixed second-order term A,, were set to zero by choosing an effective square-well diameter dsw( T) and an effective square-well range R,,(T) according to
and R,,(T)
= 2%
+ [l -
ex~(~41-‘~~~~~ - exp[-h_h)]d~
(29)
Equation (27) can be written in the form of eqn. (7), provided we make the definition
(30)
1Al=0
PAP -= N
For the compressibility of the reference system Z,, we use the square-well equation of state derived by Lee et al. (1985) ==
1+v+772-773
(l-d3
-
Z, 7&/6
[exp(GkT) - 1117 + [exp( c/2kT)
- l] 9
where n = (r/6)pd&( T). As in the previous section, we approximate the integrals in eqns. (28) and (29) with the following simple analytic expres-
139
-0.2
?=-0.4
2
-0
6
-0.6
0.0
0 +
T’=135 T’=175
0
T’=zoo
k Q
T’=3.00 T’=400
0.2
04
0.6
1.0
0.6
Fig. 5. Generalized perturbation free energies computed using a square-well reference fluid with a state-dependent diameter. The solid lines are the correlations of eqn. (34).
sions, which are accurate to within 0.05%:
(32) %v(T) a
=&w(T)
&w(T)
=
2116+
a
0.347kT/c 0.26 + kT/c
(33)
As a first approximation we set Z,, the maximum coordination number, equal to the constant 18, corresponding to h,, = 1.5. It is important to note that this approximation removes the dependence on X,, from the reference fluid equation of state. The compressibility of the square-well reference system is then given by eqns. (31) and (32). The generalized perturbation free energies PA,/N computed using the square-well reference system are shown in Fig. 5. They can be fitted to an acceptable degree of accuracy by the following quadratic function of density, used for the HCLJ fluid: *=a(&)exp(-E)p*2+[bl(&)+b2(&)l]P*
(34)
where a, b, and b, are adjustable parameters. The correlations of &4,/N from this model are the solid lines in Fig. 5. The equation of state resulting from this model for /3A,/N together with the LLS equation of state for the
140 TABLE 3 Comparison of eqn. (35) and computer simulation compressibilities Data source
a
b,
b,
AZ%
Appendix 1 Ada&i et al. (1988)
- 9.459 - 8.374
- 1.783 - 0.7127
4.581 2.947
3.03 4.18
reference system is Z=
I+17+17*-V3
(1- ?I3 +2n(&)
_
2, [exp(e/2kT) &/6 + [exp(c/XT)
exp( -g)P**+
[4(G)
1177 - I] 17 +4(&)*]p*
(35)
The adjustable parameters in this equation are determined by fitting to simulation compressibilities and are shown in Table 3. The absolute average percent deviations in compressibility, also given in Table 3, are higher than eqn. (23), but still quite reasonable. The correlations of the simulation compressibilities are shown in Fig. 6. An approximation made in eqn. (35) is that the maximum coordination number 2, is a constant. We also examined whether it was possible to improve the performance of eqn. (35) by incorporating &,-dependence into 2,. In the LLS lattice model, 2, is a discontinuous function of the width of the square-well that is calculated from geometry. Since A,, is a function of temperature in the square-well reference perturbation theory, it is desirable to represent 2, as a continuous function of Xs,. One way of doing this is
6.0
5.0
4.0
0 + 0 x v *
I
l-=1.35 T-=1.75 T'=Z.OO T-=300 T-=400 T-=600
N
0.0
0.2
0.4
0.6
0.8
1.0
PO3 Fig. 6. Correlations of computer simulation compressibilities from eqn. (35).
141 to determine 2, by fitting the LLS second virial coefficient to the exact result for the square-well fluid, from which we obtain
&) - l]/[ev( &i=)- 11 -cm =~(&,(T)3- l)[exp(
(36)
(Note that this expression for 2, is undefined at infinite temperature. However, when used in eqn. (35), 2, is multiplied by 9 a d&(T), and this product goes to zero at infinite temperature.) When this model for Z, is used in eqn. (35) and the adjustable parameters are fitted to the simulation data, the performance of the equation of state is no better than when Z, = 18. Therefore we have used Z,,, = 18 in all of our calculations. The advantage of eqn. (35) over eqn. (23) is that it can be written as a polynomial in density. This reduces the computational effort when the equation of state must be solved for the densities of two coexisting phases during a phase equilibrium calculation. Although the deviation between eqn. (35) and the simulation compressibilities is approximately twice that of eqn. (23), the use of the square-well reference system has allowed us to use a very simple model for /3A,/N, and yet develop an accurate equation of state. Other equations of state for the square-well fluid have been proposed (Aim and Nezbeda, 1983; Ponce and Renon, 1976; Lee and Chao, 1988), but none of them offer as good a compromise between simplicity and accuracy as eqn. (31). We used one of these equations of state (Lee and Chao, 1988) for the free energy of the square-well reference fluid, and computed the perturbation free energies &I,/N. The results were similar to the perturbation free energies from the LLS equation, and were fitted with eqn. (34). The resulting
Fig. 7. Phase diagram of the Lennard-Jones fluid. The symbols are the simulation data of Panagiotopoulos (1987). The solid line is the prediction of eqn. (23). The dotted line is the prediction of eqn. (35).
142 TABLE 4 Critical properties of the Lennard-Jones
fluid computed from eqns. (23) and (35)
Equation
K/c
PC03
ZC
Eqn. (23) Eqn. (35) Ada&i et al. (1988) Nicolas et al. (1979)
1.372 1.384 1.273 1.35
0.313 0.308 0.284 0.35
0.363 0.387 0.328 0.300
equation of state is no more accurate than eqn. (35), and slightly more complicated. It also cannot be written as a polynomial in density. Of special interest is the method used in deriving this equation, which is based on making systematic “corrections” to an equation of state for a simpler fluid required to fit a more realistic fluid. In this way we were able to extend the LLS square-well equation of state to the Lennard-Jones fluid by using generalized perturbation theory and incorporating temperature-dependent square-well parameters from the de Lonngi and de1 Rio theory. We have also computed the phase diagram of the Lennard-Jones fluid using eqns. (23) and (35). The values of the adjustable parameters used in both equations were those in the first lines of Table 2 and 3 respectively. A comparison of the equation-of-state phase diagrams and the simulation data of Panagiotopoulos (1987) is shown in Fig. 7. The critical properties computed from both equations are given in Table 4 and compared with the critical properties from Ada&i et al. (1988). From this table and Fig. 7 we can see that both equations of state give nearly identical critical t.emperatures and densities, and both overpredict the critical pressure. Equation (23) gives a slightly more accurate prediction of the phase diagram. APPLICATIONS
TO REAL FLUIDS OF SPHERICAL
MOLECULES
In this section we examine the utility of equations of state developed for the Lennard-Jones fluid in describing real fluids by using eqn. (23) and a modification of eqn. (35) to correlate the phase diagrams and supercritical volumes of methane, argon, krypton and xenon. For these four fluids, we considered the temperature range along the saturation line of 0.5 G T/T, G 0.95, and supercritical state points from T/T, = 1.1 to T/T, = 5.4 for methane, T/T, = 8.6 for argon, T/T, = 5.7 for krypton and T/T, = 4.14 for xenon. The sources of the experimental data were the IUPAC tables (Angus and Armstrong, 1971; Angus et al., 1978) for argon and methane respectively, and Rabinovich et al. (I988) for krypton and xenon. We first correlated the phase diagram and Pn behavior of real fluids using eqns. (23) and (35) with the values of the adjustable parameters used
143 TABLE 5 Equation of state parameters used to fit the phase diagrams of argon, krypton, xenon and methane Equation (23)
Equation (37)
a = 6.700 b, = 1.361 b,= -4.795x10-’ c=2.000
a, = - 1.325 a2 = - 3.191 b, = 6.982 b, = 1.920
X10-’
in
the last section. The values of c and 0 were adjusted to fit the experimental vapor pressures, saturated volumes, and supercritical volumes simultaneously. We found that the fit of supercritical PVT data was acceptable for both equations, but the phase diagrams were not fitted very accurately, especially in the near-critical region. Since the adjustable parameters in Tables 2 and 3 had been fitted to Lennard-Jones compressibilities in the supercritical fluid region, and no information about the phase diagram of the Lennard-Jones fluid was used to determine the parameters, this is not surprising. Also, both equations of state overpredict the critical pressure of the Lennard-Jones fluid with the parameters in Tables 2 and 3, and can be expected to perform badly near the critical point of real fluids. In order to get a more accurate fit of the phase diagram, it was necessary to adjust the parameters in both equations of state. New parameters a, b,, and b, in eqn. (23) were found by fitting the phase diagrams of the noble gases and methane. The value of the c parameter was taken from Table 2 and was not fitted to the phase diagram. The parameters given in Table 5 were used for all four fluids. Once these parameters are fitted to the phase diagram, they were used without change to predict supercritical volumes. Two modifications of eqn. (35) are necessary in order to get an accurate fit of the real fluid phase diagrams. First, the value of Z,,, was set equal to 36, twice the value used to fit Lennard-Jones data. This value of Z,,, is consistent with the value used in our previous work with the square-well fluid and with Kim et al. (1986) for more complex fluids. Second, when Z, is set equal to 36, it is also necessary to modify the temperature dependence of the p*2 term in eqn. (35) in order to obtain an accurate correlation of the phase diagram. The modified equation of state is z=
1+II+172-v3
(l-d3
-
Z, bpWk0 &/6
2(&)2]P*2+
- 11~
+ [ exp( c/2kT) [b2&
- l] 7
+b2(&)‘]p*
(37)
144 TABLE 6 Potential parameters used in eqns. (23) and (37) for methane, argon, krypton and xenon Gas
Equation
r/k(K)
a(A)
Gas
Equation
r/k(K)
a(A)
Methane
Eqn. (23) Eqn. (37) From B,
159.5 164.1 148.2
3.40 3.41 3.82
Argon
Eqn. (23) Eqn. (37) From B2
125.3 128.8 119.8
3.10 3.11 3.40
Krypton
Eqn. (23) Eqn. (37) From B,
174.6 179.5 171.0
3.31 3.32 3.60
Xenon
Eqn. (23) Eqn. (37) From B,
242.3 248.9 221.0
3.61 3.61 4.10
where 2, = 36 and a,, a2, b, and b, are adjustable parameters fitted to the phase diagram. This equation of state differs from eqn. (35) only in the value of Z, and the coefficient of the Pan term. The parameters a,, a,, b, and b, in eqn. (37) were found by fitting the phase diagrams of the noble gases and methane. Again, the parameters given in Table 5 were used for all four fluids, and they are used without change to predict supercritical volumes. The values of the Lennard-Jones potential parameters c/k and u fitted to the phase diagram are shown in Table 6. The Lennard-Jones potential parameters determined from second virial coefficients (Hirschfelder et al., 1954) are included for comparison. The values of the diameter required to fit the phase diagrams of all four fluids are substantially smaller than the values determined from second virial coefficients, while the energy is somewhat larger. In Table 7 we compare eqns. (23) and (37) with two popular engineering equations of state. Also included in the table are results from the original LLS square-well equation of state with e/k and u adjusted to fit the phase diagram, and the maximum coordination number Z,,, set equal to the constant 36. The percentage errors for the vapor pressure, saturated liquid and vapor volumes and supercritical volumes presented in Table 7 are defined for any property X as
(38) where N is the number of experimental data points. We see from Table 7 that both eqns. (23) and (37) provide quite reasonable correlations of phase behavior. Both compare favorably with the empirical Peng-Robinson (1976) and Teja-Pate1 (Pate1 and Teja, 1982) equations in correlations of the vapor pressure and saturated vapor volume, and are better for saturated liquid volumes. In particular, we point out the extremely low errors in saturated liquid volume from eqn. (23).
145 TABLE 7 Comparison of equation of state correlations of the phase diagrams and supercritical volumes of methane, argon, krypton and xenon
v VAP x
VSC%
Gas
Equation
P%
VL'Q%
Argon
Eqn. (23) Eqn. (37) LLS Peng-Robinson Teja-Pate1
1.05 1.47 3.54 0.31 0.75
0.19 1.10 3.57 9.08 3.01
1.81 0.73 ,2.41 1.26 1.20
3.39 1.44 5.37 2.79 0.71
Krypton
Eqn. (23) Eqn. (37) LLS Peng-Robinson Teja-Pate1
1.04 1.19 3.30 0.55 1.13
0.20 0.97 3.44 8.07 3.12
1.64 0.75 2.51 1.81 1.53
3.83 2.31 5.09 3.44 0.84
Xenon
Eqn. (23) Eqn. (37) LLS Peng-Robinson Teja-Pate1
0.98 1.27 3.20 1.53 2.15
0.17 0.92 3.24 7.16 3.47
1.80 0.74 1.99 2.39 2.18
3.98 2.80 5.07 3.43 1.04
Methane
Eqn. (23)
1.38 1.32 4.18 1.54 1.18
0.33 0.91 4.00 8.07 3.54
2.71 2.01 3.17 2.08 1.38
4.14 2.48 4.39 3.51 1.19
2 (37) Peng-Robinson Teja-Pate1
Our purpose in this section was not to present yet another equation of state for atomic fluids, but instead to demonstrate that the equations of state derived from generalized perturbation theory for the Lennard-Jones fluid could be used to correlate phase behavior of real fluids of spherical molecules. Although the adjustable parameters in eqn. (23) and the functional form of part of eqn. (35) had to be changed in order to fit the phase diagrams of the noble gases and methane, one should not conclude from this that the equations of state presented here are poor representations of the Lennard-Jones fluid. Evidence has existed for many years (Barker and Henderson, 1976) that the Lennard-Jones potential is a good effective potential for fluids like argon, but very little quantitative information is available concerning the effect of inaccuracy in the intermolecular potential function on the phase diagram. Also, some researchers have argued (Aim and Nezbeda, 1983) that three-body contributions to the intermolecular potential energy may account for as much as 50% of the total pressure. Clearly, intermolecular potentials more accurate than the Lennard-Jones are
146
necessary to describe the phase diagram of noble gases with a high degree of accuracy.
CONCLUSIONS
We have successfully applied generalized perturbation theory to the development of two new equations of state for the Lennard-Jones fluid. One of these equations (eqn. (23)) is more accurate than the 33-parameter Ada&i equation for the data sets considered here, even though it has only four adjustable parameters. We have shown that the choice of the hard-sphere diameter in a perturbed-hard-sphere equation of state has a small but significant effect on the density dependence of the generalized perturbation free energies. The use of the BH diameter produces generalized perturbation free energies that are almost linear functions of density, and leads to a simple van der Waals-like equation of state. However, when the correction term PA ‘ORR/N representing the deviation PAP/N from linearity is computed, the result is a complex function of density, and is difficult to fit. In contrast, when the HBH or WCA diameter is used in the calculation of PA ‘ORR/, the correction term is simpler and can be fitted with the same function used for the HCLJ fluid. Therefore, the use of either of these two diameters in the calculation of the free energy of the reference system in perturbation theory of the Lennard-Jones fluid leads to behavior of /3A,/N that is quantitatively similar to that of the HCLJ fluid. This suggests that these diameters are more realistic approximations of the “hard-core diameter” for the Lennard-Jones fluid, a conclusion that has been reached earlier by Aim and Nezbeda (1987) in a somewhat different way. The principal new results are eqns. (23) and (35), the latter of which is a straightforward, theoretically based extension of the LLS square-well equation of state to the Lennard-Jones fluid. The importance of these new equations is that their functional forms and their temperature-dependent parameters were found by considering the behavior of PAP/N in generalized perturbation theory. Simulation compressibilities were then used to determine the values of the adjustable parameters in the equations of state. Both equations are able to correlate the phase diagram and supercritical volumes of real atomic fluids, provided the adjustable parameters are fitted to the phase diagrams. While neither of these equations of state is completely predictive for Lennard-Jones or real atomic fluids, both are theoretically based, analytic, and have only a few adjustable parameters. At present, there are no statistical mechanical theories that allow the derivation of a predictive equation of state which is analytic, reasonably simple, and accurate enough for engineering purposes, even for the simple
147
fluids considered in this paper. Perturbation and integral equation theories do give accurate predictions of the PVT properties of the Lennard-Jones fluid, but require complicated numerical solutions. They also fail to give an accurate prediction of the phase diagram of either Lennard-Jones or real atomic fluids near the critical point. For non-spherical or short-chain molecular fluids such as alkanes, the theories are more complicated and less accurate. The generalized perturbation theory used here provides a common framework for developing equations of state for all of these classes of fluids.
ACKNOWLEDGMENTS
This work was supported by grant number FG-85ER13436 from the United States Department of Energy to the University of Delaware. We also thank the Pittsburgh Supercomputer Center, the University of Delaware Computing Center and the National Science Foundation for generous grants of computing time and equipment.
LIST OF SYMBOLS ACONF ACONF
A0
B,P d ECONF
&W
U UO W
YHS =0
z cs =nl
configurational Hehnholtz free energy configurational Hehnholtz free energy of the reference system generalized perturbation free energy second virial coefficient hard-core diameter configurational internal energy radial distribution function range parameter in the square-well potential intermolecular potential intermolecular potential of the reference fluid perturbation potential background correlation function of the hard-sphere fluid compressibility of the reference system compressibility of a hard-sphere fluid computed from the Carnahan-Starling equation maximum coordination number
Greek letters
P 6
l/kT, where k is the Boltzmann constant characteristic energy parameter in an intermolecular potential
148
: P P*
CT
packing fraction perturbation parameter number density reduced number density characteristic length parameter in an intermolecular
potential
REFERENCES
Ada&i, Y., Fijhara, I., Takamiya, M. and Nakanishi, K., 1988. Generalized equation of state for Lennard-Jones fluids. I. Pure fluids and simple mixture. Fluid Phase Equilibria, 39: l-38. Aim, K. and Nezbeda, I., 1983. Perturbed hard-sphere equations of state of real liquids. I. Examination of a simple equation of the second order. Fluid Phase Equilibria, 12: 235-251. Aim, K. and Nezbeda, I., 1984. Perturbed hard-sphere equations of state of real fluids. II. Effective hard-sphere diameters and residual properties. Fluid Phase Equilibtia, 17: 1-18. Aim, K. and Nezbeda, I., 1987. Perturbed hard-sphere equations of state of real fluids. III. Residual parameter up of non-polar liquids. Fluid Phase Equilibria, 34: 171-188. Allen, M.P. and Tildesley, D.J., 1987. Computer simulation of liquids. Clarendon Press, Oxford. Anderson, H.C., Weeks, J.D. and Chandler, D., 1971. Relationship between the hard-sphere fluid and fluids with realistic repulsive forces. Phys. Rev. A, 4: 1597-1607. Angus, S. and Armstrong, B. (Eds.), 1971. International thermodynamic tables of the fluid state, argon. Butterworths, London. Angus, S., Armstrong, B. and de Reuck, K.M. (Eds.), 1978. International thermodynamic tables of the fluid state-5 Methane. Pergamon Press, Oxford. Barker, J.A. and Henderson, D., 1967. Perturbation theory and equation of state for fluids. II. A successful theory of liquids. J. Chem. Phys., 47: 4714-4721. Barker, J.A. and Henderson, D., 1976. What is “liquid”? Understanding the states of matter. Rev. Mod. Phys., 48: 587-671. Beret, S. and Prausnitz. J.M., 1975. Perturbed hard-chain theory: an equation of state for fluids containing small or large molecules. AICHE J., 21: 1123-1132. de Lonngi, D.A. and de1 Rio, F., 1983. Square well perturbation theory of fluids. Mol. Phys., 48: 293-313. Dodd, L.R. and Sandler, S.I., 1989. Monte Carlo study of the Buckingham exponential-six fluid. Mol. Simul., 2: 15-32. Hansen, J.P. and McDonald, I.R., 1986. Theory of simple liquids. Academic Press, London. Hirschfelder, J.O., Curtiss, C.F. and Bird, R.B., 1954. Molecular theory of gases and liquids. John Wiley and Sons, Inc., New York. Kim, C.H., Vimalchand, P., Donohue, M.D. and Sandler, S.I., 1986. Local composition model for chainlike molecules: A new simplified version of the perturbed hard-chain theory. AICHE J., 32: 1726-1734. Lee, R.J. and Chao, K.C., 1988. Equation of state for square-well fluids. Mol. Phys., 65: 1253-1256. Lee, K.-H., Lombardo, M. and Sandler, S.I., 1985. The generalized van der Waals partition function. II. Application to the square-well fluid. Fluid Phase Equilibria, 21: 177-196.
149 Nezbeda, I. and Aim, K., 1989. On the way from theoretical calculations to practical equations of state for real fluids. Fluid Phase Equilibria, 52: 39-46. Nicolas, J.J., Gubbins, K.E., Streett, W.B. and Tildesley, D.J., 1979. Equation of state for the Lennard-Jones fluid. Mol. Phys., 37: 1429-1454. Panagiotopoulos, A.Z., 1987. Direct determination of phase coexistence properties of fluids by Monte Carlo simulation in a new ensemble. Mol. Phys., 61: 813-826. Patel, N. and Teja, A., 1982. A new cubic equation of state for fluids and fluid mixtures. Chem. Eng. Sci., 37: 463-473. Peng, D. and Robinson, D.B., 1976. A new two-constant equation of state. Ind. Eng. Chem. Fundam., 15: 59-64. Ponce, L. and Renon, H., 1976. Analytical equation for the Helmholtz free energy of a pure fluid, using the perturbation theory and a square well potential. J. Chem. Phys., 64: 638-640. Rabinovich, V.A., Vasserman, A.A., Nedostup, V.I. and Veksler, L.S., 1988. Thermophysical properties of neon, argon, krypton and xenon. Hemisphere Publishing, Berlin. Rowlinson, J.S. and Swinton, F.L., 1982. Liquids and liquid mixtures. Butterworths, London. Sandler, S.I., 1985. The general&d van der Waals partition function. I. Basic theory. Fluid Phase Equilibria, 19: 233-257. Sowers, G.M. and Sandler, S.I., 1991. Equations of state from generalized perturbation theory. I. The hard-core Lennard-Jones fluid. Fluid Phase Equilibria, 63: l-25. Verlet, L. and Weis, J., 1972. Equilibrium theory of simple liquids. Phys. Rev. A, 5: 939-952. Weeks, J.D., Chandler, D. and Anderson, H.C., 1971. Role of repulsive forces in determining the equilibrium structure of simple liquids. J. Chem. Phys., 54: 5237-5247. Zwanzig, R.W., 1954. High-temperature equation of state by a perturbation method. I. Nonpolar gases. J. Chem. Phys., 22: 1420-1426.
APPENDIX
We have used the Metropolis Monte Carlo algorithm (Allen and Tildesley, 1987) to simulate a system of 112 atoms interacting with the LermardJones potential. The simulation runs were started from a random initial configuration, averages were recorded after discarding 100 000 configurations, and the lengths of the simulation runs were 3 X lo6 configurations for pa3 G 0.3 and 1.5 X lo6 configurations for pa3 > 0.3. The compressibilities were computed from the relation N-l
N
C
C
r=l
y=r+1
rij
d”(rij)
3NkT + Z,nc
dr ‘J
where the angle brackets denote an ensemble average, and rij is the distance between atom i and atom j. The cut-off distance in the simulation is half of the box length, and long-range corrections Z,, are computed in the usual way. The compressibilities at the isotherms simulated are shown in Table Al. The compressibilities at a reduced density pa3 = 0.05 at all temperatures
150 TABLE Al Compressibilities of the LJ fluid kT/c
PO3
0.05 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
1.35
1.75
2.0
3.0
4.0
6.0
0.846 0.719 0.508 0.367 0.273 0.295 0.543 1.200 2.391 4.408
0.909 0.841 0.733 0.680 0.709 0.860 1.210 1.886 3.042 4.928
0.934 0.889 0.819 0.813 0.880 1.061 1.456 2.156 3.291 4.991
0.988 0.992 1.020 1.108 1.281 1.560 1.982 2.687 3.663 5.115
1.012 1.036 1.119 1.248 1.449 1.762 2.196 2.850 3.782 5.024
1.034 1.082 1.204 1.360 1.601 1.923 2.364 2.949 3.731 4.792
were computed from a virial expansion truncated after the first order term in density. The configurational internal energies were computed from N-l ECONF=
c r=l
N c
U(ri,)+ELRc
j=i+l
i and are shown in Table A2. The error in the simulation data is approximately 1% for the compressibility and 0.2% for the configurational internal energy.
TABLE A2 Configurational internal energies of the LJ fluid PO3
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
kT/r 1.35
1.75
2.0
3.0
4.0
6.0
-
-
-
-
-
- 0.475 - 0.934 - 1.384 - 1.781 - 2.126 - 2.377 - 2.508 -2.444 - 2.086
0.771 1.468 2.114 2.735 3.376 4.024 4.662 5.243 5.687
0.696 1.356 1.979 2.609 3.222 3.846 4.436 4.941 5.264
0.660 1.295 1.910 2.528 3.137 3.738 4.295 4.749 5.047
0.592 1.172 1.740 2.302 2.850 3.366 3.790 4.102 4.556
0.546 1.079 1.605 2.115 2.596 3.025 3.342 3.492 3.443