Molecular Simulations of CO2, H2, H2O, and H2S Gas Absorption into

Jul 28, 2015 - *Phone: +1 412 386 4406. Fax: +1 412 386 5990. E-mail: [email protected]. ... Interfacial and bulk properties of vapor-liquid equilibri...
0 downloads 10 Views 2MB Size
Article pubs.acs.org/JPCC

Molecular Simulations of CO2, H2, H2O, and H2S Gas Absorption into Hydrophobic Poly(dimethylsiloxane) (PDMS) Solvent: Solubility and Surface Tension Wei Shi,*,†,‡,§ Nicholas S. Siefert,† and Bryan D. Morreale† †

U.S. Department of Energy, National Energy Technology Laboratory, Pittsburgh, Pennsylvania 15236, United States AECOM, South Park, Pennsylvania 15129, United States § Department of Chemical and Petroleum Engineering, University of Pittsburgh, Pittsburgh, Pennsylvania 15261, United States ‡

S Supporting Information *

ABSTRACT: Henry’s law constants were calculated for H2S, CO2, H2O, and H2 gas absorption in the hydrophobic poly(dimethylsiloxane) (PDMS) solvent using an all-atom (AA) PDMS model. Calculations show that the relative gas solubility at 298 K decreases in the following order: H2S (147) > CO2 (19) ≈ H2O (15) > H2 (1). Both quantum ab initio (AI) and classical force field (FF) gas-phase calculations show that these gases interact with the PDMS molecule in the order of H2S > CO2> H2; they decrease in the same order as gas solubility. The AA PDMS model gives CO2 solubility and PDMS surface tension values close to the experimental data, with differences of 14 and 8%, respectively. In addition, by using both the all-atom and united-atom PDMS models, our simulations suggest that it is challenging to develop a solvent which both has a significantly large surface tension and exhibits large CO2 solubility at high CO2 pressure. Finally, gas absorption effects on PDMS surface tension were investigated. CO2 absorption was simulated to decrease the solvent surface tension by 3 × 10−3−4 × 10−3 N/m compared to the simulated neat PDMS solvent surface tension value of 21 × 10−3 N/m; CO2 molecules exhibit the largest concentration in the gas−liquid interface region. In contrast, H2S absorption does not decrease PDMS surface tension, which is partially due to the strong H2S−PDMS interaction compared to the CO2−PDMS interaction.

1. INTRODUCTION Chemical and physical absorption methods have been suggested for CO2 postcombustion and precombustion capture, respectively, due to low CO2 partial pressure (∼0.1 bar) in the postcombustion gas stream and high CO2 partial pressure in the precombustion capture gas stream1,2 (∼25 bar). Due to rich H2O concentration in the hot or warm high-pressure gas stream corresponding to CO2 precombustion capture, hydrophobic solvents are preferred for selective CO2 absorption over H2 and H2O.3 CO 2 and H 2 absorption in hydrophobic poly(dimethylsiloxane) (PDMS) solvents (2−50 repeating [−Si(CH3)2O−] monomers) have been experimentally investigated.3 There are several advantages for this type of polymerlike solvent, including hydrophobicity, low volatility, low viscosity, thermal stability, and commercial availability. CO2 was found to be much more soluble than H2 in this type of solvent.3 Despite the encouraging experimental results, some questions need to be addressed on a molecular level. Some examples of those questions are why are these solvents hydrophobic? Why are they CO2-selective? Considering the high CO2 partial pressure in precombustion, does the solvent volume expand significantly under this condition? Can they © XXXX American Chemical Society

absorb H2S for gas-sweetening applications? Additionally, solvent surface tension, which is one of the important factors in determining solvent foaming in absorption/desorption operations4−6 and the gas−liquid contact area in gas mass transport,7 has not been studied for the PDMS solvent. Furthermore, how does gas absorption (such as CO2 and H2S) affect the PDMS solvent surface tension? How do gas molecules behave in the gas−liquid interface region? The goal of this work is to perform molecular simulations to address the above questions. Before presenting our new work, we briefly summarize previous theoretical research. Sok et al.8 have calculated CH4 and He solubility, diffusivity, and permeability in the PDMS polymer (30 repeating [−Si(CH3)2O−] monomers) by using the united-atom PDMS model. By using an all-atom PDMS model derived from quantum ab initio calculations, Smith et al.9 have calculated the PDMS liquid (molecular weight between 310 and 1571 g/mol) X-ray structure, density, and heat of vaporization; the simulated values are in good agreement with the experimental data. Received: June 17, 2015 Revised: July 25, 2015

A

DOI: 10.1021/acs.jpcc.5b05806 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

model in the UA model. The PDMS molecule contains 177 atoms in the AA PDMS model and 69 united atoms in the UA model (Figure 3). The classical FF parameters for the UA and AA PDMS models were mainly obtained from previous work,8,9 and they are given in the Supporting Information. 2.2. Quantum Ab Initio (AI) Calculations. To evaluate the accuracy for CO2, H2, H2O, and H2S gas interactions with the PDMS molecule obtained from classical FF calculations, quantum AI gas-phase calculations have been performed for the gas−PDMS dimers by using the Gaussian 09 program.15 Note that in this work each PDMS molecule consists of 16 [−Si(CH3)2O−] monomers and contains 177 atoms when using the AA PDMS model. It would be computationally expensive to perform AI calculations for such a large PDMS molecule. Instead, a simple model of the PDMS molecule (PDMS-2) was used in AI calculations, which contains only two repeating [−Si(CH3)2O−] monomers and is terminated by −CH3 and −Si(CH3)3 groups. Optimizations were performed at the B3LYP/6-311++G(d,p) level of theory followed by single-point energy calculations using the MP2/cc-pVTZ method to compute the gas−PDMS-2 interaction. To account for the basis set superposition error, we have used counterpoise corrections in both geometry optimization and single-point energy calculations. The gas−PDMS-2 interaction was calculated to be the energy required to separate two monomers in the optimized dimer structure to an infinitely large distance while keeping the monomer structures frozen at their geometries in the optimized dimer system. 2.3. Gas Henry’s Law Constant Calculations. Although the continuous fraction component (CFC) Monte Carlo (MC) method (below) could be used to obtain the full gas absorption isotherm at different pressures, it may encounter practical difficulty when applied to certain systems. For example, if the CFC method is used to compute H2O solubility at 0.015 bar (roughly 50% water saturation pressure) and 298 K in PDMS, then about 10 000 PDMS solvent molecules are needed to obtain a statistically meaningful number of H2O molecules. This is due to small H2O solubility in PDMS and low water pressure in the gas phase. The large number of solvent molecules in turn makes the CFC MC calculation very expensive. To calculate H2O solubility in PDMS, the test-particle insertion method, which has been implemented in our in-house software,17 was used here to estimate the Henry’s law constant. The cubic simulation box contains 10 PDMS molecules, and the side length for the simulation box was set to be about 28 Å. When generating snapshots for the Henry’s law constant calculation, three independent sets of NPT molecular dynamics (MD) simulations were performed by starting from different initial configurations. Simulations were carried out at 298 K and 1 bar by using the NAMD program.18 For each set of NPT MD simulations, about 20 000 equilibrated liquid structures corresponding to the last 10 ns of the production run were saved for a later test-particle insertion calculation. Similarly, the Henry’s law constants for CO2 and H2S in PDMS were also calculated by using the test-particle insertion method. Note that it would be interesting to study the system size effect on the gas Henry’s law constant, which is not investigated here because it is out of the scope of this work. We expect that the Henry’s law constant differences among CO2, H2, H2O, and H2S obtained by using the simulation box with a side length of 28 Å would be similar to those differences corresponding to larger simulation box sizes.

Additionally, the simulated self-diffusivity and viscosity for the PDMS liquid are also in reasonable agreement with the experimental data.9 By using a similar all-atom PDMS model as used by Smith,9 Ismail et al.10 have calculated the surface tension for PDMS (20−300 repeating units) and the contact angle for H2O adsorption on the PDMS polymer surface. The literature survey indicates that there are no theoretical studies for CO2, H2, H2O, and H2S absorption in PDMS solvent. In this work, we have calculated gas Henry’s law constants and full absorption isotherms for these gases in PDMS solvent from Monte Carlo simulations. We have also calculated the surface tension values for both neat PDMS solvent and CO2−PDMS/H2S−PDMS mixtures. Two different PDMS models (all atom and united atom) were used, which allow us to investigate the important factors to determine gas solubility at low and high CO2 pressures and the important factor to determine surface tension.

2. SIMULATION DETAILS 2.1. Classical Force Field (FF) and Two PDMS Models. The classical FF potential used to simulate the PDMS solvent, CO2, H2, H2O, and H2S gases and their interactions with PDMS is given by =(r) =



kr(r − r0)2 +

bonds

+



angles

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

dihedrals N−1

+

∑ i=1

kθ(θ − θ0)2





k ψ (ψ − ψ0)2

impropers

⎧ ⎡⎛ ⎞ ⎛ σij ⎞6 ⎤ qiqj ⎫ σij ⎪ ⎪ ∑ ⎨4ϵij⎢⎜⎜ ⎟⎟ − ⎜⎜ ⎟⎟ ⎥ + ⎬ rij ⎪ ⎪ ⎢⎣⎝ rij ⎠ ⎝ rij ⎠ ⎥⎦ j=i+1 ⎩ ⎭ N

12

(1)

