Tracing the Critical Loci of Binary Fluid Mixtures Using Molecular

Aug 4, 2006 - loci for the methane-ethane system using no adjustable parameters for ... of molecular models, binary Lennard-Jones mixtures.20-25 Phase...
1 downloads 0 Views 102KB Size
Download PDF
17200

J. Phys. Chem. B 2006, 110, 17200-17206

Tracing the Critical Loci of Binary Fluid Mixtures Using Molecular Simulation Philip J. Lenart and Athanassios Z. Panagiotopoulos* Department of Chemical Engineering, Princeton UniVersity, Princeton, New Jersey 08544 ReceiVed: May 19, 2006; In Final Form: June 26, 2006

Monte Carlo simulations are used to trace the critical loci for a number of binary mixtures. In particular, we use grand canonical Monte Carlo (GCMC) simulations with histogram reweighting and mixed-field finitesize scaling to determine the mixture critical lines. Two different classes of criticality are investigated. A mixture of methane and ethane displays type I criticality, exhibiting continuous mixing between the two species across the entire composition range. A methane-water mixture shows type IIIb criticality, with a discontinuity in the critical locus. Quantitative agreement is found between simulation and experimental critical loci for the methane-ethane system using no adjustable parameters for interactions in the mixture. For the water-methane system, we investigate the effect of the combining rules for the intermolecular interaction between the two species on the mixture critical locus. We also investigate several potentials for methane: a nonpolar exponential-6, an octopolar fixed partial charge, and a polarizable fluctuating charge model. Qualitative agreement between simulations and experiments is found for all potentials, but none are able to quantitatively capture the abrupt increase in the critical temperature as methane is added to the system.

1. Introduction Molecular simulation has been used widely for the prediction of thermodynamic and transport properties of fluids and their mixtures. The vapor-liquid phase coexistence properties of binary fluid mixtures have been studied extensively using equation of state approaches and simulation techniques including Gibbs ensemble Monte Carlo1 and grand canonical Monte Carlo. Comprehensive reviews of the experimental determination and equation of state predictions of critical loci have been published.2-4 Most simulation studies of mixture phase equilibrium have been restricted to subcritical temperatures. However, there is much interest in the behavior of fluid mixtures near the critical region. Applications of supercritical fluids include separation processes5,6 and their use as media for chemical reactions.7 Binary fluid mixtures exhibit many different types of behavior in the critical region. A useful classification scheme is that of Van Konynenburg and Scott,8 who define six main classes of mixture critical behavior. Mixtures of chemically similar components such as alkanes display type I behavior, where there is continuous mixing of the fluid properties bounded by the critical points of the pure fluids. Differences in the interaction energy between the species and the molecular shapes result in behavior such as critical end points at intermediate compositions and discontinuities in the critical locus. Using the van der Waals equation of state, Van Konynenburg and Scott were able to construct all but one type of the critical loci that binary fluids exhibit experimentally.8 Several studies have been performed using more advanced equations of state for mixtures of alkanes, including the simplest one, methane-ethane.9-12 There is also interest in aqueous mixtures involving hydrocarbons in the critical region.13-19 The advantage to using an equation of state approach is that it is computationally less expensive compared to simulation techniques. Equations of state have a number of pure component and mixture parameters that must be deter* Corresponding author. E-mail: [email protected].

mined. To improve agreement with experimental results, parameters are often fitted to experimental data, reducing the predictive ability of the models. The critical region of pure fluids and fluid mixtures has also been investigated using molecular-based computer simulations. Several simulation studies have been performed for the simplest of molecular models, binary Lennard-Jones mixtures.20-25 Phase equilibrium simulations for systems with more realistic models including alkanes and water is an ongoing area of research. There are numerous alkane force fields optimized for phase coexistence properties.26-29 Fixed-charge,30 polarizable,31 and quantum mechanical32 models have been used to study the coexistence properties of pure water. Thermodynamic and structural properties of water and aqueous mixtures have been studied via simulation from subcritical to supercritical states.33-36 Although the critical point has been determined for pure fluids and for several binary mixtures, critical loci have not been studied extensively by computer simulation. A previous attempt to explicitly trace the critical locus of a binary polymer-solvent was performed using techniques similar to the ones in the present work.37,38 This mixture exhibits type III criticality with a pressure minimum and a discontinuity very close to the critical temperature of the solvent (CO2). However, the simulations did not observe the pressure minimum due to the onset of liquid-liquid immiscibility. It was shown that the value of the energy mixing parameter can alter the classification of the critical locus. A similar approach was used to predict the type I behavior of model electrolyte-salt mixtures.39 In this paper, we investigate the predictive ability of molecular simulations for the critical loci of binary mixtures. For this purpose, we study two systems of significantly different character: a nonpolar mixture of methane and ethane (type I), and a polar-nonpolar mixture of water and methane (type IIIb). In section 2, we present the potentials studied for these binary mixtures, the forms for the combining rules describing the interaction between unlike molecules, and the fluctuating charge formalism for introducing polarization into the methane poten-

