Dense fluid self-diffusion coefficient calculations using perturbation theory and molecular dynamics

Dense fluid self-diffusion coefficient calculations using perturbation theory and molecular dynamics

L. A. F. COELHO¹, J. V. OLIVEIRA¹ and F. W. TAVARES²

¹Programa de Engenharia Química, PEQ/COPPE/UFRJ,

21945-970, Rio de Janeiro - RJ, Brazil

²Departamento de Engenharia Química, Escola de Química/UFRJ,

21949-900, Rio de Janeiro - RJ, Brazil

(Received: November 17, 1998; Accepted: July 9, 1999)

Abstract - A procedure to correlate self-diffusion coefficients in dense fluids by using the perturbation theory (WCA) coupled with the smooth-hard-sphere theory is presented and tested against molecular simulations and experimental data. This simple algebraic expression correlates well the self-diffusion coefficients of carbon dioxide, ethane, propane, ethylene, and sulfur hexafluoride. We have also performed canonical ensemble molecular dynamics simulations by using the Hoover-Nosé thermostat and the mean-square displacement formula to compute self-diffusion coefficients for the reference WCA intermolecular potential. The good agreement obtained from both methods, when compared with experimental data, suggests that the smooth-effective-sphere theory is a useful procedure to correlate diffusivity of pure substances.

Keywords: Smooth-sphere theory, self-diffusion coefficient, molecular dynamics.

INTRODUCTION

Accurate diffusion coefficients of pure substances and mixtures normally play an important role in the analysis, design and optimization of separation processes in the chemical industries. As mentioned by Reed and Gubbins (1973), the kinetic theory of gases at low pressure is well established, enabling reasonable estimations of diffusion coefficients. On the other hand, at present, there is no analogous theory capable of providing satisfactory results in the estimation of these properties for dense fluids or for gases under high-pressure conditions.

There are numerous attempts to correlate or predict self-diffusion in dense fluids using theoretical and/or molecular simulation approaches (Ernst et al., 1969; Chandler, 1975; Dorfman, 1981; and Liu and Ruckenstein, 1997). A great effort to interpret self-diffusivities in dense fluids can be observed in the literature. Recently, for example, Iwai et al. (1997) used molecular dynamics to estimate transport coefficients under supercritical conditions needed to design high-pressure separation processes. An extensive review of molecular dynamics simulations and models applied to correlate self-diffusion coefficients of real substances was published by Liu et al. (1998).

Of the several theoretical approaches applied to calculate self-diffusion coefficients, the slightly modified Chapman and Enskog solution of the Boltzmann equation (Reed and Gubbins, 1973) has been used with great success. The major objective here is to test the procedure for calculating self-diffusion coefficients recently proposed by Rocha et al. (1997) by comparing systematically the results obtained with the model with those obtained with molecular dynamics simulations and with experimental data. Although a large number of molecular simulations of self-diffusion coefficients have been reported in the literature, most of them only involve hard-sphere, square-well and Lennard-Jones fluids (Liu et al., 1998). This justifies our objective to publish molecular dynamics simulations for the reference WCA intermolecular potential.

THEORY

The solution of the Boltzmann equation for hard-sphere fluid is known in the literature as Enskog’s theory for the self-diffusion coefficient (Reed and Gubbins, 1973):

(1)

where r^* = rs³ is the reduced density, r the number density, s the hard-sphere diameter, g(s) the radial distribution function at the contact point of the hard-sphere fluid, and D₀ the Chapman-Enskog low-density self-diffusion coefficient given by

(2)

By comparing the hard-sphere self-diffusion coefficients estimated from equation (1) those obtained from molecular dynamics computer simulations, with Speedy (1987) proposed an empirically more accurate equation to calculate these coefficients that better accounts for the density dependence,

(3)

In order to calculate diffusion coefficients for real (not spherical) substances, Rocha et al. (1997) used the Weeks-Chandler-Andersen (Weeks et al., 1971; Andersen et al., 1971) perturbation theory (WCA) to estimate the effective hard-sphere diameters that are dependent on temperature and density. Chandler (1975) pointed out that diffusion coefficients of non-spherical molecules differ from those of spherical ones due to the collision energy transfer from translational to vibrational and rotational energies. The author suggested a rough spherical model to account for this energy transfer and, as a consequence, a coupling factor to multiply the diffusion coefficient obtained from the smooth spherical molecular model.