where the symbols represent their conventional meanings.11 Standard Lorentz−Berthelot combining rules were used to calculate the mixed Lennard-Jones (LJ) interaction parameters. The LJ potential was switched from 10.5 to 12.0 Å. A Verlet neighbor list with a 13.5 Å radius was used. The intramolecular electrostatic and LJ interactions for atoms separated by exactly three consecutive bonds were scaled by 0.5 and were neglected for atoms separated by less than three consecutive bonds. The classical FF parameters for CO2,12 H2,13 and H2O12 were taken from previous work. For H2S, the classical FF parameters obtained by Potoff and co-workers14 were used except that the H−S bond in this work was set to be flexible; the kr force constant was obtained from quantum ab initio gasphase calculations at the MP2/aug-cc-pVQZ level of theory by using the Gaussian 09 program.15 The FF parameters for H2S were set to be ϵS = −0.5524435 kcal/mol, σS = 3.71 Å, qS = −0.252 e; ϵH = 0 kcal/mol, σH = 0 Å, qH = 0.126 e; kr,H−S = 284 kcal/mol Å2, r0,H−S = 1.34 Å; and kθ,HSH = 32.57 kcal/mol rad2, θ0,HSH = 92.5°. To compare with the experimental CO2 solubility in PDMS with an average molecular weight of 1250 g/mol,3 the PDMS molecule in this work was set to be CH3[−Si(CH3)2O−]16Si(CH3)3, which corresponds to a molecular weight of 1274.72 g/ mol, close to the experimental value of 1250 g/mol. The PDMS molecule was built by using the build/build polymers/ homopolymer module in the Materials Studio software.16 To model the PDMS molecule, both an all-atom (AA) and a united-atom (UA) model were used. The −CH3 group in the AA model was represented by a single carbon united-atom B

DOI: 10.1021/acs.jpcc.5b05806 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

⎛ ∂F ⎞ γ=⎜ ⎟ ⎝ ∂A ⎠ N , V , T F(N , V , T , A + ΔA) − F(N , V , T , A) = lim ΔA → 0 ΔA

When calculating the Henry’s law constant from the above saved structures, two million trial insertions were attempted for each snapshot. To improve the insertion efficiency, the cubelet19 and bias-removing methods17 were used. The cube length was roughly set to be ∼0.8 Å. A cubelet was marked as occupied if the distance between the center of the cubelet and any solvent atom was less than rev = s × [(σi/2) + rprobe)], where σi is the van der Waals parameter for atom i of a PDMS molecule and s is the scaling factor and was chosen to be 0.8. The radius for the probe atom, rprobe, was set to be 1.0 Å. More calculation details for the test-particle insertion method can be found in our previous work.17 The Fennel and Gezelter shift force (FGSF)12,20 method (Supporting Information) was used to compute the electrostatic interaction in the test-particle insertion procedure to speed up calculations. Note that for H2S and H2O, the O and S atoms were inserted and checked as to whether they were in the empty cubelets. The H atoms were not checked because they do not interact with solvent molecules in the model; the ϵ values in LJ potential for the H atoms of H2S and H2O were both set to be zero. For CO2, both C and two O atoms were checked as to whether they were in the empty cubelets before the CO2−solvent interaction energy was calculated. 2.4. Gas Solubility Calculations. CO2 and H2 solubilities in PDMS solvent at 298 K and different pressures were obtained from the CFC MC calculation12,13,21 by using the inhouse software. More simulation details for the CFC calculation could be found in our previous work.12,13,21 For CO2, simulations were performed at 1.25−35 bar by using a cubic solvent box; the box contains 10 PDMS molecules and has a length of about 28 Å. For H2, simulations were performed at 10−400 bar, and the cubic solvent box contains 10−100 PDMS molecules. The starting configuration for solvent molecules was obtained from the NPT MD simulation. In order to speed up calculations,12,13 the FGSF method was also used to compute the electrostatic interaction in the CFC simulation. 2.5. Surface Tension Calculations. There are mainly two types of methods to compute surface tension from molecular modeling, that is, the mechanical route22−25 and the test-area method.26 The mechanical route involves the pressure tensor calculation; the method is widely used even though it is difficult to implement due to pressure tensor calculations. The test-area method is based on the thermodynamic definition for surface tension; it is easy to implement since only potential energy is involved in calculations. The mechanical and test-area methods have been shown to give identical surface tension values.27−29 In this work, we used the test-area method to compute the surface tension, and the method is described below. For an NVT system involving two coexisting phases, such as the vapor−liquid equilibrium, the reversible work (Wrev) is equal to the Helmholtz free energy (F) change. This leads to Wrev = γ dA = dF, where γ is the surface tension, A is the interface area, and F is the total Helmholtz free energy for the two coexisting phases and the vapor−liquid interface. F is related to the partition function (Q) by F = −kbT ln Q, where kb is the Boltzmann constant and T is the temperature. The surface tension is calculated by using the following equation,

= lim − ΔA → 0

k bT Q (N , V , T , A + ΔA) ln Q (N , V , T , A ) ΔA

⎡ ⎤ ∫ exp⎣⎢ −U(Akb+T ΔA) ⎦⎥ d r ⃗ k bT = lim − ln ⎡ ⎤ ΔA → 0 ΔA ∫ exp⎣⎢ −kU(TA) ⎦⎥ d r ⃗ b

= lim − ΔA → 0

⎡ ΔU ⎤ k bT ln exp⎢ − ⎥ ΔA ⎣ k bT ⎦

A

(2)

where U is the total potential energy for the two coexisting phases (vapor + liquid) and the interface and r ⃗ denotes the xyz coordinates for all atoms in the system. Note that the average in eq 2 is based on the reference system corresponding to an interface area of A. The surface tension calculation by using the test-area method is similar to the excess chemical potential computation by using the test-particle insertion method.17,30 Furthermore, a forward (γsim,f) and backward (γsim,b) surface tension can be calculated as γsim, f =

lim



ΔA > 0, ΔA → 0

⎡ −ΔU ⎤ k bT ln exp⎢ ⎥ ΔA ⎣ k bT ⎦

A

and γsim,b =

lim

ΔA < 0, ΔA → 0



⎡ −ΔU ⎤ k bT ln exp⎢ ⎥ ΔA ⎣ k bT ⎦

A

respectively, depending on the sign of ΔA. In our work, it was found that γsim,f and γsim,b agree with each other when a small enough ΔA value was used. Similar to the Henry’s law constant calculation,17 two steps were used to calculate the surface tension. The first step is to generate configurations from MD simulation, following which a second step is performed to analyze the configuration to obtain the surface tension. For neat PDMS solvent, about 30−40 ns NPT MD simulations were first performed to obtain the PDMS liquid density and equilibrated structures. Starting from the equilibrated liquid structure and putting all molecules in image cells back to the primary cell by using the periodic boundary condition, NVT MD simulations were then performed to obtain the vapor−liquid equilibrium. Note that the simulation box sizes in the x (Lx) and y (Ly) directions obtained from the NPT MD simulation remained in the NVT MD simulation; however, the simulation box size in the z direction (Lz) was extended to a larger value in the NVT MD simulation to include two gas phases on each side of the central liquid slab. This NVT MD simulation box corresponds to two liquid−vapor interfaces and leads to ΔA = 2ΔAsim, where Asim = Lx × Ly. For neat PDMS, 500 PDMS molecules were used in the NVT MD simulation. When the UA PDMS model was used, the number of united atoms was 34 500, and the simulation box size was 124.403 Å × 124.403 Å × 165 Å. About 40 000 snapshots corresponding to the last 20 ns NVT MD production run were saved for further analysis to compute the C

DOI: 10.1021/acs.jpcc.5b05806 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

= 12.0 Å). More details to calculate γtail are given in the Supporting Information. The total surface tension, γsim,total = (γsim,f + γsim,b)/2 + γtail, was obtained by using the in-house software, which has been validated by comparing our simulated surface tension values for water and an LJ fluid with the literature simulation data (Supporting Information).

surface tension (described below). When the AA PDMS model was used, the simulation box contained 88 500 atoms, and the simulation box size was 124.519 Å × 124.519 Å × 165 Å in the NVT MD simulation. About 36 000−52 000 snapshots corresponding to the last 18−26 ns of NVT MD production runs were saved. For CO2−PDMS and H2S−PDMS gas− solvent mixtures, a similar procedure was used. In the NVT MD simulation, the box contains 500 PDMS molecules and some CO2 or H2S molecules in the liquid slab. Some CO2 or H2S molecules were also initially put in the gas phase. The surface tension for neat H2S was also calculated. In the NVT MD simulation, there are 5000 H2S molecules contained in the box (72.616 Å × 72.616 Å × 150.0 Å). All MD simulations were performed at 298 K by using the NAMD program.18 The time step was typically set to be 0.5 fs. The electrostatic interaction was calculated by using the particle mesh Ewald method.31 For the above saved snapshots, a procedure similar to the work of Vega et al.29 was used to compute the surface tension and is shown in Figure 1. Before step 1, the center of mass for

3. RESULTS AND DISCUSSION 3.1. PDMS Density. Since the solvent molar volume and free molar volume are important factors in determining gas solubility, the PDMS density was calculated at 298 K and 1 bar from NPT MD simulations. Note that the experimental PDMS density was used to tune σ parameters (Supporting Information) for the PDMS molecule. It was found that the simulated PDMS density values are very close to each other between different sets of simulations by starting from different initial configurations, using different simulation box sizes and different numbers of PDMS molecules. The AA PDMS model gives a density of 0.948 g/cm3 (Table 1), only 1.4% larger than Table 1. Simulated Density (ρ), Available Free Molar Volume Fraction (f v), Interaction Energy (E), and Surface Tension (γ) for Neat PDMS Solvent by Using Both the AllAtom (AA) and United-Atom (UA) PDMS Models at 298 K and 1 bara model

