Molecular Simulation and Regular Solution Theory Modeling of Pure

Nov 19, 2008 - Quintin R. SheridanWilliam F. SchneiderEdward J. Maginn .... Joshua D. Moon , Michelle S. Hindman , C. Heath Turner , and Jason E. Bara...
0 downloads 0 Views 185KB Size
Download PDF
16710

J. Phys. Chem. B 2008, 112, 16710–16720

Molecular Simulation and Regular Solution Theory Modeling of Pure and Mixed Gas Absorption in the Ionic Liquid 1-n-Hexyl-3-methylimidazolium Bis(Trifluoromethylsulfonyl)amide ([hmim][Tf2N]) Wei Shi and Edward J. Maginn* Department of Chemical and Biomolecular Engineering, UniVersity of Notre Dame, 182 Fitzpatrick Hall, Notre Dame, Indiana 46556-5637 ReceiVed: August 25, 2008; ReVised Manuscript ReceiVed: October 1, 2008

Monte Carlo simulations are carried out to compute the solubility of SO2, O2, and N2 in the ionic liquid 1-n-hexyl-3-methylimidazolium bis(trifluoromethylsulfonyl)amide ([hmim][Tf2N]). Simulations are also used to compute the mixed gas isotherms for the mixtures CO2/O2, SO2/N2, and CO2/SO2. The pure gas isotherms agree well with available experimental data. For the mixtures CO2/O2 and SO2/N2, the simulations predict that the mixed gas solubilities are nearly ideal, with little enhancement or competition between the two solutes. This is contrary to published experimental studies, which found that in a gas mixture, CO2 enhances the solubility of the otherwise sparingly soluble O2 and that the poorly soluble N2 inhibits the solubility of the highly soluble SO2. For the SO2/CO2 mixture, it is found that the two gases absorb independently at low pressure but compete with one another at high pressure. An energetic analysis was performed for the different solutes in the ionic liquid. The van der Waals energy between the solutes and ionic liquid was greater than the electrostatic interactions, and the van der Waals interactions were roughly the same between the solutes and the two different ions. CO2 and SO2 interact more strongly with the anion than the cation due to stronger electrostatic interactions between the solute and the anion. N2 and O2 interact weakly with the ionic liquid and show little difference in interaction energy between the cation and anion. Regular solution theory (RST) was evaluated for its ability to predict pure and mixed gas isotherms. It was found that if RST was applied in a strictly predictive mode using experimentally derived parameters, solubilities of CO2 were underpredicted by a substantial margin. 1. Introduction and Background Ionic liquids (ILs) have several properties that make them attractive for use as a separations agent, including extremely low volatility, good thermal stability, and large solubility differences among different species.1 They have been proposed for use in many applications, including extractive distillation,2 the separation of light gases,3,4 and in supported liquid membranes.5 The key to the success of all of these applications is the phase equilibrium between the species to be separated and the IL. A tremendous amount of experimental data has become available in recent years on gas and vapor solubility in ILs.6-14 Among gases, it has been shown that CO215,16 and SO217,18 are highly soluble in ILs, while other gases such as N2, O2, H2, CO, and CH4 are only sparingly soluble.10,11,19 The large differences in pure gas solubility reported in the literature obviously suggest that selectivity for important mixtures containing these gases could be high, making ILs a strong candidate for new and effective separation processes. It is important to realize, however, that the phase behavior observed for a mixture of gases cannot in general be predicted from pure gas solubility. This is unfortunate because, while single gas solubility measurements using various experimental techniques have now become almost routine, measurements of mixed gas phase equilibria are significantly more difficult, and consequently, mixture data are rare. This problem has long been appreciated in the gas adsorption community, where single gas isotherms on zeolites, carbons, and other porous materials are * To whom correspondence should be addressed. Phone: (574) 631-5687. Fax: (574) 631-8366. E-mail: [email protected].

plentiful because they can be measured easily but where there is a dearth of mixture data because the experiments are “at least an order of magnitude more complicated” than their pure gas counterparts.20 In fact, we are aware of only two studies in which mixed gas solubilities in ILs have been investigated and only one where the composition of both solutes in the two phases was measured. We note that there have been several other studies looking at solubilities of mixed hydrocarbon vapors21,22 in ILs. Hert and co-workers23 measured CO2/O2 and CO2/CH4 solubilities in the IL 1-n-hexyl-3-methylimidazolium bis(trifluoromethylsulfonyl)amide ([hmim][Tf2N]). They found that for the mixtures, the CO2 solubility decreased while the O2 and CH4 solubilities increased relative to what would be expected from pure gas-phase equilibria. It was speculated that the enhancement in solubility of the otherwise sparingly soluble O2 and CH4 was due to favorable dispersion interactions between dissolved CO2 and either O2 or CH4, while the drop in CO2 solubility was due to a reduction in free volume caused by the absorption of O2 and CH4. As discussed below, however, this group has recently reinvestigated the CO2/O2 mixture and found that their initial conclusions may be incorrect due to experimental error and/or inaccuracy.24 Huang and co-workers measured SO2 solubility in the presence of N2 in five different ILs.17 They reported that for the pure gases, SO2 is highly soluble while N2 absorption could not be detected, a finding consistent with other studies.10,18 However, when a mixture containing 10 mol % SO2 and 90 mol % N2 was exposed to the ILs, a surprising result occurred. It was observed that the amount of SO2 absorbed was drastically

10.1021/jp8075782 CCC: $40.75  2008 American Chemical Society Published on Web 11/19/2008