Chen (1983) and Liu and Ruckenstein (1997) used this theoretical approach to correlate self-diffusion coefficients of real substances. In the literature, for example Schmid et al. (1991), Chandler’s coupling factor is shown to be temperature dependent and to vary from zero to one, as expected. Although the Chandler approach is theoretically consistent, there is no model to estimate the coupling factor of real substances, and especially its dependence on temperature.

Instead, Rocha et al. (1997) included the coupling factor contribution in the smooth sphere by using perturbation theory effective hard-sphere diameters. In this approach, effective diameters are obtained to better estimate self-diffusion coefficients of pure substances. The WCA perturbation theory is presented elsewhere (Boublik et al., 1980) and will not be presented in detail here. This theory is well known for accurately describing the equilibrium thermodynamic properties of simple liquids at high densities, where repulsive forces dominate the structure of fluids. The WCA theory is based on separating rigorously the intermolecular potential into two different parts. The reference potential contribution contains only repulsive forces and the perturbation part includes the attractive contribution, as given by

(4)

(5)

(6)

where ,f^ref(r), f^p(r), and r_m are the reference (repulsive) and perturbation (attractive) potentials and the distance where the actual potential, f(r) , is at a minimum, respectively. The intermolecular potential is given by the 6-12 Lennard-Jones potential as follows:

(7)

The relationship of the reference fluid to the effective hard-sphere system is given by solving the integral equation to match the structure factor of the reference and the hard-sphere systems (Andersen et al., 1971). Souza and Ben-Amotz (1993) recently solved this integral equation and presented an accurate expression to determine the effective hard-sphere diameter, given by

(8)

where T^* = kT/e and r^* = rs³ are the reduced temperature and reduced density, respectively.

Therefore, in the Rocha et al. (1997) procedure, all the s in the self-diffusion equations (1) and (2) are replaced by the WCA effective hard-sphere diameter d(r^*, T^*), calculated by equation (8). This procedure is called here as the modified Speedy model.

MOLECULAR DYNAMICS

The molecular dynamics technique is based on solving the differential equation set formed by Newton’s equation of motion applied to all molecules. Given the molecular positions, velocities, and accelerations at time t, the new positions, velocities, and accelerations can be obtained at a later time, t + dt, with a high degree of accuracy. The continuous molecular trajectories at time t + dt are obtained here by the "velocity Verlet" algorithm (Frenkel and Smit, 1996) which consists of a Taylor expansion of positions about time t and a double Taylor expansion of velocities, one at about time t and another at about time t + (dt/2) as follows:

where r_ia, v_ia, and a_ia are the position, velocity, and acceleration, respectively, of molecule i in direction a. Molecular accelerations are evaluated directly from the Newton’s equation of motion. To calculate the thermodynamic properties in the canonical ensemble, we have used the Hoover-Nose constant temperature molecular dynamics algorithm (Allen and Tildesley, 1987; Frenkel and Smit, 1996). Thus, the equations of motion are

(12)

where m is the molecular mass, F_ia is the force on molecule i in direction a, and x is a kind of friction coefficient that appears in order to keep the temperature constant. This friction coefficient is given by the first-order differential equation (Allen and Tildesley, 1987),

(13)

where f is the number of degrees of freedom and Q is the thermal inertia parameter (in this work, we have used Q = 10^-19 J.s²/mol, as recommended in the literature, Allen and Tildesley, 1987). The temperature set point is Á and the simulated temperature, T, is calculated by the kinetic energy of the system as follows:

(14)

These equations are solved and do generate states in the canonical ensemble (Hoover, 1985 and 1986). The integration time step used is Dt = 0.005(ms²/e)^1/2. The simulations are carried out in a cubic box using the conventional periodic boundary conditions imposed in all directions and the minimum image convention. To speed up the simulation, the linked list method (Allen and Tildesley, 1987) is used to simulate 256 molecules.

