Molecular Simulations and Experimental Studies of Solubility and

J. Phys. Chem. B 2010, 114, 6531–6541

6531

Molecular Simulations and Experimental Studies of Solubility and Diffusivity for Pure and Mixed Gases of H2, CO2, and Ar Absorbed in the Ionic Liquid 1-n-Hexyl-3-methylimidazolium Bis(Trifluoromethylsulfonyl)amide ([hmim][Tf2N]) Wei Shi,*,†,‡ Dan C. Sorescu,† David R. Luebke,† Murphy J. Keller,† and Shan Wickramanayake†,‡ U.S. Department of Energy, National Energy Technology Laboratory, Pittsburgh, PennsylVania 15236 and URS Corporation, South Park, PennsylVania 15129 ReceiVed: March 2, 2010; ReVised Manuscript ReceiVed: April 7, 2010

Classical molecular dynamics and Monte Carlo simulations are used to calculate the self-diffusivity and solubility of pure and mixed CO2, H2, and Ar gases absorbed in the ionic liquid 1-n-hexyl-3-methylimidazolium bis(trifluoromethylsulfonyl)amide ([hmim][Tf2N]). Overall, the computed absorption isotherms, Henry’s law constants, and partial molar enthalpies for pure H2 agree well with the experimental data obtained by Maurer et al. [J. Chem. Eng. Data 2006, 51, 1364] and the experimental values determined in this work. However, the agreement is poor between the simulations and the experimental data by Noble et al. [Ind. Eng. Chem. Res. 2008, 47, 3453] and Costa Gomes [J. Chem. Eng. Data 2007, 52, 472] at high temperatures. The computed H2 permeability values are in good agreement with the experimental data at 313 K obtained by Luebke et al. [J. Membr. Sci. 2007, 298, 41; ibid, 2008, 322, 28], but about three times larger than the experimental value at 573 K from the same group. Our computed H2 solubilities using different H2 potential models have similar values and solute polarizations were found to have a negligible effect on the predicted gas solubilities for both the H2 and Ar. The interaction between H2 and the ionic liquid is weak, about three times smaller than between the ionic liquid and Ar and six times smaller than that of CO2 with the ionic liquid, results that are consistent with a decreasing solubility from CO2 to Ar and to H2. The molar volume of the ionic liquid was found to be the determining factor for the H2 solubility. For mixed H2 and Ar gases, the solubilities for both solutes decrease compared to the respective pure gas solubilities. For mixed gases of CO2 and H2, the solubility selectivity of CO2 over H2 decreases from about 30 at 313 K to about 3 at 573 K. For the permeability, the simulated values for CO2 in [hmim][Tf2N] are about 20-60% different than the experimental data by Luebke et al. [J. Membr. Sci. 2008, 322, 28]. 1. Introduction Currently, there is a high demand to reduce the carbon dioxide (CO2) released into the atmosphere to alleviate the greenhouse effect. In the United States, greenhouse gas emissions are primarily associated with fossil fuel combustion for energy production. Energy-related CO2 emissions represented about 82% of the total U.S. anthropogenic greenhouse gas emissions in 2006.1 Goals put forth by the Department of Energy (DOE) require 90% CO2 capture at an increase in the cost of electricity of less than 35% for postcombustion capture and less than 10% for precombustion capture.2 One method that could be used for CO2 capture in the precombustion process is based on the use of supported ionic liquid (IL) membranes (SILM)s. Ionic liquids (ILs) have many unique properties such as low vapor pressure and thermal stability. These properties ensure minimal transport media loss due to evaporation or degradation that improves membrane performance. Additionally, one can chemically modify the types of cations and anions in an IL to improve the separation of the CO2 from H2. Recently, Myers et al.3 have investigated two ionic liquids, [hmim][Tf2N] and [H2NC3H6mim] [Tf2N], encapsulated in a membrane for CO2/H2 separation. The permeability * To whom correspondence should be addressed: E-mail: shiw@netl. doe.gov. † U.S. Department of Energy. ‡ URS Corporation.

selectivity of CO2 over H2 was found to be in the range 2-20 over the temperature range 310-573 K. To achieve higher permeability selectivity for CO2 over H2 and higher permeability for CO2, it is necessary to identify better existing ionic liquids or to synthesize new ones. However, an Edisonian approach would be very time-consuming and expensive due to a large number of possible ILs. Alternatively, molecular simulations can help identify promising IL candidates for SILM applications. The CO2 is absorbed in ILs through diverse mechanisms such as physical absorption,4-17 chemical reaction,18-23 and formation of molecular complexes.24-29 In the case of physical absorption, the CO2 fills the available free molar volume of the IL30,31 and interacts preferentially with the anion rather than the cation.32,33 The CO2 solubility increases for cations with a longer alkyl side chain5 and decreases at elevated temperatures.5,13,31 Hydrogen solubility in ILs has been measured experimentally14,15,34-39 by several groups. In the case of H2 absorption in [hmim][Tf2N], significant differences occur among different experimental studies. For example, Costa Gomes15 has found that the H2 solubility decreases with increasing temperature from 283 to 343 K. In contrast, an opposite solubility-temperature relationship was observed by both Noble et al.36 and Maurer et al.37 They have observed an increase in the H2 solubility with increasing temperature. However, significant differences in the solubility values occur at temperatures higher than 313 K between the data of Noble et al.36 and Maurer et al.37 These

10.1021/jp101897b  2010 American Chemical Society Published on Web 04/23/2010

6532

J. Phys. Chem. B, Vol. 114, No. 19, 2010