Pure and Mixed Gas Absorption in [hmim][Tf2N] reduced for two of the ILs, which the authors attributed to competitive absorption of N2. Their results suggest that N2 absorbs in appreciable amounts in the mixture to compete with SO2, although as a pure gas N2 is virtually insoluble in the ionic liquid. No mechanism or explanation was proposed for this unusual behavior, but this group has begun reinvestigating these unusual results.25 Given the difficulties associated with performing mixture measurements, it is clear that methods are needed that can help predict the solubility of mixed gases in ILs. This would be true even if the mixture experiments were not so difficult to conduct since there is essentially an infinite number of potential ionic liquids and all cannot be tested experimentally. Predictive methods that give insight into what gives rise to differences in solubility are thus needed. Molecular modeling and regular solution theory (RST) are two such methods that have been used to help predict pure gas solubilities in ILs. In the present work, we compare the performance of these two methods for predicting the solubility of CO2, SO2, O2, and N2 in an IL, as well as the solubilities of CO2/O2, SO2/N2, and CO2/SO2 mixtures. Molecular modeling has been used to help understand solubility trends in ILs26,27 and to compute pure gas Henry’s law constants/chemical potentials.28-30 Maurer and co-workers computed pure gas isotherms for O2, CO2, CO, and H2 in an IL using Gibbs ensemble Monte Carlo.9,31,32 Recently, we have computed isotherms for CO2 and water in [hmim][TF2N]33 using an advanced “continuous fractional component” Monte Carlo (CFC MC) simulation procedure in the osmotic ensemble. The agreement between computed and experimental isotherms was quantitative in the latter study. Whereas mixed gas solubility experiments are significantly more difficult to execute than pure gas measurements, mixture simulations are not much more difficult to conduct than their pure gas counterparts. It is reasonable to expect that if the molecular modeling captures the physics of the gas-ionic liquids system well enough to obtain accurate pure gas isotherms, then modeling should be capable of yielding good estimates of mixed gas absorption as well. This suggests a useful strategy; simulations of pure gas solubility can be validated against (relatively) easily to obtain experimental pure gas isotherms. The validated models can then be used to predict mixed gas solubility, where experiments are much more difficulty to carry out. Such a strategy has been used successfully with adsorption in zeolites.34-36 An additional advantage is that the simulations enable one to probe local structure and energetics, thereby giving physical insight into the source of selectivity trends observed. In this work, we report the first mixed gas solubility simulations for gases in ILs. 2. Regular Solution Theory Molecular modeling is computationally intensive, especially for complex systems such as ILs, and therefore, there is an obvious advantage to having simple analytic theories that allow rapid estimates of solubility. Standard cubic equations of state have been shown to correlate gas solubility data in ILs exceptionally well.12 Likewise, quantitative structure-property relationship (QSPR) models have been developed for solubilities in ILs,37 as have activity coefficient models.21 These methods require experimental solubility data against which parameters are fit, however, and thus cannot be said to be purely “predictive”. Noble and co-workers, on the other hand, have shown that the relatively simple regular solution theory (RST) model can give qualitative and even semiquantitative predictions of gas solubility trends in ionic liquids.14,38,39 For example, trends in Henry’s law constants as a function of temperature and ionic

J. Phys. Chem. B, Vol. 112, No. 51, 2008 16711 liquid molar volume conform well to the predictions that correspond to RST. This was a remarkable finding, given that the underlying basis for RST is that the gas-ionic liquid solutions must be “regular”, that is, that the entropy change of transferring the solute from an ideal solution to the solvated system is zero.40 This is something not normally expected to be satisfied in charged molten salts. Up until now, RST has only been applied to estimate pure gas Henry’s law constants, but the theory itself can be applied to predict entire isotherms of pure and mixed gas systems. The theory can also be tested for consistency against molecular simulations since, for a given molecular model, one can compute all of the parameters that go into the RST model. By comparing the RST predictions for this model against simulations (which give an exact solution for the model), the validity of RST can be tested independent of any experimental system. Since the models used to compute properties of ILs have been repeatedly shown to yield good estimates of a wide range of IL thermodynamic properties,41 consistency among RST and simulation results implies that RST will be accurate when compared to experiment. Let subscript “1” denote a solute (gas) species, subscript “2” denote the ionic liquid solvent, and superscripts “L” and “V” denote the liquid and vapor phases, respectively. Phase equilibrium dictates that the fugacity of the solute will be equal in the liquid and vapor

f L1 ) f 01(T, P)x1γ1 ) f V1 ) Py1φV1 (T, P)

(1)

where x1 and y1 are the mole fractions of the solute in the liquid and the gas phase, respectively, f is the fugacity, and f 10(T, P) is the fugacity of the pure liquid solute at temperature T and pressure P. Given that the IL is essentially nonvolatile, we make the good assumption that y1 ) 1.0. The pure liquid fugacity is computed as f 10(T, P) ) Psatφ1V(T, Psat)exp[(V1/RT)(P - Psat)], where V1 is taken as the liquid molar volume of the solute.42 The pure solute liquid saturation pressure Psat is known experimentally for vapors. For gases that are above their critical temperature, a hypothetical pure liquid is used as the reference state, and the hypothetical saturation pressure can be estimated by extrapolating a plot of ln Psat versus 1/T. For low-solubility gases, a Henry’s Law reference fugacity is usually used, but since this requires as input an experimental Henry’s Law constant (which we seek to predict), the Raoult’s Law reference is used instead. The fugacity coefficients of the pure solute vapor at the saturation pressure and actual system pressure, φV1 (T, Psat) and φ1V(T, P), can be computed from an equation of state. This leaves the activity coefficient, γ1, as the only unknown linking the composition of the liquid phase with the temperature and pressure. RST predicts that the activity coefficient for the solute is given by