10.1021/jp0630931 CCC: $33.50 © 2006 American Chemical Society Published on Web 08/04/2006

Tracing Critical Loci of Binary Fluid Mixtures

J. Phys. Chem. B, Vol. 110, No. 34, 2006 17201

tial. The details of the simulations, including the fluctuating charge ANES procedure, are described in section 3. The predicted critical loci and a discussion of these results are presented in section 4, and conclusions are presented in section 5. 2. Intermolecular Potential Models The two mixtures studied in this work are selected based on the availability of accurate intermolecular potentials describing properties of the pure fluids. The two mixtures also have different behavior along the critical locus. The first mixture simulated is that of methane and ethane. This mixture exhibits type I criticality.40 This means that there is continuous mixing of the fluids across the entire composition range (pure methane to pure ethane). The second mixture examined is watermethane, which exhibits type IIIb behavior.41 This behavior is characterized by a discontinuity in the temperature-composition critical locus. As methane is added to water existing at its critical point, the critical temperature reaches a minimum before rapidly increasing. For all of the systems studied in this work, the dispersion interaction is modeled with the modified Buckingham exp-6 potential42 of the form

{

Udisp(r) )

[ ([

]) ( ) ]

r  6 exp R 1 R-6 R rm

-

rm r

6

for r > rmax

∞

(1)

for r e rmax

where  is the minimum potential energy, R is the repulsive steepness of the potential, and rm is the distance at which the potential is a minimum. A hard wall is introduced at rmax where the potential reaches an unphysical maximum at short distances. The characteristic size parameter σ, where the potential energy is equal to zero, is numerically determined for this potential. The interactions between different species are described by a set of combining rules for the parameters of the intermolecular potential. For the methane-ethane and water-methane mixtures, the Lorentz-Berthelot descriptions of the energy, repulsion factor, and size parameters are

ij ) xiijj

(2)

Rij ) xRiiRjj

(3)

σij ) 0.5(σii + σjj)

(4)

In addition to the Lorentz-Berthelot rules, two additional sets of combining rules were employed for the water-methane system. For the exp-6 potential, the Kong combining rules43 are

]

ijRij exp(Rij) (Rij - 6)σij

[

2σij/Rij

[

)

] [

iiRii exp(Rii) (Rii - 6)σii

[

σii/Rii

]

jjRjj exp(Rjj) (Rjj - 6)σjj

]

iiRiiσii6 jjRjjσjj6 ijRijσij6 ) Rij - 6 (Rii - 6) (Rjj - 6)

σjj/Rjj

(5)

1/2

(6)

]

σij 1 σii σjj ) + Rij 2 Rii Rjj

(7)

The unlike parameters ij, σij, and Rij are determined by solving the system of equations expressed in eqs 5-7 and are listed in Table 2. The Sadus combining rule for the molecular size was developed for modeling type III mixtures.44 The ij and Rij parameters are calculated using eqs 2 and 3. Correcting for a typographical error in the original paper, the expression for the σij parameter is

σij3 )

1 (σii + σjj)2(σii3 + σjj3)1/3 3 4x2

(8)

which represents a 2:1 geometric averaging of the Lorentz and arithmetic combining rules. The models for the alkanes were developed by Errington and Panagiotopoulos.26,28 These united-atom models consist of an exp-6 site representing each -CH3 and -CH2 group. These potentials were parametrized to match the critical point, vaporliquid coexistence densities, and the coexistence vapor pressure. The potential parameters are found in Table 1. The phase envelope for methane and ethane is presented in Figure 1 for comparison with experimental results. The water model (EP) used in this study was also developed by Errington and Panagiotopoulos.30 It has a exp-6 interaction site located on the oxygen atom, and three partial charges representing the atoms. The model was parametrized to reproduce the critical point and the vapor-liquid coexistence properties of water as seen in Figure 2. Also seen in this figure are the phase diagrams for two other water potentials, that for SPCE45 and the recently developed polarizable GCPM model.31 The EP model is superior in its ability to capture the vaporphase coexistence properties and the critical point over the SPCE model, but both suffer from underprediction of the liquid coexistence densities. The GCPM model is able to reproduce the phase envelope over an extensive temperature range, but since this potential is polarizable, simulations using this model require significantly more computational resources. We use the EP model due to its improved agreement with the experimental critical point for water over other nonpolarizable water models while requiring less computational resources. The intermolecular interaction between charge sites Uqq is described by Coulomb’s law

Uqq(rij) )

and

