A molecular dynamics study of the thermodynamics of liquid ethane

Fluid Phase Equihbria, 48 (1989) 1-10

1

Elsevier Science Publishers B.V., A m s t e r d a m - - Printed in T h e N e t h e r l a n d s

A MOLECULAR

DYNAMICS

STUDy

OF THE THERMODYNAMICS

OF LIQUID

ETHANE

ROLF LUSTIG, ALEJAN[3q3 TORO-IABBE and WILT.TAM A. STEEIF.

Department of Chemistry, The Pennsylvania State University, University Park, Pennsylvania 16802, U.S.A.

(Received September 7, 1988) ABSTRACT Lustig, R., Toro-labbe, A. and Steele, W.A., 1989. A molecular dynamics study of the thermodyr~/~ics of liquid ethane. Molecular dynamics (MD) calculations are reported in which ethane was silm~lated along four isochores spanning the entire liquid region. The molecules interact through a two-center iennard-Jones (2CIJ) pair potential derived previously in a calculation based on the thermodynamic perturbation theory due to Kohler. Pressures, residual internal energies, and the fluctuation quantities as the residual specific heat and the thermal pressure coefficient obtained from this model agree with values given by an e~pirical equation of state within the cc~bined uncertainty of both methods for all state points. Tnis suggests that the 2CIJ potential used here is an excellent effective pair potential for fluid ethane. INTRODUCTION In accurate calculations of the thermodynamic properties of real fluids from models of the intermolecular force, it is generally recessary to take explicit account of molecular shape. In the case of siaple molecules that are assumed to have no internal degrees of freedom (rigid) and negligible permanent multipole moments, the so-called n-center Iennard-Jones

(nCiJ) potential is widely used as

an efficient effective pair interaction. Of course, determination of the correct potential parameters is crucial to the success of such calculations. When these parametems are derived by fitting theoretical to experimental second virial coefficients,

it is often found that they are not capable of accounting for the

liquid state ~ e s .

Parameter determination directly from the liquid state

data is a much more suitable procedure especially when calculating along particular lines in the phase diagrams (such as the coexistence line between vapor and liquid). In principle, one could use computer simulation for this purpose, but this is expensive and not necessarily very precise if one has to adjust three or more potential parameters. However, an efficient and reliable alternative is the then~x~namic

perturbaticn theory of Weeks-chardler-Andersen

type extended to

anisotropic model molecules as su~]gested by Kohlar et al. (1979). In fact, this theory was used to determine the nCIJ potential parameters for some 20 spherical, linear (Fischer et al., 1984; Bohn et al., 1986), triangular (Lustig, 1986),

0378-3812/89/$03.50

© 1989 Elsevier Science Publishers B.V.

tetrahedral and octahedral it was

shown by molecular

derived

frc~ perturbation

(Lustig, 1987) molecules in the liquid. Very recently dynamics

(MD) simulations

that the 3CIJ potential

theory describes the thermodynamics

of real propane

correctly over the entire liquid range (Lustig and Steele, 1988). To see if this holds

true

represented

for other real

substances, we performed MD simulations

on ethane

as a 2CIJ molecule. Other recent MD studies on liquid ethane were

focussed on the density and the internal energy at zero pressure (Wojcik et al., 1982)

and

on the

specific

heat

(Fischer et al.,

1987).

To make

an overall

~ison

between the model and real ethane we present MD results here for the

pressul~,

internal

energy,

thermal

pressure

ooefficient,

specific

heat

and

adiabatic compressibility along four isochores over the entire liquid range.

MODEL AND MD ~ Q U E The methyl groups of ethane are considered to be two rigidly connected force centers. The interaction between a pair of molecules then consists of the sum of four site-site interactions which are of ~ - J o n e s

type. For a given center

of mass separation r and a given orientation of the two molecules 031 and 032 one thus has an angle-dependent potential 2 2 u(r031032) = ~E Z a=-i #=l

uij(r~8 ) , with

uIj(r~8 ) = 4E[(r--~a) 12 - (~a) 6]

where r~8 is the site-site distance, site size parameter.

The d i ~

6 is the potential well-depth and a is the

1 between the sites on the same molecule in

reduced units is L = i/a. There are currently two 2CiJ parameter sets for ethane in use. The first is as derived through MD by Wojcik et al. (1982) using densities and internal energies at zero pressure for the adjustment. dynamic