experimental results clearly indicate that it is difficult to accurately measure the solubility in ILs for very insoluble gases such as H2. Atomistic simulations represent an important tool to study ionic liquids and the progress of simulations in this field was summarized in several recent review articles.40-42 For example, transport properties such as self-diffusivity and viscosity have been calculated based on molecular dynamics simulations.43-46 Additionally, gas solubilities at infinite dilution have been computed using Monte Carlo methods such as the thermodynamic integration,47 Widom insertion,48 and expanded ensemble methods.49 Gibbs ensemble simulations have also been used to compute full gas absorption isotherms.50-52 Recently, Shi and Maginn have developed a continuous fractional component (CFC) method to accurately compute pure gas absorption isotherm.31,53,54 This method has also been extended to calculate mixed gases solubilities in ILs.33 The focus of this paper is on further analysis of the solubility and transport properties of H2, CO2, and Ar absorbed in the [hmim][Tf2N] based on atomistic simulations. As described above, there are significant discrepancies in the previous experiments for the H2 solubility. Moreover, the interaction mechanism between the solute molecules H2 and Ar and the [hmim][Tf2N] is not fully understood. The theoretical method used in this study is based on the use of CFC Monte Carlo (MC) simulations to compute the pure and mixed gases absorption in ILs and to compare the simulation results with the available experimental data. Molecular dynamics (MD) simulations were also performed to calculate the self-diffusivity of H2 and CO2 in [hmim][Tf2N]. On the basis of these data, the permeability values for CO2 and H2 were determined and compared against available experimental data. More importantly, the mechanism for gas absorption in ILs has been studied and potential IL candidates capable of efficient separation of CO2/ H2 gas mixtures have been suggested. 2. Experimental Details The analysis of gas sorption properties of ionic liquid [hmim][Tf2N] was carried out for both H2 and CO2 using a pressure decay sorption system (Setarm PCTPro -2000, Setaram Inc., Newark, CA). The ionic liquid [hmim][Tf2N] was purchased from EMD Chemicals Inc., NJ. Ultrahigh purity H2 (Matheson Trigas, NJ) was used as the sorbent gas. The ionic liquid sample was kept under vacuum for 12 h to remove dissolved gases before each experiment. The data were recorded for a pressure range from about 5 to 25 bar at temperatures 313 and 373 K. The collected data reveal the amount of gas sorption as a function of gas pressure at selected temperatures. Although the H2 sorption at low temperatures and pressures was small, the general features of the sorption graph, kinetic, and equilibrium values can be reproduced. At higher temperatures, enhanced H2 sorption enabled sorption values with better signalto-noise ratio. Nevertheless, the trend and values in all cases were reproducible. To verify our experimental setup, we have measured CO2 solubility in [hmim][Tf2N] at 333 K and 3-20 bar (Research grade 99.999% pure CO2 was purchased from Matheson Trigas, NJ). Our experimental data are in very good agreement with those from other experimental groups and simulations.31 3. Simulation Details 3.1. Classical Force Field. A classical force field has been used to simulate the ionic liquid [hmim][Tf2N], the solute molecules H2, Ar, and CO2, and the interaction between the

Shi et al. ionic liquid and the solute molecules. The interaction energy is given below

V (r) )

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

bonds

∑

angles

∑

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

dihedrals

kψ(ψ - ψ0)2 +

impropers

N-1

N

∑ ∑

{ [( ) ( ) ] σij rij

4εij

i)1 j)i+1

12

-

σij rij

6

+

qiqj rij

}

(1)

where the symbols have their conventional meaning.55 In the case of CFC MC simulations, the Lennard-Jones (LJ) parameters for the fractional molecule were scaled using the functional form proposed by van Gunstern et al.56 to aid insertions. The LJ potential was multiplied by a switching function of the form

{

Uswitch,LJ(rij) ) ULJ(rij)

rij < ron (roff - rij ) (roff + 2rij - 3ron ) 2

ULJ(rij) × 0

2 2

2

(roff2 - ron2)3

2

2

ron e rij e roff rij > roff

(2)

where ULJ(rij) denotes the full LJ interaction between atoms i and j separated by a rij distance. The ron denotes the starting of the switching region and roff indicates the cutoff distance where the LJ interaction is switched to zero. For the ionic liquid, the intramolecular electrostatic and LJ interactions for atoms separated by three consecutive bonds were scaled by 0.5 and were neglected for atoms separated by less than three consecutive bonds. Force field parameters for [hmim][Tf2N] and CO2 were taken from the previous work by Shi and Maginn.31,33 For Ar, a spherical LJ potential was used with σ ) 3.4 Å and ε/kB ) 120 K.57 For H2, two different potential models were chosen, a onecenter spherical LJ model developed by Buch58 and a two-center LJ model proposed by Cracknell.59 The LJ parameters were taken as σ ) 2.96 Å and ε/kB ) 34.2 K58 for the one-center Buch model. In the case of the two-center Cracknell H2 model, the LJ potential parameters are σ ) 2.59 Å and ε/kB ) 12.5 K, and the nominal bond length r0 between the two hydrogen atoms was set to be 0.74 Å. In the original Cracknell H2 model, H2 was taken as a rigid molecule. However, we chose a flexible model for H2 molecule in this work. This choice was motivated by the fact that a hybrid Monte Carlo (HMC)60 method was used for thermal equilibration moves and such a technique requires the use of a time reversible integrator. We note that traditional constraint dynamics methods for rigid molecules are not time reversible. A nominal force constant kb of 374 kcal/ (mol Å2) as determined at the B3LYP/6-311+G* level of theory was used for the stretching H-H bond. Standard LorentzBerthelot combining rules were used to calculate the mixed Lennard-Jones interaction parameters. Another problem studied in this paper was the effect of H2 and Ar solute polarizabilities on the gas solubility. The effect of the ionic liquid polarizability will be reported in future. Specifically, in the current work we assume that the electrostatic field on the solute molecule i due to the induced dipole on other solute molecules is small compared to the one due to the