The initial molecular positions and velocities are calculated to match with the face-centered cubic lattice positions and to obey the Maxwell-Boltzmann velocity distribution function at the desired temperature, respectively. The system is equilibrated over 4000 iteration time steps to eliminate the artificial initial configuration. After the equilibration process, the thermodynamic properties are calculated by ensemble averaging over 20 independent origins of 20000 iteration time steps. Three completely independent simulations are carried out to represent replicas of the thermodynamic properties of the "computational experiment."

Two equivalent methods are used to obtain self-diffusion coefficients by the equilibrium molecular dynamics simulation technique. One uses the mean-square displacement method (MSD) that calculates the diffusion coefficient by the Wiener-Einstein equation (Medhi, 1994),

(15)

In this equation, r_ia corresponds to the unfolded molecular position of particle i in direction a, which is accumulated before application of the periodic boundary conditions. The other method uses the velocity auto-correlation function of the Green-Kubo formula (McQuarrie, 1976)

(16)

The brackets that appear in both formulas indicate the canonical ensemble average.

RESULTS AND DISCUSSION

As mentioned by Allen and Tildesley (1987), molecular simulation can be a bridge between experiments and theoretical models as molecular simulation provides "exact" macroscopic physical properties (without simplifications) from a given intermolecular potential model. In this work, we have performed molecular simulations by using the equilibrium molecular dynamics technique to test whether or not the modified Speedy expression, Eq. (3), together with the WCA effective hard-sphere diameter, Eq. (8), gives reasonable results for the self-diffusion coefficients of the reference potential. Thus, in order to compare the transport properties obtained from molecular dynamics simulations with the theoretical WCA-Speedy results, we have adopted the reference part WCA potential, f^mf(r), the repulsive contribution, as the intermolecular potential for all simulations presented here.

In order to test the equilibrium molecular dynamics algorithm and the procedure to obtain the self-diffusion coefficients, we have performed simulations and compared the self-diffusion coefficients of the full Lennard-Jones fluids with those reported in the literature (Rowley and Painter, 1997). To test consistency in calculating the self-diffusion coefficients, we have used both methods, the mean-square displacement (MSD) and the velocity auto-correlation function, as suggested by Allen and Tildesley (1987) and Frenkel and Smit (1996). The self-diffusion coefficients obtained by both methods show excellent agreement. These results also agree with those of Rowley and Painter (1997). Here we report only the MSD method due to simplicity and good reproducibility.

In Table 1, we present the results and conditions of fitting of the self-diffusion coefficients of pure fluids using the semi-empirical Speedy formula with the effective hard-sphere diameters from the WCA perturbation theory. The Lennard-Jones parameters of distance were estimated using the Marquardt algorithm (Marquardt, 1963) and are shown in Table 2. We have found that the energy parameters are not sensitive to the experimental self-diffusivity data; therefore, they were fixed equal to those presented in the literature (Reid et al., 1987).

Table 1:
Self-diffusion coefficient average absolute deviations, obtained by the fit procedure.

Substance	T(K)	NP+	AAD%*	Data Source
CO2	323	15	16.8	Etesse et al. (1992)
Ethane	294	5	2.5	Schmid et al. (1991)
Propane	294	5	4.5	Schmid et al. (1991)
Ethylene	298	19	9.9	Arends et al. (1981)
Sulfur hexafluoride	296	6	4.8	DeZwaan and Jonas (1975)

*

;

⁺number of data points

Table 2:
Lennard-Jones size and energy parameters.

Substance	s(10-10m)*	e /k (K)+
Carbon dioxide	3.6003	195.2
Ethane	4.1047	215.7
Propane	4.6398	237.1
Ethylene	3.9671	224.0
Sulfur hexafluoride	4.7808	201.0

^*

⁺Values from the literature (Reid et al., 1987).