perturbation

coexistence

points

theory

on the

by

The seoond one was derived through thermo-

Fischer et al.

liquid

side.

Both

(1984) using

sets use the

two vapor-liquid same reduced bond

length:

Wojcik et al.:

L = 0.67, e/kB = 137.50 K, a = 3.506 A

Fischer et al.: L = 0.67, £/kB = 139.81 K, a = 3.512 A

Although these two sets are rather similar, in the parameters

may have

it was shown that small differences

rather large effects on the simulated properties

(Lustig and Steele, 1988). We have used the parameters of Fischer et al. in the sirmalations reported here. We ocmm~_nt on the other set in the discussion of the results. The MD sia~/lations were carried out by using constraint dynamics as described by Cicotti et al.

(1982). Here we use the same ccmputer program as that employed

for propane in an earlier paper (Lustig and Steele, 1988), where some details are also given. In essence, the atcmic equations of motions are solved in cartesian ooordinates

using

the

simple Verlet

algorithm.

Within this sdheme constraint

forces are taken into account to ensure that the atoms on a molecule do not fly apart. Periodic boundary conditions of the cubic simulation box and the minimum image criterion were applied to mimic an infinite system. The site-site~ potential was set to zero at a given site-site distance r c. The runs were performed at constant number of molecules N, volume V and total energy E (NVE). The time step for the numerical integration of the equations of motion was chosen such that the total energy per molecule ficant

figure.

previous

Systems

configurations

~ t u r e

(IE/N~I~O to i0) fluctuates only in the third signi-

of

108

and

256 molecules

and velocities

and were

were

usually

equilibrated

started

from

at a prechosen

for about 2000 steps, since each state point was simulated only once

we do not show error bars. However, an indication of the uncertainty of the MD points is given by their deviation frc~ smooth curves. From the figures in the next section it can be seen that this deviation is reasonably small. Details of the individual MD runs are given in the Appendix.

RESULTS Our MD results are compared to cubic interpolated values which in turn were picked

frcm an empirical

equation of state

(EOS) (Goodwin et al.,

stress this point because these semi-experimental

1982).

We

data also have some uncer-

tainty. ~he unreduced MD data points in Figures 1 to 5 were calculated using the potential parameters of Fischer et al. (1984). (Explicit values for the MD points are given in the Appendix). To cumpare our simulations with the results of other authors, values of the three highest densities

(~3

= 0.5472, 0.4681,

were studied that have also been used by Fischer et al.

0.3808)

(1987). This is parti-

cularly interesting because these authors give internal energies and their temperature derivatives that were obtained frc~ a different MD algorithm in which the separate translational and rotational equations of motion were solved using the

fifth/fourth-order

algorithm. negative

predictor-corrector

Their values than

tempez-atures,

curs due

by to

for the about

an

error

residual

ures/N£~0.07 in the

method

and

internal for

a

constant

energy

are

temperature

slightly

all three densities

long range correction.

This,

more

and all however,

cancels in the specific heat. Their corrected energies agree with ours within the error bars of the simulations

(Fischer et al.,

1988). We also note in passing

that the two different MD algorithms were checked against each other at the same state point. Pressure, energy and the mean squared force and torque are oc~mpletely consistent. Tne residual specific heat at constant volume was determined by Fischer

et

al.

(1987) by

fitting

the

internal

energies

along

isochores

to

straight lines and taking the derivative with respect to temperature. Our slopes are smaller than theirs by about 0.05R. This difference is certainly within the uncertainty of the calculations.

90.35 I

T/K 305.33 400

200 I

I

I

1

1

I

500 I

I

600 I

1

I

130

2.8-

120

2.61.6 ~

,l E

1.2p.o'3,

1.0 -

E

0.80.6W

-

60

- 50

P 4 0 M PQ

oo "

3O 2O

lil~II q

0.2-

-0.2 -

dl

.,/

£

0,4-

0-

70

d :/

14-

I0

J_....Jcp

trO I • 0.5

I 1.0

I 1.5

I 2,0

0

I 2.5 T ~/k e

I 3.0

I 3.5

I 4.0

I 4.5

-I0

FIG. i: Pressure versus temperature for fluid ethane. The line bounded by open circles (triple point, critical point) is the experimental vapor pressure curve. Dots: MD, lines: experiment.