Mixed Gases Absorbed in [hmim][Tf2N]

J. Phys. Chem. B, Vol. 114, No. 19, 2010 6533

permanent charges on solvent molecules. In this case, the general expression for the polarization energy given by61

Upol ) -

1 2

∑ ∑ µi,RindEi,Rstat

(3)

∑ Rpol|Ebistat|2

(4)

R

i

can be simplified as62

Upol ) -

1 2

i

ind In eq 3, µi,R denotes the induced dipole of the R (x, y, z) component for atom i, Estat i,R is the R component of the electrostatic field on atom i arising from permanent charges. In eq 4, Rpol is the polarizability. In this work, the polarizability values were taken to be 0.79 Å3 for H263 and 1.64 Å3 for Ar.62 We have used both the Ewald summation and the Fennel and Gezelter shift force (FGSF) method64 to calculate the electrostatic interactions. The electrostatic fields derived for the Ewald method by Fanourgakis et al.65 and for the FGSF method by Dyer et al.62 have both been implemented in our in-house software package. Additional details about the validation of the FGSF technique against the standard Ewald method to compute the electrostatic field can be found in Supporting Information. 3.2. Molecular Dynamics Simulations. Molecular dynamics simulations in NPT, NVT, and NVE ensembles were performed using the program NAMD.66 Whenever electrostatic interactions were computed, we used the particle-mesh Ewald method. Langevin dynamics with a Nose´-Hoover barostat/thermostat were applied for temperature and pressure controls. The barostat oscillation time and damping factors were set to be 1 ps, the thermostat damping factor was taken to be 5 ps-1, and the time step was 0.5-1 fs. The switching parameters in eq 2 were taken as ron ) 10.5 Å and roff ) 12.0 Å. For pure H2, a series of NPT MD simulations were performed at 313-573 K and 1-500 bar, for a system consisting of 200 H2 molecules. For pure Ar, NPT MD simulations were performed at 313-573 K and 1-300 bar for an ensemble of 200 Ar molecules. Fugacity coefficients for pure gases at different T and P were then calculated as described in the next section. For mixtures of H2/[hmim][Tf2N] and CO2/[hmim][Tf2N], NVE MD simulations were performed to calculate the selfdiffusivity (Ds) of the solute and of the cation and anion species of the solvent. In each NVE MD simulation, 2000 minimization steps were initially performed, followed by heating with a rate of 0.001-0.01 K per step from 0 K to the desired temperature. Equilibration and production were typically set to be 5 and 10 ns, respectively. The Ds values were determined from the mean square displacement in NVE production runs by

Ds ) lim tf∞

1 〈|r(t) - r(0)| 2〉 6t

(5)

In NVE MD simulations, the V value was set to be the average volume from NPT MD simulations, where N and P roughly correspond to the gas solubility calculations. 3.3. Monte Carlo simulations. 3.3.1. NPT f CFC Monte Carlo for Pure Gas Absorption. Pure gas absorption isotherms for CO2, H2, and Ar were computed using the CFC MC simulations. We used a locally developed code to perform these calculations and the computational details may be found

elsewhere.31,33,53 Here, we discuss only those issues specific to the current set of simulations. Calculations were performed with the temperature, pressure, fugacity of the solute gas, and the number of solvent molecules fixed. The LJ potential was truncated using the switching function given in eq 2 with ron ) 10.5 Å and roff ) 12.0 Å. A Verlet neighbor list with a 13.5 Å radius was used. Calculation of electrostatic interactions was done using the FGSF method64 with a cutoff distance of 12.0 Å and a damping parameter κ ) 0.2022 Å-1. All simulations were performed on systems containing 160 ion pairs. Equilibration runs of typically 2 million moves were carried out, during which various MC moves were tuned to achieve roughly 50% acceptance rates. These moves were thermal equilibration via HMC (45%); fractional solute coupling strength λ moves (50%); and system volume changes (5%), where the values in parentheses represent the fixed probabilities during the run. For the HMC calculations we used five MD steps per MC move. The time step was chosen between 0.4 and 0.8 fs to achieve about 50% acceptance in the MC move. The maximum values for volume changes were 400-1000 Å3. The CFC coupling parameter λ was changed uniformly up to a maximum value of 0.4-0.9 (depending on the state point and the system) to give an acceptance rate of about 50%. The CFC bias factors53 were adjusted during the equilibration stage to achieve as closely as possible a uniform probability distribution of λ. Optimization of these bias factors was done using the Wang-Landau updating scheme.67 After equilibration, production runs of several million moves were performed, with the exact number dictated by the convergence characteristics of the system. Block averages were used to compute the properties of interest and the corresponding uncertainties. The gas fugacity was computed as f ) Pφ, where φ is the pure gas fugacity coefficient at given T and P conditions. It was assumed that the ionic liquid was absent in the vapor phase, such that the total pressure was equal to the partial pressure of solute gas. The fugacity coefficient, the gas compressibility Zm ) PVm/kBT, and the molar volume Vm were computed from NPT MD simulations. A polynomial function of the form Zm ) 1 + BP + CP2 + DP3 was used to fit the Zm values at the same T but different P conditions. The fugacity was then computed as φ ) exp(BP + C/2P2 + D/3P3). The fugacity coefficient has also been computed using the Peng-Robinson equation of state (PReos). The φ values for H2 and Ar obtained from MD simulations were close to those obtained from the PReos. For H2 the average absolute difference in φ between MD simulations and the PReos is 2.4% while for Ar the same difference is 1.9%. For CO2, φ was computed only from the PReos to be consistent with previous simulations.31 In cases when polarization energy was included in simulations, a hard cutoff of 1.0 Å around each atom was used to exclude the overlap between the newly inserted fractional molecule and the existing atoms. In addition, the permanent charge and the polarization energy for the fractional molecule were scaled by a factor of λ5 (0 < λ < 1). 3.3.2. NPT Gibbs Ensemble CFC for Mixed Gas Absorption. To calculate mixed gases absorption in ILs, constant pressure CFC Gibbs ensemble simulations were performed using a locally developed code.33,54 The total number of each type of solute molecules, the number of solvent molecules, the temperature, and the pressure were fixed in these simulations. In all cases, 160 ion pairs were used. For the CO2/H2 mixed gases absorption, simulations were performed at 333-573 K and 50-60 bar. The numbers of CO2 and H2 molecules in the two