ρ (g/cm3)

fv

E (kJ/mol)

γ (10−3 N/m)

AA UA exp

0.948 (2) 0.941 (1) 0.935

0.118 (8) 0.143 (4) −

−159.7 (3) −299.8 (3) −

20.7−21.5 (10) 86 (2) 19

a

The interaction energy indicates the interaction between one PDMS molecule and all other PDMS molecules. Uncertainties in the last digit are given in parentheses. The experimental (exp) PDMS density32 and surface tension36,37 are also shown for comparison. The symbol “−” indicates that the experimental data are not available.

the experimental value of 0.935 g/cm3.32 The UA PDMS model gives a density of 0.941 g/cm3, also very close to the experimental data. Note that the PDMS solvent density was the only experimental input into computational simulations to tune the classical FF parameters. All other properties, such as gas solubility and surface tension, were predicted by using the tuned FF parameters. 3.2. CO2, H2, H2S, H2O and PDMS Interactions. 3.2.1. CO2, H2, H2S, and H2O Gas Interactions with PDMS-2 and PDMS. The optimized gas−PDMS-2 dimer structures obtained from AI calculations are shown in Figure 2. H2 is far away from the PDMS-2 molecule, suggesting that H2 interacts very weakly with PDMS-2. In contrast, H2O, H2S, and CO2 interact more strongly with PDMS-2 than does H2. The gas− PDMS-2 interaction energies obtained from the classical FF calculation by using the AA PDMS model were found to be close to the AI calculations (Table 2), with a difference of 0− 2.1 kJ/mol. Both AI and classical FF calculations by using the AA PDMS model indicate that gas interaction with the PDMS2 molecule decreases in the following order: H2O > H2S > CO2 > H2. Water interacts most strongly with PDMS-2 through the electrostatic (ELEC) interaction between one hydrogen atom of water and two oxygen atoms of the PDMS-2 molecule (Figure 2c). For CO2, H2, and H2S gases, the van der Waals (VDW) gas−PDMS-2 interactions predominate over the corresponding ELEC interactions.

Figure 1. Flowchart for the test-area method used to calculate surface tension by postanalyzing M snapshots obtained from the NVT molecular dynamics (MD) simulation. The NVT MD simulation involves one liquid slab in the center of a simulation box and two gas phases on each side of the liquid slab to establish liquid−gas equilibrium.