RT ln γ1 ) V1Φ22[δ1 - δ2]2

(2)

where R is the universal gas constant and δ1 and δ2 are the solubility parameters, defined as

( )

∆Uivap δi ) Vi

1/2

(3)

In eq 3, ∆Uvap is the internal energy change of vaporization.40 The volume fraction Φ2 is

16712 J. Phys. Chem. B, Vol. 112, No. 51, 2008

Φ2 )

x2V2 x1V1 + x2V2

Shi and Maginn

x1V1

)1-

(4)

(V1 - V2)x1 + V2

where V2 is the IL molar volume. On the basis of the above equations, the following equation is obtained

ln

x1 PφV1

) -ln

f 01(T, P)

-

V1[δ1 - δ2]2 RT

[

× x1V1

1-

(V1 - V2)x1 + V2

]

2

(5)

The solute mole fraction, x1, may be obtained from P, T, f 01(T, P), φ1V, V1, V2, δ1, and δ2. Hypothetical solute molar volumes V1 and solubility parameters δ1 are tabulated for many substances.42 Ionic liquid molar volumes V2 are widely available. Thus, the only unknown in eq 5 is the solubility parameter of the IL, δ2. This is a pure component property of the IL that can be determined in several ways. The solubility parameter may be determined directly from experiment if the liquid density and enthalpy of vaporization are known. Density measurements are plentiful for ILs, and three independent enthalpy of vaporization measurements have recently been reported for this IL.44-46 Two of these measurements44,45 give a nearly identical solubility parameter for this IL of 20.78 (J/cm3)0.5 at 298 K, while the third experiment46 gives a value of 22.97 (J/cm3)0.5. The solubility parameter may also be estimated using an empirical relationship involving the melting temperature38 or by treating it as an adjustable parameter that is a function of the molar volume of the IL and fitting to experimental data.39 The solubility parameter can be computed directly from a molecular simulation47 as well. Previous calculations48 for this IL have found that the solubility parameter at 293 K is 21.8 (J/cm3)0.5. As described below, the present simulations find that the value at 333 K is 20.45 (J/cm3)0.5; both are in good agreement with the available experimental data. 3. Simulation Details 3.1. Force Field. A standard force field of the following form was used for the total energy V (r)

V(r) )

∑

∑ kb(r - r0)2 + ∑ kθ(θ - θ0)2 +

bonds

angles

kχ[1 + cos(n0χ - δ0)] +

dihedrals

∑

kψ(ψ - ψ0)2 +

impropers

N-1

N

∑ ∑

i)1 j)i+1

{ [( ) ( ) ] 4ij

σij rij

12

-

σij rij

6

+

qiqj rij

}

(6)

where the symbols have their conventional meaning.49 The Lennard-Jones parameters for the fractional molecules used in the CFC MC algorithm were scaled using the functional form proposed by van Gunstern et al.50 to aid insertions. The force field parameters for SO2 were taken from previous work.51 With the exception of the harmonic bond force constants, which were computed at the B3LYP/6-311+G* level of theory, the parameters for oxygen and nitrogen were taken from the literature.52,53 Nominal force constants kb were found to be 600 kcal/(mol Å2) for oxygen and 1600 kcal/(mol Å2) for nitrogen. The nominal bond lengths r0 were taken as 1.21 Å for oxygen

and 1.0897 Å for nitrogen. The Lennard-Jones parameters for each oxygen atom were taken as σ ) 3.03 Å and /kB ) 48.3 K. For the nitrogen atoms, σ ) 3.31 Å and /kB ) 36.0 K. No partial charges were used for oxygen or nitrogen. The force field parameters for the IL were the same as those used previously by us.33 For the IL, the intramolecular electrostatic and LennardJones interactions for atoms separated by three consecutive bonds were scaled by 0.5 and were neglected for atoms closer than that. Standard Lorentz-Berthelot combining rules were used for Lennard-Jones parameters. 3.2. Single Gas Isotherms. Pure gas isotherms for SO2, O2, and N2 were computed using the CFC MC procedure in the osmotic ensemble. Details of the method may be found elsewhere;33,54 therefore, only those items specific to this set of simulations are discussed. CO2 isotherms in [hmim][Tf2N] have been computed previously33 and are used for comparison here. Simulations were run in the osmotic ensemble with the temperature and pressure (i.e., gas partial pressure) fixed. The Lennard-Jones potential was truncated using a switching function33,54 with ron ) 10.5 Å and roff ) 12.0 Å. A Verlet neighbor list with a 13.5 Å radius was used. The Fennell and Gezelter shifted force (FGSF) method was used to compute electrostatic interactions.55 The cutoff distance was taken as 12.0 Å, and the damping parameter was chosen to be κ ) 0.2022 Å-1. As shown previously,33 the FGSF method yields energies and densities in good agreement with conventional Ewald sum techniques but with less computational cost. We confirmed the accuracy of the method by performing test simulations with SO2/ N2 and CO2/O2 mixtures at room temperature and very low to intermediate densities. The differences in energy and density computed using Ewald and FGSF were small, with the best agreement coming at higher density. All simulations were performed on a system containing 160 ion pairs, using a locally developed code. Equilibration runs of 1-2 million moves were carried out, during which time the boldness of the various MC moves was tuned to achieve roughly 50% acceptance rates. The following moves were used with fixed probabilities given in parentheses: thermal equilibration via hybrid Monte Carlo56 (60%); fractional solute coupling strength (38%); system volume change (2%). For the hybrid Monte Carlo, five molecular dynamics steps with a time step of between 0.5 and 0.7 fs were used per MC move. Volume changes were uniform up to a maximum of 400-600 Å3. The CFC coupling parameter λ was changed uniformly up to a maximum value of 0.3-0.5 (depending on the state point and system) to give an acceptance rate of about 50%. The CFC bias factors54 were adjusted during this stage to achieve as closely as possible a uniform probability distribution for λ. Optimization of these bias factors was done using the Wang-Landau updating scheme.57 After equilibration, production runs of several million moves were performed, with the exact number dictated by the convergence characteristics of the system. The solute fugacity was set as f ) Pφ, where φ is the pure solute fugacity coefficient at the conditions of interest and P is the total pressure. It was assumed that no ionic liquid was present in the vapor phase, such that the total pressure was equal to the solute partial pressure. The fugacity coefficient was computed using the Peng-Robinson equation of state (PReos), and it was assumed that this value was consistent with the force field used. This is a reasonable assumption for the simple gases studied here (i.e., SO2, CO2, O2, and N2). Block averages were used to compute properties, and the standard deviation was used to estimate uncertainty.