Figure 1 shows that the pressure of the model coincides with the experimental data within the uncertainty of the simulation over the entire range of T,P (i00 to 600 K and 0 to 70 MPa). According

to

Figure

2 the residual

internal

energy

of the model

is more

negative than the EOS data by about 0.4 kJ/mol for the three highest densities at all temperatures. At the lowest density, which corresponds to about 1.5 times the critical density, the agreement between EOS and MD depends on temperature. At 370 K the sin~lation is about 0.15 kJ/mol more negative than the EOS data, whereas at 600 K the difference is about 0.4 kJ/mol. ~nis points at a difference in the specific heats as will be seen below. The ~ i n t y

for enthalpies fr~n the EOS

is given as 4 percent in the dense liquid and increases to i0 percent near the critical point (Goodwin et al., 1982).

From that we estimate the ~ i n t y

in

the residual internal energy to be 0.4 to 1.2 kJ/mol with the lower value for the

90.35

T/K 305.33 400

200

I

|

I

I

I

I

500 I

I

600 I

I

I

Cpo

p*=0.2600 / " l i ~ - ' ~ ' - - e e " ~ " - " ~ e

-8 -

p~=0.3808 ureS N---~- -~0

-I0

Ures

-12 kJ/mol p~= 0.4681 o~ee

-12

-14 ,~= 0.5472

-16

-14 tr o eea~ • -J6

I

I

I

0.5

LO

1.5

I z.0

I

I

I

2.5

3.0

3.5

1

4.0

I 4.5

-18

T

E/k~ FIG. 2: Residual internal energy versus ~ t u r e . The open circles denote experimental points on the ooexistence curve. Dots: MD, lines: experiment. highest and the higher value for the lowest density so that all values of the model agree with the semi-experimental data within the combined uncertainty of EOS and MD. Figure 3 shows

a comparison of the thermal pressure coefficients~v=(~p~T)v.

The dots are direct MD outputs obtained frc~ the fluctuation formulae (C~leung, 1977). The arrows are obtained by fitting the MD points in Figure 1 to straight lines and taking the temperature derivative. Both methods of determining Y v are completely consistent. At the two highest densities the results from the fluctuations are scattered around the experimental values. We believe that this scatter can be removed by running the simulations longer. The runs for the two lowest densities were long enough for the running averages of ~ v to converge. We conclude that there is excellent agreement between the model and the experimental thermal pressure coefficients over the entire liquid range. Figure 4 shows the residual isochoric specific heat cvreS = ~ures~T)v. Dots are again direct MD outputs frc~ the fluctuation formulae (C~eung, 1977), whereas the arzrx~ give the slopes of straight lines fitted through the isochoric MD points of Figure 2. Again, beth calculations yield consistent results. As before, we

believe that

the scatter of the results frc~ the fluctuations can be removed

T/K

200

90.35 I

I

I0

400

305.33 I

I

I

I

500

I

I

I

600 I

I

'-35

fro

5.0

"%'--.• }P*

2.5 0.5472 -

6

2.0

~v

yV o'3

M Po/K -i.5

kB

@

4

p*= 0.4681

oe~e ~ I

-l.O

C

w A

:

~

4---- p*-- 0.3808 -0.5 p* = 0.2600

C~

Cp o I

I

0.5

1.0

I

-0

I

1.5

2.0

I

I 2.5

3.0

I

I

I

3.5

40

45

T ~/k B

FIG. 3: Thermal pressure coefficient versus temperature. The open circles denote ~imental points on the coexistence curve. Dots: MD (fluctuation theory), arrows: MD (mean slope of p(T,p=const)), lines: experiment. by longer runs. qhe ovea-all ~ i n t y

of the total specific heat from the EOS

is given as 4 percent (Goodwin et al., 1982). The upper limit for the uncertainty of the ideal part of ~

(or Cv) is given as 2 percent above 325 K. Using these

values we estimate the ~ i n t y all the densities studied.

of the residual part of Cv to be about 0.2R at

Thus we conclude that the residual specific heat of

the model also reproduoes the semi-experimental data within the combined accuracy of MD and EOS for densities down to at least i. 5 times the critical density. In Figure 5 we show a quantity ~s p that is closely related to the adiabatic coefficient of bulk compressibility

~s - l = p As

(~) s = ~ T -1 + T~v2 pc~

in the study

of propane