6534


Figure 1. Simulated and experimental density versus temperature for the ionic liquid [hmim][Tf2N]. The circles represent the simulated results obtained in this study. The diamond and square symbols correspond to the experimental data from Magee et al.68 and Brennecke et al.5 The line is a linear fit of the experimental data by Magee et al. The error bars from simulations and experiments are typically smaller than the symbols and are not shown for clarity.

coexisting phases were set to be 80-360 and 210-258, respectively. In the case of H2/Ar absorption, simulations were performed at 313-573 K and 100-300 bar. The corresponding numbers for H2 and Ar molecules were set to be 60-212 and 120-240, respectively. For the CO2/Ar absorption, simulations were performed at 333-573 K and 50-70 bar. Simulations typically consisted of 2 millions equilibration moves followed by 6 million production moves. The fixed probabilities for different moves were similar to those in pure gas absorption simulations. For the λ move, the probability of choosing a particular type of solute was proportional to the concentration of that type in the combined liquid and vapor simulation boxes. The pressure was fixed and the volume of each simulation box was adjusted during the simulation. Additional details about this simulation technique can be found in reference.54 4. Results 4.1. Density of [hmim][Tf2N]. The accuracy of the force field for [hmim][Tf2N] was checked by calculations of the liquid density as a function of temperature. The corresponding calculated values are represented in Figure 1 together with the available experimental data. The simulated densities at different temperatures and ambient pressure in Figure 1 are very close to the experimental values obtained by Brennecke et al.5 and Magee et al.68 For comparison, we linearly extrapolated the experimental data obtained by Magee et al.68 at lower temperatures (T < 373 K) to the range 373-573 K. The absolute average differences of density between simulations and the experimental data by Magee68 is about 0.7% in the temperature range of 298-573 K. These findings indicate that the force field parameters of the [hmim][Tf2N] are accurate enough to predict the ionic liquid density and the corresponding molar volume. 4.2. Pure Gas Solubility. 4.2.1. Compressibility and Fugacity Coefficient. The compressibilities Zm ) PVm/kBT of H2 obtained from NPT MD simulations at 313 K using Buck58 and Cracknell59 H2 potential models are shown in Figure 2. As can be seen, the computed compressibilities are very close to the experimental data. The compressibilities at the other three temperatures 373, 473, and 573 K (results not shown), also indicate a good agreement. The absolute average differences for compressibilities between our simulations using the Cracknell H2 model59 and the experimental data are 0.4-1.2% for the 313-573 K and 1-600 bar ranges of temperatures and pressures. Similarly, in the case of the Buch H2 potential model, the corresponding absolute errors are 0.2-0.3% at 313-573 K and 1-500 bar.

Shi et al.

Figure 2. The compressibility factor obtained from NPT molecular dynamics simulations at 313 K. For comparison, the experimental data69 are also shown.

As shown in Table 1, the calculated fugacity coefficients at 313 K using the Cracknell H2 model are also very close to the experiment data.69 The largest difference between simulations and experiments is 2.8%, and the absolute average difference is 0.9%. At the other three temperatures of 373, 473, and 573 K, the absolute average differences are 0.6-0.9% for the pressure range 1-600 bar. When the Buch H2 potential model is used, the absolute average difference between simulated and experimental data for the fugacity φ is 0.2-0.5% at 313-573 K and 1-500 bar. For Ar, the computed compressibilities are also very close to the experimental data69 at 313-573 K and 1-300 bar. The absolute average differences for Zm are 0.4-0.5% and for fugacity coefficients 0.2-0.8% between simulations and experiments.69 4.2.2. H2 Solubility. The simulated H2 solubilities using different potential models are summarized in Table 2. The simulated solubility at 313 K is shown in Figure 3. In the same figure, we have also included the experimental data obtained by Maurer et al.,37 Noble et al.,36 Costa Gomes,15 and the experimental values obtained by our group. At 313 K, an experimental Henry’s law constant (H) of 1730 bar has been estimated based on a linear interpolation of the ln(H) versus 1/T data obtained by Maurer et al.37 between 293 and 413 K. As shown in Figure 3, the simulated solubilities for both potential models are very close to all four sets of experimental data at pressures below 100 bar. However, at high pressure of 200 bar, the best agreement between the calculated and experimental data is obtained for Cracknell model while larger differences were observed in the case of Buch model. In this figure, the indicated lines correspond to experimental solubilities extrapolated from lower pressure data, to values of less than 200 bar. The absolute average differences for solubilities between simulations using the Cracknell H2 potential model59 and the experimental data by Maurer et al.,37 Noble et al.,36 and Costa Gomes15 are 10, 22, and 8%, respectively. In the case of Buch H2 potential model,58 the absolute average solubility differences between simulations and the same sets of experiments are 20, 44, and 33%, respectively. The isotherms at 373 K are shown in Figure 4. Note that Noble et al.36 have not measured the H2 solubility at this temperature. Instead, they measured the solubility between 298-343 K. Hence, the Henry’s law constant 454 bar at 373 K from Noble et al.36 was estimated by extrapolating their ln(H) versus 1/T data between 298 and 343 K. As shown in Figure 4, the simulated solubilities are very close to the experimental data by Maurer et al.,37 particularly for the set of values obtained using Cracknell potential. The simulated solubilities are also