each molecule in an image cell is mapped back into the primary cell. Volumes in steps 1−3 remain the same because V = Lx × Ly × Lz = L′x × L′y × L′z; constant volume is required in the surface tension calculation of eq 2. Note that the center of mass for each molecule in step 1 is simply scaled in steps 2−3 according to the new box sizes; the internal coordinate for the molecule does not change. Consequently, the intrapotential energy values (bond, angle, dihedral, and improper) are the same between steps 1−3; only the nonbonded interaction (LJ and electrostatic interaction), Uinter, changes between steps 1−3 and contributes to the surface tension calculation. Both standard Ewald and FGSF methods have been used to calculate the electrostatic part of Uinter, and these two methods were found to give identical surface tension values. Except for γsim,f and γsim,b (Figure 1), a tail correction (γtail) to the surface tension must be added to account for the truncation of LJ interactions. Note that when calculating Uinter (Figure 1), a switch LJ interaction potential was used (ron = 10.5 Å and roff D

DOI: 10.1021/acs.jpcc.5b05806 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

predicted by the AA PDMS model (Table 2). Note that the gas interaction with the full PDMS molecule (not the simple PDMS-2 molecule) was also calculated by using both the AA and UA PDMS models. Similar to the PDMS-2 case, the UA PDMS model gives stronger CO2, H2, and H2S interactions with the PDMS molecule than with the AA PDMS model. For example, the UA PDMS model gives a 4.7 kJ/mol stronger CO2−PDMS interaction than does the AA PDMS model (Figure 3). For the H2O−PDMS interaction, the AA PDMS

Figure 2. Optimized gas−PDMS-2 dimer structures obtained from quantum ab initio gas-phase calculations at the B3LYP/6-311+ +G(d,p) level of theory. The optimized dimer structures are for (a) CO2−PDMS-2, (b) H2−PDMS-2, (c) H2O−PDMS-2, and (d) H2S− PDMS-2. Also shown are the atom−atom distances between certain atom pairs.

Table 2. Interaction Energy (INT) between CO2, H2, H2O, and H2S Gases and the PDMS-2 Molecule Obtained from Both Quantum Ab Initio (AI) and Classical Force Field (FF) Calculationsa AA INT (kJ/mol)

Figure 3. Optimized CO2−PDMS dimer structures obtained from classical force field calculations by performing 100 000 steps of dimer energy minimization, in which the PDMS structure was fixed and only the CO2 molecule was relaxed. Full van der Waals and electrostatic interactions between CO2 and PDMS were calculated without truncation. Both all-atom (AA) (a and c) and united-atom (UA) (b and d) PDMS models were used. Also shown are the interaction energies. The PDMS structures in a−d are the same except that H atoms in the AA PDMS model are removed in the UA model. CO2 positions are similar to each other between (a) and (b), (c) and (d) but different between (a) and (c), (b) and (d).

UA INT (kJ/mol)

system

AI INT (kJ/mol)

total

VDW

ELEC

total

VDW

ELEC

CO2 H2 H2O H2S

−9.4 −0.1 −23.1 −15.2

−11.5 −0.1 −22.8 −13.6

−8.5 −0.1 −5.7 −10.3

−3.0 −0.0 −17.1 −3.3

−15.6 −2.9 −18.2 −18.4

−13.5 −2.9 −6.4 −14.7

−2.1 −0.0 −11.8 −3.7

a

In the AI calculation, the gas-PDMS-2 dimer optimized structure (Figure 2) was obtained at the B3LYP/6-311++G(d,p) level of theory, following which a single-point energy calculation was performed to obtain the INT by using the MP2/cc-pVTZ theory. In the classical FF calculation, both all-atom (AA) and united-atom (UA) PDMS models were used to obtain the van der Waals (VDW) and electrostatic (ELEC) interactions. Full VDW and ELEC interactions were calculated without truncation in the classical FF calculation. The optimized dimer structure in the classical FF calculation was obtained from 100 steps of energy minimization by using the AI-optimized dimer structure as the starting configuration. Note that the internal coordinates for nonhydrogen atoms of the PDMS-2 molecule are kept the same between the AA and UA PDMS models.

model gives a 3 kJ/mol stronger energy for the ELEC interaction but a 4 kJ/mol weaker energy for the VDW interaction compared to the UA PDMS model; as a result, the UA and AA PDMS models give similar H2O−PDMS interaction energies. Note that in the AA PDMS model the atomic charge on the O atom of the PDMS molecule is 0.162e more negative compared to the values from the UA PDMS model. Additionally, the σ value of the LJ potential for the O atom in the AA PDMS model is about 0.24 Å smaller than in the UA PDMS model. These two differences between the AA and UA PDMS models lead to the H atom of the H2O molecule being brought about 0.4 Å closer to the O atom of PDMS in the H2O−PDMS dimer-optimized structure when the AA PDMS model was used; this in turn results in a stronger ELEC H2O−PDMS interaction but a weaker VDW H2O− PDMS interaction in the AA PDMS model. 3.2.2. PDMS−PDMS Interaction. The PDMS−PDMS interaction was also calculated. The interaction energy for one PDMS molecule with all other PDMS molecules in the system is defined as (Esum−Eexclude)/N, where N is the number of PDMS molecules and was set to be 100 and Esum indicates the total nonbonded interaction (VDW plus ELEC) for the system; it includes both the nonbonded interaction between different molecules and the nonbonded interaction between different atoms of the same molecule. The Esum was directly obtained from NPT MD simulation at 298 K and 1 bar. The nonbonded interaction between different atoms of the same

Similarly, the UA PDMS model predicts that H2O and H2S interact more strongly with PDMS-2 than do CO2 and H2. The UA PDMS model also predicts that the H2O−PDMS-2 ELEC interaction predominates over the VDW interaction, while for H2S, CO2, and H2, the VDW predominates over the ELEC interaction. The gas−PDMS-2 interaction energy difference between AI and the UA PDMS model is 2.8−6.2 kJ/mol, much larger than the difference between AI and the AA PDMS model. This indicates that the UA PDMS model is likely to be less accurate than the AA PDMS model to predict gas−PDMS interaction. As a final comparison, the UA PDMS model overestimates CO2, H2, and H2S gas interactions with PDMS-2 by 2.8−4.8 kJ/mol compared to the AA PDMS model (Table 2); the UA model gives a stronger VDW interaction than the AA PDMS model. In contrast, for the H2O−PDMS-2 interaction, the AA PDMS model gives a 4.6 kJ/mol stronger interaction than the UA model, mainly due to the strong ELEC interaction E

DOI: 10.1021/acs.jpcc.5b05806 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C molecule, Eexclude, accounts for the atom−atom pair interaction when the atom−atom pair is separated by three or more consecutive bonds. Eexclude was obtained by using the in-house software. To be consistent with the NPT MD simulation, the same VDW switch potential (ron = 10.5 Å and roff = 12.0 Å) was used to calculate the VDW component of Eexclude. For the ELEC part of Eexclude, the full ELEC interaction was calculated without truncation. The UA PDMS model predicts the PDMS−PDMS interaction to be −299.8 kJ/mol, much stronger than the interaction energy of −159.7 kJ/mol given by the AA PDMS model (Table 1). Additionally, it was found that the ELEC part of the PDMS−PDMS interaction is negligibly small compared to the VDW part for both the AA and UA PDMS models. In summary, the UA PDMS model gives both stronger CO2, H2, and H2S gas interactions with PDMS and also a stronger PDMS−PDMS interaction compared to the AA PDMS model. These interaction energy differences between the UA and AA PDMS models lead to different gas solubility and solvent surface tension values (below). 3.3. Henry’s Law Constants for CO2, H2, H2O, and H2S at 298 K. Since the AA PDMS model gives a more accurate gas−PDMS interaction than the UA PDMS model, most of the calculations have been done by using the AA PDMS model. The Henry’s law constants for CO2, H2, H2O, and H2S gases in PDMS solvent were calculated by using the AA PDMS model (Table 3). For CO2 and H2S, the Henry’s law constants among

(Table 2) except for H2O. Both AI and AA PDMS models indicate that H2O interacts more strongly with PDMS-2 (Table 2 and Figure 2) and PDMS (not shown here) than with H2S; however, the H2O solubility in PDMS is smaller. This is probably due to many −CH3 groups of the PDMS molecule, which may form hydrophobic pockets in the solvent phase. The hydrophobic pockets are unfavorable for H2O electrostatic interaction but favorable for H2S VDW interaction. Note that for the H2O−PDMS interaction ELEC predominates while VDW predominates for the H2S−PDMS interaction. Additionally, H2S gives a 4.6 kJ/mol stronger VDW interaction with PDMS-2 than does H2O (Table 2). As a result, the hydrophobic pockets lead to a higher H2S solubility compared to that of H2O. When comparing with experiment, the simulated solubility selectivity for CO2 over H2 is reasonably close to the experimental data. Enick et al.3 have measured CO2 and H2 solubilities in a similar PDMS 550 solvent (the average molecular weight for PDMS is 550 g/mol). The CO2 over H2 solubility selectivity at 298 K was experimentally determined to be 18 (based on the molar fraction), which is very close to our simulated value of 19 for PDMS solvent with a molecular weight of 1274.75 g/mol. Note that the PDMS molecular weight may affect the CO2 and H2 solubilities, and a more rigorous solubility and solubility selectivity comparison between simulation and experiment may be obtained by using the same or a similar PDMS molecule. By using the AA PDMS model, the H2O Henry’s law constant in PDMS at 298 K was estimated to be 8−12 bar (Table 3), which is fairly large compared to the H2O Henry’s law constants in other types of solvents. For example, the H2O Henry’s law constant at 298 K in the so-called hydrophobic [hmim][Tf2N] ionic liquid is only 0.04 ± 0.01 bar.12 The very large H2O Henry’s law constant in PDMS implies that the PDMS solvent is very hydrophobic (the larger the Henry’s law constant, the smaller the gas solubility). At room temperature and a steam partial pressure corresponding to 50% relative humidity in the gas phase, the equilibrium H2O concentration in PDMS was estimated to be 0.15% (mole fraction) and 20 ppm (weight fraction) based on the simulated H2O Henry’s law constant; the H2O solubility in PDMS is negligibly small. When the UA PDMS model was used, an even larger H2O Henry’s law constant of 17.5 ± 0.7 bar was obtained. Both the AA and UA PDMS models predict that the PDMS solvent is very hydrophobic. 3.4. CO2 and H2 Absorption Isotherms and Gas− Solvent Mixture Volumes. 3.4.1. CO2 Absorption Isotherm and Mixture Volume at 298 K. 3.4.1.1. CO2 Absorption Isotherm. CO2 solubility in PDMS was calculated by using both the AA and UA PDMS models (Figure 4 and Table 4). The simulated CO2 solubilities by using the AA PDMS model are comparable to the experimental data,3 with an absolute average difference of 14% at 15−35 bar. Although the AA PDMS model gives a 2.1 kJ/mol (Table 2) stronger CO2− PDMS-2 interaction than does the AI calculation, the AA PDMS model typically gives a smaller CO2 solubility compared to the experimental data, such as at 20 bar. One of the possible reasons for this solubility difference is that the AA PDMS model may give a smaller solvent available free volume compared to experiment. When the UA PDMS model is used, the CO2 absorption isotherm is significantly different than the experiment (Figure 4). The experiment data show that CO2 solubility increases significantly when the pressures is increased to 15−35 bar. In

Table 3. Henry’s Law Constants for CO2, H2S, H2O, and H2 Absorption in PDMS Solvent at 298 K by Using the AllAtom PDMS Modela Henry’s law constant (bar) system

RUN1

RUN 2

RUN 3

CO2 H2S H2O H2

9.0 (5) 1.06 (6) 8.4 (4) 155 (8)

7.0 (4) 0.98 (6) 11.1 (5)

8.2 (6) 1.12 (7) 12.2 (5)

a

The Henry’s law constants for CO2, H2S, and H2O were obtained from the test-particle insertion method by using three independent sets of NPT MD simulation trajectories for neat PDMS, which are indicated as RUN 1, RUN 2, and RUN 3, respectively. The Henry’s law constant for H2 was estimated from the linear regression of H2 pressure versus H2 solubility (mole fraction) below 20 bar (Table 4), not from the test-particle insertion method by using multiple NPT runs. Consequently, only one Henry’s law constant is given for H2. Uncertainties in the last digit are given in parentheses.

different runs are close to each other. For H2O, there is a large (45%) Henry’s law constant difference between RUN 1 and RUN 3, which suggests that longer NPT MD simulations are required to get a more accurate Henry’s law constant for H2O. Note that for CO2 the Henry’s law constant average for three different runs was estimated to be 8.1 ± 0.6 bar, comparable to the Henry’s law constant value of 10.8 bar, which was estimated on the basis of the simulated CO2 solubility at 1.25 bar obtained from CFC MC calculations (Table 4). The Henry’s law constant (obtained from the average of three trajectories in Table 3) suggests that gas solubility decreases in the following order: H2S (147) > CO2 (19) ≈ H2O (15) > H2 (1); the numbers in parentheses indicate gas solubility selectivity over H2. The gas solubility typically decreases in the same order as the gas−PDMS interaction F

DOI: 10.1021/acs.jpcc.5b05806 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Table 4. Simulated Gas Solubility (Mole Fraction xmolar and Mass Fraction xmass) and Mixture Molar Volume (Vmix) Obtained from the Continuous Fractional Component Monte Carlo Methoda gas

model

P (bar)

f (bar)

xmolar

xmass

Vmix (cm3/mol)

CO2 CO2 CO2 CO2 CO2 CO2 CO2 CO2 CO2 CO2 CO2 CO2 H2 H2 H2 H2 H2 H2 H2 H2 H2 H2 H2 H2

UA UA UA UA UA AA AA AA AA AA AA UA UA UA UA UA UA UA AA AA AA AA AA AA

7.004 14.0 20.0 30.0 35.0 1.25 2.5 7.004 14.0 20.0 30.0 35.0 10.0 20.0 40.0 100.0 200.0 400.0 10.0 20.0 40.0 100.0 200.0 400.0

6.7385 12.9486 17.87 25.26 28.581 1.2415 2.466 6.7385 12.9486 17.87 25.26 28.581 10.03 20.12 40.48 103.1 214.2 468.4 10.03 20.12 40.48 103.1 214.2 468.4

0.601 (4) 0.659 (5) 0.687 (4) 0.708 (2) 0.725 (2) 0.116 (8) 0.200 (14) 0.364 (5) 0.525 (7) 0.595 (4) 0.774 (4) 0.834 (4) 0.0714 (5) 0.1321 (5) 0.2250 (8) 0.455 (3) 0.592 (2) 0.715 (1) 0.070 (2) 0.129 (1) 0.225 (1) 0.400 (3) 0.538 (5) 0.680 (4)

0.0502 (8) 0.063 (1) 0.071 (1) 0.0780 (6) 0.0842 (8) 0.0048 (4) 0.0092 (7) 0.0204 (4) 0.038 (1) 0.0496 (8) 0.108 (2) 0.150 (3) 0.0001217 (9) 0.000241 (1) 0.000458 (2) 0.001366 (17) 0.002343 (18) 0.003988 (25) 0.000120 (3) 0.000235 (2) 0.000456 (3) 0.001107 (13) 0.001904 (36) 0.003422 (67)

547 (5) 471 (6) 435 (6) 408 (2) 386 (3) 1134 (11) 1040 (20) 832 (6) 637 (8) 549 (5) 338 (5) 262 (5) 1259.5 (6) 1177.7 (7) 1053 (1) 732 (4) 553 (3) 389 (2) 1255 (2) 1173 (1) 1053 (1) 795 (3) 608 (5) 423 (5)

a

Both all-atom (AA) and united-atom (UA) PDMS models were used. The numbers in parentheses are uncertainties in the last digit obtained from standard block average calculations.

It is instructive to compare CO2 solubilities obtained from the AA and UA PDMS models. Both the solvent free volume and gas−solvent interaction are two important factors in determining the gas solubility. The solvent free molar volume fraction for neat PDMS was estimated by using #U/(#O + #U), where #U and #O denote the numbers of unoccupied and occupied cubelets,17,19 respectively, which were obtained via the Henry’s law constant calculation by using the test-particle insertion method. The AA PDMS model gives a free molar volume fraction of 0.118 (Table 1) and a free molar volume of 159 ± 11 cm3/mol at 298 K whereas the UA PDMS model gives a solvent free molar volume of 194 ± 5 cm3/mol at 298 K. The UA PDMS model gives both a larger (22%) solvent free volume and a stronger (36%) (Table 2) CO2−PDMS-2 interaction compared to the AA PDMS model. As a result, the UA PDMS model gives higher CO2 solubilities than the AA PDMS model at lower CO2 pressures below 20 bar (Figure 4). In contrast, at higher CO2 pressures above 30 bar, the UA PDMS model gives a significantly smaller CO2 solubility than the AA PDMS model. This is partially due to the fact that the UA PDMS model gives a stronger PDMS−PDMS interaction compared to the AA PDMS model (Table 1), which leads to a small solvent volume expansion upon CO2 absorption (Figure 5). At 35 bar, the UA PDMS model predicts that the CO2− PDMS mixture volume increases only by 4.5% compared to the neat solvent volume; this small volume expansion leads to limited solvent free volume for additional CO2 absorption. As a result, the CO2 solubility approaches a plateau region and increases only a little at 35 bar. In contrast, the AA PDMS model gives a much larger solvent volume expansion upon CO2 absorption (Figure 5); at 35 bar, the CO2−PDMS mixture volume expansion was estimated to be 19%. This large solvent

Figure 4. Simulated CO2 absorption isotherms in PDMS solvent at 298 K by using both the all-atom (AA) (triangles) and united-atom (UA) PDMS (squares) models. Simulation error bars are smaller than the symbols and are not shown. For comparison, the experimental CO2 solubility at 298 K in a similar PDMS solvent (PDMS molecular weigh of 1250 g/mol) is also shown.3 Lines are added to guide the eyes. The (0, 0) point is also added to the experimental (solid line) and simulation data (dashed line).

contrast, the UA PDMS model predicts that the CO2 solubility increases only slightly at elevated pressures, such as above 20 bar. At lower pressures below 15 bar, the UA PDMS model overestimates the CO2 solubility by 26% compared to the experiment, which is partially due to the fact that the UA PDMS model gives a 6.2 kJ/mol stronger CO2−PDMS-2 interaction compared to the AI calculation (Table 2). On the other hand, at high CO2 pressures such as 35 bar, the UA PDMS model significantly underestimates the CO2 solubility by 44% compared to the experiment, which is partially due to the small solvent volume expansion upon CO 2 absorption predicted by the UA PDMS model (below). G

DOI: 10.1021/acs.jpcc.5b05806 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Figure 7. Simulated CO2−PDMS mixture molar volume versus CO2 solubility at 298 K, obtained from the continuous fraction component Monte Carlo simulations by using both the all-atom (AA) (triangles) and united-atom (UA) PDMS (squares) models. Lines indicate linear regressions for the simulation data.

Figure 5. Simulated CO2−PDMS mixture volume expansion relative to the neat PDMS solvent volume at 298 K by using both the all-atom (AA) (triangles) and united-atom (UA) (squares) PDMS models. The error bars are typically smaller than the symbols and are not shown. Lines are added to guide the eyes.

⎛ ∂Ṽ ⎞ ̃ + (1 − x 2) × ⎜ mix ⎟ V2 = Vmix ⎝ ∂x 2 ⎠T , P

volume expansion allows for more CO2 absorption and leads to large CO2 solubility at high CO2 pressure. The large solvent volume expansion at high CO2 pressure predicted by the AA PDMS model has also been experimentally observed for CO2 absorption in other organic solvents. For example, at 298 K and a CO2 pressure of 28 bar, the solvent volume expansion values are 18% for toluene and 26% for cyclohexanone.33 3.4.1.2. g(r) for the CO2−PDMS Mixture and CO2 Partial Molar Volume. The radial distribution function (g(r)) values for PDMS in both neat PDMS solvent and the CO2−PDMS mixture were analyzed by using the AA PDMS model (Figure 6). For neat PDMS, the g(r) data are similar to each other

where V̅ 2 is the solute partial molar volume, x2 is the solute molar fraction, and Ṽ mix is the liquid mixture (solute + solvent) molar volume. By assuming a linear relationship between Ṽ mix and x2 (Figure 7), the solute partial molar volume can be derived as the summation of the y intercept and the slope for the linear regression. By using this method, the CO2 partial molar volume in PDMS was estimated to be 59 ± 7 cm3/mol when the AA PDMS model was used. Interestingly, the UA and AA PDMS models give similar linear regressions although there are only limited data for the UA PDMS model. 3.4.2. H2 Absorption Isotherm and Mixture Volume at 298 K. H2 solubilities in PDMS by using both the UA and AA PDMS models are shown in Figure 8 and Table 4. The UA

Figure 6. Radial distribution function (g(r)) for the center of mass of the PDMS molecule at 298.2 K by using the all-atom PDMS model. The blue and black lines indicate g(r) data for neat PDMS at 1 and 35 bar, respectively. The red line is for the CO2−PDMS mixture at a CO2 pressure of 35 bar and a CO2 mass fraction of 0.162.

Figure 8. Simulated H2 solubilities in PDMS at 298 K by using both the all-atom (AA) (triangles) and united-atom (UA) PDMS (squares) models. The experimental (circles) H2 solubility data3 at 298 K in a PDMS 550 solvent (molecular weight of 550 g/mol) are also shown for comparison. The inset shows H2 solubilities at low pressures. Lines are added to guide the eyes.

between 1 and 35 bar, which implies that pressure will not significantly affect the solvent volume; simulations show that when the pressure is increased from 1 to 35 bar (without CO2 absorption), the neat PDMS solvent molar volume decreases by less than 1%. In contrast, CO2 absorption can significantly increase the solvent volume by increasing the distance between solvent molecules. For example, for neat PDMS, the first peak of g(r) greater than 1 occurs at 12.5 Å (Figure 6). Upon CO2 absorption at 35 bar, the first peak of g(r) for the PDMS molecule in the mixture is increased to a much larger value of 19 Å. Finally, the CO2−PDMS mixture molar volume versus CO2 mole fraction was computed (Figure 7 and Table 4). From thermodynamics, it can easily be shown that the solute (gas) partial molar volume in liquid can be calculated as

PDMS model gives a larger H2 solubility than does the AA PDMS model at 100−400 bar. This is also due to the fact that the UA PDMS model gives both a larger solvent free molar volume (Table 1) and a stronger H2−PDMS interaction (Table 2) than does the AA PDMS model. At 10−40 bar, the AA and UA PDMS models give similar H2 solubilities. More solvent molecules are needed to amplify the H2 solubility difference between the AA and UA PDMS models at low H2 pressures. Below 40 bar, the simulated H2 solubilities in PDMS (molecular weight of 1274.72 g/mol) are comparable to the H

DOI: 10.1021/acs.jpcc.5b05806 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Table 5. Simulated Liquid−Vapor Surface Tensions (× 10−3 N/m) by Using the Test-Area Method for Neat H2S and PDMS at 298 K, along with Surface Tension Values for CO2−PDMS and H2S−PDMS Mixtures at 298 K and Different Gas Concentrations (Mass Fraction xgas,mass)a system H2S H2S H2S H2S PDMS r1 PDMS r1 PDMS r1 PDMS r2 PDMS PDMS CO2−PDMS CO2−PDMS CO2−PDMS CO2−PDMS CO2−PDMS CO2−PDMS CO2−PDMS H2S−PDMS H2S−PDMS H2S−PDMS

PDMS model

xgas,mass

ζ

ELEC

γsim,f

γsim,b

γsim,ave

γtail

γsim,total

γexp

AA AA AA AA UA UA AA AA AA AA AA UA UA AA AA AA

1.0 1.0 1.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.04 0.04 0.04 0.018 0.008 0.065 0.065 0.005 0.005 0.005

0.0002 0.0002 0.0005 0.0005 0.0001 0.0001 0.0005 0.0005 0.0005 0.0005 0.0001 0.0001 0.0005 0.0005 0.0005 0.0005 0.0005 0.0001 0.0001 0.0005

Ewald FGSF Ewald FGSF Ewald FGSF FGSF FGSF Ewald FGSF Ewald FGSF FGSF FGSF FGSF Ewald FGSF Ewald FGSF FGSF

9.2 (1) 9.0 (1) 8.6 (1) 9.0 (1) 14.3 (9) 13.8 (9) 14.2 (10) 15.6 (10) 67 (2) 67 (2) 11.2 (17) 11.0 (16) 11.6 (17) 12.2 (8) 11.1 (11) 33 (1) 32 (1) 16.8 (13) 16.4 (13) 16.2 (14)

9.2 (1) 9.0 (1) 9.2 (2) 9.0 (2) 14.1 (8) 13.6 (8) 13.7 (7) 15.1 (12) 64 (2) 64 (2) 11.4 (16) 11.1 (16) 12.4 (14) 11.9 (7) 10.7 (10) 30 (1) 30 (1) 16.8 (14) 16.4 (13) 16.3 (14)

9.2 (1) 9.0 (1) 8.9 (2) 9.0 (2) 14.2 (9) 13.7 (9) 14.0 (10) 15.4 (12) 66 (2) 66 (2) 11.3 (17) 11.1 (16) 12.0 (17) 12.0 (8) 10.9 (11) 31 (1) 31 (1) 16.8 (14) 16.4 (13) 16.3 (14)

3.5 3.5 3.5 3.5 7.0 7.0 7.0 6.1 19.5 19.5 5.6 5.6 5.6 5.7 6.1 24.0 24.0 6.0 6.0 6.0

12.7 (1) 12.5 (1) 12.4 (2) 12.5 (2) 21.2 (9) 20.7 (9) 21.0 (10) 21.5 (12) 86 (2) 86 (2) 16.9 (17) 16.7 (16) 17.6 (17) 17.7 (8) 17.0 (11) 55 (1) 55 (1) 22.8 (14) 22.4 (13) 22.3 (14)

10.843 10.843 10.843 10.843 19.5 19.5 19.5 19.5 19.5 19.5 − − − − − − − − − −

The surface tension values for γsim,f, γsim,b, γsim,ave (γsim,ave = (γsim,f + γsim,b)/2), and γtail and the total surface tension (γsim,total = γsim,ave + γtail) are shown. Both all-atom (AA) and united-atom (UA) PDMS models and the standard Ewald (Ewald) and Fennel and Gezelter shift force (FGSF)20 methods used to compute the electrostatic (ELEC) interaction were used. For neat PDMS, two sets of NVT MD simulations (r1 and r2) were used to calculate the surface tension. Different ζ values were also used to investigate their effects on the surface tension calculation. Uncertainty in the last digit obtained from block average calculation is given in parentheses. For comparison, the experimental surface tension data (γexp) for neat H2S34 and PDMS36,37 are also included. Symbol “−” indicates that the corresponding experimental data are unavailable. a

calculated to be 12.4 × 10−3−12.7 × 10−3 N/m, about 14% different than the experimental data (10.843 × 10−3 N/m),34 which is due to the classical FF parameters for H2S. Note that when Potoff and co-workers14 developed the classical FF parameters for H2S, which were used in this work, the FF parameters reproduce the experimental vapor−liquid equilibrium density, saturation pressure, and heat of vaporization. Our simulations in this work further indicate that their H2S FF parameters also predict the H2S surface tension to be reasonably close to the experimental data. When the FGSF method is used, the γsim,f or γsim,b values agree with each other between ζ = 0.0002 and ζ = 0.0005. The γsim,f value also agrees with γsim,b for the same ζ. These results suggest that the ζ value (0.0002−0.0005) is small enough and the number of configurations (about 16 000 snap shots corresponding to 16 ns) is large enough to obtain accurate surface tension values. Actually, thermodynamic consistency requires that γsim,f and γsim,b must be equal to each other.28 When the Ewald method was used, similar findings to the FGSF method were observed. Additionally, the γsim,f and γsim,b values obtained from the Ewald method are typically in agreement with the corresponding values obtained from the faster FGSF12 method (Table 5). Finally, the coexistence bulk density for H2S at 298 K obtained from the NVT MD simulation was estimated to be 0.7452 ± 0.0002 g/cm3 for the liquid and 0.0447 ± 0.0001 g/ cm3 for the gas phase. The simulated liquid coexistence density is close to the experimental data (0.775 g/cm3),34 with a small difference of 3.8%. However, the simulated vapor coexistence density is much different (32%) than the experimental value (0.0339 g/cm3).34 Potoff et al.14 have also computed the H2S

experimental H2 solubility data3 in a similar PDMS solvent (molecular weight of 550 g/mol). Different from CO2 absorption, it was found that H2 absorption will only slightly expand the PDMS solvent volume at high H2 pressures. At 400 bar, both the AA and UA PDMS models predict that the H2−PDMS mixture volume is increased by only 0.5−1.4% compared to the neat PDMS solvent volume. H2 molecules are absorbed in the available free volume of neat PDMS solvent, and they do not expand the solvent volume as CO 2 molecules do. The partial molar volumes for H 2 absorption in PDMS were calculated to be −1 ± 7 and −22 ± 15 cm3/mol corresponding to the UA and AA PDMS models, respectively. 3.5. Surface Tension. 3.5.1. Neat H2S and PDMS at 298 K. 3.5.1.1. Neat H2S. In order to calculate the vapor−liquid interface surface tension for pure gas, the temperature needs to be lower than the gas critical temperature to establish vapor− liquid equilibrium. Note that the surface tension was calculated at 298 K in this work, which is far below the critical temperature value (373 K) for H2S34 but close to the critical temperature (304 K) for CO2.34 It is expected that H2S will exhibit an appreciable surface tension value at 298 K. In contrast, the CO2 surface tension at 298 K is very small, only about 0.5 × 10−3 N/m.34 Consequently, in this work, we report only the surface tension for H2S gas at 298 K. We have also not studied surface tensions for H2 and H2O and their absorption effects on PDMS surface tension; these two gases have low solubility values in PDMS (Table 3) at 298 K. Additionally, 298 K is far above the critical temperature (33 K) for H2.34 Surface tension values for neat H2S and PDMS are shown in Table 5. For neat H2S, the surface tension at 298 K was I

DOI: 10.1021/acs.jpcc.5b05806 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

upon CO2 absorption (Figure 5), which results in smaller CO2 solubility at high CO2 pressures (35 bar) compared to the AA PDMS model. Although more simulations or experimental data are needed, our simulation findings obtained from the UA and AA PDMS models tend to suggest that it may be challenging to develop a solvent which could exhibit both a significantly large surface tension (such as >80 × 10−3 N/m) and a large CO2 solubility at high CO2 pressures. Similar to the AA PDMS model, γsim,f agrees with γsim,b for both the FGSF and Ewald methods; the γsim,f and γsim,b values also agree with each other between the FGSF and Ewald methods. These consistencies indicate that both the ζ value (0.0005) and the length of the NVT MD simulation (40 000 snapshots for the last 20 ns production run) are appropriate for obtaining convergent surface tension values when the UA PDMS model is used. 3.5.2. CO2−PDMS and H2S−PDMS Mixtures at 298 K. In order to investigate the gas absorption effects on surface tension, CO2−PDMS and H2S−PDMS mixtures were studied, and the simulation results are shown below. 3.5.2.1. CO2−PDMS Mixture. Surface tension values for CO2−PDMS mixtures (CO2 mass fraction between 0.008− 0.065) at 298 K were calculated by using both the AA and UA PDMS models (Table 5). When the AA PDMS model was used, the CO2−PDMS mixture surface tension was calculated to be 16.7 × 10−3−17.7 × 10−3 N/m, about 3 × 10−3−5 × 10−3 N/m (20%) smaller compared to the simulated neat PDMS surface tension. In the case of using the UA PDMS model, the mixture (CO2 mass fraction of 0.065) surface tension was calculated to be 55 ± 1 × 10−3 N/m, about 30 × 10−3 N/m (35%) smaller than the simulated surface tension (86 × 10−3 N/m) for neat PDMS given by the UA PDMS model. Similar to the neat PDMS case, γsim,f agrees with γsim,b for CO2−PDMS mixtures. The γsim,f or γsim,b values also agree with each other between the FGSF and Ewald methods. In the case of using both the AA PDMS model and the Ewald method, a small ζ value (0.0001) is needed to give consistent values between γsim,f and γsim,b. In summary, both the AA and UA PDMS models show that CO2 absorption will decrease the solvent surface tension by 20−35% at a CO2 mass fraction of between 0.008 and 0.065. Similarly, it has also been experimentally determined that CO2 physical absorption in different types of solvents, such as crude oil,38 1-ethyl-3-methylimidazolium ethyl sulfate, and 1-butyl-3methylimidazolium hexafluorophosphate ionic liquids,39 will decrease the solvent surface tension. To further investigate the above surface tension decrease upon CO2 absorption, CO2 and PDMS densities were calculated by using the AA PDMS model (Figure 9). CO2 molecules are absorbed in three regions, that is, the liquid, interface, and gas regions. The interface region exhibits the largest CO2 density, followed by the liquid region; the gas phase exhibits the smallest CO2 density. Similarly, the PDMS molecules are also located in three regions. The gas phase exhibits a negligibly small PDMS density, which suggests that the PDMS solvent is nonvolatile. It is expected that CO2 molecules in the liquid and interface regions will decrease the PDMS−PDMS interaction, especially in the interface region because of significant CO2 concentration in this region. As a result, the decreased PDMS−PDMS interaction will lead to a smaller surface tension for the CO2−PDMS mixture compared to that for the neat PDMS solvent. The UA PDMS model also

coexistence density from the grand canonical MC simulation combined with the histogram-weighting technique. In their calculations, the LJ interaction was truncated at 10 Å and an analytical tail correction was applied. In our work, a switch LJ potential was used. Different methods to calculate the LJ tail interaction may affect the vapor coexistence density35 appreciably. Additionally, the size of the NVT MD simulation box in this work may affect the coexistence density due to the presence of an interface. 3.5.1.2. Neat PDMS by Using the AA PDMS Model. For neat PDMS at 298 K, the AA PDMS model gives the PDMS solvent surface tension to be 20.7 × 10−3−21.5 × 10−3 N/m (Table 5), which is 10% different than the experimental data (19.5 × 10−3 N/m36,37). By using a similar AA PDMS model, Ismail et al. have also obtained surface tension values (γsim,ave and γtail) for a similar PDMS molecule with 20 repeating [−Si(CH3)2O−] monomers,10 and their results are close to our calculations. When the FGSF method was used, the values for γsim,f and γsim,b agree with each other between ζ = 0.0001 and ζ = 0.0005; γsim,f also agrees with γsim,b for the same ζ. Additionally, two different sets of NVT MD simulations (r1 and r2 in Table 5) give surface tension values close to each other. These consistencies in surface tension calculations indicate that both ζ values and the number of configurations (about 52 000 snapshots for r1 and 36 000 snapshots for r2, which correspond to 26 and 18 ns production runs, respectively) are appropriate for obtaining reliable surface tension values. In contrast, when the Ewald method is used, γsim,f agrees with γsim,b only when a smaller ζ value (0.0001) is used (Table 5). At a larger ζ value such as 0.0005, γsim,f and γsim,b values were found not to agree with each other. The surface tension values were calculated to be γsim,f = 7.9 × 10−3 ± 0.9 × 10−3 N/m and γsim,b = 14.3 × 10−3 ± 0.8 × 10−3 N/m at ζ = 0.0005. It is interesting that γsim,f values do not agree with each other between ζ = 0.0001 and ζ = 0.0005; however, γsim,b values agree with each other. In summary, when the FGSF method was used to compute the electrostatic interaction in the surface tension calculation, γsim,f or γsim,b values agreed with each other between large (0.0005) and small (0.0001) ζ values; γsim,f also agreed with γsim,b. In contrast, when the Ewald method was used, a small ζ value, such as 0.0001, was needed to give consistent values between γsim,f and γsim,b. Additionally, the values for γsim,f or γsim,b agree with each other between the FGSF and Ewald methods when a small ζ = 0.0001 value is used. Considering that (1) the FGSF and Ewald methods give consistent surface tension values and (2) the FGSF method speeds up the electrostatic calculation compared to the Ewald method,12 we recommend using the FGSF method to compute the electrostatic interaction in the surface tension calculation. 3.5.1.3. Neat PDMS by Using the UA PDMS Model. When the UA PDMS model was used, the PDMS solvent surface tension at 298 K was calculated to be 86 × 10−3 ± 2 × 10−3 N/ m (Table 5), about 4 times larger than the experimental data. As discussed above, the UA PDMS model gives a much stronger PDMS−PDMS interaction (Table 1) compared to the AA PDMS model, which will affect the PDMS surface tension and CO2 solubility. On the one hand, the stronger PDMS− PDMS interaction predicted by the UA PDMS model leads to a 4 times larger surface tension compared to the AA PDMS model (Table 5). On the other hand, this stronger interaction for the UA model gives a smaller solvent volume expansion J

DOI: 10.1021/acs.jpcc.5b05806 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

the liquid and gas phases is much larger than the CO2 density ratio, which is consistent with the stronger H2S−PDMS interaction compared to the CO2−PDMS interaction.

4. CONCLUSIONS The AA PDMS model, as well as the AI calculation, predicts that the PDMS interaction with a gas decreases in the following order: H2O > H2S > CO2 > H2. For the H2O−PDMS interaction, the electrostatic interaction predominates. In contrast, the VDW interaction predominates for H2S, CO2, and H2 interactions with PDMS. Gas solubility follows the same order as gas interaction except for water. By using the AA PDMS model, simulations show that gas solubility (in the reverse order of the Henry’s law constant) decreases in the following order: H2S (1.05 ± 0.07 bar) > CO2 (8.1 ± 0.6 bar) ≈ H2O (10.6 ± 0.5 bar) > H2 (155 ± 8) bar. Although H2O interacts more strongly with PDMS than does H2S, its solubility in PDMS is much smaller than that of H2S. This is probably due to many −CH3 groups of the PDMS molecule, which may form hydrophobic pockets in the solvent phase. These hydrophobic regions are unfavorable for H2O electrostatic interaction but favorable for H2S VDW interactions. Consequently, PDMS is hydrophobic and exhibits a larger solubility for H2S than for H2O. Note that the CO2/H2 solubility selectivity in PDMS was estimated to be 19, reasonably close to the experimental selectivity data of 18 for CO2 over H2 in a similar PDMS solvent3 (molecular weight of 550 g/mol). When the UA and AA PDMS models are compared, the UA PDMS model gives a 22% larger solvent free volume and a 4.1 kJ/mol stronger CO2−PDMS interaction compared to the AA PDMS model. Consequently, the UA PDMS model gives a larger CO2 solubility than the AA PDMS model at CO2 pressures below 20 bar. At high CO2 pressures above 30 bar, the AA PDMS model gives a larger CO2 solubility than does the UA PDMS model, which is partially due to the large solvent volume expansion (20%) upon CO2 absorption predicted by the AA PDMS model. For neat PDMS at 298 K, when the AA PDMS model was used, the surface tension was calculated to be 20.7 × 10−3−21.5 × 10−3 N/m, 10% less than the experimental value. When the UA PDMS model was used, the surface tension was calculated to be 86 × 10−3 ± 2 × 10−3 N/m, about 4 times larger than the experimental data. The UA PDMS model gives a significantly larger surface tension (4 times) compared to that of the AA PDMS model, which is partially due to the stronger (1.9 times) PDMS−PDMS interaction predicted by the UA PDMS model. In summary, the AA PDMS model gives a larger CO2 solubility at high CO2 pressures above 30 bar but a smaller PDMS surface tension compared to that from the UA PDMS model. One may compromise between large (small) CO2 solubility and small (large) solvent surface tension. It would be challenging to develop a solvent which exhibits both a significantly large surface tension and a large CO2 solubility at high CO2 pressures. To obtain a significantly large surface tension, the solvent−solvent interaction needs to be strong, which will lead to small volume expansion upon CO 2 absorption at high CO2 pressure, and as a result a smaller CO2 solubility may be obtained. Finally, we have investigated CO2 and H2S absorption effects on PDMS surface tension. Both AA and UA PDMS models predict that CO2 absorption (mass fraction of 0.008−0.065) will decrease the solvent surface tension by 20−35%. CO2 molecules are absorbed in three regions, that is, the liquid, gas,

Figure 9. Local density distribution for CO2 and PDMS molecules in the CO2−PDMS mixture at 298 K obtained from the NVT molecular dynamics simulation by using the all-atom PDMS model. The red dashed line is for PDMS in the mixture at a CO2 mass fraction of 0.04 (xCO2 = 0.04); the red solid line is for CO2 in the mixture at xCO2 = 0.04. The blue and orange solid lines are for CO2 in the mixture at xCO2 = 0.018 and xCO2 = 0.008, respectively. For clarity, CO2 densities are multiplied by a factor of 10. Density distributions for PDMS in the mixture at xCO2 = 0.018 and xCO2 = 0.008 are similar to the dashed red line for xCO2 = 0.04 and are not shown. For comparison, the density distribution for neat PDMS (black solid line) is also shown. Two vertical dashed black arrow lines are drawn to roughly indicate CO2 absorption in three regions, that is, the liquid, interface, and gas regions at xCO2 = 0.04.

gives density distributions for the CO2−PDMS mixture that are similar to those given by the AA PDMS model. 3.5.2.2. H2S−PDMS Mixture. In contrast to the CO2 absorption effect on the PDMS surface tension, H2S absorption tends to slightly increase the solvent surface tension at small H2S concentration (Table 5) even though longer NVT MD simulations are needed to give a more accurate comparison. Note that the VDW well depth value for the S atom of H2S is about 3.5 and 10.3 times larger than the VDW well depth values for C and O atoms of CO2, respectively; this leads to a stronger H2S−PDMS interaction compared to the CO2− PDMS interaction (Table 2). Consequently, the H 2 S absorption in PDMS may slightly increase the PDMS solvent surface tension. Density distributions for the H2S−PDMS mixture (Figure 10) are similar to those for the CO2−PDMS mixture except that the H2S density peak in the interface region moves toward the liquid phase. Additionally, the H2S density ratio between

Figure 10. Local density distributions for the H2S−PDMS mixture (black lines) at 298 K obtained from NVT molecular dynamics simulation by using the all-atom PDMS model at xH2S = 0.005. For comparison, the density distributions for the CO2−PDMS mixture (red lines) at 298 K and xCO2 = 0.008 are also shown. For clarity, gas densities are multiplied by a factor of 20. Other symbols are similar to those in Figure 9. K

DOI: 10.1021/acs.jpcc.5b05806 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

of Energy’s Carbon Sequestration Program. Int. J. Greenhouse Gas Control 2008, 2, 9−20. (2) MacDowell, N.; Florin, N.; Buchard, A.; Hallett, J.; Galindo, A.; Jackson, G.; Adjiman, C. S.; Williams, C. K.; Shah, N.; Fennell, P. An Overview of CO2 Capture Technologies. Energy Environ. Sci. 2010, 3, 1645−1669. (3) Enick, R. M.; Koronaios, P.; Stevenson, C.; Warman, S.; Morsi, B.; Nulwala, H.; Luebke, D. Hydrophobic Polymeric Solvents for the Selective Absorption of CO2 from Warm Gas Streams that also Contain H2 and H2O. Energy Fuels 2013, 27, 6913−6920. (4) Thitakamol, B.; Veawab, A. Foaming Behavior in CO 2 Absorption Process Using Aqueous Solutions of Single and Blended Alkanolamines. Ind. Eng. Chem. Res. 2008, 47, 216−225. (5) Thitakamol, B.; Veawab, A. Foaming Model for CO2 Absorption Process Using Aqueous Monoethanolamine Solutions. Colloids Surf., A 2009, 349, 125−136. (6) Reynolds, A.; Verheyen, T.; Adeloju, S.; Meuleman, E.; Feron, P. Towards Commercial Scale Postcombustion Capture of CO2 with Monoethanolamine Solvent: Key Considerations for Solvent Management. Environ. Sci. Technol. 2012, 46, 3643−3654. (7) Onda, K.; Takeuchi, H.; Okumoto, Y. Mass Transfer Coefficients Between Gas and Liquid Phases in Packed Columns. J. Chem. Eng. Jpn. 1968, 1, 56−62. (8) Sok, R.; Berendsen, H.; van Gunsteren, W. Molecular-Dynamics Simulation of the Transport of Small Molecules Across a Polymer Membrane. J. Chem. Phys. 1992, 96, 4699−4704. (9) Smith, J.; Borodin, O.; Smith, G. A quantum chemistry based force field for poly(dimethylsiloxane). J. Phys. Chem. B 2004, 108, 20340−20350. (10) Ismail, A. E.; Grest, G. S.; Heine, D. R.; Stevens, M. J.; Tsige, M. Interfacial Structure and Dynamics of Siloxane Systems: PDMS-Vapor and PDMS-Water. Macromolecules 2009, 42, 3186−3194. (11) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon: Oxford, U.K., 1987. (12) Shi, W.; Maginn, E. J. Atomistic Simulation of the Absorption of Carbon Dioxide and Water in the Ionic Liquid 1-n-hexyl-3methylimidazolium bis(trifluoromethylsulfonyl) imide ([hmim][Tf2N]). J. Phys. Chem. B 2008, 112, 2045−2055. (13) Shi, W.; Sorescu, D. C.; Luebke, D. R.; Keller, M. J.; Wickramanayake, S. 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]). J. Phys. Chem. B 2010, 114, 6531−6541. (14) Kamath, G.; Lubna, N.; Potoff, J. Effect of Partial Charge Parametrization on the Fluid Phase Behavior of Hydrogen Sulfide. J. Chem. Phys. 2005, 123, 124505. (15) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A. et al. Gaussian 09, Revision A.1. Gaussian, Inc., Wallingford, CT, 2009. (16) Accelrys Materials Studio, Version v 5.5.0.0; Accelrys Inc.: San Diego, USA, 2010. (17) Shi, W.; Thompson, R. L.; Albenze, E.; Steckel, J. A.; Nulwala, H. B.; Luebke, D. R. Contribution of the Acetate Anion to CO2 Solubility in Ionic Liquids: Theoretical Method Development and Experimental Study. J. Phys. Chem. B 2014, 118, 7383−7394. (18) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kale, L.; Schulten, K. Scalable Molecular Dynamics with NAMD. J. Comput. Chem. 2005, 26, 1781− 1802. (19) Shah, J. K.; Maginn, E. J. Monte Carlo Simulations of Gas Solubility in the Ionic Liquid 1-n-butyl-3-methylimidazolium Hexafluorophosphate. J. Phys. Chem. B 2005, 109, 10395−10405. (20) Fennell, C. J.; Gezelter, J. D. Is the Ewald Summation Still Necessary? Pairwise Alternatives to the Accepted Standard for LongRange Electrostatics. J. Chem. Phys. 2006, 124, 234104.