In Tables 3, 4, 5, 6 and 7, we present the self-diffusion coefficient results for carbon dioxide, ethylene, ethane, propane, and sulfur hexafluoride obtained with the Speedy smooth effective diameter model (modified Speedy model) and from molecular dynamics simulations using the corresponding WCA reference potential parameters. In these tables, we also report the corresponding experimental data. In Figure 1, we summarize these results and compare the molecular simulation with the modified Speedy model (Figure 1a) and with experimental data (Figure 1b). The comparison shows that the equilibrium molecular dynamics results agree very well with the modified Speedy model, and also that this model gives reasonable results for the self-diffusion coefficients of real substances (Rocha et al., 1997). In Figure 2, we show the self-diffusion calculations for carbon dioxide (Figure 2a) and ethylene (Figure 2b) for a large range of densities. As mentioned previously, both the modified Speedy model and the molecular simulation technique calculate reasonably well the self-diffusion coefficients, even for high density regions.

Table 3:
Self-diffusion coefficients for carbon dioxide at 323 K. Comparison of the results obtained with the modified Speedy model and molecular dynamics simulations with experimental data (Etesse et al., 1992).

Pressure	Density	D11 x 108(m2/s)
(bar)	(g/cm3)	Experimental	Modified Speedy model	Molecular dynamics
10.34	0.0289	110.1	106.0	80.7
20.68	0.0579	54.1	52.5	49.3
33.99	0.0880	31.0	34.0	30.9
48.26	0.1088	19.6	27.2	23.8
63.91	0.1790	13.9	15.9	15.9
72.39	0.2020	12.0	13.9	14.3
84.81	0.2376	8.8	11.6	11.4
101.35	0.3454	5.2	7.5	7.6
109.90	0.4640	4.1	5.2	4.8
136.52	0.6700	2.9	3.1	2.8
206.15	0.7950	2.2	2.4	2.2
274.41	0.8530	1.9	2.1	1.9
343.36	0.8960	1.7	1.9	1.7
411.62	0.9286	1.6	1.8	1.6
482.63	0.9570	1.5	1.7	1.5

Table 4:
Self-diffusion coefficients for ethylene at 298 K. Comparison of the results obtained with the modified Speedy model and molecular dynamics simulations with experimental data (Arends et al., 1981).

Pressure	Density	D11 x 108(m2/s)
(bar)	(g/cm3)	Experimental	Modified Speedy model	Molecular dynamics
58.37	0.1261	11.36	12.40	11.24
62.18	0.1513	9.37	10.00	9.24
65.04	0.1765	7.99	8.29	7.72
67.51	0.2018	6.87	7.00	6.72
70.04	0.2270	6.21	5.99	5.52
73.24	0.2522	5.48	5.18	4.92
77.99	0.2774	4.94	4.51	4.74
85.66	0.3026	4.42	3.96	3.79
98.42	0.3279	3.96	3.48	3.30
119.37	0.3531	3.46	3.06	3.01
152.69	0.3783	3.12	2.70	2.70
203.86	0.4035	2.74	2.38	2.44
279.69	0.4287	2.41	2.09	2.21
388.53	0.4540	2.09	1.83	1.89
540.10	0.4792	1.82	1.60	1.78
745.57	0.5044	1.56	1.38	1.67
1018.81	0.5296	1.34	1.20	1.53
1374.87	0.5548	1.17	1.00	1.23
1828.47	0.5801	0.99	0.85	1.12

Table 5:
Self-diffusion coefficients for ethane at 294 K. Comparison of the results obtained with the modified Speedy model and molecular dynamics simulations with experimental data (Schmid et al., 1991).

Pressure	Density	D11 x 108(m2/s)
(bar)	(g/cm3)	Experimental	Modified Speedy model	Molecular dynamics
250	0.4320	18.7	19.8	20.9
500	0.4880	14.0	14.0	15.8
1000	0.5300	10.3	10.4	11.1
1500	0.5660	8.6	7.8	10.9
2000	0.5840	7.3	6.6	9.3

Table 6: Self-diffusion coefficients for propane at 294 K. Comparison of the results obtained with the modified Speedy model and molecular dynamics simulations with experimental data (Schmid et al., 1991).

