- Email: [email protected]
Prediction of thermodynamic properties for fluid oxygen with MD simulations
103
Fluid Phase Equilibria, 66 (1991) 103-111 Elsevier Science Publishers B.V., Amsterdam
Prediction of thermodynamic with MD simulations
properties
for fluid oxygen
Berthold Saager and Johann Fischer ’ Institut fiir Thermo- und Fluidsynamik, Ruhr-Universitiit, D-4630, Bochum 1 (Germany) (Received September 251990;
accepted in final form January 25, 1991)
ABSTRACT Saager, B. and Fischer, J. (1991). Prediction of thermodynamic with MD sirmdations. Fhdd Phase Equilibria, 66: 103-111.
properties for fluid oxygen
Pressures, internal energies and specific heats of oxygen obtained from molecular dynamics simulations on a CYBER 205 are reported. The intermolecular interaction is mode&d by a two-centre Leonard-Jones potential The three potential parameters were determined previously with a WCA-type perturbation theory by fitting to the experimental vapour-liquid coexistence curve. The agreement of the calculated with the experimental pressures is on a gas isochore within f 0.2 MPa and on two liquid isochores within f 1.7 MPa. The agreement of the simulated residual internal energies with those of the Schmidt-Wagner equation is mostly within about 1% except in the gas at a near-critical temperature. One shortcoming of the present potential model is that at high densities of p/p, = 3 the model fluid seems to become thermodynamically unstable owing to freezing problems.
INTRODUCTION
Recently, several papers reported thermodynamic properties of pure substances for the whole fluid region as obtained from computer simulations. Examples are the works for methane (Luckas et al., 1990; Saager and Fischer, 1990), ethane (Lustig et al., 1989), propane (Lustig and Steele, 1988), carbon dioxide (Luckas and Lucas, 1989) and R152a (Vega et al., 1989). The crucial point in these projects is the determination of an effective intermolecular potential model. A systematic procedure is the use of zeropressure densities and internal energies of the liquid (McDonald and Singer, ’ Author to whom correspondence 0378-3812/91/$03.50
should be addressed.
0 1991 Elsevier Science Publishers B.V. AlI rights reserved
104
1972; Singer et al., 1977). A next step along this line was to use saturated liquid densities and vapour pressures as determined from WCA-type perturbation theory (Fischer et al., 1984; Bohn et al., 1986). We should mention that in the further development we aim to overcome the limitation of perturbation theory to non-polar fluids. For that purpose we suggested recently the NpT + test particle method for the determination of phase equilibria (Miiller and Fischer, 1990) by simulations. Among the several substances for which effective intermolecular pair potentials have been determined by WCA-type perturbation theory (Fischer et al., 1984; Bohn et al., 1986; Lustig, 1986, 1987), oxygen is of special interest for simulation studies for several reasons. For testing the potential model it is important that a rather accurate equation of state is available (Schmidt and Wagner, 1985), which allows the predicted internal energies to be compared also. From a thermodynamic point of view it is remarkable that oxygen has a rather long liquid range (TJT, = 0.35) compared to the rare gases or nitrogen. Regarding the molecular modelling it is interesting that in the two-centre Lennard-Jones (2CLJ) potential, for which the parameters were obtained by fitting to the experimental vapour-liquid coexistence curve (Bohn et al., 1986), the distance between the interaction sites is only 0.71 A compared to the distance between the nuclei of 1.20 A (Gmelin, 1958). In the next section, the potential model and the details of the simulations are presented and in the final section the results for the pressure, internal energy and specific heat are compared with experimental or correlated values.
POTENTIAL MODEL AND DETAILS OF THE SIMULATIONS
A 2CLJ potential consists of two interaction sites on each molecule being a distance I apart. Between the interaction sites on different molecules, Lennard-Jones potentials with parameters E and CI are assumed. It is convenient to describe the distance I by the elongation L = l/a. The strategy to determine an appropriate set of parameters for some real fluids by using WCA-type perturbation theory (Fischer, 1980) relies on the fact that the steepness of the vapour pressure curve strongly depends on the elongation L (Fischer et al., 1984). One may start with a guessed elongation L and calculate the vapour-liquid coexistence curve for that model (Fischer et al., 1984). Then, one selects an intermediate temperature and fits the parameters c and u to the experimental vapour pressure and saturated liquid density at that temperature. In the next step, one selects a second temperature and checks whether the vapour pressure of the model agrees
105
with the experimental value there. The procedure is repeated for different values of L until the vapour pressure of the second temperature is well predicted. The parameters found for oxygen in this way (Bohn et al., 1986) are L = 0.22, c/k = 38.003 K and u = 3.2104 A. We should mention that with the 2CLJ model thus obtained for oxygen we have, besides the coexistence curve, also already calculated the second virial coefficient (Bohn et al., 1986), isochoric heat capacities in the liquid by perturbation theory (Fischer et al., 1987) and gas PVT data by simulations (Saager et al., 1990). The molecular dynamics (MD) simulations were performed in an NVT ensemble with a code that is originally due to Haile (1984) and was vectorized for the CYBER 205. The translational motions are solved in a fifth-order predictor-corrector algorithm. The rotational motions are treated by the method of quatemions using a fourth-order predictor-corrector algorithm. Periodic boundary conditions and the minimum image criterion were applied. All runs were performed with 256 particles and the cut-off radius was taken to be half of the box length. Long-range corrections were taken into account by a method due to Neumann (1987), which is the same as the later developed method of Lustig (1988). The temperature was kept constant by momentum scaling; hence, strictly speaking, we are not working in a canonical but in an isokinetic ensemble. This difference is important for the calculation of the heat capacity, which in the isokinetic ensemble cannot be obtained from the energy fluctuation formula; actually, it was computed by differentiation of the internal energies. The time step for the integration was chosen to be 0.0015 in units of aJm/E. The system was started either from an f.c.c. lattice or a previous configuration followed by an equilibration period of 4000 time steps. The production runs were extended over 15 000 to 20000 time steps. The quantities calculated directly are the pressure and the residual part of the internal energy. From plots of the running averages of these quantities, their statistical uncertainties were estimated; these will be given, together with the results, in the next section.
RESULTS AND DISCUSSION
The state points for the simulations were chosen such that the resulting pressures and isochoric heat capacities could be compared directly with experimental values. For the comparison of the pressure, simulations were performed close to one isochore in the gas, close to two isochores in the liquid and at one supercritical isotherm, which covers nearly the whole fluid region except at very high densities. The directly obtained pressures and
106 TABLE 1 Results from MD simulations: residual internal energy, pressure and their estimated uncertainties for oxygen as 2CLJ fluid (T* = kT/r, p* = PIJ~, u* = U/NC, p* = po3/c)
TV
P*
7.894114
4.031208 5.262742 6.578428 7.894114 3.368155 3.841802 4.341762 4.736468 4.999605 5.262742 5.525880 2.526116 2.789253 3.052391 3.315528
0.206180 0.241655 0.328227 0.368956 0.387309 0.415385 0.177749 0.190305 0.187085 0.184163 0.575423 0.567013 0.560236 0.570771 0.569974 0.569256 0.568579 0.697468 0.689994 0.683964 0.691963
3.578665
0.690787
3.338992 3.388992 3.488992 3.538992 4.095800 4.145800 4.245800 4.295800 2.473613 2.523613 2.623613 2.673613 2.716357 2.766357 2.866357 2.916357 1.423485 1.623485 1.556896 1.606896 1.706896 1.756896
0.572265
0.572265
0.695828
0.695828
0.809987
u* -
4.773 5.534 7.508 8.433 8.841 9.457 5.089 4.971 4,521 4.263 15.116 14.615 14.181 14.234 14.082 13.916 13.795 18.956 18.515 18.129 18.094 18.084 17.816 17.836 15.071 15.035 14.977 14.940 14.613 14.613 14.537 14.515 18.977 18.921 18.825 18.777 18.746 18.689 18.590 18.557 23.167 22.787 22.893 22.844 22.651 22.580
P*
AU*
AP*
Prod. steps
1.56 1.87 2.87 3.48 3.77 4.37 0.31 0.67 1.05 1.37 0.07 1.06 2.11 3.29 3.95 4.50 5.00 0.11 1.07 1.99 3.53 3.53 4.63 4.60 - 0.08 0.06 0.27 0.44 1.81 1.87 2.11 2.26 - 0.28 0.03 0.45 0.78 0.91 1.19 1.68 1.86 -0.86 1.26 0.67 0.93 1.94 2.34
0.025 0.025 0.025 0.025 0.025 0.025 0.050 0.040 0.040 0.040 0.020 0.020 0.020 0.020 0.020 0.020 0.020 0.015 0.015 0.015 0.015 0.015 0.015 0.015 0.015 0.015 0.015 0.015 0.020 0.020 0.020 0.020 0.015 0.015 0.015 0.015 0.020 0.020 0.020 0.020 0.035 0.035 0.035 0.035 0.035 0.035
0.015 0.02 0.03 0.03 0.03 0.03 0.01 0.01 0.01 0.01 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.07 0.07 0.07 0.07 0.10 0.10 0.10 0.10 0.17 0.17 0.17 0.17 0.17 0.17
15000 15000 20000 20000 20000 20000 15000 15000 15000 15000 20000 20000 20000 15000 15000 15000 15000 15000 20000 20000 20000 20000 20000 20000 15c!Oo 2oooo 20000 15000 20000 20000 20000 20000 15000 15000 15000 15000 20000 20000 20000 20000 15000 15000 15000 20000 20000 15000
107 TABLE 2 Comparison of pressures from MD simulations with experimental results for oxygen. The simulated pressures are given together with their estimated statistical uncertainties PMD (MPa)
TK)
Trnoll-‘)
PCXp @Pa)
153.198 200.0 250.0 300.0 128.0 146.0 165.0 180.0 190.0 200.0 210.0 96.0 106.0 116.0 126.0
8.9203 9.5504 9.3888 9.2422 28.8775 28.455 28.1153 28.644 28.604 28.568 28.534 35.0023 34.5272 34.3246 34.726
4.7775 10.7398 16.5349 21.9884 2.3713 18.4948 34.4943 52.8118 61.9364 70.8860 79.6228 3.5267 18.2560 32.4298 56.3733
136.0
34.667
72.9740 d
300.0
10.3471 12.1274 16.472 18.516 19.437 20.846
24.8740 29.9184 45.4145 55.2581 60.4637 69.4934
b ’ ’ ’ ’ = = d d d d = = = d
= ’ d d d d
4.9f0.16 10.6 f 0.16 16.7 f 0.16 21.8 f0.16 l.Of 1.0 16.9 f 1.0 33.5 f 1.0 52.2 f 1.0 62.6 f 1.0 71.3 f 1.0 79.3 f 1.0 1.8 f 1.1 17.0 f 1.1 31.6 f 1.1 56.0& 1.1 55.9f 1.1 73.5 f 1.1 73.0 f 1.1 24.8 f 0.24 29.7 f 0.32 45.6 f 0.50 55.2kO.50 59.9 f 0.50 69.2 f 0.50
- 25.41 1.75 - 1.03 0.80 1.11 1.00 0.52 0.24 - 0.24 -0.14 0.10 0.37 0.27 0.17 0.06 0.08 - 0.08 -0.00 0.27 0.60 - 0.26 0.06 0.49 0.20
a The relative difference in the density that is associated with the difference between the simulated and the experimental pressure on the basis of the Schmidt-Wagner (1985) equation of state. b From Pentermann and Wagner (1978). ’ From Weber (1970). d From Weber (1977).
residual internal energies, together with their statistical uncertainties, are given in Table 1 in reduced units for all stable simulation runs. Table 2 shows the results for the simulated pressures together with their statistical uncertainties in real units in comparison with the experimental data. For two state points, the results of two independent simulation runs are given in order to demonstrate their reproducibility. We note that generally the simulation results agree within their statistical uncertainties with the experimental data. Only on the liquid isochores at the two respectively lowest temperatures is the disagreement slightly larger. It seemed
108 TABLE 3 Comparison of residual internal energies from MD simulations with results from the Schmidt-Wagner (1985) equation of state for oxygen. (The simulated energies are given, together with their estimated statistical uncertainties) ;)mol J-‘)
(J mol-‘)
G..s
Ui%D
153.198 200.0 250.0 300.0 128.0 146.0 165.0 180.0 190.0 200.0 210.0 96.0 106.0 116.0 126.0
8.9203 9.5504 9.3888 9.2422 28.8775 28.4554 28.1153 28.644 28.604 28.568 28.534 35.0023 34.5272 34.3246 34.726
- 1801 - 1605 - 1459 - 1360 - 4781 -4615 -4464 -4466 - 4410 -4356 - 4304 - 5981 - 5839 - 5711 - 5687
136.0
34.667
- 5598
300.0
10.3471 12.1274 16.472 18.516 19.437 20.846
-
-1608f16 -1571*13 -1428f13 -1347f13 -4776& -4618* -44812 -4498& -4450f -4397k -4359f - 5990f -5850& -5728f -5717* - 5714f -5629f -5636f -1508* -1749+ -2372+ -2665k -2794f -2988*
1514 1761 2357 2636 2760 2950
(J mol-‘)
6.3 6.3 6.3 6.3 6.3 6.3 6.3 4.7 4.7 4.7 4.7 4.7 4.7 4.7 8 8 8 8 8 8
10.7 2.1 2.1 0.9 0.1 -0.1 - 0.4 -0.7 -0.9 -0.9 -1.3 -0.2 -0.2 -0.3 -0.5 - 0.5 -0.6 -0.7 0.4 0.7 -0.6 -1.1 -1.2 -1.3
a The relative difference between the simulated and the equation-of-state energies.
interesting also to express the pressure differences as differences of density, which was made on the basis of an empirical equation of state (Schmidt and Wagner, 1985). The relative density differences are also given in Table 2. We learn that, especially at the highest density, the differences are remarkably small. The largest difference occurs for a gas at a temperature close to the critical, which is simply an effect of the nearly zero slope of the isotherm there. Table 3 shows the results for the simulated residual internal energies together with their statistical uncertainties in real units in comparison with the values from the empirical equation of state of Schmidt and Wagner (1985); the state points are the same as those for which the pressures were compared. The table also contains the relative differences between the
109 TABLE 4 Comparison of isochoric heat capacities from MD simulations with experimental results (Goodwin and Weber, 1969) for oxygen C”
(-)MD
;)mol 1-l) 130.692 159.453 97.805 107.030
R
28.719 S 34.920
3.19 3.12 3.51 3.46
3.15 3.05 3.49 3.45
simulated and the equation of state values. We observe that, especially in the liquid, the agreement is very good. On the gas isochore the agreement is somewhat worse, with the largest deviation at 153 K, which is close to the critical temperature. There may be several reasons for the discrepancy. One could be that, owing to the finite number of particles (N = 256), the critical temperature of the simulated system would be higher than that of an infinite system. Then, the state point at p = 8.9 mol 1-l and T= 153.2 K could be already within the thermodynamically unstable region, where it is known that the internal energy obtained from simulations is at least particle number dependent (Heinbuch and Fischer, 1987). Further investigations are necessary in that respect and are planned, together with studies on the dependence of the coexistence curve on the particle number close to the critical point. We also calculated isochoric heat capacities by differentiation of the smoothed simulated isochoric energies as described previously (Fischer et al., 1987). The results are compared in Table 4 with experimental data (Goodwin and Weber, 1969) and the agreement is found to be within the experimental accuracy. As has been mentioned already, oxygen has a large liquid range. The data presented so far extend up to a density of p 2: 35 mol I-‘, which is at P/PC = 2.7. The triple-point density of oxygen, however, is still higher, at P/PC = 3.0. Performing simulations at that density, however, meant that the running averages mostly drifted away if the runs were made sufficiently long. Sometimes they remained stable in a liquid state for a certain time and the thermodynamic data then obtained are still quite reasonable. As an example, the isochoric heat capacity at p = 40.649 mol l-1 and T = 62.967 K was obtained from simulations as c,/R = 4.16 compared to the experimental value of 4.22. This drift of the running averages may indicate that we are already in the liquid-solid two-phase region of the model fluid. This conjecture is confirmed by the fact that the triple-point density of the LJ
110
fluid is at p/p, = 2.7 and the 2CLJ model for oxygen used here is due to its small elongation (L = 0.22) not being very acentric. A conclusion of the present work is that the simple 2CLJ model used here is able to predict the thermodynamic properties of oxygen over a large fluid region. One shortcoming is the fact that the model fluid seems to become thermodynamically unstable owing to the fact that it freezes at densities below the triple-point density of real oxygen. A proper description of the liquid-solid transition would require either a definitively more complicated effective pair potential or the inclusion of three- and higher body forces.
ACKNOWLEDGEMENTS
The authors would like to thank Professor W. Wagner (Bochum) for suggestions on this work and Dr. J.F. Ely (Boulder) for sending the report (NBSIR 77-865). The project was financially supported by the Deutsche Forschungsgemeinschaft under the priority programme “Neue Arbeitsmedien in der Energie- und Verfahrenstechnik”, AZ Fi 287/S, gefordert. Parts of this work were presented at the VDI-~e~od~~ ~olloqium 1989 in Trier.
REFERENCES Bohn, M., Lust&, R. and Fischer, J., 1986. Description of polyatomic real substances by two-center Lennard-Jones model fluids. Fluid Phase Equilibria, 25: 251-262. Fischer, J., 1980. Perturbation theory for the free energy of two-center Lennard-Jones liquids, J. Chem. Phys., 72: 5371-5377. Fischer, J., Lustig, R., Breitenfelder-Manske, H. and Lemming, W., 1984. Influence of intermol~ular potential parameters on orthobaric properties of fluids consisting of spherical and linear molecules. Mol. Phys., 52: 485-487. Fischer, J., Saager, B., Bohn, M. and Haile, J.M., 1987. Specific heat of simple liquids. Mol. Phys., 62: 1175-1185. Gmelin 1958. Handbuch der Anorganischen Chemie, Sauerstoff, System No. 3, Lfg 3. Verlag Chemie, Weinheim/Bergstrasse, 356 pp. Goodwin, R.D. and Weber, L.A., 1969. Specific heats c, of fluid oxygen from the triple point to 300 K at pressures to 350 atmospheres. J. Res. Nat. Bur. Stand., 73A: 15-24. Haile, J.M., 1984. Private communication. Heinbuch, U. and Fischer, J., 1987. On the application of Widom’s test particle method to homogenous and inhomogenous fluids. Mol. Sim., 1: 109-120. Luckas, M. and Lucas, K., 1989. Thermodynamic properties of fluid carbon dioxide from the SSR-MPA potential. Fluid Phase Equilibria, 45: 7-23. Luckas, M., Ripke, M. and Lucas, K., 1990. Prediction of the ~ermophysical properties of fluid methane from the SSR-MPA potential. Fluid Phase Eq~~b~a, 58: 35-50.
111
Lustig, R., 1986. A thermodynamic perturbation theory for non-linear multicentre LermardJones molecules with an anisotropic reference system. Mol. Phys., 59: 173-184. Lustig, R., 1987. Application of thermodynamic perturbation theory to multicentre LennardJones molecules. Results for CF,, Ccl,, neo-C,H,, and SF, as tetrahedral and octahedral models. Fluid Phase Equilibria, 32: 117-137. Lustig, R., 1988. Angle-average for the powers of the distance between two separated vectors. Mol. Phys., 65: 175-179. Lustig, R. and Steele, W.A., 1988. On the thermodynamics of liquid propane. A molecular dynamics study. Mol. Phys., 65: 475-486. Lustig, R., Toro-Labbt, A. and Steele, W.A., 1989. A molecular dynamics study of the thermodynamics of liquid ethane. Fluid Phase Equilibria, 48: l-10. McDonald, I.R. and Singer, K., 1972. An equation of state for simple liquids, Mol. Phys., 23: 29-40. Moller, D. and Fischer, J., 1990. Vapour liquid equilibrium of a pure fluid from test particle method in combination with NpT molecular dynamics simulations. Mol. Phys., 69: 463-473. Neumann, M., 1987. Private communication. Pentermann, W. and Wagner, W., 1978. New pressure-density-temperature measurements and new rational equations for the saturated liquid and vapour densities of oxygen. J. Chem. Thermodyn., 10: 1161-1172. Saager, B. and Fischer, J., 1990. Predictive power of effective intermolecular pair potentials: MD simulation results for methane up to 1000 MPa. Fluid Phase Equilibria, 57: 35-46. Saager, B., Lotfi, A., Bohn, M., Nguyen Van Nhu and Fischer, J., 1990. Prediction of gas PVT data with effective intermolecular potentials using the Haar-Shenker-Kohler equation and computer simulations. Fluid Phase Equilibria, 54: 237-246. Schmidt, R. and Wagner, W., 1985. A new form of the equation of state for pure substances and its application to oxygen. Fluid Phase Equilibria, 19: 175-200. Singer, K., Taylor, A. and Singer, J.V.L., 1977. Thermodynamic and structural properties of liquids modelled by ‘2-Lermard-Jones centres’ pair potentials. Mol. Phys., 33: 1757-1795. Vega, C., Saager, B. and Fischer, J., 1989. Molecular dynamics studies for the new refrigerant R152a with simple model potentials. Mol. Phys., 68: 1079-1093. Weber, L.A., 1970. PVT, thermodynamic and related properties of oxygen from the triple point to 300 K at pressures to 33 MN/m’. J. Res. Nat. Bur. Stand., 74A: 93-129. Weber, L.A., 1977. Thermodynamic and related properties of oxygen from the triple point to 300 K at pressures to 1000 bar. NASA Ref. Pub. 1011, NBSIR 77-865.
Fluid Phase Equilibria, 66 (1991) 103-111 Elsevier Science Publishers B.V., Amsterdam
Prediction of thermodynamic with MD simulations
properties
for fluid oxygen
Berthold Saager and Johann Fischer ’ Institut fiir Thermo- und Fluidsynamik, Ruhr-Universitiit, D-4630, Bochum 1 (Germany) (Received September 251990;
accepted in final form January 25, 1991)
ABSTRACT Saager, B. and Fischer, J. (1991). Prediction of thermodynamic with MD sirmdations. Fhdd Phase Equilibria, 66: 103-111.
properties for fluid oxygen
Pressures, internal energies and specific heats of oxygen obtained from molecular dynamics simulations on a CYBER 205 are reported. The intermolecular interaction is mode&d by a two-centre Leonard-Jones potential The three potential parameters were determined previously with a WCA-type perturbation theory by fitting to the experimental vapour-liquid coexistence curve. The agreement of the calculated with the experimental pressures is on a gas isochore within f 0.2 MPa and on two liquid isochores within f 1.7 MPa. The agreement of the simulated residual internal energies with those of the Schmidt-Wagner equation is mostly within about 1% except in the gas at a near-critical temperature. One shortcoming of the present potential model is that at high densities of p/p, = 3 the model fluid seems to become thermodynamically unstable owing to freezing problems.
INTRODUCTION
Recently, several papers reported thermodynamic properties of pure substances for the whole fluid region as obtained from computer simulations. Examples are the works for methane (Luckas et al., 1990; Saager and Fischer, 1990), ethane (Lustig et al., 1989), propane (Lustig and Steele, 1988), carbon dioxide (Luckas and Lucas, 1989) and R152a (Vega et al., 1989). The crucial point in these projects is the determination of an effective intermolecular potential model. A systematic procedure is the use of zeropressure densities and internal energies of the liquid (McDonald and Singer, ’ Author to whom correspondence 0378-3812/91/$03.50
should be addressed.
0 1991 Elsevier Science Publishers B.V. AlI rights reserved
104
1972; Singer et al., 1977). A next step along this line was to use saturated liquid densities and vapour pressures as determined from WCA-type perturbation theory (Fischer et al., 1984; Bohn et al., 1986). We should mention that in the further development we aim to overcome the limitation of perturbation theory to non-polar fluids. For that purpose we suggested recently the NpT + test particle method for the determination of phase equilibria (Miiller and Fischer, 1990) by simulations. Among the several substances for which effective intermolecular pair potentials have been determined by WCA-type perturbation theory (Fischer et al., 1984; Bohn et al., 1986; Lustig, 1986, 1987), oxygen is of special interest for simulation studies for several reasons. For testing the potential model it is important that a rather accurate equation of state is available (Schmidt and Wagner, 1985), which allows the predicted internal energies to be compared also. From a thermodynamic point of view it is remarkable that oxygen has a rather long liquid range (TJT, = 0.35) compared to the rare gases or nitrogen. Regarding the molecular modelling it is interesting that in the two-centre Lennard-Jones (2CLJ) potential, for which the parameters were obtained by fitting to the experimental vapour-liquid coexistence curve (Bohn et al., 1986), the distance between the interaction sites is only 0.71 A compared to the distance between the nuclei of 1.20 A (Gmelin, 1958). In the next section, the potential model and the details of the simulations are presented and in the final section the results for the pressure, internal energy and specific heat are compared with experimental or correlated values.
POTENTIAL MODEL AND DETAILS OF THE SIMULATIONS
A 2CLJ potential consists of two interaction sites on each molecule being a distance I apart. Between the interaction sites on different molecules, Lennard-Jones potentials with parameters E and CI are assumed. It is convenient to describe the distance I by the elongation L = l/a. The strategy to determine an appropriate set of parameters for some real fluids by using WCA-type perturbation theory (Fischer, 1980) relies on the fact that the steepness of the vapour pressure curve strongly depends on the elongation L (Fischer et al., 1984). One may start with a guessed elongation L and calculate the vapour-liquid coexistence curve for that model (Fischer et al., 1984). Then, one selects an intermediate temperature and fits the parameters c and u to the experimental vapour pressure and saturated liquid density at that temperature. In the next step, one selects a second temperature and checks whether the vapour pressure of the model agrees
105
with the experimental value there. The procedure is repeated for different values of L until the vapour pressure of the second temperature is well predicted. The parameters found for oxygen in this way (Bohn et al., 1986) are L = 0.22, c/k = 38.003 K and u = 3.2104 A. We should mention that with the 2CLJ model thus obtained for oxygen we have, besides the coexistence curve, also already calculated the second virial coefficient (Bohn et al., 1986), isochoric heat capacities in the liquid by perturbation theory (Fischer et al., 1987) and gas PVT data by simulations (Saager et al., 1990). The molecular dynamics (MD) simulations were performed in an NVT ensemble with a code that is originally due to Haile (1984) and was vectorized for the CYBER 205. The translational motions are solved in a fifth-order predictor-corrector algorithm. The rotational motions are treated by the method of quatemions using a fourth-order predictor-corrector algorithm. Periodic boundary conditions and the minimum image criterion were applied. All runs were performed with 256 particles and the cut-off radius was taken to be half of the box length. Long-range corrections were taken into account by a method due to Neumann (1987), which is the same as the later developed method of Lustig (1988). The temperature was kept constant by momentum scaling; hence, strictly speaking, we are not working in a canonical but in an isokinetic ensemble. This difference is important for the calculation of the heat capacity, which in the isokinetic ensemble cannot be obtained from the energy fluctuation formula; actually, it was computed by differentiation of the internal energies. The time step for the integration was chosen to be 0.0015 in units of aJm/E. The system was started either from an f.c.c. lattice or a previous configuration followed by an equilibration period of 4000 time steps. The production runs were extended over 15 000 to 20000 time steps. The quantities calculated directly are the pressure and the residual part of the internal energy. From plots of the running averages of these quantities, their statistical uncertainties were estimated; these will be given, together with the results, in the next section.
RESULTS AND DISCUSSION
The state points for the simulations were chosen such that the resulting pressures and isochoric heat capacities could be compared directly with experimental values. For the comparison of the pressure, simulations were performed close to one isochore in the gas, close to two isochores in the liquid and at one supercritical isotherm, which covers nearly the whole fluid region except at very high densities. The directly obtained pressures and
106 TABLE 1 Results from MD simulations: residual internal energy, pressure and their estimated uncertainties for oxygen as 2CLJ fluid (T* = kT/r, p* = PIJ~, u* = U/NC, p* = po3/c)
TV
P*
7.894114
4.031208 5.262742 6.578428 7.894114 3.368155 3.841802 4.341762 4.736468 4.999605 5.262742 5.525880 2.526116 2.789253 3.052391 3.315528
0.206180 0.241655 0.328227 0.368956 0.387309 0.415385 0.177749 0.190305 0.187085 0.184163 0.575423 0.567013 0.560236 0.570771 0.569974 0.569256 0.568579 0.697468 0.689994 0.683964 0.691963
3.578665
0.690787
3.338992 3.388992 3.488992 3.538992 4.095800 4.145800 4.245800 4.295800 2.473613 2.523613 2.623613 2.673613 2.716357 2.766357 2.866357 2.916357 1.423485 1.623485 1.556896 1.606896 1.706896 1.756896
0.572265
0.572265
0.695828
0.695828
0.809987
u* -
4.773 5.534 7.508 8.433 8.841 9.457 5.089 4.971 4,521 4.263 15.116 14.615 14.181 14.234 14.082 13.916 13.795 18.956 18.515 18.129 18.094 18.084 17.816 17.836 15.071 15.035 14.977 14.940 14.613 14.613 14.537 14.515 18.977 18.921 18.825 18.777 18.746 18.689 18.590 18.557 23.167 22.787 22.893 22.844 22.651 22.580
P*
AU*
AP*
Prod. steps
1.56 1.87 2.87 3.48 3.77 4.37 0.31 0.67 1.05 1.37 0.07 1.06 2.11 3.29 3.95 4.50 5.00 0.11 1.07 1.99 3.53 3.53 4.63 4.60 - 0.08 0.06 0.27 0.44 1.81 1.87 2.11 2.26 - 0.28 0.03 0.45 0.78 0.91 1.19 1.68 1.86 -0.86 1.26 0.67 0.93 1.94 2.34
0.025 0.025 0.025 0.025 0.025 0.025 0.050 0.040 0.040 0.040 0.020 0.020 0.020 0.020 0.020 0.020 0.020 0.015 0.015 0.015 0.015 0.015 0.015 0.015 0.015 0.015 0.015 0.015 0.020 0.020 0.020 0.020 0.015 0.015 0.015 0.015 0.020 0.020 0.020 0.020 0.035 0.035 0.035 0.035 0.035 0.035
0.015 0.02 0.03 0.03 0.03 0.03 0.01 0.01 0.01 0.01 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.07 0.07 0.07 0.07 0.10 0.10 0.10 0.10 0.17 0.17 0.17 0.17 0.17 0.17
15000 15000 20000 20000 20000 20000 15000 15000 15000 15000 20000 20000 20000 15000 15000 15000 15000 15000 20000 20000 20000 20000 20000 20000 15c!Oo 2oooo 20000 15000 20000 20000 20000 20000 15000 15000 15000 15000 20000 20000 20000 20000 15000 15000 15000 20000 20000 15000
107 TABLE 2 Comparison of pressures from MD simulations with experimental results for oxygen. The simulated pressures are given together with their estimated statistical uncertainties PMD (MPa)
TK)
Trnoll-‘)
PCXp @Pa)
153.198 200.0 250.0 300.0 128.0 146.0 165.0 180.0 190.0 200.0 210.0 96.0 106.0 116.0 126.0
8.9203 9.5504 9.3888 9.2422 28.8775 28.455 28.1153 28.644 28.604 28.568 28.534 35.0023 34.5272 34.3246 34.726
4.7775 10.7398 16.5349 21.9884 2.3713 18.4948 34.4943 52.8118 61.9364 70.8860 79.6228 3.5267 18.2560 32.4298 56.3733
136.0
34.667
72.9740 d
300.0
10.3471 12.1274 16.472 18.516 19.437 20.846
24.8740 29.9184 45.4145 55.2581 60.4637 69.4934
b ’ ’ ’ ’ = = d d d d = = = d
= ’ d d d d
4.9f0.16 10.6 f 0.16 16.7 f 0.16 21.8 f0.16 l.Of 1.0 16.9 f 1.0 33.5 f 1.0 52.2 f 1.0 62.6 f 1.0 71.3 f 1.0 79.3 f 1.0 1.8 f 1.1 17.0 f 1.1 31.6 f 1.1 56.0& 1.1 55.9f 1.1 73.5 f 1.1 73.0 f 1.1 24.8 f 0.24 29.7 f 0.32 45.6 f 0.50 55.2kO.50 59.9 f 0.50 69.2 f 0.50
- 25.41 1.75 - 1.03 0.80 1.11 1.00 0.52 0.24 - 0.24 -0.14 0.10 0.37 0.27 0.17 0.06 0.08 - 0.08 -0.00 0.27 0.60 - 0.26 0.06 0.49 0.20
a The relative difference in the density that is associated with the difference between the simulated and the experimental pressure on the basis of the Schmidt-Wagner (1985) equation of state. b From Pentermann and Wagner (1978). ’ From Weber (1970). d From Weber (1977).
residual internal energies, together with their statistical uncertainties, are given in Table 1 in reduced units for all stable simulation runs. Table 2 shows the results for the simulated pressures together with their statistical uncertainties in real units in comparison with the experimental data. For two state points, the results of two independent simulation runs are given in order to demonstrate their reproducibility. We note that generally the simulation results agree within their statistical uncertainties with the experimental data. Only on the liquid isochores at the two respectively lowest temperatures is the disagreement slightly larger. It seemed
108 TABLE 3 Comparison of residual internal energies from MD simulations with results from the Schmidt-Wagner (1985) equation of state for oxygen. (The simulated energies are given, together with their estimated statistical uncertainties) ;)mol J-‘)
(J mol-‘)
G..s
Ui%D
153.198 200.0 250.0 300.0 128.0 146.0 165.0 180.0 190.0 200.0 210.0 96.0 106.0 116.0 126.0
8.9203 9.5504 9.3888 9.2422 28.8775 28.4554 28.1153 28.644 28.604 28.568 28.534 35.0023 34.5272 34.3246 34.726
- 1801 - 1605 - 1459 - 1360 - 4781 -4615 -4464 -4466 - 4410 -4356 - 4304 - 5981 - 5839 - 5711 - 5687
136.0
34.667
- 5598
300.0
10.3471 12.1274 16.472 18.516 19.437 20.846
-
-1608f16 -1571*13 -1428f13 -1347f13 -4776& -4618* -44812 -4498& -4450f -4397k -4359f - 5990f -5850& -5728f -5717* - 5714f -5629f -5636f -1508* -1749+ -2372+ -2665k -2794f -2988*
1514 1761 2357 2636 2760 2950
(J mol-‘)
6.3 6.3 6.3 6.3 6.3 6.3 6.3 4.7 4.7 4.7 4.7 4.7 4.7 4.7 8 8 8 8 8 8
10.7 2.1 2.1 0.9 0.1 -0.1 - 0.4 -0.7 -0.9 -0.9 -1.3 -0.2 -0.2 -0.3 -0.5 - 0.5 -0.6 -0.7 0.4 0.7 -0.6 -1.1 -1.2 -1.3
a The relative difference between the simulated and the equation-of-state energies.
interesting also to express the pressure differences as differences of density, which was made on the basis of an empirical equation of state (Schmidt and Wagner, 1985). The relative density differences are also given in Table 2. We learn that, especially at the highest density, the differences are remarkably small. The largest difference occurs for a gas at a temperature close to the critical, which is simply an effect of the nearly zero slope of the isotherm there. Table 3 shows the results for the simulated residual internal energies together with their statistical uncertainties in real units in comparison with the values from the empirical equation of state of Schmidt and Wagner (1985); the state points are the same as those for which the pressures were compared. The table also contains the relative differences between the
109 TABLE 4 Comparison of isochoric heat capacities from MD simulations with experimental results (Goodwin and Weber, 1969) for oxygen C”
(-)MD
;)mol 1-l) 130.692 159.453 97.805 107.030
R
28.719 S 34.920
3.19 3.12 3.51 3.46
3.15 3.05 3.49 3.45
simulated and the equation of state values. We observe that, especially in the liquid, the agreement is very good. On the gas isochore the agreement is somewhat worse, with the largest deviation at 153 K, which is close to the critical temperature. There may be several reasons for the discrepancy. One could be that, owing to the finite number of particles (N = 256), the critical temperature of the simulated system would be higher than that of an infinite system. Then, the state point at p = 8.9 mol 1-l and T= 153.2 K could be already within the thermodynamically unstable region, where it is known that the internal energy obtained from simulations is at least particle number dependent (Heinbuch and Fischer, 1987). Further investigations are necessary in that respect and are planned, together with studies on the dependence of the coexistence curve on the particle number close to the critical point. We also calculated isochoric heat capacities by differentiation of the smoothed simulated isochoric energies as described previously (Fischer et al., 1987). The results are compared in Table 4 with experimental data (Goodwin and Weber, 1969) and the agreement is found to be within the experimental accuracy. As has been mentioned already, oxygen has a large liquid range. The data presented so far extend up to a density of p 2: 35 mol I-‘, which is at P/PC = 2.7. The triple-point density of oxygen, however, is still higher, at P/PC = 3.0. Performing simulations at that density, however, meant that the running averages mostly drifted away if the runs were made sufficiently long. Sometimes they remained stable in a liquid state for a certain time and the thermodynamic data then obtained are still quite reasonable. As an example, the isochoric heat capacity at p = 40.649 mol l-1 and T = 62.967 K was obtained from simulations as c,/R = 4.16 compared to the experimental value of 4.22. This drift of the running averages may indicate that we are already in the liquid-solid two-phase region of the model fluid. This conjecture is confirmed by the fact that the triple-point density of the LJ
110
fluid is at p/p, = 2.7 and the 2CLJ model for oxygen used here is due to its small elongation (L = 0.22) not being very acentric. A conclusion of the present work is that the simple 2CLJ model used here is able to predict the thermodynamic properties of oxygen over a large fluid region. One shortcoming is the fact that the model fluid seems to become thermodynamically unstable owing to the fact that it freezes at densities below the triple-point density of real oxygen. A proper description of the liquid-solid transition would require either a definitively more complicated effective pair potential or the inclusion of three- and higher body forces.
ACKNOWLEDGEMENTS
The authors would like to thank Professor W. Wagner (Bochum) for suggestions on this work and Dr. J.F. Ely (Boulder) for sending the report (NBSIR 77-865). The project was financially supported by the Deutsche Forschungsgemeinschaft under the priority programme “Neue Arbeitsmedien in der Energie- und Verfahrenstechnik”, AZ Fi 287/S, gefordert. Parts of this work were presented at the VDI-~e~od~~ ~olloqium 1989 in Trier.
REFERENCES Bohn, M., Lust&, R. and Fischer, J., 1986. Description of polyatomic real substances by two-center Lennard-Jones model fluids. Fluid Phase Equilibria, 25: 251-262. Fischer, J., 1980. Perturbation theory for the free energy of two-center Lennard-Jones liquids, J. Chem. Phys., 72: 5371-5377. Fischer, J., Lustig, R., Breitenfelder-Manske, H. and Lemming, W., 1984. Influence of intermol~ular potential parameters on orthobaric properties of fluids consisting of spherical and linear molecules. Mol. Phys., 52: 485-487. Fischer, J., Saager, B., Bohn, M. and Haile, J.M., 1987. Specific heat of simple liquids. Mol. Phys., 62: 1175-1185. Gmelin 1958. Handbuch der Anorganischen Chemie, Sauerstoff, System No. 3, Lfg 3. Verlag Chemie, Weinheim/Bergstrasse, 356 pp. Goodwin, R.D. and Weber, L.A., 1969. Specific heats c, of fluid oxygen from the triple point to 300 K at pressures to 350 atmospheres. J. Res. Nat. Bur. Stand., 73A: 15-24. Haile, J.M., 1984. Private communication. Heinbuch, U. and Fischer, J., 1987. On the application of Widom’s test particle method to homogenous and inhomogenous fluids. Mol. Sim., 1: 109-120. Luckas, M. and Lucas, K., 1989. Thermodynamic properties of fluid carbon dioxide from the SSR-MPA potential. Fluid Phase Equilibria, 45: 7-23. Luckas, M., Ripke, M. and Lucas, K., 1990. Prediction of the ~ermophysical properties of fluid methane from the SSR-MPA potential. Fluid Phase Eq~~b~a, 58: 35-50.
111
Lustig, R., 1986. A thermodynamic perturbation theory for non-linear multicentre LermardJones molecules with an anisotropic reference system. Mol. Phys., 59: 173-184. Lustig, R., 1987. Application of thermodynamic perturbation theory to multicentre LennardJones molecules. Results for CF,, Ccl,, neo-C,H,, and SF, as tetrahedral and octahedral models. Fluid Phase Equilibria, 32: 117-137. Lustig, R., 1988. Angle-average for the powers of the distance between two separated vectors. Mol. Phys., 65: 175-179. Lustig, R. and Steele, W.A., 1988. On the thermodynamics of liquid propane. A molecular dynamics study. Mol. Phys., 65: 475-486. Lustig, R., Toro-Labbt, A. and Steele, W.A., 1989. A molecular dynamics study of the thermodynamics of liquid ethane. Fluid Phase Equilibria, 48: l-10. McDonald, I.R. and Singer, K., 1972. An equation of state for simple liquids, Mol. Phys., 23: 29-40. Moller, D. and Fischer, J., 1990. Vapour liquid equilibrium of a pure fluid from test particle method in combination with NpT molecular dynamics simulations. Mol. Phys., 69: 463-473. Neumann, M., 1987. Private communication. Pentermann, W. and Wagner, W., 1978. New pressure-density-temperature measurements and new rational equations for the saturated liquid and vapour densities of oxygen. J. Chem. Thermodyn., 10: 1161-1172. Saager, B. and Fischer, J., 1990. Predictive power of effective intermolecular pair potentials: MD simulation results for methane up to 1000 MPa. Fluid Phase Equilibria, 57: 35-46. Saager, B., Lotfi, A., Bohn, M., Nguyen Van Nhu and Fischer, J., 1990. Prediction of gas PVT data with effective intermolecular potentials using the Haar-Shenker-Kohler equation and computer simulations. Fluid Phase Equilibria, 54: 237-246. Schmidt, R. and Wagner, W., 1985. A new form of the equation of state for pure substances and its application to oxygen. Fluid Phase Equilibria, 19: 175-200. Singer, K., Taylor, A. and Singer, J.V.L., 1977. Thermodynamic and structural properties of liquids modelled by ‘2-Lermard-Jones centres’ pair potentials. Mol. Phys., 33: 1757-1795. Vega, C., Saager, B. and Fischer, J., 1989. Molecular dynamics studies for the new refrigerant R152a with simple model potentials. Mol. Phys., 68: 1079-1093. Weber, L.A., 1970. PVT, thermodynamic and related properties of oxygen from the triple point to 300 K at pressures to 33 MN/m’. J. Res. Nat. Bur. Stand., 74A: 93-129. Weber, L.A., 1977. Thermodynamic and related properties of oxygen from the triple point to 300 K at pressures to 1000 bar. NASA Ref. Pub. 1011, NBSIR 77-865.