Pure and Mixed Gas Absorption in [hmim][Tf2N] 3.3. Mixed Gas Simulations. To model mixed gas absorption, constant pressure Gibbs ensemble simulations were conducted using a locally developed code. This code was validated by performing vapor-liquid equilibria calculations on a Lennard-Jones mixture of species “1” and “2” with σ1/σ2 ) 1, 1/2 ) 2.0 at T* ) kBT/1 ) 1.0, and P* ) Pσ13/1 ) 0.1064, as well as at T* ) 1.0 and P* ) 0.05927. Results were compared against those obtained using a traditional Gibbs ensemble Monte Carlo code supplied by Prof. Karl Johnson58 as well as against results obtained by Potoff et al.59 and found to be in excellent agreement. The total number of each type of solute molecule, the number of solvent molecules, the temperature, and the pressure were specified at the start of the simulation. In all cases, 160 ion pairs were used. For mixtures of CO2/O2, simulations were performed at 313 K and at pressures ranging from 30 to 120 bar. The total number of CO2 and O2 molecules in the two phases ranged from 200 to 260 and 140 to 200, respectively. For SO2/N2, simulations were performed at 333 K and pressures ranging from 3 to 5 bar. Note that the model saturation vapor pressure for pure SO2 was computed to be 8.6 ( 0.4 bar at 333 K,51 larger than the highest simulated pressure of 5 bar. A total of anywhere from 16 to 71 SO2 molecules was used, while 90-350 N2 molecules were simulated. For CO2/SO2, simulations were performed at 333 K and 5 and 10 bar; 330 CO2 molecules and 80 SO2 molecules were simulated. The partial pressure of SO2 in the vapor phase was much lower than the model saturation pressure. Mixture simulations typically consisted of 2 million equilibration moves and 10 million production moves. The ratio of different moves was the same as that for the pure gas isotherms. For the λ move, the probability of choosing a particular component was proportional to the concentration of that component in the two boxes. Note that unlike the single component constant volume simulations, the pressure must be specified so that the volume of each simulation box can be adjusted during the simulation. Standard acceptance rules are used for this move.60,61 4. Results 4.1. Pure Gas Solubility. 4.1.1. CO2 Solubility. In previous work33 we showed that CO2 isotherms and infinite dilution partial molar enthalpies in [hmim][Tf2N] were reproduced quantitatively with molecular simulation. To test the ability of RST to predict CO2 isotherms, the following procedure was used. Literature values42 for the CO2 molar volume, V1 ) 55 cm3/mol, and solubility parameter, δ1 ) 12.27 (J/cm3)0.5, were used without modification. The molar volume of the IL, V2, has been measured62 and found to be 338 cm3/mol at 333 K and 0.984 bar, 330 cm3/mol at 313 K and 0.987 bar, and 327 cm3/ mol at 298 K and 0.985 bar, respectively. The only remaining unknown parameter needed to predict the isotherm is the solubility parameter for the IL, δ2. As mentioned above, the solubility parameter for this IL has been measured44-46 and found to range from 20.78 to 22.97 (J/cm3)0.5. Figure 1 shows a comparison of the experimental isotherm for CO2 in [hmim][Tf2N] at 333 K and the predictions given by RST using the parameters listed above but with varying values of the IL solubility parameter. As can be observed, the RST solubilities are significantly lower than what is observed experimentally when a solubility parameter is used that is commensurate with the experimental values; at 100 bar, RST predicts that the CO2 mole fraction is 0.3, when in fact the experimental value is closer to 0.7. The RST results can be brought into closer agreement

J. Phys. Chem. B, Vol. 112, No. 51, 2008 16713

Figure 1. Comparison of experimental, simulated, and RST isotherms for CO2 absorption in [hmim][Tf2N] at 333 K. The CO2 mole fraction is given by x. The diamonds, triangles, stars, and pluses correspond to RST δ2 values of 20.50, 18.41, 17.39, and 16.36 (J/cm3)0.5, respectively. In these cases, literature values42 of V1 ) 55 cm3/mol and δ1 ) 12.27 (J/cm3)0.5 were used for CO2, and the experimental molar volume62 of 3 33 V1 ) 338 cm /mol was used for the IL. Simulation results are shown as filled circles, while the experimental data are shown as filled squares. The open circles are the RST prediction using parameters calculated directly from simulations, V1 ) 64.5 cm3/mol, δ1 ) 10.64 (J/cm3)0.5, 3 3 0.5 V2 ) 336.4 cm /mol, and δ2 ) 20.45 (J/cm ) .