[

and

1 qiqj 4π0 rij

(9)

where qi is the charge on site i, rij is the distance between the charges, and 0 is the permittivity of vacuum. For the watermethane mixture, the total energy of the system is given as the sum of the dispersion and charge interactions between all pairs of molecules, Usystem ) ∑Udisp + ∑Uqq. It is our interest to determine the effect of the polarization of methane into the water-methane system. The dipolar nature of water creates an electric field which polarizes the methane, inducing a dipole on the molecule. Since the dipole moment of water is described by a set of partial atomic charges distributed over the atomic sites, we formulated a partial-charge model for methane. This methane model is a tetrahedral arrangement of four charges on the hydrogen atoms and one charge on the central carbon, yielding an octopole moment on the molecule. The molecular geometry and partial charge values are taken

17202 J. Phys. Chem. B, Vol. 110, No. 34, 2006

Lenart and Panagiotopoulos

TABLE 1: Potential Parameters for the exp-6 Interactions of Methane, Ethane, and Water species

/kB (K)

σ (Å)

R

methanea (no charges) methane (partial charges) ethaneb waterc

160.3 157.9 129.6 159.78

3.741 3.778 3.679 3.195

15 15 16 12

a

Reference 26. b Reference 28. c Reference 30.

U(q) )

TABLE 2: Unlike Interaction Parameters of the Water-Methane Mixture for Different Combining Rules combining rules

WM/kB (K)

σWM (Å)

RWM

Lorentz-Berthelot Kong Sadus

160.04 149.25 160.04

3.468 3.518 3.475

13.42 13.86 13.42

TABLE 3: Parameters for the Octopolar and Polarizable Methane Models parameter

value

qCgp (e) -0.24 qHgp (e) +0.06 rCH (Å) 1.1 ∠HCH (deg) 109.47 χH0 - χC0 (kcal/mol e) -14.158

Table 1) is needed to obtain good agreement with the experimental phase diagram. The polarization of methane is introduced via fluctuations in the magnitudes of the partial atomic charges. The electrostatic energy of a single methane molecule with a set of charges q is given by47

parameter

value

JHH (kcal/mol e2) JCC (kcal/mol e2) JHH(rHH) (kcal/mol e2) JHC(rHC) (kcal/mol e2)

279.143 247.107 148.099 184.486

al.46

from the OPLS-AA model of Jorgensen et and the parameters are listed in Table 3. The effect of adding partial atomic charges on the vapor-liquid-phase coexistence properties of methane is seen in Figure 1. Since the electrostatic effects are of such a high order (octopolar), a minimal amount of reparametrization of the nonpolar exp-6 parameters (i

]

(10)

where Nq is the number of charge sites in the molecule, Ei0 is the static electric field, χi0 is the electronegativity, and Jii is the hardness of site i. To determine the charge magnitudes of the gas-phase molecule, a minimization of eq 10 is performed with respect to the atomic charges. Using symmetry for an isolated molecule, qC ) -4qH, one obtains the expression for the partial atomic charge of hydrogen qgp H as

qgp H )

-(χH0 - χC0)

(11)

JHH + 4JCC + 3JHH(rHH) - 8JHC(rHC)

where JHH and JCC represent the atomic hardness and JHH(rHH) and JHC(rHC) represent the intramolecular electrostatic interaction between charges on different atoms. The J values used to evaluate eq 11 were based on the fluctuating charge model of methane to study methane adsorption in zeolites using MD simulations.48 The hardness parameters Jij(rij) in this model are described by a screened Coulombic interaction. The difference in electronegativity was parametrized to reproduce the gas-phase OPLS-AA charge on the hydrogen atom (+0.06 e). The fluctuating charge parameters are presented in Table 3. For the water-polarizable methane system, the total system energy is described as

[ ∑∑[

NW+NM

Usystem ) 1 4π0

Figure 1. Phase diagrams for methane and ethane: (top) ethane exp-6 model (circles) and experimental data61 (solid line); (bottom) methane exp-6 (squares) and octopolar models (triangles) compared to experimental data61 (solid line).

Nq

∑i ∑ j>i

Udisp(rij) +

NM NqM

χm0qim +

i

m)1

1

1 4π0

NqM

Nqi

Nqj

]

im,jn

∑ Jmn(rim,jn)qimqjn

2 n)1

]

qimqjn

∑∑ m)1 n)1 r

+

- NMUgp (12)

where NW and NM are the number of water and methane molecules and Nqi is the number of charges in species i. The second line of eq 12 represents the polarization energy of the system and Ugp is the gas-phase energy of methane determine by eq 10. Details for carrying out the MC simulation using fluctuating atomic charges are provided later in the paper. 3. Simulation Details

Figure 2. Phase diagram for water. Results for the EP (square and dashed line), SPCE (triangles), and GCPM31 models (circles) compared to experimental data61 (star and solid line).