Pressure	Density	D11 x 108(m2/s)
(bar)	(g/cm3)	Experimental	Modified Speedy model	Molecular dynamics
250	0.5420	9.10	9.08	10.67
500	0.5620	7.03	7.57	10.13
1000	0.6100	5.65	5.59	6.91
1500	0.6350	4.58	4.38	5.82
2000	0.6500	3.83	3.46	6.24

Table 7:
Self-diffusion coefficients for sulfur hexafluoride at 296 K. Comparison of the result obtained with the modified Speedy model and molecular dynamics simulations with experimental data (DeZwaan and Jonas, 1975).

Pressure	Density	D11 x 109(m2/s)
(bar)	(g/cm3)	Experimental	Modified Speedy model	Molecular dynamics
n.a.	1.3350	9.40	10.16	8.54
84	1.50	7.70	7.52	6.85
161	1.60	6.60	6.16	6.00
311	1.70	5.23	4.95	5.60
574	1.80	3.83	3.89	4.42

Figure 1: Comparison of the calculated self-diffusions for carbon dioxide, ethylene, ethane, propane and sulfur hexafluoride from molecular dynamics simulations with those obtained with the modified Speedy model (a) and with experimental data (b).

Figure 2: Density dependence for the self-diffusion of carbon dioxide (a) and of ethylene (b).

CONCLUSIONS

In this work, we have used the molecular dynamics technique to test a modified Speedy model proposed by Rocha et al. (1997). This model correlates well the self-diffusion coefficients of real substances by using the smooth sphere theory with a temperature and density-dependent effective hard-sphere diameter taken from WCA perturbation theory.

To obtain real self-diffusion coefficients for the WCA reference intermolecular potential over a wide range of densities, the canonical equilibrium molecular dynamics simulations were performed. We have used the Hoover-Nosé constant temperature procedure and the "velocity Verlet" algorithm to solve Newton’s equation of motion applied to all molecules. Also, we have used the mean-square displacement formula to compute the self-diffusion coefficients. By comparing these results with the self-diffusions calculated with the modified Speedy model, we observe good agreement, suggesting that the proposed procedure to calculate pure diffusivities is theoretically consistent. The good agreement obtained by both methods, when compared with experimental data, suggests that the smooth effective sphere theory is a useful procedure to correlate diffusion coefficients of pure substances.

NOMENCLATURE

a Molecular acceleration

d Effective hard-sphere diameter

D Self-diffusion coefficient

f Number of degrees of freedom

k Boltzmann constant

Q Thermal inertia parameter

m Molecular mass

r Molecular position

T Temperature

Reduced temperature

v Molecular velocity

Greek symbols

Lennard-Jones energy parameter

Number density

Reduced density

Lennard-Jones size (diameter) parameter

Temperature set-point

Thermostat friction coefficient

ACKNOWLEDGEMENTS

This work was supported by CNPq and by PRONEX grant no. 124/96 (Brazilian Federal Government). The authors would like to thank Dr. Jaeeon Chang, Fabrício Costa, and Cláudio Dariva for their helpful comments.