with the experimental results by adjusting the value of the solubility parameter δ2. As shown in Figure 1, by reducing the solubility parameter to δ2 ) 16.36 (J/cm3)0.5, relatively good agreement is obtained at low pressure, while δ2 ) 17.39 (J/ cm3)0.5 captures the higher pressure data better. The absolute average difference in mole fractions between RST and experiment is 9.4% for δ2 ) 17.39 (J/cm3)0.5 and 13.6% for δ2 ) 16.36 (J/cm3)0.5; therefore, the former value was chosen as the “fitted” solubility parameter to be used in all subsequent RST calculations for this IL. It is interesting to note that Scovazzo and co-workers38 fit CO2 solubility data to RST and found that δ2 ) 16.4 (J/cm3)0.5 for [emim][Tf2N], which is consistent with the present results but much lower than the experimentally observed solubility parameter.44-46 Alternatively, all of the RST parameters (V1, V2, δ1, and δ2) can be computed from atomistic simulations and then used within RST to predict an isotherm. By comparing this prediction against that obtained from the molecular simulations, the internal consistency of RST can be tested. That is, the degree to which RST agrees or disagrees with the simulations in which both “use” the same set of parameters is a measure of how appropriate or inappropriate the theory is for predicting rather than correlating gas solubility in ILs. This is because the simulations provide an exact numerical solution of the phase equilibria for the model given by eq 6. Given the success with which simulations have matched experimental gas solubility in ILs,33 it is expected that the force field does a good job capturing the essential physics of gas-ionic liquid interactions. The molar volume and solubility parameter for [hmim][Tf2N] were computed by performing isothermal-isobaric Monte Carlo simulations of the IL in the gas and liquid phases. For the gasphase calculations, one pair of [hmim] and [Tf2N] ions were modeled at 333 K and a pressure of 1 × 10-5 Pa. For the liquid phase, isothermal-isobaric simulations were carried out on 160 ion pairs at 333 K and 0.984 bar. V2 and δ2 were calculated as 336.4 cm3/mol and 20.45 (J/cm3)0.5, respectively. The molar volume was nearly identical to the experimental value62 of 338 cm3/mol, while the solubility parameter was significantly higher than that obtained by fitting RST to the experimental data but agreed very well with experimentally obtained values. For CO2, the molar volume and solubility parameter were computed by conducting gas-phase simulations at 333 K and 4

16714 J. Phys. Chem. B, Vol. 112, No. 51, 2008

Shi and Maginn TABLE 1: The Simulated Mole Fraction x and Molar Mixture Volume Vmix at Different Temperatures T, Pressures P, and Fugacities f for Pure Gas Absorption in the Ionic Liquid [hmim][Tf2N]a

Figure 2. Comparison of CO2 isotherms in [hmim][Tf2N] at 313 and 298 K predicted from RST and the corresponding experimental data. Literature42 RST parameters were used for CO2 (V1 ) 55 cm3/mol and δ1 ) 12.27 (J/cm3)0.5). The previously fit solubility parameter for the IL, δ2 ) 17.39 (J/cm3)0.5, was used along with the experimental62 IL molar volume at 333 K (V2 ) 338 cm3). The temperature-dependent experimental molar volumes (V2 ) 327 cm3 at 298 K and V2 ) 330 cm3 at 313 K) were also used. The temperature dependence of the molar volume has little effect on the predicted isotherm. The lines are guides to the eye.

bar and liquid-phase simulations at 333 K and 200 bar. This latter state point was intended to mimic the hypothetical “liquid” CO2. The values of V1 and δ1 were found to be 64.5 cm3/mol and 10.64 (J/cm3)0.5 at 333 K, respectively, which are similar to those suggested in the literature for 298 K42 (i.e., 55 cm3/ mol and 12.27 (J/cm3)0.5). As shown in Figure 1 by the open circles, the RST prediction using these calculated parameters significantly underpredicts the solubility and does not agree even qualitatively with the simulation results. On the other hand, the simulated isotherm is in near perfect agreement with the experimental data. This suggests that while RST is able to correlate experimental data reasonably well, there are limitations as to its pure predictive capabilities; some of the physics responsible for the high solubility of CO2 in ILs (despite the large difference between ionic liquid and CO2 solubility parameters) is not captured by the RST model but appears to be accounted for in the simulations. This can only be compensated for by treating the IL solubility parameter as an adjustable parameter within RST. The ability of RST to predict CO2 isotherms at temperatures other than 333 K was tested by using the fitted solubility parameter δ2 ) 17.39 (J/cm3)0.5 and the experimental molar volume of the ILs. Figure 2 shows a comparison between RST predictions and experimental results.62 At 313 K, the absolute average difference in x1 between the RST model and experiment is 19.1%, while at 298 K, it is 18.6%. The shape of the isotherms at 298 and 313 K is not captured particularly well, particularly at high pressure. On the basis of these results, it appears that the RST model is able to match CO2 solubility to within about 10-20% when a single adjustable parameter (the solubility parameter) is fit to an experimental isotherm. The shapes of the isotherms are not captured very well, however, especially at high pressures. We therefore conclude that RST is most useful for correlating low-pressure solubility data in the Henry’s Law regime and for making limited predictions of Henry’s Law constants. It is less useful for making pure predictions of isotherms. This is to be expected, given the simplifying assumptions upon which it is based. RST is still a useful tool, considering the ease with which it can be applied. 4.1.2. O2 Solubility. Experimental data for O2 solubility in [hmim][Tf2N] at 298 K was published by Hert and co-workers23 in 2005. The Henry’s Law constant, determined by fitting a straight line forced through the origin to the data, was found to