Grand canonical Monte Carlo (GCMC) simulations were performed for the methane-ethane and water-methane systems. In this ensemble, the chemical potentials of both species µ1 and µ2, the system volume V, and the temperature T are held constant. During the course of the simulation, fluctuations in the system density occur due to the insertion and removal of particles in the simulation volume. Particles in the simulation box are also displaced and rotated. The mix of moves for all simulations was 35% translation, 15% rotation, and the balance attempted insertions or deletions selected with equal probability. For the methane-ethane, a simulation box of volume V ) 15000 Å3 was used, while a smaller box V ) 5000 Å3 was used for the water-methane systems. These volumes were selected

Tracing Critical Loci of Binary Fluid Mixtures

J. Phys. Chem. B, Vol. 110, No. 34, 2006 17203

Figure 3. Ordering operator distribution PL(x) for a methane-ethane mixture (xM ∼ 0.51) at T ) 260 K and V ) 15 000 Å (circles) as compared to the three-dimensional Ising distribution54 (solid line).

because they represent the system size for which the ethane and water potentials were parametrized. Lennard-Jones interactions were truncated at half of the simulation box length, and a tail correction was implemented through the algorithm of Theodorou and Suter.49 The long-range electrostatic interactions were calculated using the point charge form of the Ewald summation50,51 with insulating boundary conditions (s ) 1). A total of 276 k-vectors were used in the reciprocal-space calculation and a real-space dampening factor κ of 6.0 for the water-methane mixture. Histogram reweighting techniques52,53 were used to determine the vapor and liquid compositions and pressure at coexistence for the binary mixtures using the observations of the particle number for each species and the system energy from a limited number of simulations. The value of the grand partition function Ξ was calculated from the histogram observations and is known to within an arbitrary constant. The pressure of the system was determined using its thermodynamic relation to Ξ

βPV ) ln Ξ + c

(13)

where β is the inverse temperature. Using the ideal gas law βPV ) N, a plot of ln Ξ vs N is used to determine the intercept c under ideal gas conditions (F f 0). To determine the critical points, we used the mixed-field finite-size scaling methods of Wilding and Bruce54 assuming that the fluid mixture follows Ising criticality. Figure 3 presents a match of our simulation results for the methane-ethane mixture to the universal order parameter curve. From this match, we were able to obtain values for the critical pressure and the critical composition at a fixed temperature. The statistical uncertainties for the constructed critical loci were obtained using three different random number seeds to generate three independent critical point estimates. Results were averaged, and the statistical uncertainty was calculated as the standard deviation between these three sets of runs. The critical locus presented for each system is for the single volume size listed above. To test system-size effects, the volume of the methane-ethane mixture was doubled, increasing the box length about 25%. Little change was seen in the temperature-composition critical locus, but an increase of ∼5% in the pressure-temperature critical locus was seen (results not presented). This discrepancy is attributed to the error between the experimental and simulated critical pressures of the pure fluids. Reparametrization of the pure fluid potentials for the larger system volume would bring the predicted results in agreement with the initial simulation results. Fluctuating Charge ANES. The incorporation of polarization effects via the adiabatic nuclear and electronic sampling (ANES)