(Lustig and Steele,

expression by using for C v on the r.h.s. This pseudo

coefficient

is directly

1988),

(cvreS+2.5R)

accessible

~s p differs

from this

instead of (cvreS+cvid).

through a fluctuation

(C~eung, 1977). We report these results for the sake of

formula

convpleteness only. It is

difficult to assign an uncertainty to ~s p from the SOS so that it is difficult to cc~mnent on the cor~0arison to MD. However, the feature of Figure 5 is essen-

90.35 2.4

T/K 30533 400

200

tr~

l

t

I

I

I

I

500 I

600

I

I

I

2.22.01.8-

"p.V,=-0.5472

1.6

cres N ks

1,4

-?_2

1.2

1.0 0.8 0.6 0.4 0.2 0

I 1.0

I

0.5

I 1.5

I 2.0

I 2.5

p~= 0.2600 I I I 3.0 3.5 4.0

I 4.5

T elk s

FIG. 4: Residual specific heat versus tesi0erature. The open circles denote experimental points on the coexistence curve. Dots: MD (fluctuation theory), a~zuws: MD (mean slope of ures(T,p=oonst)), lines: experiment. tially the same as we found for propane

(Lustig and Steele,

1988).

The MD re-

sults are qualitatively in agreement with the EOS, which means that ~ ses with

decreasing

density

and

decreases

with

~ i n g

increa-

temperature

along

isochores. However, quantitatively the MD points are smaller than the EOS data by a factor of about 1/3 at all isochores considered.

DISCUSSION We have shown earlier (Lustig and Steele, 1988) that it is possible to obtain agreement between simulation and experiment

in the p,T-diagram

fpr propane by

interpolating linearly between two parameter sets 6 and a that are not too far apart. In the case of ethane a linear interpolation between the parameter sets of Fischer et al.

(1984)

(used here)

and Wojcik et al.

(1982) should result in

essentially the same pressure isochores. Due to the smaller E of the latter set the corresponding MD internal energies are closer to the EOS data, which in fact were used for the adjustment.

Frc~ our results on propane we conclude that all

the other properties considered here are not significantly affected by changing between either of the two sets. ~hus,

in the scope of our estimations for the

uncertainty of the EOS data, any linear interpolation between these two sets

90.:35

0.40

T/K 305.:33 400

200

I

I

I

I

500

600

I

I

9 8

0.35

-7

0.30

6 0.25

p~= 0.2600

P

- 5 ~r3

0.20

Bs (I0 3 M Pa)

4

0.15 -3

0.10 -

p~=0.3808t I

~ "

--2

-

0.05 0

-I I

0.5

L-~I

•

1.0

I

I

1.5

2.0

I

2.5 T

I

3.0

I

3.5

I

4.0

I

4.5

0

Elk a

FIG. 5: "Adiabatic ooefficient of bulk ocmvpressibility" versus ~ t u r e . The experimental lines are recalculated frc~ EOS data (see text). The open circle~ denote experimental points on the coexistence curve. Dots: MD (fluctuation theory), lines: experiment. should be equally suitable for the calculation of the thermodynamic properties in the dense liquid. However, we can pick an optimum set from this range by fitting the second virial coefficients. Bohn et al. (1986) have shown that the parameter set of Fischer et al. (1984) yields theoretical second virial coefficients that are within

the

experimental

uncertainty throughout the measured temperature

range. Now, a reduction of £/kB below about 139 K moves the theoretical second virial coefficients outside the experimental ttncertainty for all temperatures, towards less negative values. virial

coefficient

parameter set of

B/c:n3mol-I

For example, at 200 K the e x p e r ~ t a l is

-410~-i0, the

theoretical

ones

second with

the

Fischer et al. are -400 and with the set of Wojcik et al. about

-384. ~nus, E/kB =139.81 K can be considered as the best value for quantitative studies. To ~ i z e

the results, the t h ~ c

quantities of the 2CIJ model used

here are in excellent agreement with the semi-experimental data of real ethane, including the second virial coefficient (Bohn et al., 1986), but excluding ~ . q~lis suggests that three-body interactions have no significant effect upon the thermodynamics of this dense liquid, qhis statement is backed up by our detailed analysis of three-body effects in liquid propane (Lustig and Steele, 1988).

