Macromolecules 1993,26, 2655-2662
2666
Microscopic Equations-of-State for Hydrocarbon Fluids: Effect of Attractions and Comparison with Polyethylene Experiments John G.Curro,'J Arun YethirajJ Kenneth S. SchweizerJ John D.McCoy,# and Kevin G. Honnellll Sandia National Laboratories, Albuquerque, New Mexico 87185,Department of Materials Science & Engineering and Materials Research Laboratory, University of Illinois, Urbana, Illinois 61801,Materials & Metallurgical Engineering Department, New Mexico Institute of Mining & Technology, Socorro, New Mexico 87801,and Phillips Petroleum Company, Bartlesville, Oklahoma 74006 Received October 26,1992; Revised Manuscript Received February 2,1993 ABSTRACT Microscopic equations-of-state are developed for n-alkanes and polyethylene based on the polymer reference interaction sitemodel (PRISM)integral equation theory and a generalizedFlory approach. The molecules are modeled as a series of overlapping spheres (methylene groups) with constant bond length and bond angles; internal rotations are accounted for by the rotational isomeric state approximation. The interaction between siteson differentmoleculesis taken to be of the Lennard-Jonesform. The thermodynamic properties of the fluid are obtained via standard perturbation theory in which the potential is divided into a repulsive reference system and an attractive perturbation. The reference system is approximated by a hard-core repulsion in which the hard-sphere diameter d ( T ) is estimated for polyethylene from wide-angle X-ray scattering experiments. The PRISM theory is used to calculate the hard-sphere chain contribution to the equation-of-stateby three different thermodynamic routes: (1)integrating the compressibility, (2) evaluating the density profile at a hard wall, and (3) using a hard-sphere 'charging" method analogous to the virial approach in monatomic liquids. The generalized Flory dimer (GFD) theory is used to obtain a fourth equation-of-statefor the hard-sphere chains. The attractive perturbation is treated with first-order perturbation theory, making use of the radial distributionfunctiongo(r)of the reference system. The various equations-of-statepresented differ in the route to the hard-chain pressure;PRISM is used in all cases to treat the attractions. Excellent agreement for the equation-of-stateis found between the hybrid GFD/PRISM calculationsand moleculardynamicssimulationsof n-butane and experimental pressure-volume-temperature (PVT) measurements on polyethylene melts. The compressibility and charging routes predict pressures which are too low and too high, respectively, for polyethylene. The wall route yields pressures in good agreement at experimental densities but predicts a melt which is too compressible.
I. Introduction A microscopic description of the equation-of-state of polymer liquids is a well-known problem which has been investigated by many workers.'* Early polymer equationof-state theories have necessarily consisted of lattice, cell, or hole model descriptions. Although approaches such as those due to Simha and SomcynskyZ or Sanchez and Lacombe3work quite well when compared to experimental data, they are based on qualitatively useful but quantitatively undefined parameters such as free volume and number of external degrees of freedom. It is not obvious how these parameters are related to the molecular architecture of the actual polymer chains. More recently, equation-of-statetheories for tangent hard-sphere chainsp7 and square-wellchains8have been developed in continuous space. The idealized nature of these models does not allow them to reliably predict actual experimental data, however, computer simulations can be used to test their validity. In the present investigation, we will focus on calculating the equation-of-state of polyethylene and alkane liquids in continuous space. In order to compare with experimental pressure-volume-temperature (PVT) data, we model the chains with a realistic backbone structure (rotationalisomeric state model) interacting with realistic intermolecular potentials (Lennard-Jones). Recently, integral equation techniques have been developed to calculate the intermolecular structure of polymer melhg-13 These methods are a generalization to +
Sandia National Laboratories.
* University of Illinois.
New Mexico Institute of Mining & Technology.
11 Phillips Petroleum Co.
polymers of the reference interaction site model (RISM theory) of Chandler and Andersen.IP16 This new theory, which we call polymer RISM or PRISM theory, has also been expanded to allow both the intramolecular and intermolecular structure of homopolymer liquids to be determined in a self-consistentmanner.17 Extensionshave also been made of the PRISM theory to binary blend~.I8-~~ Good agreement of the PRISM theory with molecular d y n a m i c ~ ~and ~ J Monte ~ 8 ~ ~Carlo26simulations was found for the radial distribution function of hard-core chain liquids. Comparisons have also been made between the structure factor for polyethylene melts obtained from the PRISM theory and X-ray scattering experiments of Narten and H a b e n s c h u ~ s Remarkably . ~ ~ ~ ~ ~ good agreement was found by modeling the polyethylene as overlapping hardsphere chains. On the basis of these comparisons between simulation and experiment, it seems that the PRISM theory adequately describes the intermolecular structure of homopolymer liquids. The level of accuracy for polymers appears to be about the same as corresponding RISM calculations for small, rigid molecules as far as intermolecular structure or packing is concerned. A natural question to ask is, can the PRISM approach also be employed to calculate thermodynamic properties of polymer melts? In earlier work7 Schweizer and Curro used PRISM theory to compute the equation-of-statefor freelyjointed chains using a virial route to the pressure in which a superposition approximation was made for the threebody correlation function. The results were inconclusive, however, since the superposition approximation was later shown to be i n a c c ~ r a t e . ~ ~ ~ ~ ~
0024-929719312226-2655$04.00/0 0 1993 American Chemical Society
2666 Curro et al.
In a recent publication31 (designated as paper 1) we reexamined this question by computing the equation-ofstate of alkane chain liquids interacting with hard-core repulsions between intermolecularsitea (methylenegroups). In this study we calculated the equation-of-state from PRISM theory by three different routes and also from the generalizedFlory dimel.4~~ (GFD)theory,suitably extended to overlapping hard-sphere chains. Unfortunately, since the PRISM theory is an approximate theory for intermolecular structure and the pressure is known to be an extremely sensitivefunction ofg(r),significantly different results for the hard-chain pressure are found from the three routes. In the present investigation we continue this equationof-state theme by considering polymer chains interacting with a Lennard-Jones potential. Our strategy will be to use our previous hard-sphere chain model as a reference system and then to include the effect of attractive interactions as a perturbation.32 The hard-core diameter of the reference system is then chosen in an optimum fashion to match the continuous soft-core repulsion of the Lennard-Jonespotmtial. These calculationscan thus be compared with PVT experiments on polyethylene melts. It should be mentioned that recently Yethiraj and Hallsb calculated the equation-of-stateof square-wellchain fluids of n = 4,8, and 16units by using a GFD/PRISM approach. They found good agreement with computer simulations at meltlike densities. An obvious difficulty in making comparisons to experimental data, in contrast to simulations, is the ambiguity that exists regardingthe interatomicpotentials. Assuming a Lennard-Jonesform, a clean comparison with experiment is obscured by two adjustable parameters. In the present investigation, we eliminate one of these parameters by choosing the hard-core diameter to give agreement with the first peak in the experimental structure factor of polyethylene. We have some feeling for the magnitude of the remaining adjustable parameter, the depth of the attractivewell, based on known values for short alkanes.48 We will begin by briefly reviewing the PRISM theory and its application to computing the hard-core polymer equation-of-~tate.3~ Our strategy for using these hardcore results, together with perturbation theory, to model Lennard-Jonesinteractions between methylene groups is also presented in section 11. Numericalresults for butane are presented in section 111 and compared with the simulations of Ryckaert and B e l l e m a ~ w . ~Calculations ~!~~ for polyethylene by the various routes are also presented and critically compared with experimental PVTdata. The paper concludes with a brief discussion in section IV. 11. Theory The PRISM theoryg-13 is an extension to flexible polymers of the reference interaction site model of Chandler and AndersenlP16 and has been discussed in detail previously.lOJ1 In this approach, each polymer is taken to consist of a chain of spherically symmetric sites; intermolecular interactions occur only between these sites. The intersite direct correlation function C(r)is related to the intersite radial distribution function g(r) through a generalized Ornstein-Zernike equation. For long chains we can neglect end effecta; consequently the correlations between any pair of intermolecular sites are taken to be the sameregardless of their position along the chain. After the appropriate preaveraging over the intramolecular structure of the polymer molecule, the generalized Orn-
Macromolecules, Vol. 26, No. 11, 1993
stein-Zernike equation becomes
h(k) = &k) C(k) &k) + p,&k) &k) h(k) = &k) c(k)&k)
(1) where the carets denote Fourier transformation with respect to wavevector k, h(r) is the total correlation function defined as h(r) = g(r) - 1,and pm is the monomer or site density. All the information about the intramolecular structure of the polymer molecules enters the theory through the function o(r)
where the sum is extended over all intramolecular pairs of N repeat units per chain, and wij(r) is the probability density that sites i and j are a-distance r apart. &k) is the structure factor defined as w(k) + pmh(k).In general, the intramolecular probability function w(r) and the intermolecular correlations Cfr) andg(r) would have to be determined self-consistently. Although we have recently developed methods for dealing with this difficult selfconsistentcal~ulation,~~ for simplicity this self-consistency calculation will be avoided in the present paper. Instead, we make use of the Flory ideality hypothesis35 which appears to be an excellent approximation in dense homopolymer melts based on SANS measurements and computer simulations. With this approximation, we can then calculate w(r)from a separate single-chaincalculation in which long-range, excluded-volume and many-body, medium-induced interactions are omitted. The Flory ideality hypothesis is expected to improve17in the case of real polymer chains such as polyethylenewhich, in contrast to freelyjointed chains, have some stiffnesswhich prevents local overlap of intramolecular chain segments. Introduction of the o(r)function appropriate to the particular polymer structure intothe generalizedOmsteinZernike equation above leads to a relation between intermolecular correlation functions C(r) and g(r). An additional closure relation is needed in order to uniquely solve for these functions. For hard-core interactions, we have the exact expression h(r) = -1, r < d (3a) resulting from the impenetrability of hard spheres of diameter d. On the basis of an analogy with the PercusYevick (PY) theory for hard Chandler and Andersen suggested1”16 the atomic-like closure
C(r)z 0, r > d (3b) taken to be approximatelyvalid outside the hard core. We have found this PY type closure to be accurate for the intermolecular structure of hard-sphere polymer chains by comparison with MD17924725 and Monte Carlo%simulations, as well as, with experimental X-ray scattering 1-3 experiments on polyethylene l i q ~ i d s . ~ ~Equations l~8 then form a nonlinear integral equation for the intermolecularg(r) of a homopolymer melt. Thisintegral equation can be solved numerically by either a variational method35 or an iterative approach. Complicationsarise if the interaction potential contains an attractive branch. In the case of polymer blends, use of the atomic-like MSA closure for attractions leads to a of the nonclassical molecular weight critical temperature at variance with recent simulations3’j and SANS e~periments.~’Recent theoretical studies% suggest that atomic-like closures are not appropriate for attractive interactions in polymers and new molecular-
Macromolecules, Vol. 26, No. 11, 1993 12
I"
0
,
"
"
Microscopic Equations-of-State for Hydrocarbon Fluids 2657 3
"
0.5
1
1.5
2 kl
2.5
3
3.5
Figure 1. Intramolecular structure factor ;(k) for n-butane
(dashed curve) (d = 3.77 A,T = 300 K)and polyethylene (solid curve) (d = 3.90A,T = 430 K,N = 6429)plotted in Kratky form.
based closures have been proposed. In the present work there is no difficulty in using the hard-core closure (eq 3b), since only the repulsive portion of the intermolecular potential is analyzed using integral equation theory; the contribution to the pressure due to the attractive tail of the potential will be accounted for by perturbation theory. In paper 1 we modeled the methylene groups along the chain backbone of polyethylene as overlapping hard spheres of diameter d. The backbone structure of polyethylene was calculated using a rotational isomeric state approach39with bond length I, bond angle 8, and torsional Totational angles 4. The intramolecular structure factor o(k),the Fourier transform of w(r) defined in eq 2, cannot be computed analytically for the rotational isomeric state model. Inpur previous studies on polyethylene27j28J1we estimated w(k)by treating the terms in eq 2 corresponding to the local structure exactly and approximatingthe longrange contribution^.^^ The short-rangestructureis treated explicity by enumerating the rotational isomeric states for distances li - jl I5 along the bactbone. Thus the pentane effect is included exactly. The oi,(k) corresponding to distances li - jl 1 6 are approximately calculated from a semiflexible chain model, subject to the condition that the second and fourth momenta of r , agree with the rotational isomeric state model. In the rotational isomeric state calculation the bond length was taken as 1.54A, the bond angle as 112O, and gauche states of i12Oo relative to trans were used for polyethylene. A trans/gauche energy of 500 cal/mol and an additional gauche+/gauche- en9rgy of 2000 cal/mol was used. In the butane calculation, w(k) was evaluated exactly for the rotational potential, bond length, and bond angles employed by Ryckaert and B e l l e m a n ~ ~in3 >their ~ ~ simulations. The results for the continuouspotential of butane were indistinguishablefrom arotational isomeric state calculation with a trans/gauche energy of 700 cal/mol. A plot of the intramolecular structure factor is shown in Figure 1for both polyethylene and butane. Note that, in polyethylene, there appears a well-defined intermeciate scaling regime in the range 0.1 < kl < 0.4 for which wi,(k) l/(k1)2 expected for long chains. For butane no such scaling regime is observed. For the calculation of the equation-of-state, we expect all wavevector contributions to w(k) to be important. The intermolecular radial distribution function has been obtained in paper 1by solving eqs 1-3 using the intramolecular structure factors for polyethylene and butane. Representativeg(r)'sare shown in Figure 2 for polyethylene melts at three different densities. Having solved for the intermolecular structure, it is in principle possible to then obtain the thermodynamic properties of the melt. Because the PRISM theory for the structure is only approximate, the calculated ther-
-
0 4 " 0
4
"
4
"
"
"
8
"
12
"
"
'
16
20
(1,
Figure 2. Intermolecular radial distribution function g(r) computed for polyethylene (d = 3.90A, T = 430 K,N = 6429) using PRISM theory. The solid curve is for density pmd3= 2.2, the short dashed curve for pmd3= 2.0, and the dash-dot curve for pmd3 = 1.8.
modynamic functions depend on the particular statistical mechanical route used. In paper 1 three different routes were employed to extract the hard-core contribution to the pressure from the PRISM results the compressibility, charging, and wall methods. In the compressibility route, one simply integrates the i s o t h e r d compressibility KT from zero density to the density of intarest. The campressibility is related to the structure factor at zero wave vector according to PmkBTKT = &O). The pressure can then be seen to be (4)
where p = l / k ~ T . The analog of the virial route for polyatomicliquidswas put forth by ChandleP and consists of gradually turning on the hard-sphere diameter, starting from a point, as a charging parameter X varies from 0 to 1. This charging procedure gives the Helmholtz free energy for a homopolymer fluid as (5)
where A and A0 are the free energies of the hard-core fluid and an ideal gas of chains, respectively, and gWuf+) is the contact value for g(r) for a liquid comprised of interaction sites of hard-core diameter Ad. The pressure is then obtained from eq 5 by numerical differentiation with respect to volume. In the monomeric fluid limit, eq 5 lea& to the familiar virial formula for the pressure. Yet a third route to the pressure was utilized recently by Dickman and Halla for extractingthe pressure from Monte Carlo simulations and used by Yethiraj and Hall41 in conjunction with PRISM theory. In this wall approach the pressure is simply related to the density of chain sites at a hard wall P / P , = g,(O) (6) whereg,(z) is the totalsite density profile in the direction perpendicular to the wall, the wall is impenetrable to the centers of the chain sites. Using the g,(z) predicted by PRISM theory, eq 6 was found to be accurate at high densities41but breaks down at low densities. In paper 1 we extended the generalized Flory dimer (GFD) theory415of Hall, Dickman, and Honnellfor tangent hard-sphere chains to the case of overlappinghard spheres, a model appropriate for alkanes, polyethylene, and other chain molecules. In this approach the pressure of a liquid composed of hard chains is approximately mapped onto a fluid of hard dimers andmonomers. The compressibility factor Z(N,$ b P / P m of an N-mer fluid of monomer density pm or packing fraction q is approximatelyrelated to the corresponding monomer and dimer compressibility
Macromolecules, Vol. 26,No.11, 1993
2658 Curro et al.
factors: NZ(N,q)r 2(1+ YN)Z(2,4- YNz(l,d (7a) where YN is related to Ue(N),the volume excluded by an N-mer to a monomer, averaged over all conformations of the N-mer, i.e.
u,(N) - u,(2) (7b) ue(2) - u,(l) The Carnahan-Starling42 and Tilde~ley-Streett~~ equations are used for Z(1,q) and 2(2,q), respectively. In paper 1we showed that the identification of the monomer and dimer "reference" fluids whose compressibility factors appear in eq 7a is not unique. One way of obtaining the reference monomer and dimer fluids is to merely relax all the bonding constraints in the polymer fluid (GFD-A); since the packing fraction is kept constant, this results in a decrease in the site density of the monomer and dimer fluids from that of the polymer fluid. The second way of obtaining the referenoe monomer and dimer fluids is to require that both the packing fraction and the site density of the monomer and dimer fluids remain the same as in the polymer fluid (GFD-B); this results in the monomer and dimer molecules having different hard-core diameters than the polymer molecules. In paper 1it was concluded that the GFD-B method was more appropriate (and more accurate) than the GFD-A and will therefore be employed in the present investigation. In order to describe polyethylene, we represent the CH2 groups as single sites interacting with a united atom potential of the Lemard-Jones form YN
=
u(r) = 4t[(u/r)12- (u/rT] (8) where u and Q are the usual Lennard-Jonesparameters. In it is customary standard perturbation theory of to divide the potential into a reference part udr) and a perturbative part u&). It is assumed that the intermolecular radial distribution function go@)can be calculated for the reference system with potential uo(r). A free energy expansion is then performed about the reference system, assuming that the perturbative potential u,(r) is small. To first order the Helmholtz free energy per unit volume is given by A = A,
+ z1p 2 J g 0 ( r ) u,(r) di + ...
(9)
Two well-known perturbation approaches that have been applied to molecular fluids are the BarkeI-Henderson (BH) the0ry4~9~~ and the Weeks-Chandler-Andersen (WCA) These two approaches differ in their choice of u&) and u&). For the BH theory, the reference and perturbative potentials are taken to be the repulsive and attractive branches of the Lennard-Jones potential, respectively. u,(r) = u(r); r I Q
=O;
r>a
u,(r) = 0; r I u = v(r); r > u The WCA theory uses the separation
(10)
for which the perturbation varies more slowly than the BH perturbation over the range of r near the maximum in go@) for monatomic fluids. Because of the correlation hole effect12in polymeric fluids, the first peak is significantly smaller (see Figure 2 for polyethylene) and, dependingon the density and chain length, may disappear altogether. As a consequence,the significant improvement in WCA over BH perturbation theory demonstrated in monatomic fluids may not carry over to polymer liquids. Ineqe lOand 11it canbeseen thatthereferencepotential udr) is a continuous, soft-core repulsion. It is customary to use a hard-core potential for the reference syatem and then to choose the hard-core diameter d in an optimum fashion to mimic the continuous repulsion. For the BH theory, this optimum hard core depends only on temperature. d = Jom{l - exp[-Ou,(r)l) dr
(12)
In the WCA theory d is calculated using a "blip" function technique. For the case of the Percus-Yevick, hard-sphere closure in eq 3b, d is chosen to satisfy
hdrzexp[-Ouo(r)lCo(~;d) dr + Jdmr2(1 - exp[-Ou,(r)l)g, (r;d)dr = 0 (13) where COand go are the direct correlation and radial distribution functions for the hard-core system with hardcore diameter d. Since both of these correlation functions are state-dependent functions, d in the WCA theory depends on both temperature and density, in contrast to the BH theory where d is only temperature dependent. As a result, the BH theory is much easier to implement than the WCA theory. In practice, we have found the difference in the d's predicted by eqs 12 and 13 to be quite small (-0.5 %) for polyethylene melts. It should be mentioned that eq 12can be to be a first-order approximation to the WCA expression in eq 13. 111. Results and Discussion
A. Butane. The first of our calculations that we will discuss is on n-butane. Ryckaert and Belleman~~a3~ performed MD simulations for three different states of n-butane near the coexistence curve. In these simulations the temperature and molar volume were fixed and the compressibility factor @l'/pm was calculated. For butane the monomer density Pm is 4 times the molecular density. In the simulations butane is modeled as a Lennard-Jones fluid with the u and t parameters chosen to match experimental data. Since the pressures chosen by Ryckaert and Bellemans to study in their simulations are near zero and the slope of the pressure volume curves is very large, a direct comparison between the simulated and calculated compressibility factors is not very meaningful. A more useful comparison, which we attempt to make here, is on the basis of the molar volume. In Table I we show the calculated molar volumes for butane corresponding to the compressibility factors found in the simulations of Ryckaert and Bellemans. These calculations were performed using the three routes to the pressure: wall, compressibility (Comp.), and generalized
Macromolecules, Vol. 26, No. 11, 1993
Microscopic Equations-of-State for Hydrocarbon Fluids 2659
Table 1. Comparison of Molar Volumes Obtained by Various Routes (Together with Barker-Henderson Perturbation Theory for the Attractive Interactions) with
4.6 4.6
'
'
-
BH Theory
l
the Simulations of n-Butane by Ryckaert and B e l l e m a n ~ ~ ~ J ~
P (cmYmol) T(K) (BP/~,)sim. RBsim. 199.9O 0.0855 86.7 291.5b 99.7 0.035 274.0' -0.138 99.5 IJ
C ~ K 83.4 147 109
w . all ComD. GFD-B 62.4 64.4 84.7 77.4 70.7 97.9 69.1 67.2 91.6
0
3.0
elk = 72 K, IJ = 3.923 A, d = 3.717 A, 1 = 1.53 A. elk = 72 K, = 3.923 A, d = 3.668 A, 1 = 1.53 A. elk = 84 K, IJ = 3.920 A, d =
3.6
'"""
30
0
40
3.694 A, 1 = 1.53 A.
Flory dimer (GFD-B). The hard-sphere diameter d of the reference system was optimized using BH theory according to eq 12 using the Lennard-Jones parameters employed in the Ryckaertand Bellemans simulations. The attractive contribution to the pressure was determined from BH perturbation theory using eqs 8-10. The radial distribution function go@)of the reference system in eqs 9 and 12 was computed from PRISM theory for the optimized hard-core system. The different theoretical predictions in Table I (columns 4-7) differ only in the route employed to obtain the hard-chain pressure. The calculations of the charging, compressibility, and wall routes were carried out using PRISM theory for the hard-sphere contribution to the pressure, as described above and in paper 1. The calculations for the GFD-B route are described in detail in paper 1. The GFD-B route is actually a hybrid method that uses GFD computationsfor the hard-sphere reference system and PRISM to treat the attractions. It can be seen from Table I that at T = 291.5apd 274.0 K the predicted molar volumes followthe pattern V(Chrg.) > V(GFD) > P ( W d ) > P(Comp.), the same trend seen for the hard-sphere chains at T = 300 K shown in Figure 9 of paper 1 at fixed pressure. At the lower temperature of 199.9 K, however, the ordering of the molar volumes becomes inverted: V(GFD) > V(Chrg.) > V(Comp.) > V(Wal1). This same ordering is also seen in Figure 9 of paper 1at very high densities or packing fractions. This inversion in the ordering of the purely hard-sphere chains shifts to lower densities as the temperature is lowered. For the hard-sphere chain calculation,temperature enters the problem indirectly through the optimum hard-core diameter d ( T ) and the intramolecular structure function
"
"
"
"
"
"
50 dkg
"
'
60
70
''
60
(OK)
Figure 3. Lennard-Jones parameters: u versus elks. The solid curve is for polyethylene ( N = 6429) calculated at fixed hardcore diameter d = 3.90 A and temperature T = 430 K,according to Barker-Henderson theory using eq 12 in the text. The filled circles correspondt o parametersfound by Lopez-Rodriguez and co-worker~.~~ The filled square represents the value found by Ryckaert and B e l l e m a n ~ ~for e ~butane. ~ The X corresponds to the parameters used in Figure 7. 5 4 E
Q
g
3 2 1
0 1
14
18
22
26
3
P$'
Figure 4. Hard-core equation-of-statefor polyethylene (5" = 430 K,N = 6429) computed by various routes.
Lennard-Jones parameters of a series of n-alkanes by comparing theoretical and experimental second virial coefficients. The theoretical second virial coefficientswere evaluated using a two-chain Monte Carlo procedure. In this comparison, the authors allowed the well depths of methyl and methylene groups to vary independently. Interestingly, they found the c parameter reflecting attractive interactions between CH2 groups to be significantly smaller than corresponding interactions between CH3moieties, suggesting that the Ryckaertand Bellemans potential would not be expectedto apply to polyethylene. &I. These considerations illustrate the difficulties of unamThe GFD-B calculations agree very well with the biguously comparing theoretical and experimental equasimulation molar volumes (Wacceptable Lennard-Jonesparameters u and E through the chose their Lennard-Jones parameters to match experiBarker-Henderson relation (eq 12) or the WCA theor (eq 13). A plot of Q versus t corresponding to d = 3.90 mental data at three states, it is tempting to apply these and satisfying eq 12 is shown in Figure 3 for polyethylene same values toother alkane chains at other thermodynamic at 430 K, By choosing d from scattering experiments in conditions. Unfortunately, the Lennard-Jones potential is known to be inadequate for alkanes and p~lyethylene.~~ this manner, we have effectively reduced the number of adjustable parameters from two to one. As a result, the apparent Lennard-Jones parameters Figure 4 depicta the equation-of-state of overlapping depend on chain length and thermodynamic state. Recently, Lopez-Rodriguez and co-worker~~~ obtained the hard-sphere chains for polyethylene at 430 K (d = 3.90 A)
x
Macromolecules, Vol. 26, No. 11, 1993
2660 Curro et al. ....,
0.6 I
I
I
I
I
.
.,.
.
.
I
.
,
r/k,=40DK BH
-
0.5
I
I
. .
" ' ' I . ' , ' . ' ' I I I I
i
I
04
I I 0
"
04
/r
,
.
.
.
I
.
"
05
06
.
.I .
I ! ,k/kB=72'1p
I
.
07
,
.
,....,....,....,', 08
09
1
11
r
I
' I . , ,
7
I
1 12
06
07
Figure 5. Equation-of-state computed for polyethylene (d = 3.90 A, T = 430 K, N = 6429) using PRISM theory with BarkerHenderson theory. The solid curves correspond to the charging route. The dashed lines correspond to the compressibilityroute. The solid circles are the experimental PVT data of Olabisi and Simha.50 obtained in paper 1using the methods discussed above. Note the wide variation in the results, depending on the route used to extract the pressure. The charging and compressibility routes bracket the GFD and wall results, in contrast to butane (see Figure 9 in paper l),where the GFD/charging and wall/compressibility results cross at high densities. The hard-sphere chain equations-of-stateby the various routes form the basis of the reference system for the BH or WCA perturbation theory in order to treat the full Lennard-Jones potential. Because the BH approach is much easier to implement (since d depends only on temperature) than the corresponding WCA method and because the BH theory has been shown4 to work well for molecular fluids, most of the present calculations utilize BH theory. Figure 5 displays these equation-of-state calculations at various values of the attractive well depth for the charging and compreesibility routes. The points are from experimental PVT measurementss0of Olabisi and Simha in the melt phase (p, = [23.288/d31p,d3, where d is in angstroms). The charging route, shown by the solid curves for t / k g = 40 and T = 72 K, clearly gives values of the compressibility factor in excess of the experiment at the experimental densities. Agreement with experiment would improve if r / k g was increased beyond the Ryckaert and Bellemans value of 72 K. Such an increase,however, would not seem to be justified in view of the investigation of Lopez-Rodriguez and co-workers,48suggesting that the interaction between methylene groups is much lees than that between methyl groups. The compressibility route, shown by the dashed linea in Figure 5 for the same well depths, gives compressibility factors which are much smaller (and negative) for the experimental densities. Lowering the well depth much below 40K is not consistent with the Lopez-Rodriguez study either, as can be seen in Figure 3. Furthermore, if €/kgis lowered to attempt to fit the experimental data, the slope of the theoretical curve will be too small. Thus, we conclude that neither the charging nor compressibility routes give satisfactory agreement with experimental PVT data on polyethylene melts at 430 K. Figure 6 shows the calculation of the equation-of-state from the wall route. It can be seen that, for € / k B = 40 K, reasonable agreementis found with experiment, although the slope of the theoretical curve is somewhat low. The corresponding calculation using WCA theory for the
09
08
P, (alcm3)
1
11
12
P, (glcm3)
Figure 6. Equation-of-state computed for polyethylene (d = 3.90 A, T = 430 K, N = 6429) using PRISM theory via the wall
route. The solid (dashed) curve uses BH perturbation theory with elk = 40 K (72 K). For comparison, the dotted curve was calculated from WCA perturbation theory at elk = 40 K. The solid circles are the experimentalPVTdataof Olabisi and Simha.50 06
0.4 E
0'5
e t
%
1
I""
" "
"
'
'"1
1 I
I
i
0.2
tt 06
'
I
~Jk,=38.7"K
0.3
01
'
1 I"
i 07
1
I
I
08
09
1
11
12
P, W m 3 )
Figure 7. Equation-of-state computed for polyethylene (d = 3.90 A, T = 430 K, N = 6429) using PRISM theory via the hybrid
GFD-B approach discussed in the text. Barker-Henderson perturbationtheory ia used. The solid circles are the experimental PVT data of Olabisi and Simha."
attractive interactions is also shown for comparison. Figure 7 depicts the hybrid GFD-Bmethod for d k B = 38.7 K. In this case it can be seen that the agreement with experiment at T = 430 K is almost quantitative. Furthermore, this combination of u and e indicated by the X in Figure 3 seems completely consistent with the LopezRodriguez investigation. Similar conclusions can be reached by examining the theoretical calculations at T = 500 K shown in Figure 8. These calculations were performed assuming the hardcore diameter d at 500 K remains the same (3.90 A) as at 430 K. The well depth of elkg = 38.7 K that fits the experiments quantitatively at 430 K leads to density predictions which are slightly high via the GFD route at 500 K. This indicates that the theoretical thermal expansion coefficient is lower than that seen in the Olabisi and Simha data. If we assume that the V-T curve is linear over the temperature range of 430-500 K, then we can estimate the thermal expansion coefficients obtained for the various routes to the equation-of-state. These estimates are shown in Table 11. Note that the wall method leads to thermal expansions which are too high for polyethylene; the GFD-B route, on the other hand, is slightly lower than the Olabisi-Simha data. Unfortunately, X-ray scattering experiments have not been performed at 500 K. On the basis of previous e ~ p e r i e n c e , ~ ~
Macromolecules, Vol. 26, No. 11, 1993
Microscopic Equations-of-State for Hydrocarbon Fluids 2661
1
0.8
2
0.6
0.4
02
/'
0
,
06
.
"
I
'
07
.;/ '
I
I
/ , , , , , , , , , . ~ , , , , , '
'
08
'
'
'
09
1
11
12
P, (g/cm3)
Figure 8. Equation-of-state computed for polyethylene (d = 3.90 A, T = 500 K, N = 6429)using PRISM theory and BarkerHenderson perturbation theory. The solid circles are the experimental PVT data of Olabisi and Simha.50
Table 11. Thermal Expansion Coefficients Estimated for the Various Routes to the Equation-of-State for Polyethylene Melts T (K) 430-500 430-500 430-500
route wall GFD-B exptl (Olabisi & Simhaso)
a (x1W K-l)
19.0 6.4 8.5
we expect that the hard-core diameter d ( T ) which matches the fiist peak in the structure factor would increase slightly with temperature. Such a shift in d would increase the thermal expansion coefficient and improve the agreement with the experiments of Olabisi and Simha at 500 K. It should be pointed out that, from the BH or WCA theories, one would expect the hard-core diameter to decrease with increasing temperature, in contradiction to what we anticipate from X-ray scattering experiments. A possible explanation51for this reverse trend is that we represent the CH2 groups with a united atom potential. One might anticipate that better interlockingof the hydrogens would occur as the temperature is lowered, thereby lowering the effective hard core.
IV. Conclusions The equation-of-stateof n-butane and polyethylenehas been calculated from microscopic theory by four independent routes. The charging route appears to give pressureswhich are too high, and the compressibility route leads to pressures which are too low. This is not surprising in view of the findings from the purely hard-sphere chain calculationsstudied in paper 1which shows similar trends. For polyethylene, the wall route gives much better agreementwith experiment,however,the melt is calculated to be too compressiblerelative to experiment. An inherent difficulty in using PRISM theory (or any integral equation theory) for the reference system equation-of-state is traceable to its approximate nature. Since integral equation theories of atomic and molecular liquids are not exact, they are not thermodynamically consistent. Thus differences in the equation-of-state result from the different methods for extracting the pressure. As seen in paper 1 these differences are large in the case of polyethylene hard-sphere chains. The hybrid GFD-B method, which uses the hard-chain GFD equation-of-state for the reference system and the radial distribution function from PRISM theory, gives remarkable agreement with both MD simulationsof butane and PVT experiments on polyethylene melts. Such an
approach not only takes advantage of the GFD theory which has been demonstrated to work well for hard-chain systems but also employs the PRISM theory in a mode which is known to work well for calculating the structure of the polymer liquid. Of course,the main difficultyin making an unambiguous comparison between theory and experiment for the equation-of-state rests on our uncertain knowledge of the correct interparticle potential. Although the LennardJones parameters that lead to agreement between hybrid GFD theory and experiment seem reasonable, a definitive comparison is difficult. Such uncertainties would disappear if realistic simulations of long-chain alkane systems at high density were available.
Acknowledgment. Work at Sandia National Laboratories was supported by the U.S. Department of Energy under Contract DE-AC047DP00789. Work at the University of Illinois was sponsored by the Division of Materials Sciences, Office of Basic Energy Sciences, US. Department of Energy. The authors thank N. J. Curro, Rice University, for performing some of the numerical equation-of-state calculations. References and Notes (1) Curro, J. G. J. Macromol. Sci., Reu. Macromol. Chem. 1974, C l l , 321. (2) Simha, R.; Somcynsky, T. Macromolecules 1969,2, 342. (3) Sanchez, I. C.; Lacombe, R. H. Macromolecules 1978,11,1145. (4) Dickman, R.; Hall, C. K. J. Chem. Phys. 1986,85,4108. (5) Honnell, K. G.; Hall, C. K. J. Chem. Phys. 1989,90,1841. (6) Wertheim, M. S. J. Chem. Phys. 1987,87,7323. Chapman, W. G.; Jackson, G.; Gubbins, K. E. Mol. Phys. 1988,65,1057. Chiew, P. C. Mol. Phys. 1990, 70, 129. Boublik, T.; Vega, C.; DiazPena, M. J. Chem. Phys. 1990,93,730. (7) Schweizer, K. S.; Curro, J. G. J. Chem. Phys. 1988,89,3342; J. Chem. Phys. 1988,89,3350. (8) (a) Yethiraj, A.; Hall, C. K. J. Chem. Phys. 1991,95,8494; (b) J. Chem. Phys. 1991,95, 1999. (9) Schweizer, K. S.; Curro, J. G. Phys. Reu. Lett. 1987, 58, 246. (10) Curro, J. G.; Schweizer, K. S. Macromolecules 1987,20,1928. (11) Curro, J. G.; Schweizer, K. S. J. Chem. Phys. 1987, 87, 1842. (12) Schweizer, K. S.; Curro, J. G. Macromolecules 1988,21,3070. (13) Schweizer, K. S.; Curro, J. G. Macromolecules 1988,21,3082. (14) Chandler, D.; Andersen, H. C. J. Chem. Phys. 1972,57,1930. (15) Chandler,D. InStudies inStatistical Mechanics VZZe Montroll,
E. W., Lebowitz, J. L., Eds.; North-Holland Amsterdam, The
Netherlands, 1982; p 274 and references cited therein. (16) Chandler, D. Chem. Phys. Lett. 1987, 140, 108. (17) Schweizer, K. S.; Honnell, K. G.; Curro, J. G. J. Chem. Phys. 1992,96, 3211. (18) Schweizer, K. S.; Curro, J. G. Phys. Rev. Lett. 1988, 60, 809. (19) Schweizer, K. S.; Curro, J. G. J. Chem. Phys. 1989,91, 5059. (20) Curro, J. G.; Schweizer, K. S. Macromolecules 1990,23, 1402. (21) Schweizer, K. S.; Curro, J. G. Chem. Phys. 1990, 149, 105. (22) Schweizer, K. S.; Curro, J. G. J. Chem. Phys. 1991,94, 3986. (23) Curro, J. G.; Schweizer. K. S. Macromolecules 1992.24.6737. (24) Curro, J. G.; Schweizer;K. S; Grest, G. S.; Kremer, K: J. Chem. Phys. 1989,89, 1357. (25) Honnell, K. G.; Curro, J. G.; Schweizer, K. S. Macromolecules 1990.23, 3496. (26) Yethiraj, A.; Schweizer, K. S. J. Chem. Phys. 1992, 97, 1455. Yethiraj, A.; Hall, C. K. J. Chem. Phys. 1992, 96, 797. (27) Honnell, K. G.; McCoy, J. D.; Curro, J. G.; Schweizer, K. S.; Narten, A.; Habenschuss, A. J. Chem. Phys. 1991,94,4659. G.;McCoy,J.D.;Curro, (28) Narten,A.;Habenschuss,A.;Honnell,K. J. G.;Schweizer, K. S. J. Chem. Soc., Faraday Trans. 1992,88, 1791. (29) Yethiraj, A.; Hall, C. K.; Honnell, K. G. J. Chem. Phys. 1990, 93, 4553. (30) Gao, J.; Weiner, J. H. J. Chem. Phys. 1989, 91, 3168. (31) Yethiraj, A,; Curro, J. G.; Schweizer, K. S.; McCoy, J. D. J. Chem. Phys. 1993,98,1635. (32) Hansen, J. P.; McDonald, I. R. Theory of Simple Liquids, 2nd ed.; Academic Press: London, 1986. (33) Rvckaert. J. P.: Bellemane. A. Chem. Phvs. Lett. 1976.30.123. (34) Ryckaert; J. Pi;Bellemans, A. Faraday" Discws., Chem.' Soc. 1978, 66, 95.
Macromolecules, Vol. 26, No. 11, 1993
2662 Curro et al. (35) Lowden, L. J.; Chandler, D. J. Chem. Phys. 1974, 61, 5228; 1973,59,6587; 1975,62,4246. (36) Deutch, H. P.; Binder, K. Europhys. Lett. 1992,17,697. (37) Gehlsen, M. D.; R o d a l e , J. H.; Bates, F. s.; wignall, D.; Hansen, L.; Almadal, K. Phys. Reu. Lett. 1992,68,2452. (38) Yethiraj, A.; Schweizer, K. S. J. Chem. Phys. 1992,97, 5927. J. D.; Honnell, K. G,; Curro, J. G.;Sehweimr, K. (39) MCCOY, Honeycutt, J. D. Macromolecules 1992,25,4905. (40) Dickman, R.; Hall, C. K. J. Chem. Phys. 1988,89,3168. (41) Yethiraj, A.; Hall, C. K. J. Chem. Phys. 1991,95,3749. (42) Carnahan, N. F.; Starling, K. E. J. Chem. Phys. 1969,51,635. (43) Tildesey, D. J.; Streett, W.B. Mol. Phys. 1980,41,341. (44) Barker, J. A,; Henderson, D. J. Chem. Phys. 1967,47, 4714; Annu. Rev. Phys. Chem. 1972,23,439 Rev. Mod. phys. 1976, 48,587.
s,;
(45) Lombardero,M.;Abascal, J. L. F.;Lago,S. Mol.Phys. 1981,42, 999. Abascal, J. L. F.; Martin, M.; Lombardero, M.; Vasquez, J.; Banon, A.; Santamaria, J. J. Chem. Phys. 1985,82, 2445. (46) Weeks, J. D.; Chandler, D.; Andereen, H. C. J. Chem. Phys. 1971,54, 5237. (47) Pant, p. v. K.; Boyd, R. H. Mac~omolecules19929% 494. (48) Lopez-Rodriguez,A.; Vega, C.; Freire, J. J.; Lago, S. Mol. Phys. 1991, 73, 691. (49) McCoy, J. D.; Honnell, K. G.; Schweizer, K. S.; Curro, J. G. J. Chem.Phys. 1991,95,9348. Slonimskii, G. L.; Askadskii,A. A.; Kitaigorodskii, A. I. Polym., Sci. USSR 1970, 12, 556.
(50)Olabiei, 0.; Simha, R. Macromolecules 1975,8,206. (51) Honnell, K. G.; et al., to be submitted to J. Chem. Phys.