TABLE 1: The Simulated Molar Volume Vm, the Compressibility Zm, and the Fugacity Coefficient O Obtained from NPT MD Simulations at 313 K Using the Cracknell H2 Potential Model59a simulated (Cracknell model)59 T (K)

P (bar)

Vm (cm3/mol)

313 313 313 313 313 313 313 313 313 313

1 25 50 100 150 200 250 300 400 500

26065 (65) 1055 (3) 534 (1) 274.4 (0.3) 187.0 (0.4) 143.8 (0.2) 117.7 (0.2) 100.5 (0.2) 78.9 (0.1) 65.9 (0.1)

exp

Zm 1.001 1.013 1.025 1.054 1.077 1.104 1.130 1.158 1.212 1.265

(0.002) (0.003) (0.002) (0.001) (0.002) (0.002) (0.002) (0.002) (0.002) (0.002)

φ 1.001 1.013 1.026 1.053 1.080 1.109 1.138 1.169 1.232 1.299

(0.001) (0.001) (0.003) (0.004) (0.004) (0.005) (0.006) (0.006) (0.006) (0.007)

Vm (cm3/mol)

Zm

φ

26056 1056.5 535.8 275.5 188.9 145.6 119.7 102.5 80.9 68.0

1.0 1.014 1.029 1.058 1.088 1.118 1.149 1.181 1.243 1.306

1.0 1.007 1.021 1.051 1.082 1.115 1.148 1.183 1.257 1.336

a The experimental values69 are also included for comparison. Note that the experimental Zm and φ values were computed using the experimental Vm data. The uncertainties for the simulated values are given in parentheses.

TABLE 2: Simulated Molar Fractions x and Mixture Molar Volumes Vmix at Different Temperatures T, Pressures P, and Fugacities f for Pure H2 Absorption in the Ionic Liquid [hmim][Tf2N]a H2 model

T (K)

P (bar)

f (bar)

x

Vmix (cm3/mol)

C C C C C C C C C C C C C C C C C C C C C C B B B B B B B B B B B B B B B

313 313 313 313 313 313 313 373 373 373 373 373 373 373 473 473 473 573 573 573 573 573 313 313 313 313 313 313 373 373 373 373 573 573 573 573 573

25 50 100 200 300 400 500 25 50 100 200 300 400 500 50 100 200 25 50 100 200 300 50 100 200 300 400 500 50 100 200 300 25 50 100 200 300

25.325 51.3 105.3 221.8 350.7 492.8 649.5 25.3 51.2 104.8 219.4 344.7 481.2 630 50.65 103.1 215 25.125 50.6 102.8 213.2 330.6 51.4 105.7 223.8 356.1 504.0 669.5 51.25 105.2 221.6 350.4 25.225 50.95 103.9 215.6 335.4

0.0175 (9) 0.033 (1) 0.063 (3) 0.126 (3) 0.165 (4) 0.209 (3) 0.231 (4) 0.0174 (7) 0.034 (1) 0.064 (2) 0.127 (3) 0.144 (4) 0.190 (6) 0.221 (4) 0.045 (1) 0.075 (1) 0.153 (3) 0.0285 (5) 0.054 (1) 0.111 (2) 0.194 (2) 0.230 (6) 0.024 (1) 0.050 (1) 0.081 (3) 0.135 (3) 0.147 (3) 0.195 (2) 0.0230 (6) 0.048 (1) 0.086 (3) 0.131 (4) 0.0265 (3) 0.051 (1) 0.106 (2) 0.186 (3) 0.247 (4)

324.9 (4) 318.7 (3) 309.1 (9) 289.8 (12) 275.2 (13) 260.7 (9) 252.2 (12) 336.4 (3) 329.9 (3) 320.3 (7) 298.2 (9) 290.9 (11) 275.1 (18) 263.9 (14) 354.3 (3) 340.0 (6) 313.4 (10) 390.1 (5) 379.1 (6) 360.1 (5) 325.7 (7) 303.4 (16) 320.7 (3) 312.3 (4) 301.0 (7) 284.2 (10) 278.9 (10) 263.7 (8) 331.9 (3) 324.0 (4) 310.6 (8) 295.9 (11) 391.4 (5) 382.5 (5) 363.5 (7) 330.9 (9) 304.2 (8)

a The Cracknell (C)59 and Buch (B)58 H2 potential models were used. The numbers in parentheses are uncertainties in the last digit obtained from standard block averages.

close to the experimental values obtained by us. In contrast, simulations predict solubilities about 300-500% smaller than the extrapolated experimental data by Noble et al.36 This may be partly due to our extrapolation of their experimental Henry’s law values from low to high temperatures. Additional simulations were also performed at 473 and 573 K and a similar good

