J. Phys. Chem. C 2007, 111, 16013-16020
16013
Predictions of Transport Properties in Gaseous Mixtures of Sulfur Hexafluoride and Nitrogen† Aurelio Olivet and Lourdes F. Vega*,‡ Institut de Cie` ncia de Materials de Barcelona, Consejo Superior de InVestigaciones Cientı´ficas (ICMAB-CSIC), Campus de la UAB, 08193 Bellaterra, Spain ReceiVed: June 4, 2007; In Final Form: August 12, 2007
We present and discuss here results concerning the implementation of molecular dynamics simulations with the Green-Kubo formalism to predict transport coefficients in SF6/N2 mixtures under conditions of practical interest. SF6 was modeled using a flexible molecular force field recently proposed [Olivet, A.; Vega, L. F. J. Chem. Phys. 2007, 126, 144502], while the force field of Galassi and Tildesley was used for N2 modeling [Galassi, G.; Tildesley, D. J. Mol. Simul. 1994, 13, 11]. The influence of mixture composition on the mutual diffusion coefficients and shear viscosities was investigated by a series of molecular dynamic simulations of SF6/N2 mixtures with different compositions at 300 K and 1 MPa. Temperature and pressure dependencies of the mutual diffusion coefficient were also investigated for a fixed composition of SF6/N2, covering the temperature range from 260 to 340 K and two pressure values (1 and 2 MPa). Simulations with these anisotropic force fields led to mutual diffusion coefficients that diminish when the SF6 concentration in the mixture is increased. This behavior differs from the hypothesis of no composition dependence assumed for several estimation methods available in the literature for gaseous binary mixtures at low pressures. Such a composition dependence is explained on the basis of the deviations from ideality of SF6/N2 mixtures. The effect of the molecular flexibility considered in the SF6 molecular model is analyzed by comparing the transport coefficients estimated with the aforementioned force fields to those estimated with simple Lennard-Jones potentials. Differences between predictions with the anisotropic force fields and the simpler Lennard-Jones potentials were more significant in the case of the shear viscosity; while viscosities estimated with the Lennard-Jones potentials showed no dependence on composition, the estimations based on the anisotropic force fields exhibited a maximum value for the viscosity, which is the same behavior described by an empirical method. Results presented here are pure predictions, thus covering part of the lack of information existing about the transport coefficients of these mixtures, which is of relevance for several industrial processes.
1. Introduction Molecular simulation tools have acquired a degree of maturity at present times up to the point of being capable of making quantitative predictions of physicochemical properties of systems that are currently used in the industry. Although by no means will they replace experimental data, which is always needed under certain conditions and for the final design, there are situations in which molecular simulation predictions, arising from well-tested intermolecular potentials and models, can cover a gap in the data required for practical applications. The success in the use of these techniques for practical applications is due to several concomitant factors. Among them are the high speed of present computers, the development of simulation algorithms to tackle different problems, and the development of more refined force fields, in which more details of the physics of the system are kept, searching for quantitative predictions. The results presented here are a combined effect of some of these factors. They serve as a validating test for the recently developed refined force field for sulfur hexafluoride (SF6), and calculations with this refined model are possible thanks to the computer †
Part of the “Keith E. Gubbins Festschrift”. * To whom correspondence should be addressed. E-mail: lvega@ icmab.es. Phone: +34 935 929 950. Fax: +34 935 929 951. ‡ Present address: MATGAS-Air Products. Campus de la UAB, 08193 Bellaterra, Spain.
power available today. The aim of this work is to use simulations to predict the transport behavior of mixtures of sulfur hexafluoride and nitrogen in order to test the behavior of these mixtures (of practical interest) for which experimental data are not available and only empirical correlations are used. Nowadays, SF6 represents the best insulating gas available in the electrical industry. Its excellent dielectric and arcquenching properties make it suitable for use in currentinterruption equipment and numerous gas-insulated equipment, such as circuit breakers, transmission lines, transformers, and electrical substations.1,2 Furthermore, SF6 is used in several other applications because of its other useful properties like being nontoxic, thermally stable, and chemically inert. In spite of these favorable properties, SF6 is an efficient infrared absorber, and due to its high chemical stability, its atmospheric degradation occurs very slowly. Hence, efficient methods are under development for handling and recovering SF6 after industrial usage, and at the same time, alternative gases and/or gaseous mixtures for the same final applications are also being investigated.2,3-8 Among the many gases and gaseous mixtures studied in the literature, a widespread acceptance exists about the suitability of SF6/N2 mixtures as an alternative for pure SF6 in different applications. SF6/N2 mixtures are cheaper than pure SF6 and, under the appropriate pressure conditions, can perform equally
10.1021/jp0743222 CCC: $37.00 © 2007 American Chemical Society Published on Web 10/09/2007
16014 J. Phys. Chem. C, Vol. 111, No. 43, 2007 well for both electrical insulation and arc interruption.2 An important issue concerning the use of these mixtures is the SF6 recovery. Due to the difference in the boiling points of SF6 and N2, their separation by liquefaction would require an excessive energy consumption to cool and pressurize the mixture. At room temperature (293 K), pure SF6 liquefies at 2 MPa, while a mixture containing 10% of SF6 liquefies at 20 MPa.5 Recently, different methods to separate SF6/N2 mixtures have been discussed and experimentally tested, most of them using adsorbents. Toyoda and co-workers5 designed a prototype equipment based on the pressure swing adsorption (PSA) technique, which used the molecular sieving effect of a Ca-A type zeolite to absorb N2 but not SF6. Yamamoto and coworkers6 proposed a system with two polymer membranes connected in cascade. They analyzed the influence of different operating variables, such as the gas feeding pressure and the membrane temperature, on the gaseous mixture separation. Murase and co-workers7 developed a different sort of zeolite filter, without making use of the PSA technique, but employing a Na-X type zeolite with a high selectivity to adsorb SF6. Shiojiri and co-workers8 separated SF6 from gaseous mixtures containing N2 using a porous Vycor glass membrane, discussing the influence of surface diffusion and capillary condensation on the membrane operation. The performance of the aforementioned processes is determined by the operating variables at the continuous level and the different phenomena occurring at the atomistic level. The design of these processes requires the knowledge of several of the physical properties of the mixture to be separated, including transport properties such as the diffusion coefficient and the shear viscosity. In order to obtain these properties, theoretical methods describing bulk transport coefficients in binary gaseous mixtures are widely used, and several empirical correlations are available in the literature.9,10 These methods include low and high pressures, mixtures of nonpolar gases, mixtures containing polar gases, and binary mixtures at supercritical conditions. In the case of SF6/N2 mixtures, the validity of these methods has been poorly tested because of the scarce number of experimental data available for these mixtures.11-13 In particular, the extended review on gaseous diffusion coefficients presented by Marrero and Mason11 contains only a few data for SF6/N2 mixtures, specifically, some experimental measurements at atmospheric pressure between 328 and 473 K and some predicted values at the same pressure and higher temperatures. Additional data were reported by Kestin and co-workers,12 who measured the viscosities at atmospheric pressure of two SF6/N2 mixtures with N2 mole fractions of 0.3357 and 0.8780, covering a temperature range from 298 to 483 K. Moreover, these authors reported binary diffusion coefficients derived from the mixture viscosity data, using the Chapman-Enskog solution of the Boltzmann equation, supplemented by an extended law of corresponding states. Finally, Haghighi and co-workers13 used a similar approach to determine both the viscosity and the diffusion coefficients for equimolar SF6/N2 mixtures, also at atmospheric pressure but in a wider temperature range from 273 to 1273 K. Since these are the only experimental data we could find, the fact is that there is a lack of information concerning transport coefficients in SF6/N2 mixtures at low pressure, presumably because these mixtures are assumed to behave near ideality at these conditions and, hence, the calculations rely on the application of existing estimation procedures. From the modeling point of view, although SF6 and N2 have been successfully modeled as spherical molecules in particular applications, it is known that the molecular flexibility of SF6 plays a key role in
Olivet and Vega some of its properties; indeed, deviations from ideality for specific properties could arise from this molecular flexibility and, also, as a consequence of the different molecular sizes of N2 and SF6. One of the purposes of the present work is to test the influence of the molecular flexibility on the transport properties of the mixture as obtained by molecular simulations. An additional goal is to use these results to compare with the existing correlations for these transport properties, in order to validate them. This work is a research effort toward relating the macroscopic behavior of SF6 with its microscopic behavior, from a general perspective.14-16 In this context, and because of the aforementioned applications, the main goal of the present work is to provide predictions on bulk transport properties (diffusion coefficients and shear viscosities) of different gaseous SF6/N2 mixtures at low pressure from a molecular perspective. An additional goal is to further check the validity of the recently developed force field for SF6 simulations.16 Molecular dynamics (MD) simulations have been carried out for SF6/N2 mixtures of different compositions at 300 K and 1 MPa. The influence of the temperature and pressure on the diffusion coefficient has been investigated for a mixture containing 25% in volume of SF6 in the temperature range 260-340 K for pressure values of 1 and 2 MPa. Although the number of published work concerning the simulation of transport properties of pure fluids is very large compared to those of binary mixtures, some work can be found on the behavior of binary mixtures using both simple Lennard-Jones potentials17-22 and anisotropic models for complex and/or polar molecules.23-28 Of special relevance to the present work is the evaluation of the force field influence on the transport properties of pure and multicomponent systems performed by Dysthe and co-workers.29,30 They evaluated the available force fields for linear and branched alkanes, finding evidence for the influence of flexible terms, such as the torsion potential included in the alkanes force fields, on the accurate prediction of transport properties. The rest of the paper is organized as follows. The formalism used to calculate the diffusion coefficients and the shear viscosity are explained in sections 2 and 3, respectively. Molecular force fields used for modeling SF6 and N2 are described in section 4. Relevant details concerning the simulation work are explained in section 5 followed by the main results and discussions in section 6. Section 7 provides the concluding remarks and guidelines for future work. 2. Formalism To Calculate the Diffusion Coefficients Molecular diffusion has to be understood as the net transport of mass within a single-phase pure compound, or a mixture, in the absence of a mixing effect produced by mechanical means or by convection.9 The driving force in a molecular diffusion process can be an external force field acting over the molecular species or a gradient in the following variables: pressure, temperature, or concentration. This work deals only with the diffusion process due to a concentration gradient in isothermal and isobaric systems. Under these conditions, species A and B in a binary mixture will diffuse among themselves with net molar diffusion fluxes proportional to their respective composition gradients, as indicated by Fick’s first law:
JM A ) -cDAB∇xA
(1)
JM B ) -cDBA∇xB
(2)
Transport Properties in Gaseous Mixtures
J. Phys. Chem. C, Vol. 111, No. 43, 2007 16015
where c refers to the total molar concentration and ∇xA and ∇xB represent the gradient of molar fractions for both species. The coefficients of proportionality labeled as DAB and DBA are defined as mutual diffusivities or mutual diffusion coefficients. In a binary system, since ∇xA + ∇xB ) 0, DAB ) DBA. Thus, the molecular diffusion process in a binary system can be described by the unique mutual diffusion coefficient DAB. When a pure component is studied, the coefficient appearing in the Fick’s first law is defined as the self-diffusion coefficient, DAA. This represents a measure of the compound mobility into itself. In the special case, where A and B exhibit very similar molecular properties, the self-diffusion coefficient of pure compounds will be approximately equal to the mutual diffusion coefficient, DAB. If the mobility of compound A is greater than that of compound B, the self-diffusion coefficients indicate the lower and upper bounds of the mutual diffusion coefficient: DAA e DAB e DBB. Another type of diffusion coefficient studied in mixtures is the tracer diffusion coefficient, DA′B. This one is related to the mobility of component A in the presence of a second component B. At the dilute limit of A in B, the tracer diffusion coefficient of A approaches its mutual diffusion coefficient at infinite dilution: DA′B f D0AB. Contrarily, at the pure component limit of A, its tracer diffusion coefficient approaches its self-diffusion coefficient: DA′B f DAA. The three types of diffusion coefficients defined above can be calculated in the framework of molecular simulations. Zhou and Miller31 derived the Green-Kubo formulas to calculate the mutual diffusion coefficients in n-component isothermal and isobaric systems. In their development, the diffusivity matrix containing the diffusion coefficients is represented as the product of a kinetic matrix and a thermodynamic matrix. The first one is expressed in terms of the time integrals of the flux correlations, and the second one is written as a function of the spatial integrals of the radial distribution functions. In a binary system constituting NA particles of the species A and NB particles of the species B, the Green-Kubo formula used to calculate the mutual diffusion coefficients in the barycentric reference frame takes the following form:
DAB )
Q 3NxAxB
∫0∞ 〈J(t + t0)‚J(t)〉 dt
(3)
where N is the total number of particles in the system, with N ) NA + NB. The molar fraction of species A and B are xA and xB, and J(t) is referred to the relative velocity vector, which is defined as a function of the center of mass particles’ velocities (vi) and the molar fractions of both species as follows: NA
J(t) ) xB
1 1 + xAxBF(ΓAA + ΓBB - 2ΓAB)
(4)
Equations 1 and 2 indicate that the flux of diffusing species is proportional to the concentration gradient. However, diffusion can be affected by the intermolecular interactions occurring at the microscopic level, introducing complex composition dependence in the diffusion flux. The variable Q in eq 3 is defined as a thermodynamic factor which must be included to correct such effects. The Q factor is related to the compositional derivative of chemical potential by Q ) (xA/RT)(∂µA/∂xA). This factor can be calculated in molecular dynamics simulations by the Kirkwood and Buff expression, which relates Q with the radial distribution functions, as follows:32
(5)
where the number density F is equal to N/V, with V the volume of the simulated system, and the coefficient Γij represents the spatial integral of the radial distribution function (gij) for the pair of species ij as follows:
∫0∞ [gij(r) - 1]4πr2 dr
Γij )
(6)
In the case of ideal binary mixtures, the thermodynamic factor Q is very close to unity. If the Q factor is equal to 1, the mutual diffusion coefficient defined by eq 3 is equivalent to the Maxwell-Stefan diffusion coefficient.26 Combining eqs 3 and 4 the mutual diffusion coefficient can be expressed as
DAB )
Q 3NxAxB
〈(
∫0∞
NA
xB
∑ i)1
NB
vi (t + t0) - xA
(
xB
NA
∑ i)1
)
vj(t + t0) ‚ ∑ j)1 NB
vi (t) - xA
vj(t) ∑ j)1
)〉
dt (7)
Arranging this expression, it is possible to express the mutual diffusion coefficient as a linear combination of the self-diffusion coefficients of species A and B, which are calculated as the integrals of the velocity autocorrelation functions (VACF) of both species, and a distinct diffusion coefficient DdAB, which contains the velocity crossed correlation functions (VCCF). Under certain conditions, these crossed terms and the thermodynamic factor Q could be considered to have a slight influence over the mutual diffusion coefficients and, hence, are ignored. In the calculations reported here, none of these simplifications were done; all terms involved in eq 3 were included in the estimation of the mutual diffusion coefficients. 3. Formalism To Calculate the Shear Viscosity Shear viscosity accounts for the resistance of a fluid to move when submitted to a shear stress. It represents the proportional constant relating the flux of momentum with the gradient of velocity. The Green-Kubo equation used to calculate the shear viscosity is equal to the time integral of the pressure tensor autocorrelation function, as indicated below:33
η)
〈
V
〉
∞ Pxy(t - t0)Pxy(t0) ∫ ∑ 0 3k T x y and indicates the inclusion of pressure tensor elements xy, xz, and yz. Since the pressure tensor autocorrelation function is a collective function, it cannot be statistically improved by the number of particles in the system.34 The elements of the pressure tensor are calculated during the simulation by the following expression:
Pxy )
1 V
[∑ j
mjVxjVyj +
1 2
]
rxijfyij ∑ i*j
(9)
16016 J. Phys. Chem. C, Vol. 111, No. 43, 2007
Olivet and Vega
TABLE 1: SF6 and N2 Force Field Parametersa molecular weight
S (g/mol) F (g/mol) N (g/mol) r0 (Å) Kr (kJ/(mol Å2)) θ0 (degree) Kθ (kJ/mol) σ (Å) (kJ/mol) σ (Å) (kJ/mol)
S-F stretching F-S-F bending F-F Lennard-Jones N-N Lennard-Jones a
32.0660 18.9984 14.0067 1.565 693.48 90 307.36 2.7692 0.60807 3.3100 0.31013
Values taken from refs 16 and 36.
m being the mass of the particle j, V its velocity, and r and f represent the distance and the force between the two particles, respectively. Subscripts x and y are referred to space components. 4. Molecular Force Fields We explain next the molecular potentials used for both the compounds in the mixture. SF6 has been described by an accurate flexible model recently proposed by us.16 This force field was obtained by simultaneously fitting both vapor-liquid equilibria and transport properties to available experimental data. The suitability of this force field for the SF6 simulations was demonstrated by an exhaustive evaluation regarding several equilibrium, interfacial, and transport properties over a wide range of conditions. Given the accuracy shown as compared to several of the experimental properties, and the special feature of considering flexibility, responsible for some of the properties exhibited by this compound, we have decided to use this force field in the present work. The force field defining the SF6 molecule contains two parts: a LJ potential quantifying the F-F intermolecular interactions, and a second part dealing with the intramolecular interactions related to molecular flexibility. This second part includes 6 S-F bond stretches and 12 bending deformations, 1 for each F-S-F angle within the SF6 molecule. The SF6 force field can be expressed by the following equations:
U ) Uinter + Uintra M-1 M
Uinter ) Uintra ) M
[
∑ ∑ ∑ ∑ 4 m)1 n)1 i∈m j∈n
(10)
[( ) ( ) ] σ
rij
12
-
σ
6
(11)
rij
]
∑ ∑ Kr(r - r0)2 + F-S-F∑angles Kθ(θ - θ0)2 m)1 S-F bonds
(12)
In these equations, M is the number of SF6 molecules in the simulated system and rij is the distance between the fluorine atom i belonging to molecule m and the fluorine atom j belonging to molecule n. Since SF6 has an octahedral molecular geometry, the average F-S-F angle (σ0) is equal to 90°. The average length for the S-F bond (r0) was taken as 1.565 Å, the same value used in a rigid model of SF6 proposed by Pawley.35 The rest of parameters appearing in the force field, i.e., the length parameter σ, the well depth , the stretching constant Kr, and the bending constant Kθ were adjusted by the multivariable optimization procedure described in ref 16. The values used for these parameters are provided in Table 1. To be consistent with the force field used for SF6, a force field without charged sites was also used for N2, specifically, the one proposed by Galassi and Tildesley.36 This force field uses a rigid dumbbell representation for N2 molecules, where
N-N intermolecular interactions are quantified by a LJ potential, like the one in eq 11 (parameters are also provided in Table 1). Sokhan and co-workers37 studied transport properties of N2 in single-walled carbon nanotubes using two different force fields, the one of Galassi and Tildesley and a second one including charged sites, finding essentially the same results for both of them. To check the validity of the Galassi-Tildesley model for our purposes, MD simulations were carried out to calculate two points over the vapor-liquid coexistence curve at 70 and 98 K. The saturated vapor densities and the vapor pressures obtained by simulations presented a relative error around 10%, with respect to the experimental data reported in ref 40, while in the case of the saturated liquid densities this error was less than 1%. Although these results do not represent an exhaustive evaluation of the force field, they showed that reasonable good estimations for N2 could be obtained. The interactions between SF6 and N2 molecules were calculated following the standard Lorentz-Berthelot combining rules. In order to compare the performance of the abovementioned refined force fields with that of the standard LJ force fields, additional simulations were run using the LJ parameters for SF6 and N2 taken from the work of Kestin et al., as reported in refs 38 and 39. 5. Simulation Details Molecular simulations were performed with the molecular dynamics simulation package DL_POLY.41,42 All simulations were carried out with a system of 512 molecules, where the respective number of SF6 and N2 molecules were defined by the desired composition of the mixture. Transport coefficients were calculated at 300 K and 1 MPa for varying mixture compositions starting from pure N2 to pure SF6, including 10, 25, 50, 75, and 90% of SF6 by volume in these mixtures. The temperature dependence of the mutual diffusion coefficient was also analyzed at a fixed 25% SF6 content by volume for temperatures ranging from 260 to 340 K and at the pressure values of 1 and 2 MPa. A time step equal to 0.001 ps and a cutoff radii equal to 19 Å were used in all the simulations. Transport coefficients were calculated by means of the Green-Kubo formalism explained in sections 2 and 3. This approach requires the knowledge of the center of the mass velocities to estimate the diffusion coefficients and the knowledge of the pressure tensor elements to estimate the viscosities. MD simulations in the NVE ensemble were performed to record these data. Since transport coefficients were required at specific conditions of temperature and pressure, NPT MD simulations were performed first. These simulations ensured the density corresponding to the required temperature and pressure values. As starting configurations for the NPT simulations, SF6 and N2 molecules were placed on a regular lattice into the simulation box. The Nose`-Hoover thermostat and the Hoover barostat43 implemented in the DL_POLY code were used for these simulations. After some preliminary runs, the thermostat and barostat constants were fixed at 0.01 and 0.2 ps, respectively. The number of simulation steps for NPT was equal to 10 × 106 time steps. For the averaging stage performed through the NVE simulations, at least 20 × 106 time steps were run. The number of time steps in the NVE simulations was chosen to ensure enough statistics to determine transport coefficients. Velocities at each 500 time steps and pressure tensor elements at each 10 steps were recorded. Preliminary simulations were also performed to verify the appropriateness of these time resolutions for transport coefficients estimations.
Transport Properties in Gaseous Mixtures
Figure 1. Normalized VACF used to calculate the tracer diffusion coefficients of SF6 (dashed line) and N2 (dotted dashed line), and normalized RVCF used to calculate the mutual diffusion coefficient (solid line). This data corresponds to a SF6/N2 mixture with an SF6 concentration of 50% by volume. The inset in the figure depicts the running integral of the autocorrelation functions (see text for simulations details).
Figure 2. Radial distribution functions for a mixture with an SF6 concentration of 50% by volume: dashed line, SF6-SF6; dotted dashed line, N2-N2; solid line, SF6-N2. The inset in the figure shows fluctuations of the radial distribution functions (see text for details).
6. Results and Discussion The molecular systems studied along this work correspond to low density gaseous states. Under these conditions very long simulations are required to obtain reliable statistics in the estimation of transport coefficients. The typical time of the simulations reported here was over the 20 ns. Figure 1 depict an example of the normalized VACFs and the relative velocity correlation functions (RVCF) used to calculate the tracer and mutual diffusion coefficients, respectively. This figure corresponds to a SF6/N2 mixture with a SF6 concentration of 50% by volume. Clearly, the VACFs are more accurate than the RVCFs. The latter is a collective property that cannot be improved by the number of molecules in the simulated system. The VACFs of SF6 and N2 decay to zero around 150 and 100 ps, respectively, and in both cases remain at zero without appreciable fluctuations. In contrast, the RVCFs in Figure 1 present important fluctuations after decaying to zero. Those fluctuations produce rather inaccurate time integral functions, which prohibit a direct reading of the mutual diffusion coefficient. Here, following the previous work, the integral over time of both VACF and RVCF was adjusted to an empirical function
J. Phys. Chem. C, Vol. 111, No. 43, 2007 16017
Figure 3. Diffusion coefficients obtained from MD simulations for the anisotropic force fields as a function of the mixture composition: circles, mutual diffusion coefficients (DSF6-N2); down triangles, SF6 tracer diffusion coefficients (DSF6′-N2); up triangles, N2 tracer diffusion coefficients (DN2′-SF6); squares, self-diffusion coefficients (DSF6 and DN2). These results correspond to SF6/N2 mixtures at 300 K and 1 MPa. Mutual diffusion coefficients predicted with the correlation of Fuller et al. are indicated by the dashed line. Dotted lines are guides for the eyes.
Figure 4. Mutual diffusion coefficients from MD simulations as a function of the mixture composition: circles, simulations with the SF6 flexible model of Olivet and Vega and the N2 force field of Galassi and Tildesley; open squares, simulations done with simple LennardJones force fields for both SF6 and N2. Results correspond to SF6/N2 mixtures at 300 K and 1 MPa.
of the form DAB(t0) ) Kτ[1 - exp(-t0/τ)].44 This function represents the time integral of an exponential decaying autocorrelation function, in which K and τ stand for fitting parameters obtained by the least-squares method, applying a weighted factor equal to (1/t0)2. This factor assigned lower weight to the data correlated at long delay times, for which poor statistics is expected. As it was explained in section 2, the thermodynamic factor Q can be calculated by integrating radial distribution functions with respect to distance, combining eqs 5 and 6. Since the radial distribution functions tend toward unity as the distance increases, the term [gAB(r) - 1] in eq 6 must tend toward zero at larger r values. In practice, this term varies close to zero due to the small statistical fluctuations exhibited by the radial distribution functions. Therefore, when the term [gAB(r) - 1]4πr2 is calculated, small fluctuations are magnified as r increases, making the calculation of the Γij coefficients inaccurate.28 When radial distribution functions are calculated, a compromise exists between the distance resolution and the number of simulation
16018 J. Phys. Chem. C, Vol. 111, No. 43, 2007
Olivet and Vega
TABLE 2: Mutual Diffusion Coefficient, Tracer Diffusion Coefficient, Thermodynamic Factor Q, and Shear Viscosity Calculated from MD Simulations at 300 K and 1 MPa using the Olivet-Vega Flexible Force Field for SF6 and the Force Field of Galassi and Tildesley for N2a % vol SF6
MD T (K)
MD P (MPa)
0 10 25 50 75 90 100
300 ( 1 299 ( 3 299 ( 3 303 ( 3 302 ( 3 301 ( 3 297 ( 3
0.99 ( 0.03 0.95 ( 0.30 0.96 ( 0.46 1.10 ( 0.80 1.04 ( 0.96 1.11 ( 1.14 1.00 ( 1.07
a
Q factor
106 × DSF6-N2 (m2/s)
0.9913 0.9831 0.9709 0.9797 0.9905
2.151 ( 0.013 1.01 ( 0.03 1.05 ( 0.03 0.88 ( 0.04 0.83 ( 0.04 0.70 ( 0.10 0.298 ( 0.006
106 × DSF6′-N2 (m2/s) 0.852 ( 0.006 0.694 ( 0.004 0.442 ( 0.003 0.3589 ( 0.0011 0.301 ( 0.003
106 × DN2′-SF6 (m2/s)
η (µPa s)
1.95 ( 0.02 1.687 ( 0.006 1.169 ( 0.017 0.997 ( 0.018 0.867 ( 0.015
17.36 ( 0.64 17.52 ( 0.75 17.92 ( 1.00 16.76 ( 1.00 15.42 ( 0.77 14.64 ( 0.80 14.73 ( 0.44
Self-diffusion coefficients and shear viscosities for pure SF6 and N2 are also included.
TABLE 3: Mutual Diffusion Coefficient, Thermodynamic Factor Q, and Shear Viscosity Calculated from MD Simulations around 300 K and 1 MPa using Simple Lennard-Jones Force Fields for Both SF6 and N2 % vol SF6 10 50 90
MD T (K)
MD P (MPa)
301 ( 1 0.95 ( 0.30 303 ( 3 1.10 ( 0.80 301 ( 3 1.11 ( 1.14
106 × DSF6-N2 Q factor (m2/s) 0.9950 0.9709 0.9905
0.90 ( 0.11 0.87 ( 0.02 0.81 ( 0.02
η (µPa s) 16.1 ( 0.9 16.4 ( 0.9 16.2 ( 0.9
steps. Statistical fluctuations are supposed to increase when the distance resolution becomes smaller and to diminish for longer simulations. Nevertheless, fluctuations are difficult to avoid, in particular for the case of a flexible molecule, where the center of mass is continuously changing. Figure 2 shows an example of the radial distribution functions obtained in this work. In the inset of this figure the aforementioned statistical fluctuations can be observed. In previous simulation works, some authors have calculated the thermodynamic factor Q from the activity coefficient of the mixture.19,27 The approach taken here was to adjust the radial distribution function to a smooth function, in order to obtain the appropriate tendency to reach unity at long distances. Thus, the term [gAB(r) - 1]4πr2 was numerically integrated, and the Γij coefficients directly obtained from the simulations were used to estimate the thermodynamic factor Q. Tracer and mutual diffusion coefficients as a function of the mixture composition are depicted in Figure 3. Mutual diffusion coefficients were calculated from MD simulation results using eq 3, so they include the correction introduced by the Q factor. Diffusion coefficients and the Q values are listed in Table 2. Self-diffusion coefficients for pure SF6 and N2 were also calculated from MD simulations, as shown in Figure 3 and listed in Table 2. At the temperature and pressure conditions investigated (300 K and 1 MPa), both SF6 and N2 exist in gaseous state and their behavior is nearly ideal. However, the estimated N2 self-diffusion coefficient is approximately 7 times that of SF6. This difference is a consequence of the smaller size and lower molecular weight of the N2 molecule. When both gases coexist in a mixture, their tracer diffusion coefficients decrease as the SF6 content in the mixture increases. Mutual diffusion coefficients illustrated in Figure 3 are between the tracer diffusion coefficients of both species, and they diminish as the SF6 concentration increases, but with a smaller slope than in the case of the tracer diffusion coefficients. In order to validate the order of magnitude of the coefficients predicted here, the mutual diffusion coefficients obtained from MD simulations were compared with the estimations made with the Fuller et al. correlation.45 This is a simple correlation based on the Chapman-Enskog equation, whose parameters were adjusted by a regression analysis of several experimental data. A review presented in ref 9 recommends its use with gaseous
mixtures of nonpolar gases at low pressures. Fuller et al. correlation defines the mutual diffusion coefficient as
DAB )
0.00143T1.75
1/3 2 [ ∑V)1/3 A + (∑V)B ]
P(MAB)1/2 (
(13)
In this equation, T is the temperature in K, P is the pressure in bar, and MAB ) 2[(1/MA) + (1/MB)]-1, with MA and MB the molecular weights of species A and B in g/mol. The factor (ΣV) is calculated for each component in the mixture by summing certain atomic diffusion volumes. For SF6 and N2 this factor is equal to 71.3 and 18.5, respectively. For the SF6/N2 mixture at 300 K and 1 MPa, this correlation predicts a mutual diffusion coefficient equal to 0.977 × 10-6 m2/s. This value has been depicted in Figure 3 along with all the composition ranges, demonstrating that the order of magnitude of the mutual diffusion coefficients estimated from MD simulations is in good agreement with the empirical predictive method. The correlation of Fuller et al., similar to other estimation methods based on the diffusion theory for gaseous binary mixtures at low pressure, establishes that the mutual diffusion coefficients are essentially independent of mixture composition.9 This assumption is valid for mixtures of ideal gases with spherical molecules, for which molecular interactions can be defined only by binary collisions and quantum effects are not relevant. For some gaseous mixtures, small variations from 0 to 5% have been observed in experimental measurements at atmospheric pressure.11 Such variations are supposed to increase when the pressure increases and the mixtures deviate from ideality. The composition dependence of the mutual diffusion coefficient observed in Figure 3 may be explained on the basis of different factors. First, it could be a consequence of the molecular geometry of SF6 and N2. Since these molecules are not entirely spherical and have different sizes, their mixtures behave slightly far from ideal. Estimations of the Q factor listed in Table 2 confirm that under the studied conditions, SF6/N2 mixtures slightly differ from the ideal behavior. A second factor which could explain the composition dependence observed in the diffusion coefficients is the presence of molecular flexibility in the SF6 force field. In order to analyze this last hypothesis, a different set of MD simulations was carried out using simple LJ potentials for both SF6 and N2.38,39 These simulations were run for three different compositions, 10, 50, and 90%, with the same simulation details explained in section 5 and with the standard Lorentz-Berthelot combining rules for crossed intermolecular interactions. Figure 4 depicts the mutual diffusion coefficients obtained with both sets of models, the anisotropic force fields described in section 4 (the same data depicted in Figure 3), and the LJ potentials. The coefficients estimated in the last case have been listed in Table 3. In simulations with the LJ models, the mutual diffusion
Transport Properties in Gaseous Mixtures
J. Phys. Chem. C, Vol. 111, No. 43, 2007 16019
Figure 5. Mutual diffusion coefficients obtained with MD simulations versus temperature for a SF6/N2 mixture with an SF6 concentration of 25% by volume: circles, simulations at 1 MPa; squares, simulations at 2 MPa; dashed lines, estimations with the Fuller et al. correlation.
TABLE 4: Mutual Diffusion Coefficient and Thermodynamic Factor Q Calculated from MD Simulations for a Mixture with an SF6 Concentration of 25% by Volume MD P (MPa)
MD T (K)
Q factor
106 × DSF6-N2 (m2/s)
1.08 ( 0.64 0.95 ( 0.52 0.96 ( 0.48 1.10 ( 0.53 0.99 ( 0.45 1.94 ( 1.25 2.10 ( 1.40 2.16 ( 1.25 2.06 ( 1.09 2.01 ( 0.99
263 ( 3 285 ( 3 299 ( 3 322 ( 4 339 ( 4 264 ( 3 282 ( 3 301 ( 3 325 ( 4 341 ( 4
0.9663 0.9817 0.9794 0.9839 0.9883 0.9263 0.9251 1.0466 0.9693 0.9786
0.72 ( 0.02 0.97 ( 0.04 1.05 ( 0.08 1.01 ( 0.02 1.159 ( 0.015 0.38 ( 0.03 0.34 ( 0.02 0.51 ( 0.02 0.58 ( 0.02 0.58 ( 0.04
coefficient also decreases as the SF6 composition increases; nonetheless, there is a weaker dependence on composition than in the case of using anisotropic force fields. Although this result is not conclusive to quantify the effect of the flexibility in the SF6 force field on the estimation of the diffusion coefficients of the mixture, it reveals the importance of considering this feature when calculating transport properties. The temperature dependence of the mutual diffusion coefficient was also studied. With this aim, MD simulations were performed with a mixture of fixed SF6 composition (25% in volume) at two different pressure values (1 and 2 MPa) in the temperature range 260-340 K. Results from these simulations are depicted in Figure 5 and listed in Table 4. For comparison purposes, the values estimated with the correlation of Fuller et al. were also included in the figure. Again, a good agreement was found between the simulation results and the empirical estimation method. As expected, the mutual diffusion coefficient increases with temperature and decreases with pressure. Finally, the shear viscosities obtained from MD simulations of the anisotropic models, the simple LJ models, and estimation methods, at 300 K and 1 MPa, are shown in Figure 6. As for the mutual diffusion coefficients, several techniques are available to estimate the viscosity in pure gases and gaseous mixtures based on either the Chapman-Enskog theory or the law of corresponding states. Reference 9 presents a detailed review about this matter, recommending the use of Lucas’ method with gaseous mixtures at low pressures. This method calculates pseudocritical and other mixture properties, considering the composition and pressure dependencies on the mixture viscosity. Under the conditions investigated (300 K and 1 MPa), Lucas’ method predicts a maximum shear viscosity value around the
Figure 6. Shear viscosities obtained from MD simulations as a function of the SF6/N2 mixture composition at 300 K and 1 MPa: circles, simulations with the SF6 flexible model of Olivet and Vega and the N2 force field of Galassi and Tildesley; open squares, simulations done with simple Lennard-Jones force fields for both SF6 and N2; dashed line, shear viscosities predicted with the method of Lucas. Dotted lines are for guidance to the eyes.
SF6 concentration of 25%, in accordance with the behavior obtained from MD simulations (see Figure 6 and Table 2). Overall, the behavior between the viscosities estimated with the method of Lucas and the ones predicted by using the anisotropic force fields is consistent, while both differ by an average relative error of approximately 12%. On the contrary, results obtained from MD simulations with the simple LJ potentials, also presented in Figure 6 and Table 3, do not show the composition dependence predicted by the other two methods. Again, this last observation highlights the importance of considering anisotropic force fields, and especially flexible molecular features, when calculating transport properties of fluids. 7. Concluding Remarks We have presented and discussed here the predictions on the transport properties of the SF6/N2 mixtures as obtained by MD simulations with the Green-Kubo formalism. Three models have been used for this purpose: anisotropic force fields for both components, simple LJ force fields for the same components, and available empirical correlations to estimate the transport properties of gaseous mixtures at low pressures. Under the studied conditions, the empirical estimation methods available in the literature assume that the mutual diffusion coefficient do not depend on the mixture composition. In this work, the correlation of Fuller et al. was compared with the simulation results. In spite of yielding the same order of magnitude, the diffusion coefficients estimated from simulations exhibited a composition dependence, in which this property decreases as the SF6 content in the mixture increases. This behavior was explained on the basis of the molecular features of SF6 and N2. These molecules are not entirely spherical, and in the case of the SF6, the molecular flexibility strongly influences the behavior of transport properties. In fact, MD simulations with simple LJ models for the same fluids showed weaker composition dependence than that seen with anisotropic models. Differences between predictions with the anisotropic force fields and those with the simpler LJ potentials were more significant with the shear viscosity. While the viscosities estimated with the LJ potentials showed no dependence on composition, the estimations based on the anisotropic force fields exhibited a maximum value of the viscosity at a given
16020 J. Phys. Chem. C, Vol. 111, No. 43, 2007 composition. This behavior agreed with the predictions obtained by applying the Lucas empirical method for mixture viscosities. The results of this work reveal that molecular simulations can be used as a predictive tool for determining the transport coefficients in SF6/N2 gaseous mixtures. Moreover, further refinement on empirical transport coefficients estimation methods seems to be necessary, as well as obtaining more experimental data related to the transport coefficients of these mixtures. Regarding the SF6 force fields used in this work, the results presented here constitute a new proof of the advantage of fitting the force field parameters using both equilibria and transport properties. Such an approach seems to bring a greater transferability to the force field parameters over a wide range of conditions and properties. Currently, work is in progress on estimating the transport coefficients of SF6/N2 mixtures in confined media in order to study the separation processes of these compounds. Acknowledgment. This work is dedicated to Professor Keith E. Gubbins for his outstanding contribution to the application of molecular simulation tools to predict the behavior of experimental systems. Financial support in part has been provided by the Spanish Government (Project CTQ2005-00296/ PPQ), the Catalan Government (2005SGR-00288), Ormazabal Corporate Technology, and MATGAS 2000 A.I.E. Computational time provided by the Super Computing Center of Catalonia (CESCA) is also gratefully acknowledged. References and Notes (1) Christophorou, L. G.; Van Brunt, R. J. IEEE Trans. Electr. Insul. 1995, 2, 952. (2) Christophorou, L. G.; Olthoff, J. K.; Green, D. S. Gases for Electrical Insulation and Arc Interruption: Possible Present and Future AlternatiVe to Pure SF6; C 13.46:1425; NIST Technical Note 1425; U.S. Department of Commerce: Washington, DC, 1997. (3) Niemeyer, L. A Systematic Search for Insulation Gases and their Environmental Evaluation. In Gaseous Dielectrics VIII, 8th International Symposium on Gaseous Dielectrics, Virginia Beach, Va, 1998; Christophorou, L. G., Olthoff, J. K., Eds.; Plenum Press: New York, 1998; p 459. (4) Shiojiri, K.; Yamasaki, A.; Fujii, M.; Kiyono, F.; Yanagisawa, Y. EnViron. Prog. 2006, 25, 218. (5) Toyoda, M.; Murase, H.; Imai, T.; Naotsuka, H.; Kobayashi, A.; Takano, K.; Ohkuma, K. IEEE Trans. Power DeliVery 2003, 18, 442. (6) Yamamoto, O.; Takuma, T.; Kinouchi, M. IEEE Trans. Electr. Insul. Mag. 2002, 18, 32. (7) Murase, H.; Imai, T.; Inohara, T.; Toyoda, M. IEEE Trans. Dielectr. Electr. Insul. 2004, 11, 166. (8) Shiojiri, K.; Yanagisawa, Y.; Yamasaki, A.; Kiyono, F. J. Membr. Sci. 2006, 282, 442. (9) Poling, B. E.; Prausnitz, J. M.; O’Connell, J. P. The Properties of Gases and Liquids, 5th ed.; McGraw-Hill: New York, 2004.
Olivet and Vega (10) Perry’s Chemical Engineers’ Handbook, 8th ed.; Perry, R. H., Green, D. W., Eds.; Chemical Engineering Series; McGraw-Hill: New York, 1998. (11) Marrero, T. R.; Mason, E. A. J. Phys. Chem. Ref. Data 1972, 1, 3. (12) Kestin, J.; Khalifa, H. E.; Ro, S. T.; Wakeham, W. A. Physica 1977, 88A, 242. (13) Haghighi, B.; Djavanmardi, A. H.; Papari, M. M.; Najafi, M. J. Theor. Comput. Chem. 2004, 3, 69. (14) Olivet, A.; Duque, D.; Vega, L. F. J. Chem. Phys. 2005, 123, 194508. (15) Olivet, A.; Duque, D.; Vega, L. F. J. Appl. Phys. 2007, 101, 023308. (16) Olivet, A.; Vega, L. F. J. Chem. Phys. 2007, 126, 144502. (17) Jolly, D. L.; Bearman, R. J. Mol. Phys. 1980, 41, 137. (18) Schoen, M.; Hoheisel, C. Mol. Phys. 1984, 52, 33. (19) Stoker, J. M.; Rowley, R. L. J. Chem. Phys. 1989, 91, 3670. (20) Zhou, Y.; Miller, G. H. Phys. ReV. E 1996, 53, 1587. (21) Fernandez, G. A.; Vrabec, J.; Hasse, H. Int. J. Thermophys. 2004, 25, 175. (22) Fernandez, G. A.; Vrabec, J.; Hasse, H. Int. J. Thermophys. 2005, 26, 1389. (23) van de Ven-Lucassen, I. M.; Vlugt, T. J.; van der Zanden, A. J.; Kerkhof, P. J. Mol. Phys. 1998, 94, 495. (24) Dysthe, D. K.; Fuchs, A. H.; Rousseau, B. Int. J. Thermophys. 1998, 19, 437. (25) Rivera, J. L.; Alejandre, J. Fluid Phase Equilib. 2001, 185, 389. (26) Zhou, Z.; Todd, B. D.; Travis, K. P.; Sadus, R. J. J. Chem. Phys. 2005, 123, 054505. (27) Zhang, C.; Yang, X. Fluid Phase Equilib. 2005, 231, 1. (28) Zhang, L.; Wang, Q.; Liu, Y. C.; Zhang, L. Z. J. Chem. Phys. 2006, 125, 104502. (29) Dysthe, D. K.; Fuchs, A. H.; Rousseau, B.; Durandeau, M. J. Chem. Phys. 1999, 110, 4060. (30) Dysthe, D. K.; Fuchs, A. H.; Rousseau, B. J. Chem. Phys. 2000, 112, 7581. (31) Zhou, Y.; Miller, G. H. J. Phys. Chem. 1996, 100, 5516. (32) Kirkwood, J. G.; Buff, F. P. J. Chem. Phys. 1951, 19, 774. (33) Rappaport, D. C. The Art of Molecular Dynamics Simulations; Cambridge University Press: Cambridge, U.K., 1995; p 114. (34) Haile, J. M. Molecular Dynamics Simulation; Wiley: New York, 1992; p 277. (35) Pawley, G. S. Mol. Phys. 1981, 43, 1321. (36) Galassi, G.; Tildesley, D. J. Mol. Simul. 1994, 13, 11. (37) Sokhan, V. P.; Nicholson, D.; Quirke, N. J. Chem. Phys. 2004, 120, 3855. (38) Hellemans, J. M.; Kestin, J.; Ro, S. T. J. Chem. Phys. 1972, 57, 4038. (39) Aziz, R. A.; Slaman, M. J.; Taylor, W. L.; Hurly, J. J. J. Chem. Phys. 1991, 94, 1034. (40) Nowak, P.; Kleinrahm, R.; Wagner, W. J. Chem. Thermodyn. 1997, 29, 1157. (41) www.ccp5.ac.uk/DL_POLY/. (42) Smith, W.; Yong, C. W.; Rodger, P. M. Mol. Simul. 2002, 28, 385. (43) Nose`, S. J. Chem. Phys. 1984, 81, 511. (44) Rey-Castro, C.; Vega, L. F. J. Phys. Chem. B 2006, 110, 16157. (45) Fuller, E. N.; Ensley, K.; Gidding, J. C. J. Phys. Chem. 1969, 73, 3679.