algorithm has been described for single-component and fluid mixtures.55-58 For the polarizable methane molecules, the electronic states are characterized by the atomic charge magnitudes. The main element of this algorithm is that nuclear moves are made at a system temperature T, and for each nuclear move a number of electronic moves are made at an “electronic” temperature Telec , T. During the course of the Monte Carlo simulation, for each nuclear move (particle translation, rotation, insertion, or deletion depending on the ensemble) there is a number of perturbations Nfluc made in the charge magnitude on different sites in the system. Individual electronic moves are accepted/rejected based on standard MC acceptance criteria using Telec and then the overall nuclear + electronic moves are accepted/rejected based on T. If a nuclear move is rejected, then all nuclear and electronic coordinates must be returned to their values before the attempted nuclear move was made. We were able to exploit many features of the water-methane system when incorporating polarization into the methane potential. First, the density of the system, consisting primarily of water, is much less than that of liquid water at ambient temperature (compare 0.322 g/cm3 at 647 K to 0.997 g/cm3 at 298 K). This results in a decreased electric field at the methane site in critical water than in liquid water at a lower temperature, decreasing the number of fluctuation steps Nfluc needed to equilibrate the charges in the molecule. The procedure for making a change in the magnitude of the methane charges is as follows. First, a random charge site i is selected among the five charges in the methane molecule. A perturbation in the charge ∆q is generated as ∆q ) qmaxξ, where qmax is the maximum charge perturbation magnitude and ξ is a random number on the interval [-1,1]. The selected site charge is then adjusted as qi ) qi + ∆q. To maintain charge neutrality in the molecule, the four other charge sites j are adjusted as qj ) qj - 0.25∆q. All of these steps are performed during each charge fluctuation step. For ANES-MC, a number of parameters needed to be determined. These include the electronic temperature Telec, the number of electronic moves Nfluc made during each MC step, the maximum charge perturbation qmax, and the maximum displacement of molecules made during translational moves. In previous studies of water56,57 and a methanol-methane mixture,58 Telec is chosen to be 5 K. The same value was used in this study, as Telec is much less than the temperature range of interest (∼650 K). We also selected qmax to be ∼10% of the gas-phase charge magnitude of the hydrogen atom (qmax ) 0.005). Convergence of the system energy for a mixture of water and polarizable methane was found for Nfluc ) 15. Finally, the maximum displacement of the species during translation moves was adjusted until 30 of displacement moves were accepted during the equilibration period of the simulation. 4. Results and Discussion Methane-Ethane Mixture. The Tc-xM projection of the critical locus is presented in Figure 4. The type I behavior of this mixture is seen in the continuous curve that links the critical temperature of ethane (302.5 K) to that of methane (195.5 K) over the entire composition range. As seen, there is very good quantitative agreement with the experimental results. Figure 5 presents the critical locus for the methane-ethane mixture in the Pc-Tc plane. Once again, very good agreement is found between the simulation using the Lorentz-Berthelot combining rules and the experimental results. A pressure maximum is predicted at a mixture composition of xM ≈ 0.57 using a quadratic fit to the simulation results.

17204 J. Phys. Chem. B, Vol. 110, No. 34, 2006

Lenart and Panagiotopoulos

Figure 4. Temperature-composition critical locus for a methaneethane mixture: simulation results (open circles) and experimental data40 (filled circles). Simulation statistical uncertainties are smaller than the symbol size. Line is drawn through the experimental data to guide the eye.

Figure 6. Temperature-composition critical locus for a watermethane mixture for different combining rules. Simulation results for Lorentz-Berthelot (triangles), Kong (squares), and Sadus rules (diamonds) compared to experimental data41 (filled circles). Simulation statistical uncertainties are smaller than the symbol size. Lines are drawn to guide the eye.

Figure 5. Pressure-temperature critical locus for a methane-ethane mixture. Symbols and line are the same as those given in Figure 4. Simulation statistical uncertainties are smaller than the symbol size.

Figure 7. Pressure-temperature critical locus for a water-methane mixture for different combining rules. Symbols are the same as those given in Figure 6. Simulation statistical uncertainties that are smaller than the symbol size are not shown. Lines are drawn to guide the eye.

Due to the similarity of the potential parameters, the unlike interaction parameters from other combining rules may provide the same result. Applying the Kong combining rules of eqs 5-7, the results are found to differ by less than 0.4% for the energy parameter and less than 0.1% for the size and repulsive parameters. These small differences are expected to have a negligible effect on the results of Figures 4 and 5. The equivalence of results obtained using different combining rules has been shown in the simulation of an alkane mixture at subcritical temperatures.43 Water-Methane Mixture. Now that we are confident that our method allows us to trace the critical locus for a binary fluid mixture, we studied the phase behavior of the more complex water-methane system. In this type IIIb mixture, the critical locus originating from the critical point of methane ends in an upper critical end point. However, due to a lack of experimental data in this region, we restricted our study to the water-rich region of the phase diagram for comparison with existing experimental data. The type IIIb behavior of the watermethane mixture is seen in the experimental results in Figure 6. As methane is added to the system, there is a decrease in the critical temperature of the system until it reaches a minimum around xW ≈ 0.73. Further introduction of water results in a sharp increase in the critical temperature. Our initial investigation of the water-methane system was to determine the effect of the combining rules used for unlike intermolecular interactions on the critical locus. For a simple model for methane (a single exp-6 interaction site), a comparison of the Lorentz-Berthelot, Kong, and Sadus combining rules on the Tc-xW critical locus is presented in Figure 6. The

differences between the experimental and simulation results for pure water are due to discrepancies in the parametrization of the water model, which we used without alteration. All combining rules used show the characteristics of a type IIIb fluid: a decrease of Tc upon addition of methane, a minimum in the critical locus, and an increase in Tc with subsequent addition of methane. The differences are seen in the magnitude of these changes. The Lorentz-Berthelot combining rules predict a critical locus that reaches its minimum around T ) 627 K and xW ) 0.76. In contrast, the minimum in the critical locus using the Kong combining rules is shifted to much higher temperature and water mole fraction (T ) 638 K and xW ) 0.83). We attribute this to the fact that the WM parameter for the unlike interaction energy is greater for the Lorentz-Berthelot combining rules as seen in Table 2. This produces a more attractive potential between the water and methane species, resulting in an increased lowering of the critical temperature as methane is added to the system. Since the Kong combining rules result in a smaller interaction energy, the decrease in the critical temperature is not as large. The Sadus rule was developed for nonpolar mixtures which exhibit type III criticality in attempts to account for the nonspherical molecular shape of chain alkanes. The critical locus predicted by the Sadus combining rule is almost equivalent to that using Lorentz-Berthelot. This result is expected since the size parameter using the Sadus rule is only 0.2% larger. The choice of the combining rules on the Pc-Tc projection of the critical locus is seen in Figure 7. The addition of methane in the system results in an order of magnitude increase in the

