Fluid Phase Equilibria 155 Ž1999. 75–83
Second virial coefficients of Lennard–Jones chains Y.C. Chiew
a,b,)
, V. Sabesan
a
a
b
Department of Chemical and Biochemical Engineering, Rutgers UniÕersity, Piscataway, NJ 08854, USA Department of Chemical Engineering, National UniÕersity of Singapore, Kent Ridge Crescent, Singapore 119260, Singapore Received 1 April 1998; accepted 5 November 1998
Abstract The second virial coefficients B2 of Lennard–Jones chain fluids were calculated through Monte Carlo integration as a function of chain length m Žup to 48 segments. and temperature. We found that at a fixed temperature the second virial coefficient decreases with chain length. At low temperatures, the virial coefficient changes sign from positive to negative as m increases. The simulation data also provide an estimate for the theta temperature TQ at which the attractive and repulsive interactions cancel each other for dilute solutions. It is found that the theta temperature TQ for Lennard–Jones chains with m ) 32 is 4.65 independent of chain length m. A comparison of simulated values of B2 with those evaluated from two different perturbation theories for chain fluid shows that these approximate theories underestimate the second virial coefficients of Lennard–Jones chains. q 1999 Elsevier Science B.V. All rights reserved. Keywords: Statistical thermodynamics; Chain fluids; Molecular simulation; Perturbation theories
1. Introduction Recently there has been considerable theoretical effort directed at modeling the equilibrium volumetric Ž P–Õ–T . and phase behavior of polymers and chain-like fluids. Several equations of state based on well-defined chain molecular models have been developed. Among these molecular models is the Lennard–Jones chains ŽLJC. . In this model, each chain molecule is assumed to consist of m linearly bonded segments that are allowed to rotate around its adjacent neighbors freely. Segments in these molecules interact through the Lennard–Jones potential energy. Despite the simplicity of these models, they do capture effects of chain connectivity, repulsion and attraction between chain
)
Corresponding author. Tel.: q65-874-8044; Fax: q65-779-1936; E-mail:
[email protected]
0378-3812r99r$ - see front matter q 1999 Elsevier Science B.V. All rights reserved. PII: S 0 3 7 8 - 3 8 1 2 Ž 9 8 . 0 0 4 5 7 - 9
Y.C. Chiew, V. Sabesanr Fluid Phase Equilibria 155 (1999) 75–83
76
segments. Also, it offers the advantage of being continuous and combining harsh repulsion with short-range attraction. Banaszak et al. w1x used a first order thermodynamic perturbation theory to calculate the Helmholtz free energy and compressibility factor of Lennard–Jones chains. O’Lenick and Chiew w8x calculated these properties for Lennard–Jones chains based on a variational perturbation theory. In these two approaches, the theories yield the ideal gas limit at zero density; however, they only approximately capture the low density P–Õ–T behavior as they do not give exact values of the second virial coefficient B2 of Lennard–Jones chain fluids. The purpose of this paper is to examine and calculate the second virial coefficient B2 of Lennard–Jones chains through Monte Carlo integration. Monte Carlo integration permits us to compute ‘exact’ values of B2 as a function of temperature T and chain length m. The second virial coefficient becomes equivalent to the osmotic second virial coefficient of a polymer solution if one interprets the segment–segment interactions as solvent-mediated potential of mean force. The theta temperature can be interpreted as the temperature at which B2 vanishes. Earlier Monte Carlo studies have calculated the second virial coefficients of athermal hard-sphere chains w11x and square-well chains w2,5x. In this work, short Lennard–Jones chains with chain length ranges from m s 2 to m s 48 and reduced temperature T ) s kTr´ ranges from 4 to 8 are examined. This paper is organized as follows. In Section 2 we briefly describe the Monte Carlo integration method used. Section 3 presents results obtained from this study. Concluding remarks are reported in Section 4.
2. Monte Carlo method The virial expansion for the compressibility factor Z ' Prrc kT of a fluid can be written as Z s 1 q B2 rc q B3 rc3 q PPP
Ž1.
where rc refers to the number of chains per unit volume, B2 and B3 represent the second and third virial coefficients, respectively, P is pressure, T is absolute temperature, and k denotes Boltzmann’s constant. The second virial coefficient B2 can be calculated from Refs. w10,5x: B2 s
1 2
yU1Ž V 1 .r kT yU2 Ž V 2 .r kT
He
e
12 Ž R 12 .r kT
Ž 1 y eyU
. d R12 d V 1d V 2rn
Ž2.
where ns HeyU1Ž V 1.r kT d V 1 HeyU2Ž V 2 .rkT d V 2 , U1Ž V 1 ., U2 Ž V 2 . are intra-molecular energies of the two single polymer molecules, U12 Ž R 12 . is the intermolecular potential between the two chains. The integrals Hd V 1 and Hd V 2 are taken over all internal degrees of freedom of the two polymers, and the integral Hd R 12 is taken over all separations between the centers of mass of the two chains. To calculate the second virial coefficient B2 , single polymer chain conformations V 1 and V 2 must first be generated before the integral given in Eq. Ž2. can be evaluated. Chain conformations were generated using the pivot algorithm w5x. In this algorithm, the ‘trial’ conformation of a chain is generated as follows: a segment on a chain, of a fixed length m, is randomly selected as the pivot point. The short end of the chain is then rigidly rotated through a randomly selected angle to a trial position. The intra-chain energy U trial Ž V 1 . of polymer chain for the trial conformation is then calculated. The trial conformation will be accepted as the new conformation with probability
Y.C. Chiew, V. Sabesanr Fluid Phase Equilibria 155 (1999) 75–83
77
minwexpŽ U trial y U old .rkT, 1x. This is in essence the Metropolis algorithm and allows us to sample low energy chain conformations. In our work, a new chain conformation sample is obtained after a series of 20 pivot rotation attempts. When new conformations of two chains were generated, their respective centers of mass Ž X 1 and X 2 . and radii of gyration Ž R g1 and R g2 . were calculated. The centers of mass of these two chains were initially placed at a separation R 12 s 10 R g2 apart at a randomly chosen mutual orientation V 1 – V 2 . The integral given by Eq. Ž 2. was then numerically evaluated as a function R 12 , by moving polymer 1 towards polymer 2 at intervals D R 12 s 0.1 s , to yield a value for B2 . Note that the integration over the intra-molecular energies and their associated Boltzmann weight are accounted for in the way the chain conformations are generated through the pivot algorithm. Without generating new chain conformations for both chains, a new mutual chain orientation V 1 – V 2 between the two chains was obtained by randomly and rigidly rotating chain 1 Žvia rotation of three Euler angles. while holding the orientation of chain 2 fixed. With this new V 1 – V 2 , ten additional values of B2 were computed through integration over R 12 . This procedure reduces the number of new chain conformations Ž V 1, V 2 . that have to be generated and proves to be an efficient method for generating B2 values. Note that the above procedure Ž namely, generating new V 1 – V 2 by holding V 2 fixed. introduces very little inaccuracy and error in the computed value of B2 since only a small number of relative orientations V 1 – V 2 are sampled before completely new conformations are generated Ž V 1, V 2 .. 3. Results and discussion The variation of the reduced second virial coefficient B2rm2s 3 of short Lennard–Jones chains with chain length m at T ) s kTr´ s 4.0 is shown in Fig. 1. The second virial coefficients of Lennard–Jones spheres are included in the figure. It can be seen that as chain length increases from m s 2 to 4, B2rm2s 3 changes from positive to negative suggesting that the effective interaction between two chains becomes increasingly attractive. Similar plots for B2rm2s 3 vs. m at higher temperatures, i.e., T ) s 5 and 8 are displayed in Figs. 2 and 3, respectively. In both cases, the second
Fig. 1. The reduced second virial coefficient B2 r m2s 3 vs. chain length m at reduced temperature T ) s kTr ´ s 4.0.
78
Y.C. Chiew, V. Sabesanr Fluid Phase Equilibria 155 (1999) 75–83
Fig. 2. The reduced second virial coefficient B2 r m2s 3 vs. chain length m at reduced temperature T ) s kTr ´ s 5.0.
virial coefficients B2rm2s 3 are found to be positive indicating that the integrated interactions between Lennard–Jones chains are effectively repulsive at these two temperatures. The second virial coefficients of Lennard–Jones chains B2rm2 vs. chain length m for different values of reduced temperature T ) s ´rkT are shown in Table 1. Fig. 4 shows the second virial coefficient B2rm2s 3 plotted as a function of reduced segment–segment energy ´rkT or the reciprocal of dimensionless temperature T ) for different chain length. The second virial coefficient B2rm2s 3 for Lennard–Jones spheres Ž m s 1., calculated using the correlation proposed by Koutras et al. w6x, is included in the figure and represented by the solid curve. The theta temperature, TQ) , which is defined as the temperature at which the second virial coefficient vanishes can be obtained from the figure. At the theta temperature the repulsive and attractive
Fig. 3. The reduced second virial coefficient B2 r m2s 3 vs. chain length m at reduced temperature T ) s kT r ´ s8.0.
Y.C. Chiew, V. Sabesanr Fluid Phase Equilibria 155 (1999) 75–83
79
Table 1 Second virial coefficients of Lennard–Jones Chains B2 r m 2 vs. chain length m for different values of reduced temperature T ) s ´ r kT T)
m
B2 r s 3
B2 r s 3 m 2
4.0
2 4 8 16 24 32 40 48 2 4 8 32 48 2 4 8 16 32 48
0.098 y1.166 y6.411 y28.630 y64.770 y120.520 y201.300 y290.490 1.050 2.150 4.650 46.480 103.000 2.270 6.150 18.400 58.100 200.620 436.400
0.0244 y0.0729 y0.1002 y0.1118 y0.1124 y0.1177 y0.1258 y0.1261 0.2625 0.1344 0.0727 0.0454 0.0447 0.5675 0.3844 0.2875 0.2270 0.1959 0.1894
5.0
8.0
interactions between two chain molecules cancel each other. From Fig. 4 it can be seen that the theta temperature for Lennard–Jones spheres occurs when 1rTQ) ( 0.29. The theta temperature for Lennard–Jones chains can easily be obtained by interpolating the simulation data. The two curves that go through the m s 8 and m s 48 data were best fits of the data using a second order polynomial. It is seen that for m s 8, the reciprocal of the theta temperature 1rTQ) ( 0.22; while for m s 32 and m s 48, 1rTQ) ( 0.215. These results suggest that TQ) increases from a value of approximately 3.45 for m s 1 to a value of approximately 4.65 for m s 48. Specifically, we find that TQ) increases for small m and seems to become insensitive to chain length for m ) 32. The value of 1rTQ) ( 0.215 is lower than the value of 0.32 for square-well chains w5x. These authors also found that the theta temperature is independent of chain length for moderate and long square-well chains. 3.1. Comparison with chain equations of state Several equation of state models for Lennard Jones chains have been reported in the literature. Among them include the thermodynamic perturbation theory Ž TPT. by Banaszak et al. w1x and the variational theory by O’Lenick and Chiew w8x. Both theories are approximate and are based on different perturbation approaches. In the TPT formulation and following Banaszak et al. w1x, the equation of state of Lennard–Jones chain molecules may be written as: Z s mZ LJS q Z chain
Ž3.
Y.C. Chiew, V. Sabesanr Fluid Phase Equilibria 155 (1999) 75–83
80
Fig. 4. Variation of the reduced second virial coefficient B2 r m2s 3 with 1r T ) for ms1, 8, 16, 32 and 48. The ms1 data are taken from Ref. w1x. The theta temperature corresponds to the temperature at which the second virial coefficient B2 r m 2s 3 vanishes.
where Z LJS Ž r ) , T ) . is the compressibility factor of the non-bonded Lennard–Jones spheres while Z chain is given by
E Z chain s Ž 1 y m . 1 q r )
Er )
ln g LJS Ž s , r ) ,T ) .
Ž4.
Here, r ) is the segment density, T ) is the dimensionless temperature, and g LJS Ž s , r ) , T ) . refers to the radial distribution function of the Lennard–Jones sphere evaluated at r s s . Several versions of the TPT theory for Lennard–Jones chains have been reported in the literature; they differed in the expressions used and approximations made for Z LJS Ž r ) , T ) . and g LJS Ž s , r ) , T ) .. The compressibility factor Z LJS Ž r ) , T ) . of Lennard–Jones spheres is well represented by the Nicolas BWR equation w7x. To second order in r ) , the Nicolas equation gives Z LJS s 1 q
a1 T)
r ) q PPP
Ž5.
where
a1 T
)
s x1 q
x2
'T
)
q
x3 T
)
q
x4 T
)2
q
x5 T )3
.
Ž6.
The empirical parameters are x 1 s 0.862309, x 2 s 2.976219, x 3 s y8.402230, x 4 s 0.105414, and x 5 s y0.856458. The radial distribution function g LJS Ž s , r ) , T ) . can be approximated by g LJS Ž s , r ) ,T ) . ( g HS Ž d; r d 3 . exp Ž yw Ž r . rkT .
Ž7.
Y.C. Chiew, V. Sabesanr Fluid Phase Equilibria 155 (1999) 75–83
81
Fig. 5. The second virial coefficient B2 r m 2s 3 is plotted as a function of chain length m for varying temperature T ). Symbols represent Monte Carlo simulation results while curves are values computed from the TPT theory, Eq. Ž9..
where g HS Ž d; r d 3 . is the contact radial distribution function of a hard-sphere system with a temperature dependent diameter d that is evaluated using the Barker–Henderson perturbation theory, i.e., 0.3837 q 1.068rT )
d
s
s 0.4293 q 1rT )
.
Ž8.
If we combined Eqs. Ž4. – Ž7. with Eq. Ž3., and expand the compressibility factor in density, we found that the second virial coefficient can be written as B2 m 2s 3
a1
s
q
3
Ž 1 y m . 5p d
T)
m
12
ž / s
.
Ž9.
Fig. 5 shows values of reduced second virial coefficient obtained from this study and those calculated from the TPT theory outlined above. It is seen that the TPT theory consistently underestimates the second virial coefficient. In the variational theory of O’Lenick and Chiew w8x, the compressibility factor of the chain fluids can be written as P
rc kT
P
s
ref
pert
P
ž / ž / rc kT
q
Ž 10a.
rc kT
where P
ž / rc kT
ref
sm
1qhqh2yh3
Ž1 y h .
3
y Ž m y 1.
1 q hr2
Ž1 y h .
2
,
Ž 10b.
Y.C. Chiew, V. Sabesanr Fluid Phase Equilibria 155 (1999) 75–83
82
and pert
P
ž /
12 mh s
rc kT
T a )
6
JAŽ h ,m . q
12 mh T a )
6
ž
1
a6
/
y 1 J B Ž h ,m . .
Ž 10c.
Here, JAŽ h ,m . ' IA q h
E IA
ž / Eh
and J B Ž h ,m . ' I B q h m
E IB
ž / Eh
Ž 10d. m
where IA and I B are evaluated from `
ž
IA s
H0
IB s
H0
`
4
4
z 12 4
z 12
y
z6
/
HSC g inter Ž z ; r )a 3 ,m . z 12 dz
HSC g inter Ž z ; r )a 3 ,m . z 2 dz .
Ž 11a. Ž 11b.
HSC Ž . In Eqs. Ž10a. and Ž 10b. the function g inter z refers to the inter-chain segment–segment correlation function which is computed from the Percus–Yevick approximation w3,4x. The values of integrals IA and I B , and the optimal diameter were numerically computed as outlined in Refs. w8,9x for different values of chain density and temperature. These numerical results were used in combination with Eqs. Ž10a. , Ž10b. and Ž10c. to calculate the second virial coefficients. Fig. 6 shows values of B2rm2s 3 calculated from the variational theory and Monte Carlo simulations. Again it is found that, like the
Fig. 6. The second virial coefficient B2 r m 2s 3 is plotted as a function of chain length m for varying temperature T ). Symbols represent Monte Carlo simulation results while curves are values computed from the variational theory, Eqs. Ž10a., Ž10b. and Ž10c..
Y.C. Chiew, V. Sabesanr Fluid Phase Equilibria 155 (1999) 75–83
83
TPT theory, the variational theory underestimates B2rm2s 3. The discrepancies between the simulation data and those calculated from these theories can be attributed to the failure of these theories to take into account intra-chain interactions correctly.
4. Summary The second virial coefficients of Lennard–Jones chain molecules were obtained through Monte Carlo integration. We found that, at a fixed temperature, the second virial coefficient decreases with chain length. The dimensionless theta temperature TQ) s kTQ r´ for the system is found to be 4.65 in the long-chain limit. Further, a comparison of the simulation data with values obtained from the TPT and variational theory shows that both theories underestimate B2 . Acknowledgements This research is partially supported by the Petroleum Research Fund administered by the American Chemical Society.
References w1x w2x w3x w4x w5x w6x w7x w8x w9x
M. Banaszak, Y.C. Chiew, R. O’Lenick, M. Radosz, J. Chem. Phys. 100 Ž1994. 3803–3807. C.P. Bokis, M.D. Donohue, C.K. Hall, Ind. Eng. Chem. Res. 33 Ž1994. 146–150. Y.C. Chiew, J. Chem. Phys. 93 Ž1990. 5067–5074. Y.C. Chiew, Mol. Phys. 73 Ž1991. 359–373. J. Dautenhahn, C.K. Hall, Macromolecules 27 Ž1994. 5399–5412. N.K. Koutras, V.I. Harismiadis, D.P. Tassios, Fluid Phase Equilibria 77 Ž1992. 13–38. J.J. Nicolas, K.E. Gubins, W.B. Streett, D.J. Tildesley, Mol. Phys. 37 Ž1979. 1429–1454. R. O’Lenick, Y.C. Chiew, Mol. Phys. 85 Ž1995. 257–269. Von Solms, N., O’Lenick, R., Chiew, Y.C., 1998. Lennard–Jones chain mixtures: variational theory and Monte Carlo simulation results. Molecular Physics, in press. w10x Yamakawa, H., 1971. Modern Theory of Polymer Solutions. Harper and Row, New York. w11x A. Yethiraj, K.G. Honnell, C.K. Hall, Macromolecules 25 Ž1992. 3979–3983.