and interface regions. The interface region exhibits the largest CO2 density, followed by the liquid region, and the gas phase exhibits the smallest CO2 density. The presence of CO2 in both the liquid and interface regions will decrease the PDMS− PDMS interaction, which will decrease the PDMS surface tension. For H2S−PDMS, H2S absorption slightly increases the PDMS surface tension, which is partially due to the stronger H2S−PDMS interaction compared to the CO2−PDMS interaction. Similar to CO2, H2S molecules are absorbed in three regions. Although H2S exhibits the highest density at the interface, which is followed by the liquid and gas phases, its density peak in the interface region moves away from the gas phase and is brought closer to the liquid region compared to CO2.



ASSOCIATED CONTENT

* Supporting Information S

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.5b05806. Classical force field parameters for the AA and UA PDMS models, validation of the FGSF method against the standard Ewald method to compute electrostatic interactions for PDMS, tail correction (γtail) formula for the surface tension for a complex PDMS molecule and gas−PDMS mixture, and validation of the in-house software to compute surface tensions. (PDF)



AUTHOR INFORMATION

Corresponding Author

*Phone: +1 412 386 4406. Fax: +1 412 386 5990. E-mail: [email protected]. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS We thank Bob Enick, Peter Koronaios, David Luebke, Hunaid Nulwala, Fan Shi, Janice Steckel, and David Hopkinson for helpful discussions. This technical effort was performed in support of the National Energy Technology Laboratory’s ongoing research in computational chemistry under RES contract DE-FE0004000. This project was funded by the Department of Energy, National Energy Technology Laboratory, an agency of the United States Government, through a support contract with URS Energy & Construction, Inc. Neither the United States Government nor any agency thereof nor any of their employees nor URS Energy & Construction, Inc. nor any of their employees makes any warranty, expressed or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed or represents that its use would not infringe on privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof. The views and opinions of the authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof.