Figure 3. The simulated and experimental isotherms for H2 absorption in [hmim][Tf2N] at 313 K. The circles and diamonds represent simulated results using the Cracknell59 and Buch58 H2 potential models. The error bars from simulations are typically smaller than the symbols and are not shown in the figure. The cross symbols correspond to the experimental data obtained by us. The red solid line, the blue dotdashed line, and the green dotted line correspond to the experimental data obtained by Maurer et al.,37 Noble et al.,36 and Costa Gomes.15 For improved clarity the absorption isotherms at low pressures are also provided in the inset figure.

Figure 4. The simulated and experimental isotherms for H2 absorption in [hmim][Tf2N] at 373 K. The symbols and legends are the same as in Figure 3 except the triangle symbols which represent the experimental data by Maurer et al.37

agreement with extrapolated Maurer et al.37 data was obtained. The absolute average differences were 10-18%. Our predicted solubilities at 473 and 573 K are significantly different from the extrapolated experimental data obtained by Noble et al.36

6536


Figure 5. ln(H) versus 1/T from simulations and experiments for H2 absorption in [hmim][Tf2N]. The simulated Henry’s law constants were obtained using the Cracknell H2 potential model59 (filled circles). The experimental Henry’s law constants are from Maurer et al.37 (triangles), Noble et al.36 (diamonds), and Costa Gomes15 (squares). The solid line is a linear fitting of the simulation values. The error bars for both simulations and experiments are also shown.