solute

T (K)

P (bar)

f (bar)

x

Vmix (cm3/mol)

O2 O2 O2 O2 O2 O2 O2 O2 SO2 SO2 SO2 SO2 SO2 SO2 SO2 SO2 SO2 N2 N2 N2 N2 N2 N2 N2 N2 N2 N2 N2 N2 N2 N2

373 373 373 373 313 313 313 313 333 333 333 333 333 298 298 298 298 373 373 373 373 333 333 333 333 298 298 298 298 298 298

25 50 100 200 12.5 25 50 100 0.25 0.5 1.0 2.0 3.0 0.0375 0.075 0.15 0.4 100 200 300 500 50 100 200 500 20 40 80 120 200 500

24.8275 49.36 97.82 194.04 12.385 24.5525 48.3 93.88 0.24935 0.49745 0.9899 1.9592 2.9082 0.03748 0.074925 0.149685 0.39776 101.3 207.8 323.1 591 49.76 99.6 202.4 572 19.8360 39.4 78.0240 116.4720 194.76 544.5

0.041 (1) 0.062 (3) 0.121 (4) 0.194 (6) 0.022 (2) 0.054 (3) 0.074 (4) 0.147 (3) 0.085 (2) 0.136 (5) 0.211 (2) 0.382 (1) 0.531 (1) 0.057 (1) 0.070 (1) 0.127 (2) 0.319 (1) 0.063 (4) 0.115 (4) 0.182 (4) 0.209 (4) 0.039 (1) 0.079 (4) 0.138 (3) 0.198 (3) 0.034 (2) 0.050 (2) 0.070 (2) 0.086 (4) 0.130 (6) 0.233 (3)

330.9 (4) 322.7 (11) 304.6 (10) 279.6 (18) 321.7 (6) 314.1 (9) 305.5 (12) 284.5 (9) 306.8 (6) 292.7 (11) 269.6 (8) 222.5 (4) 180.7 (3) 309.7 (4) 304.7 (5) 288.2 (5) 234.4 (4) 320.6 (11) 304.1 (11) 284.3 (15) 271.9 (14) 320.5 (3) 307.8 (10) 289.2 (7) 269.4 (10) 318.1 (6) 312.0 (5) 304.6 (7) 299.5 (12) 285.2 (18) 253.5 (7)

a Uncertainties are computed as standard deviations of block averages; the numbers in parentheses are the uncertainties in the last digit.

