1222
J. Phys. Chem. 1994,98, 1222-1231
Simulation of the Temperature Dependence of Mechanical Properties of Polyethylene Daniel J. Lacks and Gregory C. Rutledge' Department of Chemical Engineering, Massachusetts Institute of Technology, Cambridge, Massachusetts 01239 Received: June 24, 1993; In Final Form: November 9, 1993"
Simulations of crystalline polyethylene at finite temperature are carried out using a molecular mechanics force field for the interatomic potential and quasi-harmonic lattice dynamics for the vibrational free energy. T h e thermal expansion coefficient, determined by direct minimization of the free energy, is in excellent agreement with experimental results for temperatures u p to 250 K and remains in reasonable agreement with experiment throughout the range of temperatures for which experimental results exist. The axial Young modulus a t 300 K is found to be 280 GPa, which agrees with the results of Raman scattering experiments with corrections for the effects of interlamellar coupling. The vac and Vbc Poisson ratios a r e found to increase substantially with temperature, whereas the Vab and Vba Poisson ratios are relatively temperature independent. T h e Gruneisen parameter is calculated as a function of temperature and is in agreement with the experimental values. The results of calculations using classical mechanics and quantum mechanics for the vibrational energy are compared to assess the importance of quantum effects as a function of temperature. Quantum effects are largest for the c-axis thermal expansion coefficient, which is in error by a factor of 2 a t room temperature when classical mechanics is used.
I. Introduction The properties of semicrystalline polymer materials are poorly understood, even though they are important and widely used materials. The difficulties arise from the fact that fully crystalline samples are impossible to synthesize, and the properties of the material depend on the extent of the crystallization and the orientation of the crystallites, as well as other factors. Investigators have attempted to determine the intrinsic crystalline properties by carefully designing the experiments to probe only the crystalline part of the sample. Nevertheless, coupling via continuity of chains between amorphous and crystalline phases is unavoidable. For example, the axial Young modulus of crystalline polyethylene is not known to better than 50% accuracy, even though it has been intensively studied for almost 30 years-careful experiments by many investigators find different values for the modulus, depending on the experiments used.' Due to the inherent difficulty in determining the intrinsic properties of polymer materials from experiment, computer simulation, which does not suffer from the experimental difficulties involving the amorphous regions, can yield much information and insight. Materials can be modeled as collections of atoms that are in thermal motion under the influence of the bonded and nonbonded interactions. The simulation can be considered in two parts: the determination of the interactions between the atoms and the simulation of the thermal (vibrational) motion of the atoms. Recently, Karasawa, Dasgupta, and Goddard (KDG) introduced an accurate force field describing the interatomic interactions in polyethylene.2 In their work, KDG calculated properties of the static lattice only and did not completely incoporate thermal effects. For example, KDG used only the potential energy, and not the free energy, in the determination of the elastic constants a t finite temperature. Sorensen et al. also developed an accurate force field for polyethylene and carried out static lattice calculation^.^ The results of the calculations by KDG and Sorensen et al. are very similar. The present paper focuses on the simulation of the thermal motion of the atoms to determine the temperature dependence of the properties of polyethylene. The simulation of the dynamics of the atoms (or more rigorously thenuclei), and the associated energy and free energy, is a difficult problem. In a correct treatment, the calculation should be @
Abstract published in Advance ACS Abstracts, January 1, 1994.
0022-365419412098-1222%04.50/0
quantum mechanical, as quantum mechanical effects are expected to beimportant when thevibrationalfrequencies v > k T / h (=200 cm-l a t room temperature), and in polyethylene 75% of the vibrational frequencies are over 700 cm-l. A classical calculation will generally overestimate thermal effects, as the vibrational modes become thermally excited independent of the vibrational frequency in classical mechanics, whereas in quantum mechanics the modes become excited only when the thermal energy can contribute a quantum of energy hv. The nature of the "motion" of the ground-state quantum vibration and the classical vibration also differs. In the classical vibration, the particle oscillates back and forth, moving fastest in the middle and slowest at the ends of its path (when it turns around), leading to a probability distribution which peaks at the ends of the path and is lowest in the middle. In contrast, in the ground-state quantum vibration, the probability distribution peaks in the middle and is lowest at the ends, which is incompatible with an interpretation of a particle moving back and forth in a harmonic potential. The exact quantum mechanical solution requires methods, such as path integral Monte Carlo, which are currently only practical for small systems of atoms with simple (central force) interacti~ns.~ However, within the quasi-harmonic approximation the solution can be efficiently ~ b t a i n e d .The ~ quasi-harmonic approximation neglects terms higher than second order in the Taylor expansion of the potential energy for the set of lattice parameters and reduces the many-body problem to many one-body problems.6 The quasiharmonic approximations is different from the harmonic approximation in that the lattice parameters can vary with temperature, causing the harmonic potential energy surface and the vibrational frequencies to vary with temperature. Therefore, anharmonic effects due to volume changes are incorporated in the quasi-harmonic model. The anharmonic effects which result from the fact that the expectation value for the position of the excited states of an anharmonic oscillator differs from that of the ground state, and the anharmonic coupling of the vibrational modes, are neglected.' To test the quasi-harmonic approximation, experiments which measured the change in vibrational frequencies with temperature and pressure were carred out for ionic crystals* and for p~lyethylene.~These experiments showed that for the modes examined the frequency variations were related to the volume variations in a nearly identical way for both the temperature and pressure experiments, indicating that the 0 1994 American Chemical Society
Mechanical Properties of Polyethylene dominant anharmonic effects are due to volume variations (which are included in the quasi-harmonic approximation) and that the thermally induced anharmonic coupling effects (which are not included in the quasi-harmonic approximation) are relatively small. We examined the accuracy of the present variable-volume quasi-harmonic method for the determination of equilibrium structures and properties elsewhere by comparing the quasiharmonic results to exact results obtained by Monte Carlo simulation and found that properties are accurately calculated up to two-thirds the melting temperature for the Lennard-Jones solid.I0 The calculation of the vibrational motion for the complete, anharmonic potential energy surface can only realistically be carried out with classical mechanics. The argument can be made that at room temperature the thermal behavior of polyethylene is controlled by the low-frequency modes which correspond to interchain vibrations, and since quantum effects are not important for these low-frequency modes at room temperature, classical mechanics should be sufficiently accurate for the simulation. This approach assumes that inappropriate treatment of the highfrequency vibrations will not significantly affect the results. A classical simulation will clearly deteriorate at low temperatures. With classical mechanics, the simulation can be carried out using molecular dynamics or Monte Carlo method^.^ Besides neglecting all quantum effects, these methods converge very slowly for lattice parameters in constant-pressure simulations,ll making it extremely difficult to simulate properties such as thermal expansion coefficients and elastic constants. Also, it is not straightforward toobtain statistical properties, such as the freeenergy andentropy, with these methods, and several lengthy simulations are necessary to obtain such properties for a single t e m p e r a t ~ r e . ~ 11. Computational Method
The Born-Oppenheimer approximation separates the fully quantum mechanical calculation of the electronic and nuclear energy into a calculation of the electronic energy of the system a t fixed nuclear coordinates and a calculation of the vibrational energy of the nuclei in the potential energy surface generated by the electronic energy.I2 a. Electronic Energy. For the study of thermal properties, it is convenient to parametrize the potential enegy surface generated by the electronic energy with a molecular mechanics force field. Note that the electronic energy includes all bonding and nonbonding interactions. The parametrization can be carried out by fitting the force field parameters to quantum mechanical calculations, experimental data, or a combination of the two. In the present calculations, we use the force field developed by KDG.2 This force field was obtained by fitting the force field parameters to a b initio calculations of butane, experimental spectroscopic data for butane, and experimental spectroscopic, structural, and energetic data for crystalline polyethylene. The KDG force field was designed to provide accurate vibrational frequencies as well as accurate structural and energetic results. For the static polyethylene lattice, the KDG force field calculates the lattice parameters and dissociation energy to within about 1% of the experimental values, and the average deviation between the calculated and observed infrared and Raman vibrational frequencies is approximately 25 cm-1.2 In our computer program, the long-range Coulombic and dispersion interactions are summed efficiently with Ewald sum techniques6J3to an accuracy of 0.0001 kcal/mol. The repulsive nonbonded interactions are summed in a direct space sum to an accuracy of 0.0001 kcal/mol. Analytical first and second derivatives of the electronic energy with respect to the atomic coordinates are calculated for use in the geometry optimizations and the calculation of vibrational frequencies. b. Vibrational Energy. The quasi-harmonic approximation, which we now review, is used for the vibrational motion of the
The Journal of Physical Chemistry, Vol. 98, No. 4, 1994 1223 nuclei. Complete derivations and discussions of this approximation are found in the l i t e r a t ~ r e . ~The . ~ potential energy surface V(q) in terms of the positions of the Nt atoms in the crystal, q = (q1,42, ..., 4 3 ~ ~can 1 , be expanded in a Taylor series about the equilibrium positions, 90,
The term V(q0) is an additive constant that represents the electronic energy of the system a t the equilibrium positions of the atoms and can be disregarded for now. Also, at the equilibrium positions of the atoms, the term (aV/aqi)l,,,, = 0. In the quasiharmonic approximation, the terms higher than second order are neglected, and so the potential energy surface is represented by
For a crystal, with lattice vectors 1 = m l a l + m2a2 + m3a3, where ai are the unit cell parameters and mi are integers, the periodicity of the crystal enforces periodicity in the displacements
qi - qio = uie&.I,
(3) where ui is the amplitude of the displacement of coordinate i, k is the wavevector describing the displacement, and li is the lattice vector for the atom associated with coordinate i. The potential energy for the displacement described by wavevector k can then be expressed in terms of the displacements of the Nuatoms in a unit cell (rather than the Nt atoms in the crystal),
3Nu 3Nu
(4)
+
where the sum over the lattice vectors I (I = m l a l + m2a2 m3a3) refers to a triple sum over the integers ml, 1122, m3, from --OD to -OD. It is convenient to switch to the mass-weighted coordinates, for which
and
where mi is the mass of the atom associated with coordinate i. The dynamics of the nuclei in a crystal is thus reduced to coupled harmonic oscillations. By diagonalizing the matrixfl,(k), a new set of coordinates, the normal-mode coordinates f ( k ) , can be obtained as linear combinations of the coordinates u' such that the potential energy is represented by 13Nu
v ( F , ~ )= ;Cgi(k)
EO)*
(6b)
i= 1
With the normal-mode coordinates, the problem reduces to
1224
Lacks and Rutledge
The Journal of Physical Chemistry, Vol. 98, No. 4, 1994
uncoupled one-dimensional harmonic oscillators, for which the solutions are known in both quantum mechanics and classical mechanics. The Hamiltonian for the crystal can be written as
where Pi;@)is the momentum associated with coordinate &(k). The number of wavevectors k, Nk, is equal to the number of unit cells in the crystal. For an infinite crystal, there are an infinite number of wavevectors. A unique set of wavevectors is defined by the Brillouin zone (BZ), for which ?r/ai < ki < ?r/ai,i = 1-3, for an orthorhombic lattice with lattice parameters ai. The vibrational freuqency oi(k) associated with coordinate &(k) is given by
The partition function for the vibrational motion of the crystal is the product of the partition functions for the uncoupled, onedimensional oscillators. The quantum and classical partition functions are given by Nk
3Nu
1
For an infinite crysal, k is a continuous variable, and the sum over k can be replaced by an integral. The Helmholtz free energy per unit cell, A, is obtained by 1 A,, = - -kT
In Q ,, =
Nk
3Nu 1 kTJBzdk c - h o i ( k ) ; 2
+ In
1
A,, = - -kT Nk
In Q,, = JBZdk
(lob)
The above derivation is for a given set of lattice parameters, and thevibrational frequencies oi(k)and consequently the Helmholtz free energy vary with the lattice parameters. To carry out the integrations over the Brillioun zone, we used Gauss-Legendrequadrature with a 4 X 4 X 4 mesh of integration points. For the k3 dimension (parallel to the chain axis), only points with k3 > 0 were used, as 02(-k) = u2(k).6 The Brillouin zone integrations are accurate to within 0.002 kcal/mol. We also performed some calculations with a 6 X 6 X 6 mesh, and the results werevery similar to the results with the smaller integration. c. Crystal Free Energy. Most experiments are performed at constant temperature and stress (or pressure), and so the equilibrium structure will be that which minimizes the Gibbs free energy. The Gibbs free energy, G(T,u),as a function of
temperature T and.stress u is C( T,a) = V(a)
+ A( T,a) - T
~
~
u
~ 1) ~ (1
ij
where V(a) is the electronic energy for the crystal with lattice parameters a (with the atoms in the minimum-energy positions for this set of lattice parameters), A(T,a) is the Helmholtz free energy for the vibrational motion, T O is the crystal volume in the absence of strain, ui, is the ij component of stress, and ei, is the ij component of the strain defined by the lattice parameters a. In the present simulations we applied the stresses ull, u22, and u33; the corresponding strains are defined by14
where ai are the lattice parameters and ai, are the lattice parameters at minimum Gibbs free energy in the absence of imposed stress. Equation 11 assumes that the excited electronic states are much higher in energy than the ground state, so that -kT In Qelectronic = W. We obtained our equilibrium structures, at a given temperature and stress, in the following manner. Because we use the parameters of KDG, and KDG found the static lattice to be orthorhombic,2we constrain the lattice to be orthorhombic. We also consider the two chains of the unit cell to be conformationally equivalent, as is observed experimentally-the space group for polyethyleneis Pnam.Is Thevariables used todescribe theatomic positions in the unit cell are 3N - 4 Cartesian coordinates of the atoms in the first chain (where N is the number of atoms in the chain per unit cell, three coordinates fix the unit cell origin, and a fourth coordinate fixes the setting angle with respect to the x-z plane so that the setting angle can be varied explicitly), the setting angle for the first chain, and the setting angle and three translational offset coordinates for the second chain. For a given set of lattice parameters, the electronic energy is minimized with respect to the variables describing the atomic positions in the unit cell, using analytical derivatives and a variable metric minimization algorithm.16 After minimization of these internal coordinates, the numerical integration of the vibrational free energy over the Brillioun zone is carried out: for each k point, analytical second derivatives, and from these the vibrational frequencies and the vibrational free energy, are calculated. The total Gibbs free energy is determined with eq 11, for this set of lattice parameters. The lattice parameters are then varied to minimize the Gibbs free energy, and the lattice parameters which yield the minimum Gibbs free energy are taken as the equilibrium lattice parameters for the given temperature and stress. By minimizing free energy with respect to atom coordinates and lattice parameters, this method does not require any experimental data for thermal expansion.
III. Results and Discussion
a. Effects of Zero-Point Energy. Previous theoretical investigations of polyethylene determined the equilibrium crystal structure at “zero temperature” by minimizing the potential energy (electronic energy) with respect to the crystallographic parame t e r ~ This . ~ ~ procedure ~ would be correct if classical mechanics were valid for the vibrations of the nuclei. However, at zero temperature, the vibrational motion must be treated quantum mechanically. In a quantum calculation, the total energy at zero temperature consists of the electronic energy plus the zero-point ho for a vibrational vibrational energy, which is equal to mode of frequency u.12 The total energy, including the zeropoint energy, must be minimized to give the zero temperature structure. As the vibrational frequencies vary with crystal structure, the quantum result for the zero temperature structure will differ from the classical static structure.
E
~
~
The Journal of Physical Chemistry, Vol. 98, No. 4, 1994 1225
Mechanical Properties of Polyethylene
TABLE 1: Lattice Parameters at Zero Temperature experiment (T= 4 K)a static lattice (classical) quantum quantum, lattice modes only quantum, CH2 group modes only quantum, CH stretch only Avitabile et al., ref 22.
(a>
a(&
b(4
c(A)
7.121 7.192 7.331 7.265 7.230 1.232
4.851 4.801 4.812 4.821 4.806 4.845
2.548 2.546 2.558 2.544 2.556 2.547
The results for the zero-temperature calculations are given in Table 1. The inclusion of the zero-point energy leads to an expansion of the lattice by almost 4%. The nonbonded force field parameters were obtained by KDG by fitting the structure, dissociation energy, and lattice mode frequencies of static lattice calculations of polyethylene to experimental data a t low temperatures. (The dissociation energy calculation included zeropoint energy, but for the minimum-energy static structure.2) The present results show that to obtain highly accurate force fields, the quantum zero temperature structure, rather than the static structure, should be used to fit the force field parameters. We note, however, that the results for the internal parameters, including the bond lengths, bond angles, and the setting angle, are very similar to those of the static lattice, which have been reported by KDG.Z The factors underlying the expansion due to the zero-point energy (the “zero-point expansion”) are different than those that cause thermal expansion. For thermal expansion, only modes that are thermally excited will influence the structure. For the zero-point expansion, in contrast, all modes contribute equally. The effects of the different vibrational modes on the zero-point energy were examined in order to understand the causes of the lattice expansion. The vibrational modes can be separated into three groups: the lattice modes, which have frequencies less than 600 cm-’; the ‘CH2 group modes” (rock, wag, twist, and scissor modes and the skeletal modes), which have frequencies between 700 and 1500 cm-*; and the CH stretching modes, which have frequencies greater than 2800 cm-I. To understand the effects of the different groups of vibrational modes, quantum calculations were carried out in which the energy from only one of the groups of vibrational modes was included. The results are given in Table 1 and show that all three groups of modes contribute to the expansion along the a axis. Along the b axis, only the lattice modes and the CH stretch modes lead to the expansion, while the CH2 group modes do not have much effect. The most interesting case is along the c axis, where the CH1 group modes lead to expansion, while the lattice modes, in contrast, lead to slight contraction. (The CH stretch modes do not have an effect along the c axis because the C H bonds are directed perpendicular to the c axis.) As discussed below, the c lattice parameter contracts with increasing temperature, and this negative thermal expansion occurs because only the lattice modes, and not the CHI group modes, are significantly excited a t temperatures below polyethylene’s melting point. The changes we obtain for the infrared and Raman frequencies with the unit cell volume are in qualitative agreement with experiments9J7J* and other c a l c ~ l a t i o n s . We ~ ~ note that the experimental changes in the frequency of the symmetric CH stretch due to volume changes alone are obscured by Fermi resonance interactions, which are neglected in a quasi-harmonic model. However, the frequency of the asymmetric C H stretch, which is of a different symmetry than the symmetric C H stretch and does not have significant Fermi resonance interactions, has been found experimentally to decrease with increasing unit cell volume,*8in agreement with our calculations. b. Thermal Expansion. Fairly significant differences are obtained by different experiments for the lattice parameters of polyethylene. Davis et ~ 1showed . ~ that ~ the lattice parameters vary with the lamella thickness and the related extent of
7 7
7.1
100
0
200
400
300
T (K)
5.00
R
I
4.95
4.90
4.85
4.80
0
2.555 2.555
100
200
300
400
f”--, -
2.550
-
I-
2.545
2.540‘
0
’
I 100
200
300
400
T (K)
Figure 1. Lattice parameters as a function of temperature. The filled
circles are the quantum calculations, the filled squares are the classical calculations,the open squares are the experimentalresults of Davis et al. (ref 20), the open triangles are the experimental results of Avitabile et al. (ref 22), and the open circles are the experimental resultsof Kobayashi and Keller (ref 21). crystallinity. The variation becomes more pronounced at higher temperatures. In the comparisons of the resultsof our simulations with experiment, we tried to use experimental results for samples with highcrystallinity and thick lamella, as the present calculations are for fully crystalline polyethylene with infinitely thick lamella. The thermal expansion simulations were carried out by minimizing the free energy at zero stress, with respect to the lattice parameters. The results for thechangeinlattice parameters with temperature are given in Figure 1. For the c lattice parameter, the quantum results are in good agreement with the experimental results of Kobayashi and KellerZ1(c = 2.552 A a t room temperature) but do not agree quite as well with the results of Davis et al. and Avitabile et ~ 1 (c < . 2.55 ~ ~A). Kobayashi and Keller comment that they found c < 2.55 A for unannealed samples, but upon annealing, the c lattice parameter becomes 2.552.21 Therefore, the disagreement with the experiments of
1226 The Journal of Physical Chemistry, Vol. 98, No. 4, 1994
Lacks and Rutledge (a)
h
2.00e-4 v
1.00e-4 ou
0
100
200
300
400
T (K)
Figure 2. Volume as a function of temperature. The filled circles are the quantum calculations,the filled squares are the classical calculations, the open squares are the experimental results of Davis et al. (ref 20), and the open triangles are the experimental results of Avitabile et al. (ref 22).
0
50
100
150
200
250
350
300
T (K)
(b) 1.2Ce-4 I
Davis et al. and Avitabile et al. may be due to complications in the experiments from the noncrystalline parts of the samples. The volume expansion is shown in Figure 2. The classical results, and not the quantum results, give very good agreement with experiment for the volume because the nonbonded force field parameters were fit to the classical zero-temperature structure, rather than the quantum structure. The results for the thermal expansion coefficients are shown in Figures 3 and 4. The thermal coefficients are defined as
I
/
1.OOe-4
4.0%-5
7
1
0
50
"
100
'
"
'
150
"
200
'
'
250
~
300
350
T (K)
(c) where the ai are the lattice parameters. The ai were obtained by taking the numerical derivatives of the equilibrium lattice parameters with respect to temperature. A salient feature of these plots is that the quantum a goes to zero as the temperature goes to zero, while the classical a does not. The calculated aa and ab are not in very good agreement with the experimental results of Davis et aLzo Unfortunately, there are no other highquality experimental results available for comparison. The calculated a,,on the other hand, are in virtually perfect agreement with the experimental results of White and Choyz3 and in good The (Ybulk are in agreement with the results of Davis et excellent agreement with the results of Engeln et to temperatures over 250 K and remain in fairly good agreement throughout the range of temperatures in which experimental results are available. The experiments of Engeln et al. are expected to be the most accurate available, as the sample in their experiments was 98% crystalline and had a lamella thickness of 4400 A.24 The sample of Davis et al., in contrast, had a lamella thickness of only 385 A.20 To assess the accuracy of our simulations, we note that a,and (Ybulk are in excellent agreement with themost recent experimental results-the results of White and Choy were published in 1984, and those of Engeln et al. were published in 1985; all of the other results were published prior to 1975. The poor results for a. at first glance appear to be due to the breakdown of the quasihormonic approximation at high temperatures, as the experimental a, increase more rapidly at the higher temperatures than do the calculated a,. However, the situation for a b is reversed-the calculated a b increase more rapidly than the experimental a b . In Figure 5 we plotted (11,b = a, a b , and the calculated a o b are in reasonable agreement with the results of Davis et al. It appears, therefore, that there are compensating errors for a, and a b either in the experiments of Davis et al. or in the present calculations. The significance of the quantum effects at room temperature can be surmised from our simulations. The difference between the quantum and the classical result for the volume is 2%, and
+
-5.000-6 h
c
-1.00e-5 n o on0
-1.50e-5.
-2.00e-5
.
-2.50e-5' 0
'
'
50
100
150
200
250
300
1 I
350
T (K)
Figure 3. Anisotropic thermal expansion coefficients as a function of temperature. The filled circles are the quantum calculations, the filled square are the classical calculations, the open squares are the experimental results of Davis et al. (ref 20), the open triangles are the experimental results of White and Choy (ref 23), and the open circle is the experimental res ult of Kobayashi and Keller (ref 21).
the differences for the a, b, and c lattice parameters are 0.7%, I%, and 0.4%. From our analysis presented in Table 1, we see that these differences correspond roughly to the effects of the zero-point energy from the CH2 group mode and the CH stretch mode frequencies. Since these modes have high vibrational frequencies, the differences between the quantum and classical results will persist until temperatures of several thousand kelvin. The classical a. and ab are smaller than the quantum values by about lo%, but the classical acis smaller than the quantum ac by almost a factor of 2. The importance of the quantum effects for acis discussed below in terms of the Gruneisen parameters. c. Elastic Constants. As mentioned in the Introduction, the results of different experiments for the value of the axial Young modulus for crystalline polyethylene are not in agreement. One of the experimental methods uses X-ray diffraction to measure the change in the lattice parameters of the crystalline regions of
The Journal of Physical Chemistry, Vol. 98, No. 4, 1994 1227
Mechanical Properties of Polyethylene 1.2Oe-4 1.m-4
0
A 0
50
100
150
250
200
300
350
T (K)
Figure 4. Bulk thermal expansion coefficientsas a function oftemperature. The filled circles are the quantum calculations, the filled squares are the classical calculations, the open squares are the experimental results of Davis et al. (ref 20), and the open triangles are the experimentalresults of Engeln et al. (ref 24).
5.m-51
o.m+oY. 0
J "
50
'
100
"
150
"
200
'
250
.
'
300
'
350
T (K) Figure 5. Thermal expansion coefficients for the direction perpendicular to the chain axis as a function of temperature. The filled circles are the quantum calculations and the open squares are the experimental results of Davis et a/. (ref 20).
the sample under stresses which are applied to the sample as a whole. Most of these experiments lead to an axial modulus of 235 GPa at room temperature and 255 GPa at 115 K,1925$26 although one of these experiments reports a modulus of 266 GPa at room t e m p e r a t ~ r e .This ~ ~ method, however, assumes that the actual stresses on the crystalline regions are equal to the stress applied to the sample as a whole. (The Young modulus for the crystalline region is determined by the stress and the strain in the crystalline region; the strain in the crystalline region is measured directly by X-ray diffraction.) Proponents of this method claim that the stress on the crystalline regions is indeed equal to the macroscopic stress, as the experiments are carried out for samples with different microstructures and the same results are obtained. Raman scattering experiments on crystallinealkanes have been used to find the axial modulus.25 The vibrational frequencies of the "accordian" modes, which correspond to the longitudinal acousticvibrations in the limit of the infinite alkane chain length, are measured. The ratio of the vibrational frequencies to the "wavevector" (taken to be equal to 2u/h, where his the wavelength for the vibrational mode and is assumed to be the length of the chain alkane for the accordian modes) for these modes gives a speed of sound, ui. The modulus Cii (in Voigt notation) is then determined from the relation
cii= v;p where p is the density. The Raman scattering experiments lead to an axial modulus of 358 GPa.28 Strobl and E ~ k e land , ~ ~later Kobayashi et al.,30pointed out that this analysis assumes that there is no vibrational coupling between the lamella. In the absence of this coupling, the ratio of the frequencies of the first three Raman-active modes (Le., with 1, 3, and 5 nodes) should be 1:3:5. However, this ratio is not strictly observed, and Strobl and Eckel argue that the discrepancy indicates that there is a
significant coupling between the lamella.29 The magnitude of the interlamellar coupling can be surmised by analyzing the deviation of the vibrational frequencies from the ratio 1:3:5. To determine the intrinsic modulus (i.e., the modulus for infinite chain molecules in a perfect crystal), corrections for the interlamellar coupling were included, and the Raman scattering experiments indicate the intrinsic axial modulus to be 280-290 GPa.28,29The modulus found in these experiments (using eq 14) is the stiffness modulus, which is similar to, but not equal to, the Young modulus. Inelastic neutron scattering experiments are used to determine the modulus by using the slope of the measured longitudinal acoustic mode dispersion curve in the limit as the wavevector goes to zero to find the speed of sound in the crystalline regions of the sample.31 The modulus is then obtained from eq 14. These experiments find the axial modulus to be 329 GPa at room t e m p e r a t ~ r e . ~The ] effects of interlamellar coupling on the modulus obtained by this method must be considered, although Fanconi and Rabolt argue that these effects are minor.32 The macroscopic axial modulus (i.e., the axial modulus for the entire sample and not just the crystalline regions) is smaller than the crystalline modulus because the amorphous regions have a much smaller axial modulus than do the crystalline regions. Elasticity experiments have been carried out on highlycrystalline samples. The highest macroscopic axial modulus observed in stress-strain experiments is 232 GPa.33 It is not clear whether this result corresponds to the intrinsic modulus or whether the modulus would be still higher in a more perfectly crystalline sample. Another experiment, using a dynamic tester, leads to a macroscopic axial modulus of 262 GPa at I1 K, and when corrections are included to compensate for systematic errors in the technique, a modulus of 288 GPa is obtained.34 The elastic stiffness moduli Cij are defined by (using Voigt notation)14 3
ui =
pij ej j= I
where ui is the i component of the stress, and e, is t h e j component of the strain. The isothermal stiffness moduli, Fi,,can be defined asI4
&]
CTij= 1 aZA ' i "j
T.rtziJ
where A is the Helmholtz free energy, Vis the unstrained volume and superscript T indicates isothermal conditions. The definition of strain used in the modulus calculations is given in eq 12. This definition of strain is different than that used in the calculation of the thermal expansion coefficientsand Gruneisen parameters (eqs 13a and 22a are based on the definition ei = Aai/aio). These two definitions do not lead to any significant difference in the calculated properties for the small strains considered in the present calculations. We chose to use the two different definitions of strain in order to be consistent with previously published results for the different properties. The elastic compliance moduli Si, are defined by14 3
ti
= csijuj j= 1
Note that the compliance moduli correspond to constant-stress
Lacks and Rutledge
1228 The Journal of Physical Chemistry, Vol. 98, No. 4, 1994
0
240
t
LL"
0
I
TABLE 2: Isothermal Mechanical Properties of Polyethylene Obtained from Quantum Calculations OK 100K 200K 300K 12.6 12.4 316 6.5 2.1 4.3 10.5 9.5 318 0.105 0"
I 100
200
300
400
T (K)
Figure 6. Axial Young modulus as a function of temperature. The filled circles are the quantum calculations, the filled squares are the classical calculation,the large open squareis the inelasticneutron scattering result (ref 31), the small open squares are X-ray diffraction results (ref l), the small open circle is an X-ray diffraction result (ref 27), the large open circle is the Raman scattering result (ref 28), the open triangles are the Raman scattering results with correctionsfor the interlamellar coupling (refs 29 and 30), and the + is the dynamic tester result (ref 34).
0"
00 00
1 .606
11.5 11.5 314 6.3 2.4 4.6 9.5 8.8 312 0.1 14 9.49 4.48 -0.706 4.42 1.17
10.2 10.4 303 5.4 3.3 5.1 7.9 7.8 301 0.132 12.5 5.93 -0.910 5.82 0.66
8.8 8.8 290 4.3 4.5 5.8 7.5 6.9 283 0.149 14.3 8.45 -1.12 7.21 0.48
400K7.2 7.7 270 3.0 7.0 6.7 5.7 6.1 259 0.202
18.5 12.1 -1.71 9.63 0.33
Set equal to the theoretical value of thermal expansion at zero temperature. At T = 10 K. I
situations, while the stiffness moduli correspond to constantstrain situations. The Young moduli, Ei, are the inverses of the diagonal compliance moduli,
Ei = l/Sii and the Poisson ratios are defined as the ratios of the off-diagonal compliances to the diagonal compliances, vij = -sijlsjj
The Poisson ratios represent the magnitude of the coupling of the lattice coordinates, in which a stress along one direction leads to a strain in another direction. To calculate the Young moduli and Poisson ratios, we performed a series of simulations at constant temperature and several different stresses to determine the change in equilibrium lattice parameters as a function of stress. This procedure, which directly mimics the X-ray scattering stress-strain experiments, gives the compliance moduli as in eq 17. Along the chain axis stresses of -2, -1, 1, and 2 GPa were applied, which led to strains of up to 0.5%. Perpendicular to the chain axis stresses of -0.06, -0.03, 0.03, and 0.06 GPa were applied, which also led to strains of up to 0.5%. By applying both negative and positive stresses, the problems due to higher-order, nonlinear elastic effects are minimized. Since the simulations are at constant temperature, the compliance moduli obtained are the isothermal compliance moduli. The results for the axial Young modulus are shown in Figure 6. Our results are in excellent agreement with the corrected Raman scattering results of Strobl and E ~ k e and l ~ Kobayashi ~ et aL30 Since we showed above that our simulations are accurate for thermal expansion (especially along the c axis), and since thermal expansion can be considered as the elastic response to a thermal stress (see discussion of the Gruneisen parameters below), we feel that our results are strong evidence that the Young modulus is approximately 280 GPa at room temperature and that the experimental analyses of Strobl and Eckel and Kobayashi et al. are correct. The zero-temperature Young modulus is found to be 317 GPa. The difference between the classical and the quantum result for the Young modulus is 7 GPa at room temperature. We noted above that the Raman scattering and neutron scattering experiments actually measure the stiffness modulus, rather than the Young modulus. To understand the extent of the difference between the values of these two moduli, we give the results for the stiffness moduli in Table 2, obtained by evaluating eq 15 numerically. The stiffness moduli can also be obtained by inverting the compliance matrix (as the lattice is orthorhombic, the compliance and stiffness matrices are block diagonal and the
0
100
200
300
400
T (K)
Figure 7. Isothermal Stiffnessmoduli as a functionof temperature. The the filled triangles are C22, and the open circle is filled circles are 91, the experimental CIIof Twistleton et al. (ref 35). shear components of these matrices do not have to be considered)-numerical errors cause the stiffness moduli calculated by the two methods to differ by an average of less than 1 GPa, which gives an indication of the precision of our elasticity calculations. At low temperatures, the axial Young modulus and the axial stiffness modulus are virtually identical, while at higher temperatures the axial stiffness modulus is slightly larger than the Young modulus. In this way only a small portion of the difference between the Raman and neutron scattering experiments on the one hand and the X-ray diffraction experiments on the other is explained. The results for the transverse stiffness moduli from the quantum calculation are shown in Figure 7. and c T 2 2 haveverysimilar values throughout the temperature range. The calculated cTl modulus is in good agreement with the inelastic neutron scattering result of Twisleton et The Poisson ratios are plotted in Figure 8. The vue and Vbc Poisson ratios, and in particular vue, increase substantially with temperature. The veUand V& are very small and are less than 0.02 for all temperatures. The Vu6 and Vb" remain relatively constant with temperature. The off-diagonal stiffness moduli are given in Table 2. The c T , 2 modulus decreases with temperature, while the c T 1 3 and c T 2 3 moduli increase with temperature for temperatures up to the melting point. KDG calculated the Young modulus at finite temperatures by calculating the second derivatives of the potential energy, rather than the free energy, for the experimentally determined structures at the various temperatures.2 We note that we can reproduce the results of KDG when we carry out our stiffness modulus calculations by evaluating eq 15 numerically, but considering only the potential energy rather than the freeenergy. Considering only the potential energy leads to an axial Young modulus of 337 GPa at T = 4 K, which decreases to 318 GPa at T = 303 K.2The
Mechanical Properties of Polyethylene
The Journal of Physical Chemistry, Vol. 98, No. 4, 1994 1229
in,
"-1
I
volume. The heat capacity for a vibrational mode with frequency w is obtained quantum mechanically as
>
and the total heat capacity at constant volume is 3N..
V."
0
200
100
300
400
c,
T (K)
Figure8. Poisson's ratios as a function of temperature. The filled circles are voc, the filled triangles are v k , the filled squares are V#b, and the open squares are v b . The we. and Ucb are less than 0.02 for all temperatures. These are results of the quantum calculation.
'
6
9 x
0.18
.
(23b) SBZdk
The y can also be defined in terms of the vibrational frequencies in the quasi-harmonic appro~imation,2~
R
0.20 -
h
I
=
yi = -
c,
where w is the vibrational frequency. In the present calculations we obtain t h e y by taking the numerical derivatives of the entropy with the respect to the lattice parameters, as in eq 22a. The entropy is calculated by
0.16 -
0
200
100
300
400
T (K)
Figure 9. Isothermal bulk compressibility as a function of temperature. These are results of the quantum calculation.
inclusion of the vibrational free energy, however, lowers the Young modulus by 20 GPa at T = 4 K and by 40 GPa at T = 303 K. Also, calculations which consider only the potential energy find the C13 and C23 moduli to decrease with temperature,2 whereas calculations which include the vibrational free energy show that these moduli increase with temperature. However, it should be emphasized that the breakdown of the quasi-harmonic approximation for temperatures above roughly two-thirds the melting temprature limits the accuracy of elastic constants calculated above 300 K. The moduli we have presented above are isothermal moduli. The values of the adiabatic moduli relative to the isothermal moduli are discussed below in terms of the Gruneisen parameters. d. Bulk Compressibility. The isothermal bulk compressibility, x, is defined as
X=#glT where P is the pressure and V is the volume. An alternative expression for x is
j=l
j=1
The results for x,which were obtained from eq 21, are shown in Figure 9. e. Gruoeisen Parameters. The Gruneisen parameters y are defined as36 yi=--4'
as
C, dai
where S is the entropy and C,is the heat capacity at constant
1
-
(25)
where Avib is the vibrational Helmholtz free energy and &b is the vibrational energy, given by Evib,qm
--
I
JBZdk
The utility of the Gruneisen parameters is for the description of thermal expansion, which can be expressed a d 4
where V is the volume and STij is the isothermal compliance. Since ai is a strain per degree of temperature, the quantity (C,/ V)yj represents a thermally-induced stress per degree of temperature. The term Cu/V is an energy density per degree of temperature, and so the Gruneisen parameter yj can be interpreted as the thermally-induced stress per unit of thermal energy. With eq 27, the thermal strain along an axis can be considered as the elastic response to the thermal stress along that axis, plus the elastic response to the thermal stresses along the other axes. Therefore, the accuracy of our calculated CY (for which there is much experimental data for comparison) indicates that our calculated STij(for which there is little experimental data for comparison) are also accurate. The results for yi are shown in Figure 10. The yoand yb are positive, while ycis negative: the entropy increases with increasing a and b lattice parameters, while the entropy decreases with increasing c. Equation 24 shows that y represents the change in vibrational frequency with strain, but only of those vibrational frequencies that contribute to the heat capacity (i.e., have a significant probability of being excited). As discussed above in section a, the lattice mode frequencies show a net decrease with
Lacks and Rutledge
1230 The Journal of Physical Chemistry, Vol. 98, No. 4, 1994 30 -
4++--1
25
2
-
20 -
0
T (K) Figure 10. AnisotropicGruneisen parameters as a function of temperature. The filled circles are 71,the filled triangles are 72,and the filled squares are 73. These are results of the quantum calculation.
1
*I
tA
-0
100
200
300
400
T (K)
Figure 11. Bulk Gruneisen parameter as a function of temperature. The filled circles are thequantum result, theopen squares are the experimental results of Engeln et al. (ref 24), and the open trianglesare the experimental results of Gibbons (ref 36).
c-axis strain, while the CH2 group mode frequencies show an increase. Because the lattice mode frequencies are less than 600 cm-l, while the CHI group mode frequencies are greater than 700 cm-I, only the lattice modes are significantly excited a t temperatures less than 400K. Therefore, ycis determined almost exclusively by the lattice modes. We now explain why the classical acis significantly smaller than the quantum ac (see Figure IC). In classical mechanics, which corresponds to a high-temperature limit of the quantum mechanics, all vibrational modes contribute the same amount to the heat capacity. Therefore, in the classical calculation, the CHI group mode frequencies contribute to the sum in eq 24 with terms of the opposite sign to those of the lattice modes, leading to a decreased Gruneisen coefficient and therefore a decreased magnitude of thermal expansion relative to the quantum calculation. The reason that the classical calculation leads to negative thermal expansion while the zero-point energy leads to positive expansion over the static (classical) calculation, even though in both cases all of the vibrational modes contribute fully, is that the driving force for thermal expansion is to increase the entropy, which classically is proportional to W - I , while the driving force for the zero-point expansion is to decrease the zero-point energy, which is proportional tow. Therefore, in the classical calculation (or at the high-temperature limit of a quantum calculation) the thermal expansion is influenced more by lower frequencies (e.g., the lattice modes), while the zero-point expansion is influenced more by higher frequencies (e.g., CH2 group modes). The results for ybull; are shown in Figure 1 1. The agreement with the available experimental results is very good. (The experimental -,'bulk were not measured directly but calculated by other investigators from the measured thermal expansion coefficient, bulk compressibility, and heat capacity.) We see that the unusually low ybulk for polyethylene is due to the partial cancellation of the positive yi perpendicular to the chain axis with the negative yi along the chain axis.
100
200
300
400
Figwe 12. Heat capacityat constant pressure as a functionof temperature. The filled squares are the results of our quantum calculation, and the open circles are the experimental results of Gaur and Wundcrlich (ref 37). The heat capacity obtained from a classical calculation will be greater than 74 J/(mol K) for all temperatures.
To check the consistency of our results, we used eq 27 to calculate ai and compared these values of ai with the values of ai obtained above by taking the numerical derivatives with respect to temperature of the equilibrium lattice parameter. The results for ai calculated from the two methods are very similar, which shows that our calculations of the various properties are consistent and indicates that the numerical errors in our procedures are relatively small. The Gruneisen parameters can be used to obtain the adiabatic elastic constants, Gi,, from the isothermal constants, CTu,l4
For all the moduli, the difference between the adiabatic and isothermal moduli is found to less than 1 GPa at all temperatures. f. Heat Capacity. The expression for the heat capacity at constant volume, Cv,was given above in eq 23. The heat capacity a t constant pressure, C,, can be obtained byI4 3
cp= C"[1 + T&Yi] i= 1 The heat capacities are shown in Figure 12 and compared with the experimental C, of Gaur and W~nderlich.~'The agreement between our results and experiment is excellent to over 350 K. Our results are very similar to those of KDG, but we note that KDG used experimental data to set the finite temperature lattice parameters and to obtain C, from C,, whereas our simulations use only the force field as input. A calculation using classical mechanics would give C, greater than 74 J/(mol K) for all temperatures.
IV. Conclusions The structural and mechanical properties of polyethylene were simulated at temperatures ranging from 0 to 400 K. The major
approximations involved in the simulations were a molecular mechanics force field for the interatomic interaction energies and the quasi-harmonic approximation for the vibrational free energy. Experimental results are most abundant and accurate for the thermal expansion coefficients, and so this property allows the best assessment of the quality of our results. The thermal expansion coefficients are calculated very accurately to temperatures over 250 K and remain fairly accurate to at least 330 K, which is the highest temperature for which experimental results exist. The decreasing agreement of the calculated thermal expansion coefficients above 250 K is most likely due to the worsening of the quasi-harmonic approximation, although it could also be due
Mechanical Properties of Polyethylene to the decrease in accuracy of the KDG force field as the system moves further from the potential energy minimum, inaccuracies in the experimental results which are expected to increase with temperature, or a combination of these reasons. Studies on rare gas systems found that the quasi-harmonic approximation is accurate to temperatures of about two-thirds of the melting t e m p e r a t ~ r e . ' ~Since , ~ ~ the melting temperature of polyethylene is 41 1 K, our results appear to be in agreement with these estimates. The axial Young modulus is calculated to be 280 GPa a t room temperature, which agrees very well with the results of Raman scattering experiments that have been corrected for interlamellar coupling. We expect that the value for the axial modulus is accurate because the overall thermal expansion and the thermal expansion in the direction along the chain axis were found to be very accurate, and the thermal expansion coefficients and the elastic moduli are related parameters (see eq27). Thecalculations of Sorensen et al., which use a force field totally independent of the KDG force field, find a static lattice axial modulus very similar to that of KDG.3 We take this agreement as further evidence of the accuracy of the KDG force field and thus our simulations. The very good agreement with experiment of the thermal expansion coefficient a t lower temperatures indicates that the KDG force field is accurate, a t least near the potential energy minima. Especially encouraging is the excellent agreement with experiment of the ac,as this represents the thermal expansion along the chain axis, and the bonded interactions determined by the KDG force field are the most important along this axis. The KDG force field could be improved further by using the quantum zero-temperature structure, rather than the classical zerotemperature structure, to fit the nonbonded parameters. The results of calculations using classical mechanics and quantum mechanics for the vibrational energy were comapred to assess the importance of quantum effects as a function of temperature. Quantum mechanics is absolutely necessary for the simulation of thermal expansion coefficients at temperatures and a b less than 150 K. At room temperature, the classical CY. are smaller than the quantum values by about lo%, but the classical aCis smaller than the quantum a,by almost a factor of 2. For the axial Young modulus, the quantum result is lower than the classical result by about 20 GPa at zero temperature and 7 GPa at room temperature. Finally, we note that empirical methods (such as Einstein, Debye, and Tarasov models) have been used to calculate volume thermal expansion, bulk Gruneisen parameters, and vibrational density of states.24*39q40 These methods, however, rely on experimental heat capacity data for the system being studied, whereas our molecular methods require only the interatomic force field as input and can therefore be applied to phases for which highquality heat capacity data are not available. Our methods also hold the advantage that the directional thermal expansion and Gruneisen parameters can be calculated, whereas the empirical methods only allow the determination of the bulk thermal expansion and Gruneisen parameters. Molecular methods can thus lead to improved understanding of the connection between the material properties and the molecular structure.
The Journal of Physical Chemistry, Vol. 98, No. 4, 1994 1231 Acknowledgment. The authors are grateful to Texaco for providing the financial support for this work through the TexacoMangelsdorf Career Development Chair to G.C.R. References and Notes (1) Nakamea, K.; Nishino, T.; Ohkubo, H. J . Macromol. Sci., Phys. 1991, B30, 1 and references therein. (2) Karasawa, N.; Dasgupta, S.; Goddard, W. A. J . Phys. Chem. 1991, 95,2260. (3) Sorensen, R. A.; Liau, W. B.; Kesner, L.; Boyd,R. H. Macromolecules 1988,21,200. (4) Allen, M. P.; Tildesley,D. J. Computer Simulation ofLiquids; Oxford University Press: Oxford, 1987. (5) Boyer, L. L. Phys. Rev. B 1981,23, 3673. (6) Born, M.; Huang, K. Dynamical Theory of Crystal Lnttices; Oxford University Press: Oxford, 1954. Venkataraman, G.; Feldkamp, L. A.; Sahni, V. C. Dynamics of Perfect Crystals; MIT Press: Cambridge, 1975. (7) Leibfried, G.; Ludwig, W. In SolidState Physics; Seitz, F., Turnbull, D.. Eds.: Academic Press: New York. 1961: Vol. 12. o 275. ' (8) Taylor, J. A. Haque, M. S.; Potts, J. E. Pagk,rJrB.; Walker, C. T. Phys. Rev. B 1975,12,5969. (9) Wu, C. K. J. Polym. Sci: Polym. Phys. Ed. 1974,12,2493. (10) Lacks, D. J.; Rutledge, G. C., manuscript in preparation. (1 1) Ray, J. R. Comout. Phvs. R ~ D1988.8. . 109. SDrik. M: ImDrev. . . R. W.f Kiein, M. L. Phys.-Rev. B*1984,'29,4368.(12) K a r p l u s , M . ; P o r t e r , R. N . A t o m s and M o l e c u l e s ; Benjamin/Cummings: Menlo Park, 1970. (13) Karasawa, N.; Goddard, W. A. J . Phys. Chem. 1989,93,7320. (14) Weiner, J. H. Statistical Mechanics of Elasticity; John Wiley & Sons: New York, 1983. (15) Tadokoro, H. Structure of Crystalline Polymers; Krieger Publishing: Malabar, 1990. (16) Press, W. H.; Flannery, B. P.; Teukolsky, S.A.; Vetterling, W. T. 1
.
.
Numerical Recipes: The Art of Scientific Computing Cambridge University Press: Cambridge, 1986. (17) Wu, C. K.; Nicol, M. J . Chem. Phys. 1973,58,5150. (18) Wu, C. K.; Nicol, M. Chem. Phys. Lett. 1973,18, 83. (19) Tashiro, K.; Minami, S.; Wu, G.; Kobayashi, M. J . Polym. Sci.: Polym. Phys. Ed. 1992,30, 1143. (20) Davis, G. T.; Eby, R. K.; Colson, J. P. J . Appl. Phys. 1970,41,4316. (21) Kobayashi, Y.;Keller, A. Polymer 1970, 11, 114. (22) Avitabile, G.; Napolitano, R.; Pirozzi, B.; Rouse, K. D.; Thomas, H. W., Wills, B. T. M. J. Polym. Si., Polym. Lett. Ed. 1975,13,351. (23) White, G. K.; Choy, C. L.J. Polym. Sci.: Polym. Phys. Ed. 1984,
22,835. (24) Engeln, I.; Meissner, M.; Pape, H. E. Polymer 1985,26, 364. (25) Matsuo, M.; Sawatari, C. Macromolecules 1986,19,2036. (26) Sakurada, I.; Nukushina, Y.; Ito, T. J. Polym. Sci. 1962,57, 651. (27) Clements, J.; Jakeways, R.; Ward, I. M.; Longman, G. W. Polymer 1979,20,295. (28) Shauffele, R. F.; Shimanouchi, T. J. Chem. Phys. 1967,47, 3605. (29) Strobl, G. R.; Eckel, R. J. Polym. Sci.: Polym. Phys. Ed. 1976,14, 913. (30) Kobayashi, M.; Sakagami, K.; Tadokoro, H. J . Chem. Phys. 1983, 78, 6391. (31) Holliday,L.; White,J. W.PureAppl.Chem.1971,26,545. Feldkamp, L. A.; Venkaterman, G.; King, J. S. Neutron InelasticScattering Proc. Symp. Copenhagen; IAEA: Vienna, 1968;Vol. 2,p 159. (32) Fanconi, B.; Rabold, J. F. J . Polym. Sci.: Polym. Phys. Ed. 1985, 23, 1201. (33) Kunugi, T.; Oomori, S.;Mikami, S. Polymer 1988,29,814. (34) Barham, P. J.; Keller, A. J . Polym. Sci., Polym. Lett. Ed. 1979,17, 591. (35) Twisleton, J. F.; White, J. W.; Reynolds, P. A. Polymer 1982,23, 578. (36) Gibbons, T. G. J . Chem. Phys. 1974,60,1094. (37) Gaur, U.; Wunderlich, B. J . Phys. Chem. Ref. Data 1981,10,119. (38) Lutsko, J. F.; Wolf, D.; Yip, S.J . Chem. Phys. 1988,88, 6525. (39) Wunderlich, B. J. Chem. Phys. 1962,37, 1207. (40) Shen, M.; Hansen, W. N.; Romo, P. C. J . Chem. Phys. 1969,51,425.