CONCLUSION MD calculations on a 2CIJ model were carried out to simulate ethane over the entire liquid region. All the thermodynamic quantities considered (except ~ ) are reproduced within the combined uncertainty of MD and semi-experimental data. We note that the 2CLJ model used here was not derived through MD but by a quite different approach based on a thermodynamic perturbation theory of Weeks-Chandler-Andersen type

(Fischer et al.,

1984).

In this

and a previous paper on

propane (Lustig and Steele, 1988), we have shown that this perturbation theory gives the correct nCIJ model for the two molecules studied. Our results demonstrate once more that the nCIJ potential can be an excellent effective pair potential in fluid phases of silmple molecules. As in the case of propane, the effects of three-body interactions on the thermodynamics of dense liquid ethane must be very small.

A C K N ~ This work was supported by the National Science Foundation Grant CHE 8619063. Ra~NCES Bohn, M., Lustig, R. and Fischer, J., 1986. Description of polyatomic real substances by two-center lennard-Jones model fluids. Fluid Phase Equilibria, 25: 251-262. Cheung, P.S.Y., 1977. On the calculation of specific heats, thermal pressure coefficients and ccmloressibilities in molecular dynamics simulations. Molec. Phys., 33: 519-526. Cicotti, G., Ferrario, M. and Ryckaert, J.-P., 1982. Molecular dynamics of rigid systems in cartesian coordinates. A general formulation. Molec. Phys., 47: 1253-1264. Fischer, J., Lustig, R., Breitenfelder-Manske, H. and Leam/r~, W., 1984. Influence of intermolecular potential parameters on orthobaric properties of fluids consisting of spherical and linear molecules. Molec. Phys., 52: 485-497. Fischer, J., Saager, B., Bohn, M., Oelschlager, H. and Haile, J.M., 1987. Specific heat of s i d l e liquids. Molec. Phys., 62: 1175-1185. Fischer, J., Saager, B., Bohn, M., Oelsclllager, H. and Haile, J.M., 1988. Private ~m,~nication. Goodwin, R.D., Roder, H.M. and Straty, G.C., 1982. Thermophysical properties of ethane from 90 to 600 K at pressures to 700 bar. Nat. Bur. Stand. Technical Note 684, Boulder. Kohler, F., Quirke, N. and Perram, J.W., 1979. Perturbation theory with a hard dlanbbell reference system. J. C~e~. Phys., 71: 4128-4131. iustig, R., 1986. A thermodynamic perturbation theory for non-linear multicentre ~ - J o n e s molecules with an anisotropic r e f ~ system. MOlec. Phys., 59: 173-194. Lustig, R., 1987. Application of thermodynamic perturbation theory to multicentre L ~ - J o n e s molecules. Results for CF4, CCI 4 , neo-C5H!2 and SF 6 as tetrahedral and octahedral models. Fluid Phase Equilibria, 32: 117-137. Lustig, R. and Steele, W.A., 1988. On the thermodynamics of liquid propane. A molecular dynamics study. Molec. Phys., 65: 475-486. Wojcik, M., Gubbins, K.E. and Powles, J.G., 1982. ~he thermodynamics of s~tric two-centre ~ - J o n e s liquids. Molec. Phys., 45: 1209-1225.

10

% ~ 0 0 0 0

0

0

0

0

~

0

0 0 ~ 0 0

,c: ~

0

~

0

0 0 0 0 0

0

0

.? 0 0 0 o 0 . . . . .

~

0

0

0

° . . .

0

0

0

0

0

0

0

, ° . , , ,

°

0

0

°

°

0 .

°

:Z ~n

0

0

~

° ° ° . . .

.q 0 ~--I

:>t

,-4

~

0 0 ~ 0 0 , , , , .

0

I

l

l

l

l

-M -.4 ~ ~

0

~

0

0

~

~

0

~

~

0

~

0 ,-I ~

~

0

~

~n 0

~ ~ . . .

I

I

I

I

I

I

~'~ :>

~

0

o

~

dSdSd

. ° 0 . ° °

0~

. , . . o ,

~

~

0

~ 0

o_,

0

°

,

0 ~ 0 0 ~

I I I I I I I I I

!

0 0 %0 C~

.

~

I I I I I I

° 0 ° ° °

°

, ° ° 0 ° °

ZZZZZ I I I I I

0

°

. . ° ,

CO 0 O0 CG

~--I O0 ~0 0

c~ L~

c~