Tracing Critical Loci of Binary Fluid Mixtures

Figure 8. Temperature-composition critical locus for a water-octopolar/ polarizable methane mixture. Results for the octopolar (triangles) and polarizable methane models (squares), equation of state41 (open circles), and experimental data41 (filled circles). Simulation statistical uncertainties that are smaller than the symbol size are not shown. Lines are drawn to guide the eye.

experimental system pressure for the highest concentration of methane measured. In the high-pressure range, the repulsive interaction between molecules dominates. Once again, the discrepancy between the simulation and experimental results for the critical pressure of pure water is due to the parametrization of the potential (183 bar for the model vs 220.6 bar experimentally).30 As methane is introduced into the system, all of the investigated combining rules underpredict the critical pressure of the mixture. The Lorentz-Berthelot and Sadus combining rules are more commensurate with the pressures measured experimentally, but the Kong combining rules suffer from an underprediction of over 50% for the highest pressure simulated. We then studied the role of methane polarization on the prediction of the critical locus. The first test was to see the effect of using a fixed-charge octopolar methane potential in the water-methane mixture. Figure 8 shows that the nonpolarizable methane model results in a critical locus that underpredicts the experimental results for the entire range of interest in the TcxW plane. The Lorentz-Berthelot combining rules are used for the unlike pair interactions, and this model predicts interactions that are too attractive between the water and methane molecules. This is reasonable due to the added electrostatic interactions that occur between the molecules. Then, as more methane is added to the system, the repulsive interaction between molecules becomes more prevalent, resulting in the increase in the critical temperature. Although the critical temperatures are underestimated, the minimum in the critical locus is located near the experimental value (xW ∼ 0.76). With the addition of polarization to the methane model, we see that the critical locus is even further underpredicted, but the composition of the temperature minimum is in agreement with experiment. As methane is polarized, the induced dipole interacts favorably with the water as well as other methane molecules. Also shown for comparison is the predicted result from an equation of state based on a variation of the Heilig-Franck equation, which includes a square-well attractive interaction, hard-sphere repulsion,41 and two adjustable parameters to optimize agreement with experiment. It is seen in Figure 8 that the equation of state result also underpredicts the critical locus and overpredicts the location of the temperature minimum (xW ∼ 0.82), indicating that the mixture square-well parameter may be too attractive and the potential not repulsive enough. Figure 9 presents the Pc-Tc critical locus for the wateroctopolar/polarizable methane mixtures. In contrast to what is seen in Figure 7 for the methane model with no electrostatic

J. Phys. Chem. B, Vol. 110, No. 34, 2006 17205

Figure 9. Pressure-temperature critical locus for a water-octopolar/ polarizable methane mixture. Symbols are the same as those given in Figure 8. Simulation statistical uncertainties that are smaller than the symbol size are not shown. Lines are drawn to guide the eye.

moment, the octopolar model predicts a critical locus that is similar to the result from the equation of state. These attractive potentials predict critical temperatures that are lower than experimentally seen but retain features that are similar to the experimental results; the pressure range is more commensurate with the experiments than that seen in Figure 7 for the Kong combining rules. Figure 9 also shows that the pressure predictions using the octopolar methane model are similar to those from the equation of state, although the behaviors of the temperature-composition critical loci (Figure 8) differ. 5. Conclusions In this paper, we have shown that estimates for the critical loci of binary mixtures can be obtained using a combination of grand canonical Monte Carlo simulations, histogram reweighting, and mixed-field finite-size scaling techniques to estimate the critical parameters. For a methane-ethane mixture, we are able to quantitatively predict the experimental results across the entire composition range using the conventional LorentzBerthelot combining rules for the interaction parameters between methane and ethane. This is true for both the Tc-xM and PcTc projections of the critical locus. For a water-methane mixture, we are able to get a qualitative description of the type IIIb critical behavior but are not able to get a quantitative prediction of the experimental properties. Using various combining rules for the interaction between water and methane, we were able to determine that the LorentzBerthelot and Sadus descriptions provide a better prediction of the experimental trend. The introduction of an octopolar electrostatic moment to the methane as well as polarization resulted in potentials that were too attractive and underestimated the critical temperature for the entire composition range studied. The choice of combining rules is important for describing the interaction between different molecular species. Studies by Delhommelle and Millie´59 and Al-Matar and Rockstraw60 have shown that the energy parameter is overpredicted using the Lorentz-Berthelot geometric average and that the Kong rules provide better agreement with experiment. For the methaneethane mixture, there is a negligible difference between the parameters determined using either the Lorentz-Berthelot or Kong combining rules. However, the choice of the combining rule is very important for the water-methane system. We have found better agreement using the more attractive geometric averaging of the energy parameters than that of the Kong combining rule. This was also found in a simulation study of other polar-nonpolar mixtures and the more attractive interaction energy was seen as compensating for the lack of polarization