REFERENCES

(1) Figueroa, J. D.; Fout, T.; Plasynski, S.; McIlvried, H.; Srivastava, R. D. Advancesn in CO2 Capture Technology - The US Department L

DOI: 10.1021/acs.jpcc.5b05806 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C (21) Shi, W.; Maginn, E. J. Continuous Fractional Component Monte Carlo: An Adaptive Biasing Method for Open System Atomistic Simulations. J. Chem. Theory Comput. 2007, 3, 1451−1463. (22) Rowlinson, J. S.; Widom, B. Molecular Theory of Capillarity; Clarendon: Oxford, U.K., 1982. (23) Kirkwood, J.; Buff, F. The Statistical Mechanical Theory Of Surface Tension. J. Chem. Phys. 1949, 17, 338−343. (24) Irving, J.; Kirkwood, J. The Statistical Mechanical Theory Of Transport Processes 0.4. The Equations Of Hydrodynamics. J. Chem. Phys. 1950, 18, 817−829. (25) Walton, J.; Tildesley, D.; Rowlinson, J.; Henderson, J. The Pressure Tensor at The Planar Surface of A Liquid. Mol. Phys. 1983, 48, 1357−1368. (26) Gloor, G.; Jackson, G.; Blas, F. Test-area Simulation Method for the Direct Determination of the Interfacial Tension of Systems with Continuous or Discontinuous Potentials. J. Chem. Phys. 2005, 123, 134703. (27) Ibergay, C.; Ghoufi, A.; Goujon, F.; Ungerer, P.; Boutin, A.; Rousseau, B.; Malfreyt, P. Molecular Simulations of the n-alkane Liquid-Vapor Interface: Interfacial Properties and Their Long Range Corrections. Phys. Rev. E 2007, 75, 051602. (28) Biscay, F.; Ghoufi, A.; Lachet, V.; Malfreyt, P. Calculation of the Surface Tension of Cyclic and Aromatic Hydrocarbons from Monte Carlo Simulations Using an Anisotropic United Atom Model (AUA). Phys. Chem. Chem. Phys. 2009, 11, 6132−6147. (29) Vega, C.; de Miguel, E. Surface Tension of the Most Popular Models of Water by Using the Test-Area Simulation Method. J. Chem. Phys. 2007, 126, 154707. (30) Frenkel, D.; Smit, B. Understanding Molecular Simulation; Academic Press: San Diego, CA, 2002. (31) Essmann, U.; Perera, L.; Berkowitz, M.; Darden, T.; Lee, H.; Pedersen, L. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103, 8577−8593. (32) Gelest Inc. http://shop.gelest.com/Product.aspx?catnum/DMST11&Index/0&TotalCount/1, accessed 2015-05-01. (33) Chang, C. Volume Expansion Coefficients and Activitycoefficients of High-pressure Carbon-dioxide Dissolution in Organic Liquids at 298-K. J. Chem. Eng. Jpn. 1992, 25, 164−170. (34) NIST Web Book. http://webbook.nist.gov/cgi/cbook. cgi?Contrib/LMF, Accessed 2015-04-02. (35) Rosch, T. R.; Maginn, E. J. Letter to the Editor. AIChE J. 2011, 57, 1100−1101. (36) Sauer, B.; Dee, G. Molecular-weight and TemperatureDependence of Polymer Surface-Tension - Comparison of Experiment with Theory. Macromolecules 1991, 24, 2124−2126. (37) Dee, G.; Sauer, B. The Surface-Tension of Polymer Blends Theory and Experiment. Macromolecules 1993, 26, 2771−2778. (38) Yang, D.; Tontiwachwuthikul, P.; Gu, Y. Interfacial Tensions of the Crude Oil Plus Reservoir Brine + CO2 Systems at Pressures up to 31 MPa and Temperatures of 27 Degrees C and 58 Degrees C. J. Chem. Eng. Data 2005, 50, 1242−1249. (39) Reyes, G.; Cartes, M.; Rey-Castro, C.; Segura, H.; Mejia, A. Surface Tension of 1-Ethyl-3-methylimidazolium Ethyl Sulfate or 1Butyl-3-methylimidazolium Hexafluorophosphate with Argon and Carbon Dioxide. J. Chem. Eng. Data 2013, 58, 1203−1211.

M

DOI: 10.1021/acs.jpcc.5b05806 J. Phys. Chem. C XXXX, XXX, XXX−XXX