• Allen, M.P. and Tildesley, D.J., Computer Simulation of Liquids, Oxford Science Publications (1987).
• Andersen, H.C.; Weeks, J.D. and Chandler, D., Relationship Between the Hard Sphere Fluid and Fluids with Realistic Repulsive Forces, Phys. Rev., A4:1597-1607 (1971).
• Arends, B.; Prins, K.O. and Trappeniers, N.J., Self-Diffusion in Gaseous and Liquid Ethylene, Physica A, 107:307-318 (1981).
• Boublik, T.; Nezbeda, I. and Hlavaty, K., Statistical Thermodynamics of Simple Liquids and Their Mixtures, Elsevier Scientific Publishing Company (1980).
• Chandler, D., Rough Hard Sphere Theory of the Self-Diffusion Constant for Molecular Liquids, J. Chem. Phys., 62:1358-1363 (1975).
• Chen, S.H., A Rough-Hard-Sphere Theory for Diffusion in Supercritical Carbon Dioxide, Chem. Eng. Sci., 38:655-660 (1983).
• DeZwaan, J. and Jonas, J., Density and Temperature Effects on Motional Dynamics of SF₆ in the Supercritical Dense Fluid Region, J. Phys. Chem., 63:4606-4612 (1975).
• Dorfman, J.R., Advances and Challenges in the Kinetic Theory of Gases, Physica A 106:77-101 (1981).
• Ernst, M.H.; Haines, L.K. and Dorfman, J.R., Reviews of Modern Physics, 41:296-315 (1969).
• Etesse, P.; Zega, J.A. and Kobayashi, R., High Pressure Nuclear Magnetic Resonance Measurement of Spin-Lattice Relaxation and Self-Diffusion in Carbon Dioxide, J. Chem. Phys., 97:2022-2029 (1992).
• Frenkel, D. and Smit, B., Understanding Molecular Simulation: From Algorithms to Applications. Academic Press, New York (1996).
• Hoover, W.G., Canonical Dynamics: Equilibrium Phase-Space Distribution, Phys. Rev. A, 34:2499-2503 (1985).
• Hoover, W.G., Constant Pressure Equation of Motion, Phys. Rev. A, 34:2499-2503 (1986).
• Iwai, Y.; Higashi, H.; Uchida, H. and Arai, Y., Molecular Dynamics Simulation of Diffusion Coefficients of Naphthalene and 2-Naphthol in Supercritical Carbon Dioxide, Fluid Phase Equilibria, 127:251-261 (1997).
• Liu, H. and Ruckenstein, E., Predicting the Diffusion Coefficients in Supercritical Fluids, Ind. Eng. Chem. Res., 36:888-895 (1997).
• Liu, H.; Silva, C.M. and Macedo, E.A., Unified Approach to the Self-Diffusion Coefficients of Dense Fluids over Wide Ranges of Temperature and Pressure Hard-Sphere, Square-Well, Lennard-Jones and Real Substances, Chem. Eng. Sci., 53:(13), 2403-2422 (1998).
• Marquardt, D.W., An Algorithm for Least Squares Estimation of Nonlinear Parameters, J. Soc. Ind. Appl. Math, 11: 431-440. (1963).
• McQuarrie, D.A., Statistical Mechanics, Harper Collins Publishers, New York (1976).
• Medhi, J., Stochastic Processes, John Wiley & Sons Inc., New York (1994).
• Reed, T.M. and Gubbins, K.E., Applied Statistical Mechanics, McGraw-Hill Chemical Engineering Series, New York (1973).
• Reid, R.C.; Prausnitz, J.M. and Poling, B.E, The Properties of Gases and Liquids, 4th Edition, McGraw-Hill, New York (1987).
• Rocha, S.R.P.; Oliveira, J.V. and Rajagopal, K., An Evaluation of Density Corrections for Estimating Diffusivities in Liquids and Liquid Mixtures, Chem. Eng. Sci., 52:1097-1109 (1997).
• Rowley, R.L. and Painter, M.M., Diffution and Viscosity Equations of State for a Lennard-Jones Fluid Obtained from Molecular Dynamics Simulations, International Journal of Thermophysics, 18, 1109-1121, (1997).
• Schmid, A.G.; Wappmann, S.; Has, M. and Ludemann, H.D., Self-diffusion in the Compressed Fluid Lower Alkanes: Methane, Ethane, and Propane, J. Chem. Phys., 94:5643-5649 (1991).
• Souza, L.E.S. and Ben-Amotz, D., Optimized Perturbed Hard Sphere Expressions for the Structure and Thermodynamics of Lennard-Jones Fluids, Molec. Phys., 78:137-149 (1993).
• Speedy, R.J., Diffusion in the Hard Sphere Fluid, Molec. Phys., 62:509-515 (1987).
• Weeks, J.D.; Chandler, D. and Andersen, H.C., Role of Repulsive Forces in Determining the Equilibrium Structure of Simple Liquids, J. Chem. Phys., 54:5237-5247 (1971).

Publication Dates

• Publication in this collection
  16 Dec 1999
• Date of issue
  Sept 1999

History

• Accepted
  09 July 1999
• Received
  17 Nov 1998