17206 J. Phys. Chem. B, Vol. 110, No. 34, 2006 in the underlying potentials.43 By incorporating the polarization of methane into our mixture, we have seen that the predicted results follow a trend of being too attractive; however the Kong energy parameter may not bring the results into experimental agreement. By adding the explicit polarization of water into the mixture, we may find that the Kong combining rules would bring the simulated results into closer agreement with experiment, but the computational expense in performing these simulations is not yet feasible to test this hypothesis. It is important to note that the unlike interaction potentials were parametrized using only the results from the different combining rules selected. We employed no adjustable parameters in attempts to try and obtain a quantitative fit with the experimental data, unlike the equation of state results. It is possible that the addition of an empirically determined parameter to describe the ij and σij parameters may improve agreement with experimental data, at the cost of rendering the models nonpredictive. Acknowledgment. Funding was provided by the National Science Foundation through a Graduate Research Fellowship to P.J.L. Additional support was provided from the Department of Energy, Office of Basic Energy Sciences (Grant DE-FG0201ER15121). References and Notes (1) Panagiotopoulos, A. Z. Mol. Phys. 1987, 61, 813. (2) Hicks, C. P.; Young, C. L. Chem. ReV. 1975, 75, 119. (3) Sadus, R. J. AIChE J. 1994, 40, 1376. (4) Mu¨ller, E. A.; Gubbins, K. E. Ind. Eng. Chem. Res. 2001, 40, 2193. (5) Eckert, C. A.; Knutson, B. L.; Debenedetti, P. G. Nature 1996, 383, 313. (6) Blanchard, L. A.; Brennecke, J. F. Ind. Eng. Chem. Res. 2001, 40, 287. (7) Savage, P. E.; Gopalan, S.; Mizan, T. I.; Martino, C. J.; Brock, E. E. AIChE J. 1995, 41, 1723. (8) Van Konynenburg, P. H.; Scott, R. L. Philos. Trans. Soc. London, Ser. A 1980, 298, 495. (9) Kiselev, S. B.; Rainwater, J. C. Fluid Phase Equilib. 1997, 141, 129. (10) Van Nhu, N.; Kohler, F. Fluid Phase Equilib. 1995, 105, 153. (11) Polishuk, I.; Wisniak, J.; Segura, H. Fluid Phase Equilib. 1999, 164, 13. (12) Llovell, F.; Vega, L. F. J. Phys. Chem. B 2006, 110, 1350. (13) McCabe, C.; Galindo, A.; Cummings, P. T. J. Phys. Chem. B 2003, 107, 12307. (14) Polishuk, I.; Wisniak, J.; Segura, H.; Yelash, L. V.; Kraska, T. Fluid Phase Equilib. 2000, 172, 1. (15) Matubayasi, N.; Nakahara, M. J. Chem. Phys. 2000, 112, 8089. (16) Neichel, M.; Franck, E. U. J. Supercrit. Fluids 1996, 9, 69. (17) Galindo, A.; Whitehead, P. J.; Jackson, G.; Burgess, A. N. J. Phys. Chem. 1996, 100, 6781. (18) Duan, Z.; Møller, N.; Weare, J. H. Geochim. Cosmochim. Acta 1992, 56, 2619. (19) Anderko, A.; Pitzer, K. S. AIChE J. 1991, 37, 1379. (20) Singer, J. V. L.; Singer, K. Mol. Phys. 1972, 24, 357. (21) Nakanishi, K.; Toukubo, K. J. Chem. Phys. 1979, 70, 5848. (22) Harismiadis, V. I.; Koutras, N. K.; Tassios, D. P.; Panagiotopoulos, A. Z. Fluid Phase Equilib. 1991, 65, 1. (23) Johnson, J. K.; Gubbins, K. E. Mol. Phys. 1992, 77, 1033.