Partial molar enthalpy has also been computed. In this case, the Henry’s law constant (H) was calculated from a linear fit of the fugacity f to the solubility x at low pressures. The simulated H values were estimated to be 1552 ( 66 bar at 313 K, 1505 ( 31 bar at 373 K, 1126 ( 25 bar at 473 K, and 936 ( 34 bar at 573 K, respectively. The simulated and experimental H values are shown in Figure 5. For convenience, the experimental H constants having units of MPa/(mol · kg-1) was converted to bar units. Note that a larger H constant corresponds to smaller solubility. When the slope of ln(H) versus 1/T is positive, this indicates that H decreases when the temperature is increased, that is, the solubility increases with temperature. Alternatively, when the slope is negative, the corresponding solubility decreases with temperature increase. The infinite partial molar enthalpy was calculated from the linear fitting of ln(H) versus 1/T (solid line in Figure 5) and a simulated value of 3.0 ( 0.6 kJ/mol was determined. This value is close to the experimental one of 4.09 ( 0.12 kJ/mol obtained by Maurer et al.,37 but it is significantly different from the experimental one of 16 ( 8 kJ/mol obtained by Noble et al.36 Both the calculated molar enthalpy and the two experimental values obtained by Maurer and Noble are positive, implying that the H2 solubility in [hmim][Tf2N] increases with temperature. Our own experimental data also support this trend. In contrast, as shown in Figure 5, Costa Gomes15 found experimentally that the H2 solubility decreases when the temperature is increased. This finding is also contradicted by the other experimental measurements.39 Raeissi and Peters have measured experimentally the H2 solubility in similar ionic liquids, that is, [bmim][Tf2N] and [emim][Tf2N], and they have found that H2 solubility39 increases with temperature over 333-453 K. This is in agreement with our simulations and the three sets of experiments by Maurer,37 Noble,36 and our group. Another property of interest is the enthalpy of mixing. This was calculated as ∆Hen ) Hen,mixture - xHen,H2 - (1 - x)Hen,IL, where Hen,mixture is the enthalpy of the mixture, and Hen,H2 and Hen,IL are the enthalpies of pure H2 and of pure ionic liquid. The Hen,mixture value was obtained from gas solubility calculations, while in the case of the pure gas and of the ionic liquid, the corresponding enthalpies were obtained from NPT MD simulations. As shown in Figure 6, the ∆Hen values are very small and both negative and positive values were determined.

Shi et al.

Figure 6. Simulated enthalpy of mixing versus solubility for H2 absorption in [hmim][Tf2N] using the Cracknell H2 potential model.59 The error bars from simulations are also shown.

Figure 7. (a) The computed mixture molar volume for H2 absorption in [hmim][Tf2N] at 313-573 K using the Cracknell H2 potential model.59 The lines represent linear fitting of the simulation data. (b) The variation of the volume expansion upon H2 absorption relative to the pure [hmim][Tf2N] as function of H2 concentration. Also shown are error bars from simulations.

These results suggest that the interaction between H2 molecules and the [hmim][Tf2N] is very weak. The molar volume versus concentration for H2 absorption in [hmim][Tf2N] is shown in Figure 7a. The partial molar volume was computed based on a linear fitting of the simulated data and the corresponding values are -4.9 ( -4.2 cm3/mol at 313 K, -12.6 ( 2.8 cm3/mol at 373 K, -1 ( 30 cm3/mol at 473 K, and -13 ( 20 cm3/mol at 573 K. The partial molar volumes are negative. These findings are in contrast to CO2 absorption in [hmim][Tf2N], when a positive partial molar volume31 was found. In Figure 7b, we show the volume expansion as a function of H2 concentration. Most of the ionic liquid volume expansion values are negative. These results suggest that H2 is very sparsely soluble in [hmim][Tf2N] and high pressures are required to achieve appreciable H2 solubility. At these high pressures, a small compression of the system takes place leading to a negative volume expansion and negative partial molar volumes. The free molar volume of the ionic liquid and the interaction between the gas solute and the ionic liquid are two important factors that determine the gases solubility. As mentioned above, the interaction between the H2 and the [hmim][Tf2N] is very weak. This implies that the free molar volume of the ionic liquid is an important parameter which determines the H2 solubility. To analyze this hypothesis further, we present the variation of the molar volume of the ionic liquid versus the Henry’s law


Figure 8. The calculated Henry’s law constants versus the molar volume for H2 absorption in [hmim][Tf2N] using the Cracknell H259 model at 313, 373, 473, and 573 K. The simulated pure ionic liquid molar volumes are used in this figure. The error bars from simulations are also shown.

constant for H2 absorption in [hmim][Tf2N] at different temperatures in Figure 8. As can be seen, the molar volume for [hmim][Tf2N] is indeed highly correlated with solubility. The Henry’s law constant decreases linearly with the molar volume. Note that for data in Figure 8, the molar volume instead of the molar free volume for the ionic liquid was used. This is because these two quantities are highly correlated and calculation of the free molar volume of the ionic liquid is nontrivial. When the temperature is increased, the molar volume of the ionic liquid expands, leading to high H2 solubilities and low H values. This simple relationship between the molar volume and the solubility can also be used to interpret the H2 solubility in other ionic liquids. For example, the H2 solubility in ionic liquids such as [bmim][PF6] and [bmim][BF4] have been experimentally determined14,35,38 and larger values have been obtained for the former case. This is partly due to a larger molar volume of 205 cm3/mol for [bmim][PF6] than the value of 190 cm3/mol for [bmim][BF4], at ambient conditions. To compare the H2 solubility difference between the [hmim][Tf2N] and [bmim][PF6], solubility calculations for H2 in [bmim][PF6] have also been performed. The Henry’s law constant was calculated to be 4439 ( 129 bar at 313 K, a value very close to the experimental one of 4152 ( 46 bar as determined by Maurer.38 The H value for H2 absorption in [bmim][PF6] is about 2.9 times larger than that for H2 absorption in [hmim][Tf2N]. This result is also consistent with a smaller molar volume of [bmim][PF6] than [hmim][Tf2N]. The H2 polarizability effect on the solubility has also been evaluated. For simplicity, the one-center Buch H2 potential model was used in this case. For pure H2, the polarizability does not affect the thermodynamics and the transport properties because there is no permanent charge in the Buch H2 potential model. As a result, there are no contributions to the total energy and total force due to the H2 polarizability, as shown in eq 4. However, in the case of H2 absorption in the ionic liquid, the H2 polarizability will contribute to both the energy and the electrostatic force. The H2 solubilities at 313 K with and without polarizability contributions included are shown in Figure 9. As can be seen the H2 polarizability effect on H2 solubility is very small. A similar behavior was also observed at two other temperatures of 373 and 573 K. 4.2.3. CO2 and Ar Solubilities. The CO2 solubility in [hmim][Tf2N] at 573 K was computed and the corresponding H was estimated to be 335 ( 7 bar. This value together with


Figure 9. The simulated H2 solubility in [hmim][Tf2N] at 313 K with (squares) and without (circles) polarizability using the Buch H258 potential model. The error bars from simulations are also shown.

Figure 10. A van’t Hoff plot of the computed and experimental Henry’s law constants for CO2 absorption in [hmim][Tf2N]. The experimental data (triangles) are from Shiflett and Yokozeki.13 The circle symbols indicate the simulation data from Shi and Maginn.31 The square symbols represent the simulation results obtained in this work. The line is a linear fitting of the experimental data.

previous results obtained by Shi and Maginn31 at 333, 373, and 423 K are indicated in Figure 10 and Table 3. Both previous calculated data at low temperatures and current simulations at high temperature agree with the experimental data or with the extrapolated ones obtained by Shiflett et al.13 For Ar, the simulated solubilities are summarized in Table 3 and the corresponding absorption isotherm at 313 K is shown in Figure 11. The H values for Ar were computed to be 449 ( 24 bar at 313 K, 622 ( 27 bar at 373 K, and 608 ( 32 bar at 573 K, respectively. As shown in Figure 11 and in Table 2 and Table 3, Ar exhibits solubilities about 1.5-3 times higher than H2 at 313-573 K. These results are similar to the experimental solubilities for Ar and H2 in other ionic liquids. For example, Jacquemin et al.14,35 have found experimentally that Ar exhibits about 1.6-4.2 higher solubility than H2 in both [bmim][PF6] and [bmim][BF4] systems. 4.2.4. Energy Analysis. Another parameter of interest is the interaction energy between the solute and the IL. This quantity has been determined by summing directly the interactions between the solute and IL in the primary cell and in the 124 neighboring cells. The configurations in these calculations were obtained from NPT MD simulations where the T, P, and solubility values were close to those used in the Monte Carlo solubility calculations. The interaction energy is shown in Figure

6538


Shi et al.

TABLE 3: Simulated Molar Fractions x, and Mixture Molar Volumes Vmix at Different Temperatures T, Pressures P, and Fugacities f for Pure CO2 and Ar Absorbed in the Ionic Liquid [hmim][Tf2N]a gas

T (K)

P (bar)

f (bar)

CO2 CO2 CO2 Ar Ar Ar Ar Ar Ar Ar Ar Ar Ar Ar

573 573 573 313 313 313 373 373 373 373 573 573 573 573

8 17.15 83.6 25 50 100 25 50 100 200 25 50 100 200

7.976 17.04024 81.29264 24.75 49.1 96.7 24.95 49.85 99.6 200.4 25.15 50.6 102.5 211.2

Vmix (cm3/mol)

x 0.031 0.052 0.245 0.060 0.109 0.180 0.037 0.080 0.137 0.207 0.045 0.083 0.146 0.287

(2) (3) (6) (2) (4) (4) (3) (2) (3) (3) (2) (3) (2) (5)

393.2 (10) 383.7 (11) 319.2 (23) 310.6 (6) 296.3 (12) 274.1 (13) 329.1 (8) 317.0 (9) 297.9 (12) 273.8 (11) 387.2 (7) 372.3 (7) 346.4 (9) 299.0 (13)

a The numbers in parentheses are uncertainties in the last digit obtained from standard block averages.

Figure 11. The simulated isotherms for Ar and H2 (using the Cracknell H2 potential model) absorption in [hmim][Tf2N] at 313 K. The squares are for Ar and the circles are for H2. Also shown are the error bars from simulations.

Figure 12. The interaction energy for each solute molecule due to all the other solute molecules (SS), cations (SC), and anions (SA). Results are from averages at low solute solubilities. The mole fraction is about 0.12 for H2, 0.13 for Ar, and 0.15 for CO2. The temperatures are 313 K for Ar and H2 and 333 K for CO2. Note that the results for CO2 are taken from the previous simulations33 and they are included here for comparison. For H2, the Cracknell potential model is used.59

12. It was found that the CO2 interaction with the IL is about two times stronger than the one for Ar/IL interaction, and about six times stronger than the one for H2/IL. These interaction

energies are consistent with the calculated solubilities, that is, CO2 exhibits the largest solubility while the smallest is obtained for H2. 4.3. Mixed Gas Solubilities. 4.3.1. H2/Ar Solubility. The solubilities for mixed gases of H2 and Ar in [hmim][Tf2N] were computed at 313-573 K and 100-300 bar. The results are shown in Table 4. The selectivity for Ar over H2 was computed as (xAr/xH2)/(yAr/yH2), where x and y are the mole fractions in the liquid and gas phases. Simulations give selectivity values of 1.4-3.3, very close to the Henry’s law constant ratio values of 1.5-3.5 between the pure H2 and pure Ar. The enhancement factor, defined as the solubility of the solute in the mixed gas system divided by the pure gas solubility, has also been computed at the same temperature and fugacity. The solubility was found to decrease by about 10-35% for H2 and 2-20% for Ar due to the presence of each other solute molecule. These enhancement factors correspond to solubilities of 0.026-0.122 for H2 and 0.086-0.179 for Ar. Typically, the supported ionic liquid membrane experiments are performed at ambient pressures.3,70 This implies very low solubilities for both the H2 and Ar in the ionic liquid. Under these conditions, the solubilities of these two gases are not significantly influenced by one another. If a membrane experiment is performed at somewhat higher pressures corresponding to high solubilities for H2 and Ar, one needs to consider the effect on the solubility due to the presence of the two solute gases. 4.3.2. CO2/H2 Solubility. Simulated solubility results for mixed gases of CO2/H2 in [hmim][Tf2N] are shown in Table 5. The selectivity of CO2 over H2 at 573 K is close to H ratio between the pure H2 and CO2. Interestingly, the solubility selectivity decreases significantly from about 30 at 333 K to about 3 at 573 K. This is due to opposite trends, that is, increase of H2 solubility and decrease in CO2 solubility when the temperature is increased. This decreased selectivity with increased temperature partly contributes to the decrease of the permeability selectivity at elevated temperatures in membrane separation. These results are fully consistent to experimental data where the permeability selectivity3 for CO2 over H2 has been found to decrease from about 9 at 310 K to about 1.3 at 573 K. 4.3.3. CO2/Ar Solubility. Simulated solubility results for mixed gases of CO2 and Ar in [hmim][Tf2N] are shown in Table 6. The selectivity of CO2 over Ar decreases from about 9 at 333 K to about 2 at 573 K. The selectivity value at 573 K is very close to the Henry’s law constant ratio of 1.8 between the pure Ar and CO2. 4.4. Self-Diffusivity Coefficients. The self-diffusivity was calculated as (1/6)[d(∆r2)/dt] from NVE MD simulations using a procedure similar to that developed by Maginn et al.44 In this case, the trajectories were translated to the same origin before the mean-squared displacement ∆r2 was calculated. The corresponding results are shown in Table 7. We note that to get reliable estimates for the self-diffusivity values from MD simulations, one needs to perform NVE MD simulations rather than the Langevin NVT MD simulations. We find that the Ds values for the [hmim][Tf2N] from Langevin NVT MD simulations are typically about 2-4 times smaller than the values from the NVE microcanonical MD simulations. At 313-573 K, both the solute and the ionic liquid typically exhibit diffusive dynamics with β ) d(log∆r2)/dt values close to 1. Hence, simulations are long enough to obtain reliable Ds values at these temperatures. At all temperatures, it was found that the cation diffuses faster than the anion as shown in Table



TABLE 4: The Equilibrium Mole Fraction of x in the Liquid and y in the Vapor for Mixed Gases of H2/Ar Absorbed in [hmim][Tf2N]a y T (K)

P (bar)

313 313 313 573 573

100 200 300 200 300

x Ar

H2 0.492 0.490 0.488 0.479 0.480

(2) (2) (1) (1) (1)

0.508 0.510 0.512 0.521 0.520

EF

H2 (2) (2) (1) (1) (1)

0.026 0.041 0.066 0.097 0.122

Ar (2) (2) (2) (2) (2)

0.086 0.141 0.156 0.150 0.179

(2) (5) (3) (2) (2)

selectivity

HH2/HAr

H2

Ar

3.2 3.3 2.3 1.4 1.4

3.5 3.5 3.5 1.5 1.5

0.79 0.65 0.68 0.89 0.8

0.80 0.80

Molecular Simulations and Experimental Studies of Solubility and

Recommend Documents