be 966 ( 21 bar, indicating very low oxygen solubility. Recently, this group has re-examined these data and carried out additional measurements.43 They now find that O2 is actually much more soluble in this IL than previously reported; their revised Henry’s Law constants are 515 ( 10 bar at 298 K and 561 ( 18 at 313 K.43 The computed solubility of O2 in [hmim][Tf2N] at 313 and 373 K and various pressures is given in Table 1 and compared against the original “old”23 and revised “new”43experimental data in Figure 3. The simulated Henry’s Law constants are 623 ( 20 and 790 ( 28 bar at 313 and 373 K, respectively (see Table 2). Note that there are no experimental data available at 313 K with a fugacity larger than 10 bar. Hence, the experimental solubility at a fugacity larger than 10 bar was estimated by linear extrapolation of the low-pressure data, designated by the line in Figure 3. The simulations at 313 K show an absolute average difference with the experiments of 12.4%, consistent with the ∼10% difference between experimental and computed Henry’s Law constants at 313 K. Both the experiments and simulations show a very small or slightly negative enthalpy of dissolution for oxygen. RST predictions of O2 solubility were made by using the same experimental pure component molar volume and fitted solubility parameter for [hmim][Tf2N] as that used for CO2 solubility calculations (i.e., V2 ) 338 cm3/mol, δ2 ) 17.39 (J/cm3)0.5), while literature values42 were used for O2 (V1 ) 33.0 cm3/mol and δ1 ) 8.18 (J/cm3)0.5). The RST calculations at 313 K were performed at conditions shown in Table 1. Figure 3 shows that

Pure and Mixed Gas Absorption in [hmim][Tf2N]

Figure 3. A log-log comparison of experimental, simulated, and RST isotherms of O2 in [hmim][Tf2N]. The “new” experimental data43 is a revised estimate of the original “old” data.23 The solid line is a linear fit of the low-fugacity “new” experimental data at 313 K.

J. Phys. Chem. B, Vol. 112, No. 51, 2008 16715

Figure 5. Comparison of experimental18 and simulated isotherms for SO2 in [hmim][Tf2N].

TABLE 2: The Simulated and Experimental Henry’s Law Constants for O2, N2, and SO2 Absorption in [hmim][Tf2N] at Different Temperaturesa Henry’s law constant (bar) solute

temperature (K)

O2 O2 SO2 SO2 N2 N2 N2

313 373 333 298 373 333 298

simulation 623 (20) 790 (28) 3.6 (0.5) 0.9 (0.3) 1772 (51) 1410 (42) 1138 (47)

experiment 561 (18)43 4.09 (0.06)18 1.64 (0.01)18 848 (27)63

a

Estimates were made by fitting the low-pressure results to a line forced to pass through the origin. The standard error from the fitting is shown in the parentheses. Also shown are experimental values for the Henry’s law constants. The symbol - means the experimental values are not available.

Figure 4. Comparison of simulated and RST isotherms for N2 in [hmim][Tf2N] (main figure) and simulated and experimental63 isotherms for N2 in [hmim][Tf2N] (inset). The solid line in the inset is a linear fit to the experiment data.

the RST results predict significantly lower solubilities than what is observed experimentally, even using the fitted solubility parameter; the absolute average difference in x1 between the RST and the revised experimental data is 35.5%. 4.1.3. N2 Solubility. Figure 4 shows the results for N2 absorption from simulation, experiment, and RST. Simulation results are also shown in Table 1, with Henry’s Law constants given in Table 2. A linear fit to the experimental data63 at fugacities below 80 bar results in an estimated Henry’s Law constant at 298 K of 848 ( 27 bar. Given the difficulty in measuring the solubility of gases such as N2 gravimetrically, the uncertainty in the experimental Henry’s Law constant is likely underestimated. The simulations yield Henry’s Law

Figure 6. The same as Figure 5 except that the fugacity is scaled by the saturation pressure.

constants of 1138 ( 47, 1410 ( 42, and 1772 ( 51 bar at 298, 333, and 373 K, respectively. The simulated Henry’s Law constant at 298 K is 34% larger than the experimental value, and the absolute average difference in mole fractions between simulation and experiment is 25.6%. The RST predictions again used the same IL parameters (V2 ) 327 cm3/mol at 298 K and V2 ) 338 cm3/mol at 333 K, δ2 ) 17.39 (J/cm3)0.5), while for N2, values of V1 ) 32.4 cm3/mol and δ1 ) 5.28 (J/cm3)0.5 were taken from the literature.42 The absolute average difference between RST predictions and experiment is 68.2%. 4.1.4. SO2 Solubility. Table 1 shows the simulated solubility for SO2, while Figure 5 shows the comparison with experiment. The agreement is quite good. Calculated Henry’s law constants at 298 and 333 K are 0.9 ( 0.3 and 3.6 ( 0.5 bar, while the experimental values are18 1.64 ( 0.01 bar at 298 K and 4.09 ( 0.06 bar at 333 K. It was shown experimentally that there is very little temperature dependence to the SO2 isotherm if the vapor fugacity is reduced by the saturation pressure. Figure 6 shows a comparison of the computed and measured solubility of SO2 as a function of reduced fugacity. Note that the simulation results were reduced by the saturation pressure calculated for the SO2 model, which were found to be51 8.6 ( 0.4 bar at 333 K and 2.1 ( 1.3 bar at 298 K. These values are slightly different from the experimental values64 of 11.038 bar at 333 K and 3.9304 bar at 298 K. Both the simulations and experiments fall on a single curve. The partial molar volume for SO2 absorption in [hmim][Tf2N] at 333 K was estimated from a linear fit to the simulated mixture molar volumes. The resulting value of 21 ( 10 cm3/mol is between the values of 40.7 ( 1.0 cm3/mol for CO2 at 333 K and 6 ( 6 cm3/mol for H2O at 298 K in the same IL.33 Since small partial molar

16716 J. Phys. Chem. B, Vol. 112, No. 51, 2008

Shi and Maginn

volumes tend to correlate with stronger interactions between the solute and solvent, these results suggest that SO2 has a binding affinity intermediate between relatively weakly associating CO2 and strongly associating H2O. No attempt was made to model SO2 solubility with RST due to the scarcity of tabulated solubility parameter data. 4.2. Mixed Gas Solubility. 4.2.1. CO2/O2 Solubility. Mixtures of CO2 and O2 were simulated at 313 K and various total pressures and partial pressures of the two gases. A typical result from a simulation showing the number of solute molecules in each phase is shown in Figure 7, where the values include “whole” and “fractional” molecules given by the value of N and λ, respectively. The last 5 million simulation steps were averaged to obtain compositions. Table 3 shows the results of simulations along with RST predictions. The selectivity for CO2 over O2 was computed as (xCO2/xO2)/(yCO2/yO2), where x and y refer to the liquid and vapor phase mole fractions, respectively. The simulations predict that the selectivity ranges from 8 to 20; experimental selectivities have not been reported. To apply RST to mixtures, the gas phase was modeled using the PReos. The activity coefficients are given by

RT ln γi ) Vi[δi - δ¯ ]2

(7)

j ) ∑j Φjδj, j ) 1, 2, 3 denote CO2, O2, and the IL, where δ respectively, and Φj ) xjVj/∑k xkVk. RST predicts selectivities of between 5 and 9, and the absolute average difference in x for CO2 between RST and simulation is 40.0%. For O2, the absolute average difference in x between RST and the simulation is 29.4%. Experimentally, Hert et al.23 (i.e., the “old” experimental data referred to above) found that CO2 absorption enhances the solubility of O2 in [hmim][Tf2N]. To quantify this, they computed an enhancement factor (EF) for O2, which is defined as the oxygen solubility in the mixed gas system divided by the oxygen solubility in the pure O2 system at the same O2 fugacity. An EF greater than unity indicates that O2 solubility increases over what would be expected from simple ideal mixing. They found that EFs for O2 ranged from 1.7 to 4.9 for the CO2/O2/[hmim][Tf2N] system. Upon reanalyzing their results and the uncertainty of the measurements, however, this group has subsequently found that the EFs are likely near unity and that most likely there is little enhancement in oxygen solubility.24 We have computed the EF from simulations in the following manner. The solubility for pure O2 absorption in the IL

Figure 7. The number of species (N + λ) for each component in the different phases obtained during the simulation of the CO2/O2 mixture at 313 K and 30 bar. The number of O2 molecules in the IL phase has been shifted upward by 50 for clarity.

corresponding to the fugacity in the gas mixture was estimated from interpolation of any two data sets in Figure 3, between which the fugacity of the O2 in the gas mixture lies. EFs were then computed by comparing the mixture solubilities with the pure gas solubilities. Interestingly, EFs from the simulations were always less than unity, ranging from 0.35 and 0.75. Thus, the simulations indicate that CO2 actually slightly reduces the solubility of O2. These findings are more in line with the reanalyzed results from the Brennecke group than the original mixture experimental results. For example, taking results in Table 3 at 313 K, 50 bar of pressure (fugacity of O2 equal to 32.92 bar), and yO2 ) 0.676, the simulation predicts that xO2 ) 0.039. This is very close the “new” experimental value at nominally the same conditions; at 313 K, 52.88 bar, and an O2 fugacity of 34.84 bar, the experiments43 indicate that xO2 ) 0.038. Likewise, at 65 bar and yO2 ) 0.677, which corresponds to a fugacity of 42.6 bar for O2, the simulations predict that xO2 ) 0.052, while the experiments indicate that at 72.95 bar and a fugacity of 42.99 bar, xO2 ) 0.050. In the “old” experiments by Hert et al.,23 it was postulated that the enhancement of O2 solubility was due to CO2 cosolubilizing the sparingly soluble O2 via favorable dispersion interactions. In that study, the concentration of CO2 in the IL phase was always quite small, ranging from a mole fraction of 0.03 to 0.09. This is because low total pressures ranging from 3 to 16 bar were used. To test whether or not dispersion interactions between CO2 and O2 would be significant under the experimental conditions, isothermal-isobaric molecular dynamics simulations were conducted in which the mole fractions of CO2 and O2 were both set to 0.051 at a temperature of 313 K. These conditions are comparable to those of the Hert et al. experiments. The interaction energy between CO2 and O2 was found to be extremely small, suggesting that it is unlikely that CO2 could interact with O2 to any appreciable extent, and certainly not enough to enhance O2 solubility by a significant amount. This, combined with the revised experimental data, suggests that CO2 does not enhance the solubility of O2 in this IL under the present conditions. 4.2.2. SO2/N2 Solubility. Mixture simulations for SO2 and N2 were carried out at 333 K, with results shown in Table 4. These calculations were motivated by the unusual finding reported by Huang et al.17 They found that pure SO2 is highly soluble in 1-n-butyl-3-methylimidazolium tetrafluoroborate [bmim][BF4], 1-n-butyl-3-methylimidazolium bis(trifluoromethylsulfonyl)amide ([bmim][Tf2N]), and three guanidiniumbased ILs. They also found that pure N2 is virtually insoluble in all five ionic liquids. In a gas mixture containing 10 mol % SO2 and 90 mol % N2, however, the experiments suggest that N2 inhibits SO2 absorption for the imidazolium-based ILs but does not alter the SO2 solubility in the guanidinium-based ILs to a great extent. Using mixtures and conditions that are similar to those used by Huang et al., we find that N2 does not absorb to any appreciable extent from the mixture and that SO2 absorption is not significantly different from what is observed with pure SO2. Since so little N2 absorbed in the liquid under these conditions, it was not possible to compute a selectivity or enhancement factor. The inhibition of SO2 absorption by N2 can be quantified by comparing the pure gas Henry’s law constant with the j SO , defined as mixture Henry’s law constant, H 2