Lenart and Panagiotopoulos (24) Potoff, J. J.; Panagiotopoulos, A. Z. J. Chem. Phys. 1998, 109, 10914. (25) Mej×c1a, A.; Pa`mies, J. C.; Duque, D.; Segura, H.; Vega, L. F. J. Chem. Phys. 2005, 123, 034505. (26) Errington, J. R.; Panagiotopoulos, A. Z. J. Chem. Phys. 1998, 109, 1093. (27) Chen, B.; Siepmann, J. I. J. Phys. Chem. B 1999, 103, 5370. (28) Errington, J. R.; Panagiotopoulos, A. Z. J. Phys. Chem. B 1999, 103, 6314. (29) Ungerer, P.; Beauvais, C.; Delhommelle, J.; Boutin, A.; Rousseau, B.; Fuchs, A. H. J. Chem. Phys. 2000, 112, 5499. (30) Errington, J. R.; Panagiotopoulos, A. Z. J. Phys. Chem. B 1998, 102, 7470. (31) Paricaud, P.; Pøedota, M.; Chialvo, A. A.; Cummings, P. T. J. Chem. Phys. 2005, 122, 244511. (32) McGrath, M. J.; Siepmann, J. I.; Kuo, I. F. W.; Mundy, C. J.; VandeVondele, J.; Hutter, J.; Mohamed, F.; Krack, M. J. Phys. Chem. A 2006, 110, 640. (33) Lin, C. L.; Wood, R. H. J. Phys. Chem. 1996, 100, 16399. (34) Errington, J. R.; Boulougouris, G. C.; Economou, I. G.; Panagiotopoulos, A. Z.; Theodorou, D. N. J. Phys. Chem. B 1998, 102, 8865. (35) Jedlovszky, P.; Vallauri, R.; Richardi, J. J. Phys.: Condens. Matter 2000, 12, A115. (36) Konrad, O.; Lankau, T. J. Phys. Chem. B 2005, 109, 23596. (37) Virnau, P.; Mu¨ller, M.; MacDowell, L. G.; Binder, K. J. Chem. Phys. 2004, 121, 2169. (38) Binder, K.; Mu¨ller, M.; Virnau, P.; MacDowell, L. G. AdV. Polym. Sci. 2005, 173, 1. (39) Hynninen, A. P.; Dijkstra, M.; Panagiotopoulos, A. Z. J. Chem. Phys. 2005, 123, 084903. (40) Bloomer, O. T.; Gami, D. C.; Parent, J. D. Inst. Gas Technol. Res. Bull. 1953, 22, 1. (41) Shmonov, V. M.; Sadus, R. J.; Franck, E. U. J. Phys. Chem. 1993, 97, 9054. (42) Buckingham, R. A. Proc. R. Soc. London, Ser. A 1938, 168, 264. (43) Potoff, J. J.; Errington, J. R.; Panagiotopoulos, A. Z. Mol. Phys. 1999, 97, 1073. (44) Sadus, R. J. J. Phys. Chem. 1993, 97, 1985. (45) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269. (46) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. J. Am. Chem. Soc. 1996, 118, 11225. (47) Rick, S. W.; Stuart, S. J. ReV. Comput. Chem. 2002, 18, 89. (48) Smirnov, K. S.; van de Graaf, B. J. Chem. Soc., Faraday Trans. 1996, 92, 2475. (49) Theodorou, D. N.; Suter, U. W. J. Chem. Phys. 1985, 82, 955. (50) Frenkel, D.; Smit, B. Understanding Molecular Simulation: From Algorithms to Applications; Academic Press: San Diego, 2002. (51) Nymand, T. M.; Linse, P. J. Chem. Phys. 2000, 112, 6152. (52) Ferrenberg, A. M.; Swendsen, R. H. Phys. ReV. Lett. 1988, 61, 2635. (53) Ferrenberg, A. M.; Swendsen, R. H. Phys. ReV. Lett. 1989, 63, 1195. (54) Wilding, N. B.; Bruce, A. D. J. Phys.: Condens. Matter 1992, 4, 3087. (55) Chen, B.; Siepmann, J. I. Theor. Chem. Acc. 1999, 103, 87. (56) Chen, B.; Potoff, J. J.; Siepmann, J. I. J. Phys. Chem. B 2000, 104, 2378. (57) Chen, B.; Xing, J.; Siepmann, J. I. J. Phys. Chem. B 2000, 104, 2391. (58) Lenart, P. J.; Panagiotopoulos, A. Z. Ind. Eng. Chem. Res. 2006, DOI 10.1021/ie051302i. (59) Delhommelle, J.; Millie´, P. Mol. Phys. 2001, 99, 619. (60) Al-Matar, A. K.; Rockstraw, D. A. J. Comput. Chem. 2004, 25, 660. (61) NIST Chemistry WebBook. http://webbook.nist.gov/chemistry (accessed 05/04/06), NIST Standard Reference Database Number 69, June 2005 Release.