404
Ind. Eng. Chem. Res. 2003, 42, 404-407
GENERAL RESEARCH Molecular Dynamics Calculation of the Diffusivity of Sodium Chloride in Steam Allan H. Harvey Physical and Chemical Properties Division, National Institute of Standards and Technology, Boulder, Colorado 80305-3328
Raymond D. Mountain* Physical and Chemical Properties Division, National Institute of Standards and Technology, Gaithersburg, Maryland 20899-8380
Molecular dynamics simulation is used to calculate the diffusivity of sodium chloride at infinite dilution in high-temperature steam at seven state points representing conditions of interest for the study of impurity deposition in steam turbines. Water is modeled with an existing four-site, polarizable potential, and the ions are modeled as charged, polarizable Lennard-Jones spheres. At the conditions studied, the sodium and chloride ions diffuse together; the simulations determine the diffusivity of this ion pair to within approximately 20%. The results can be fitted satisfactorily to the form of a simple kinetic-theory expression, allowing for extrapolation to lower densities that are less amenable to simulation. At turbine conditions, no experimental diffusivity data are available, and existing semiempirical estimation techniques are unreliable, so the results from this work provide the best values available for this industrially important property. 1. Introduction A large fraction of worldwide electricity generation involves the generation of steam (either directly in a boiler or in a heat recovery system following a combustion turbine) and passage of the steam through turbines. Because of the potential for damage to the turbine from corrosion or other consequences of deposited impurities, enormous efforts are made to purify the water and eliminate the possibility of deposition in the turbine sections of plants. Costs associated with impurities in water and steam for power generation have been estimated at $200 million per year in the U.S. alone.1 The impurities that cause the most serious problems in steam turbines are sodium chloride, sodium hydroxide, and silica. Greater efficiency and reliability can be achieved if the solubility and precipitation behavior of these solutes in steam is better understood. One of us (A.H.H.) was involved in a study of available solubility data for important salts and minerals in steam, with the goal of improving steam purity recommendations.2,3 One finding of this study was that the solubility data (especially for NaCl) indicated that precipitation problems should exist at impurity levels lower than one would expect on the basis of turbine operating experience. This led Bellows4 to suggest, on the basis of a mass-transfer analysis, that the deposition of solids in turbines at typical conditions could be mass-transfer-limited. His analysis was handicapped by the total lack of data for the diffusivity of the solutes in steam; he relied on a hard-sphere calculation with estimated diameters for * Corresponding author. E-mail:
[email protected]. Tel.: 301-975-2484. Fax: 301-869-4020. 10.1021/ie020562m
the species. Such an approach cannot be expected to be quantitatively accurate for systems with strong specific interactions, such as water with salt and mineral species. Experimental measurement of these diffusivities is difficult. The only remotely relevant measurements known to us are those of Butenhoff et al.5 for the diffusion of NaNO3 in high-temperature, high-pressure water. Although these data are at conditions relevant for supercritical water oxidation and other hydrothermal reactions, the densities are at least an order of magnitude greater than those relevant for steam turbines. These data do help illustrate the point that better data and models are needed for these systems; the simple diffusivity calculation used by Bellows4 deviates from the data of Butenhoff et al. by approximately a factor of 6.6 Alternatively, diffusivities can be calculated from molecular dynamics simulations with an accuracy limited only by computer resources and the quality of the intermolecular potentials. In this case, the lack of experimental data and the inapplicability of existing models means that modest errors (on the order of 10% or even 50%) in the diffusivity due to imperfection in the potentials would still provide a significant improvement over existing knowledge. Of the three most important solutes mentioned above, NaCl is the simplest; in addition, there exist realistic potentials for systems of H2O molecules with Na+ and Cl- ions. We therefore chose to simulate the diffusion of NaCl in supercritical steam. There have been several studies of the mobility of infinitely dilute aqueous Na+ and Clions (simulated without the counterion present) at supercritical conditions for densities as low as 200
This article not subject to U.S. Copyright. Published 2003 by the American Chemical Society Published on Web 12/10/2002
Ind. Eng. Chem. Res., Vol. 42, No. 2, 2003 405 Table 1. Interaction Parameters site
(kJ‚mol-1)
σ (nm)
q
R (nm3)
O H M Na Cl
0.764 0 0 0.544 0.418
0.325 0 0 0.235 0.415
0 0.519 -1.038 1 -1
0 0 0.00144 0.00013 0.00337
kg‚m-3.7,8 The only previous study of NaCl diffusion in supercritical H2O is that of Koneshan and Rasaiah.9 The lowest density they considered was 175 kg‚m-3, which is substantially greater than the densities of interest in steam turbines, 30 kg‚m-3 or less. Although application to steam turbines was the original motivation for this work, the diffusion of sodium chloride and other salts in steam is also of importance in supercritical water oxidation (SCWO), an emerging technology for the destruction of hazardous wastes. Solid precipitation and deposition is a major problem for SCWO systems; better understanding of solute diffusion would help in the design of this environmental technology. 2. Models and Simulation Details The simulated system consists of 214 water molecules, 1 Na+ ion, and 1 Cl- ion. These objects are placed in a cubic simulation cell with periodic boundary conditions in all three dimensions. The size of the cell is set so that the desired density of the mixture is realized. The equations of motion for the particles are integrated using an iterated version of the Beeman algorithm10-12 with a time step of 1 fs. The water molecules are taken to be rigid objects, and the orientational degrees of freedom are described using quaternions.13-15 The desired temperature is maintained using thermostats as formulated by Martyna et al.16 The interactions between the water molecules are described using a four-site, polarizable model that has been found to provide reliable representations of the structure and energy of small clusters of water molecules.17 The hydrogen sites (H-sites) are located 0.095 72 nm from the oxygen site (O-site) and the H-O-H bond angle is 104.5°. The M-site, where the polarizability is placed, is on the bisector of the H-O-H angle and is located 0.0215 nm from the O-site in the direction of the H-sites. The O-site has a Lennard-Jones interaction. The charges on the M- and H-sites result in a permanent dipole moment equal to the moment of an isolated water molecule. The ions are treated as polarizable, charged particles that also interact via a Lennard-Jones interaction.18,19 The ion-water and unlike-ion Lennard-Jones parameters are obtained using the Lorentz-Berthelot combining rules.20 The potential parameters are listed in Table 1. When this model was used for dense fluid systems, the interactions, including the Coulomb interactions, were truncated at 0.9 nm.17 This approximation is plausible for a dense system where the fluid is effectively uniform for such separations of molecules. This is not the case in steam, so we carefully examined the interactions for larger separations. We determined the largest deviations of the interaction energy from zero for water-water and water-ion pairs for several separation distances. This was done by fixing the center-ofmass separations and examining several thousand possible relative orientations of the pair. We found that the maximum deviation of the potential energy, includ-
ing polarization energy, was less than 10 K (0.08 kJ‚mol-1) when the separation was greater than 1.6 nm. The maximum deviation from zero of the force on a molecule at this separation is at least an order of magnitude smaller than the individual site-site forces. For smaller separations, the maximum deviations of the potential energy and the force from zero are considerably larger. We truncated the interactions at 1.6 nm, thus avoiding the use of the computationally more expensive Ewald summation treatment of the longrange Coulomb interactions,21 as the deviations from zero of the forces on the molecules are acceptably small. We emphasize that this truncation is acceptable because of the cancellation of terms that occurs in the absence of large ion-ion separations and because of the relatively large size of the simulation cells, 6.32 nm for the highest density. The polarization-based interactions were generated by first solving for the induced dipole moments selfconsistently and then determining the charge-dipole and dipole-dipole interactions. The induced moments, b µi on molecule (or ion) i, are obtained as solutions to the set of equations
b µ i ) RiE B 0i + Ri
uik b µk ∑ k*i
(1)
where uik is the dipole tensor,22 Ri is the polarizability at site i, and E B 0i is the electric field at site i due to the charges on the other molecules (or ions). Six iterations of these equations is sufficient to obtain a self-consistent solution.23 The diffusivities, Ds, are obtained as the limiting slopes of the mean-square displacements of the particles as a function of time
Ds )
d 1 lim 〈|r b(t) - b r (0)|2〉 6 tf∞ dt
(2)
where the angular brackets indicate an average over time origins and over the number of particles.24 In these simulations, the long-time limit is 100 ps. There are 20 time origins, so the duration of a simulation is 2 ns. This choice of limiting time is based on our observation that the long-time-limit behavior for the lowest density considered is reached after about 30-50 ps. The initial conditions for these simulations are obtained from an earlier study of ions in supercritical water.25 The coordinates were scaled so that the low densities needed to model steam were realized. These scaled configurations are then equilibrated by allowing the system to evolve for at least 100 ps before any data are taken. These relatively long equilibration times are needed so that the ion pairs have the opportunity to coalesce. We find that this always happens and that the ions remain close together at the densities and temperatures of interest. The temperature, T, density, F, and pressure, p, of the simulated states are listed in Table 2. These states are a compromise between the range of temperatures and densities encountered in steam turbine operation and the problems involved in simulating very low density fluids. As the density of the fluid decreases, the “mean free time” between collisions increases, and therefore, the necessary duration of a simulation also increases. The lowest density considered here, 7.6 kg‚m-3 (or about 8 × 10-3 times the density of liquid water), is about the limit for reasonable run
406
Ind. Eng. Chem. Res., Vol. 42, No. 2, 2003
Table 2. Temperature, Density, Pressure, and State Label for the Thermodynamic States of Steam Examined Here Ta (K)
Fa (kg‚m-3)
pb (MPa)
label
573 673 673 673 773 773 773
13.1 25.6 13.1 7.6 25.6 13.1 7.6
3.17 7.14 3.86 2.29 8.51 4.51 2.66
1 2 3 4 5 6 7
a The temperature and density are input quantities for the simulation. b The pressures listed are not computed in the simulation, but instead are calculated as a function of the input temperature and density from ref 26.
times. As will be shown below, lower densities are not required as this density is effectively in the regime where kinetic theory provides an adequate description of diffusion once the effective collision parameter has been determined. These simulations provide an effective value for the collision parameter. As a test of the water model, we generated the diffusivity of steam at 673 K and 100 kg‚m-3 to compare with an experimental value.27 The simulation result, 340 × 10-9 m2‚s-1, is a bit larger than the experimental value of 290 × 10-9 m2‚s-1. The uncertainty in the experimental value is on the order of 10%. This suggests that our self-diffusion coefficients might be correspondingly high. If that is the case, our calculations provide an upper bound on the diffusivity. As a further test, we used molecular dynamics to estimate the diffusivity of NaNO3 in water at 673 K and 600 kg‚m-3 where an experimental value for the mutual diffusion coefficient has been determined. At the low concentrations involved, the mutual diffusion coefficient is effectively equal to the diffusivity obtained from the mean square displacement. A model potential for NaNO3 has been reported in the literature and is used here.28 Because this model is of the Buckingham (exp-6) form, we used a similar potential for water.29 These models do not include explict polarizability. When Kong combining rules30,31 for the unlike interactions are used, we find that the diffusivities for water and the ions are 57 × 10-9 and 24 × 10-9 m2‚s-1 respectively. These values are in good agreement with the experimental values of 57 × 10-9 m2‚s-1 for water27 and 21 × 10-9 m2‚s-1 for the ions.5 This exercise indicates that reliable estimates for diffusivities are possible using molecular dynamics simulations. 3. Hard-Sphere Model In the absence of simulation or experimental data, Bellows4 resorted to a hard-sphere model to estimate diffusion coefficients in his study of deposition in turbines. Because there are some misprints in his paper, we describe kinetic theory calculations here before comparing them with our simulation results. For a binary mixture of hard spheres, the diffusivity at low and moderate densities such as those considered in this study can be calculated with good accuracy from the kinetic theory derived by Chapman and Enskog.32 The diffusivity is given by 1/2 3 (πkBT/M12) D12 ) 8 nπσ 2 12
(3)
Table 3. Diffusivities of the Ions for the Seven States Described in Table 2, along with the Hydration Numbers, 〈N〉, and the Variances of the Hydration Numbers, σN, of the Ions label
DNaCla
〈N〉
σN
DCE
1 2 3 4 5 6 7
360 218 363 679 258 430 770
1.6 2.4 1.4 1.0 2.0 1.3 0.9
0.9 1.1 0.9 0.8 1.0 0.9 0.7
1180 655 1280 2200 701 1370 2360
a Diffusion coefficients are in units of 10-9 m2‚s-1. Simulation values are uncertain at the 20-30% level.
where M12 is the reduced mass defined by M12 ) 2[(1/M1) + (1/M2)]-1, M1 and M2 are the molecular masses of the two species, n is the number density of molecules in the mixture, and σ12 is the binary collision diameter, which is the arithmetic mean of the diameters of the two species. Bellows4 took the diameter of water molecules to be 0.272 nm. NaCl in steam is believed to be associated with approximately four water molecules, so he used an estimated diameter of 0.525 nm for this hydrated species, along with a mass corresponding to that of one NaCl molecule plus four H2O molecules. We used these parameters in eq 3 to produce the hardsphere values, DCE, reported in the next section. We also performed a simulation of a binary mixture with these parameters (using r-12 spheres rather than perfectly hard spheres) and verified that the Chapman-Enskog diffusivities and our simulations were consistent to within a few percent. 4. Results Two sets of quantities are extracted from the simulated molecule trajectories during the simulations. The first is the mean-square displacement over time intervals, and the second is the number of water molecules that are “solvating” the ions. A water molecule is counted if the distance between the Na+ ion and the oxygen site or the distance between the Cl- ion and a hydrogen site is 0.47 nm or less. This choice of solvation distance includes both “first” and “second” neighbor molecules.25 The diffusivities for NaCl using Chapman-Enskog theory with Bellows’ diameters,4 DCE, and from the simulations, DNaCl, are listed in Table 3 for the seven states listed in Table 2. The hydration of the ions is also summarized in Table 3. The average number of water molecules, 〈N〉, that satisfy the solvation condition is listed along with the rms deviation from the mean, σN. These numbers are based on samples taken every 0.1 ps (20 000 samples). The number of solvating water molecules varies from 0 to 2, except for states 2 and 5, where the number varies from 1 to 3. 5. Discussion It is interesting to note that the simple hard-sphere calculation consistently overpredicts the simulated diffusivities by roughly a factor of 3. This can be attributed to the binary collisions in this system being very different from hard-sphere collisions. In more rigorous kinetic theory,32 the denominator in eq 3 has an additional factor Ω(1,1)* known as a “collision integral”, which incorporates detailed scattering calculations for binary intermolecular collisions and which has a value
Ind. Eng. Chem. Res., Vol. 42, No. 2, 2003 407
of 1 for hard spheres. Although Ω(1,1)* is intractable for our anisotropic system, we can estimate it by determining an “effective temperature” for diffusion. We make a Boltzmann-weighted average over all orientations of the water-ion pair interaction energy as a function of separation and find the minimum of this averaged interaction (in units of the Boltzmann’s constant times the temperature, kBT) to be an attractive interaction with a strength of 7550 K at a separation of 0.265 nm. This provides a characteristic temperature scale for the diffusion process. Our simulation conditions therefore have a “reduced” temperature on the order of 0.1. For the simpler Lennard-Jones potential, Ω(1,1)* at a reduced temperature of 0.1 has a value of approximately 4.33 This suggests that the collision integral is the main explanation for the factor of 3 by which the hard-sphere calculations differ from the results of our simulations. The diffusivities of NaCl for seven states of steam have been generated. It is reasonable to ask whether this is sufficient information to estimate the diffusivity for other states of steam. Kinetic theory32 states that the diffusivity at low and moderate densities is given as
D ) AT1/2/n
(4)
where T is the temperature in Kelvin, n is the number density in units of m-3, and A incorporates information about the sizes of the molecules and the nature of their collisions. A is weakly temperature-dependent for real fluids because of the temperature dependence of the collision integral, but for our purpose of obtaining an approximate expression for steam turbine conditions, this dependence can be neglected. Using the simulation results for DNaCl, we find that A is 6.9 × 1019 m-1‚s-1 K-1/2 when averaged over the seven states. Using this value for A reproduces the simulation results within (15% for all states. This estimate for A does not take into account the differing degrees of hydration of the NaCl pair, as it is not obvious how that might be done. The main result of this study is the expression for DNaCl in terms of the quantity A, the temperature, and the number density. This expression should provide adequate estimates for the diffusivity needed for turbine deposition studies. Acknowledgment We thank Dr. J. C. Bellows for bringing this problem to our attention and for helpful discussions. Literature Cited (1) Jonas, O. Cost of “Corrosion and Scale” in U.S. Utilities. In Proceedings, 1985 Fossil Plant Water Chemistry Symposium; Electric Power Research Institute: Palo Alto, CA, 1986. (2) Harvey, A. H.; Bellows, J. C. Evaluation and Correlation of Steam Solubility Data for Salts and Minerals of Interest in the Power Industry, NIST Technical Note 1387; NIST: Gaithersburg, MD, 1997. (3) Bellows, J. C.; Harvey, A. H. Steam Solubilities for Combustion Turbine Steam Cooling. Int. J. Thermophys. 1999, 20, 197. (4) Bellows, J. C. Deposition Rates of Impurities from Steam. In Steam, Water, and Hydrothermal Systems: Physics and Chemistry Meeting the Needs of Industry; NRC Press: Ottawa, Canada, 2000; pp 883-890. (5) Butenhoff, T. J.; Boemans, M. G. E.; Buelow, S. J. Mass Diffusion Coefficients and Thermal Diffusivity in Concentrated Hydrothermal NaNO3 Solutions. J. Phys. Chem. 1996, 100, 5982. (6) Bellows, J. C. Siemens Westinghouse Power Corp., Orlando, FL. Personal communication, 2001.
(7) Lee, S. H.; Cummings, P. T.; Simonson, J. M.; Mesmer, R. E. Molecular dynamics simulation of the limiting conductance of NaCl in supercritical water. Chem. Phys. Lett. 1998, 293, 289. (8) Rasaiah, J. C.; Noworyta, J. P.; Koneshan, S. Structure of Aqueous Solutions of Ions and Neutral Solutes at Infinite Dilution at a Supercritical Temperature of 683 K. J. Am. Chem. Soc. 2000, 122, 11182. (9) Koneshan, S.; Rasaiah, J. C. Computer simulation studies of aqueous sodium chloride solutions at 298 K and 683 K. J. Chem. Phys. 2000, 113, 8125. (10) Schofield, P. Computer simulation studies of the liquid state. Comput. Phys. Commun. 1973, 5, 17. (11) Beeman, D. Some multistep methods for use in molecular dynamics simulations. J. Comput. Phys. 1976, 20, 130. (12) Mountain, R. D.; Brown, A. C. Molecular dynamics study of the liquid and plastic phases of neopentane. J. Chem. Phys. 1985, 82, 4236. (13) Evans, D. J.; Murad, S. Singularity free algorithm for molecular dynamics simulation of rigid polyatomics. Mol. Phys. 1977, 34, 327. (14) Sonnenschein, R. An Improved Algorithm for Molecular Dynamics Simulation of Rigid Molecules. J. Comput. Phys. 1985, 59, 347. (15) Rapaport, D. C. Molecular Dynamics Simulation Using Quaternions. J. Comput. Phys. 1985, 60, 306. (16) Martyna, G. J.; Klein, M. L.; Tuckerman, M. Nose´-Hoover chains: The canonical ensemble via continuous dynamics. J. Chem. Phys. 1992, 97, 2635. (17) Dang, L. X.; Chang, T.-M. Molecular dynamics study of water clusters, liquid, and liquid-vapor interface of water with many-body potentials. J. Chem. Phys. 1997, 106, 8149. (18) Smith, D. E.; Dang, L. X. Computer simulations of NaCl association in polarizable water. J. Chem. Phys. 1994, 100, 3757. (19) Pyper, N. C.; Pike, C. G.; Edwards, P. P. The polarizabilities of species present in ionic solutions. Mol. Phys. 1992, 76, 353. (20) Hansen, J. P.; McDonald, I. R. Theory of Simple Liquids; Academic Press: New York, 1986; p 179. (21) Sangster, M. J. L.; Dixon, M. Interionic potentials in alkali halides and their use in simulations of the molten salts. Adv. Phys. 1976, 25, 247. (22) Vesley, F. J. N-particle Dynamics of Polarizable Stockmayer-type Molecules. J. Comput. Phys. 1977, 24, 361. (23) Mountain, R. D. Comparison of a fixed-charge and a polarizable water model. J. Chem. Phys. 1995, 103, 3084. (24) Frenkel, D.; Smit, B. Understanding Molecular Simulation, From Algorithms to Applications, 2nd ed.; Academic Press: New York, 2002; pp 100-102. (25) Mountain, R. D. Molecular dynamics study of ions in water. In Steam, Water, and Hydrothermal Systems: Physics and Chemistry Meeting the Needs of Industry; NRC Press: Ottawa, Canada, 2000; pp 509-516. (26) Harvey, A. H.; Peskin, A. P.; Klein, S. A. NIST/ASME Steam Properties; NIST Standard Reference Database 10, version 2.2; National Institute of Standards and Technology: Gaithersburg, MD, 2000. (27) Lamb, W. J.; Hoffmann, G. A.; Jonas, J. Self-diffusion in compressed supercritical water. J. Chem. Phys. 1981, 74, 6875. (28) Lyden-Bell, R. M.; Ferario, M.; McDonald, I. R.; Salje, E. A molecular dynamics study of orientational disordering in crystalline sodium nitrate. J. Phys.: Condens. Matter 1989, 1, 6523. (29) Errington, J. R.; Panagiotopoulos, A. Z. A Fixed Point Charge Model for Water Optimized to the Vapor-Liquid Coexistence Properties. J. Phys. Chem. B 1998, 102, 7470. (30) Kong, C. L. Combining rules for intermolecular potential parameters. II. Rules for the Lennard-Jones (12-6) potential and the Morse potential. J. Chem. Phys. 1973, 59, 2464. (31) Potoff, J. J.; Errington, J. R.; Panagiotopoulos, A. Z. Molecular simulation of phase equilibria for mixtures of polar and nonpolar components. Mol. Phys. 1999, 97, 1073. (32) Hirschfelder, J. O.; Curtiss, C. F.; Bird, R. B. Molecular Theory of Gases and Liquids; John Wiley & Sons: New York, 1954. (33) Klein, M.; Smith, F. J. Tables of Collision Integrals for the (m,6) Potential Function for 10 values of m. J. Res. NBS 1968, 72A, 359.
Received for review July 29, 2002 Revised manuscript received November 5, 2002 Accepted November 5, 2002 IE020562M