j SO ) lim H 2 Pf0

jf SO 2 xSO2

(8)

Pure and Mixed Gas Absorption in [hmim][Tf2N]

J. Phys. Chem. B, Vol. 112, No. 51, 2008 16717

TABLE 3: The Simulated Equilibrium Mole Fractions for the Vapor (y) and Liquid (x) for the CO2/O2 Mixture in [hmim][Tf2N] at 313 K and a Series of Different Total Pressures P, along with the Mole Fractions for the Species in the IL Phase Predicted from Regular Solution Theory (RST), the Selectivity Values for CO2 over O2 Both from Simulations and RST, And the Simulated Enhancement Factor (EF) for O2 Absorptiona y

x (simulation)

x (RST)

selectivity

P (bar)

CO2

O2

CO2

O2

CO2

O2

simulation

RST

EF O2

90 120 50 30 65 50

0.630(5) 0.658 (5) 0.516 (2) 0.529 (1) 0.323 (4) 0.324 (3)

0.370 (5) 0.342 (5) 0.484 (2) 0.471 (1) 0.677 (4) 0.676 (3)

0.591 (5) 0.572 (5) 0.415 (2) 0.397 (1) 0.384 (4) 0.385 (3)

0.022 (5) 0.034 (5) 0.027 (2) 0.018 (1) 0.052 (4) 0.039 (3)

0.448 0.521 0.252 0.173 0.189 0.158

0.0431 0.056 0.029 0.017 0.050 0.039

15.8 8.7 14.4 19.6 15.5 20.6

6.1 4.8 8.2 9.1 7.9 8.5

0.35 0.49 0.52 0.69 0.75 0.64

a Uncertainties were computed as standard deviations of block averages from simulations; the numbers in parentheses are the uncertainties in the last digit.

TABLE 4: The Equilibrium Mole Fraction of the Vapor (y) and Liquid (x) for the SO2/N2 Mixture in [hmim][Tf2N] along with the Fugacity Coefficients for Components in the Gas Mixture Computed from the Peng-Robinson Equation of Statea Φigas

y

x

P (bar)

SO2

N2

SO2

N2

SO2

N2

MH (bar)

3 3 5

0.075 (4) 0.115 (1) 0.095 (1)

0.925 (4) 0.885 (1) 0.905 (1)

0.984 0.9828 0.9725

0.9996 0.9998 0.9995

0.046 (2) 0.056 (2) 0.101 (3)