Molecular Simulations of the Thermophysical ... - ACS Publications

Article pubs.acs.org/JPCC

Molecular Simulations of the Thermophysical Properties of Polyethylene Glycol Siloxane (PEGS) Solvent for Precombustion CO2 Capture Wei Shi,*,†,‡,¶ Nicholas S. Siefert,† Hseen O. Baled,† Janice A. Steckel,† and David P. Hopkinson† †

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: The thermophysical properties for neat polyethylene glycol siloxane solvent (PEGS) along with CO2, H2, H2O, and H2S gas absorption in PEGS at 298−373 K were investigated via molecular simulations. The predicted neat PEGS density, heat capacity, surface tension, and CO2 and H2 solubilities in PEGS solvent agree reasonably well with the experimental data, with typical differences of 0.8−20%, while the predicted PEGS solvent viscosity is 1.7−2.5 times larger than the experimental data. Gas solubility in PEGS at 298 K decreases in the following order, H2O (31000) > H2S (230) > CO2 (33) > H2 (1), which follows the same order as the gas−PEGS interaction. In contrast, gas diffusivity in PEGS at 298 K decreases in an opposite way, H2 (1) > CO2 (0.22) ≈ H2S (0.12) > H2O (0.018). The numbers in parentheses are the corresponding values relative to H2. Compared to the widely studied poly(dimethylsiloxane) (PDMS) solvent, PEGS is more hydrophilic due to its stronger interaction with H2O and fewer branched −CH3 groups, which in turn leads to fewer hydrophobic pockets. The CO2/ H2 solubility selectivity in PEGS is larger than that in PDMS due to a stronger interaction with CO2 in PEGS. Finally, it was found that CO2 absorption in PEGS could significantly improve the CO2−PEGS solution dynamics by 5−6 times, resulting in a decrease in solution viscosity and increase in diffusivity. These CO2 absorption effects are due to solution volume expansion upon CO2 absorption compared to the neat PEGS solvent volume and the possibility that CO2 acts as a “lubricant” to decrease the solvent−solvent interaction.

1. INTRODUCTION Countries across the globe are acting together to reduce greenhouse gas emission because of growing evidence of harm from rising temperatures and rising sea levels.1,2 Chemical and physical absorption methods have been suggested for CO2 postcombustion and precombustion capture, respectively, from coal-based power generation emissions due to low CO2 partial pressure (∼0.1 bar) in flue gas and high CO2 partial pressure (∼25 bar) in fuel gas. Due to rich H2O concentration in the hot or warm gas stream corresponding to CO2 precombustion capture, hydrophobic solvents are preferred for selective CO2 absorption over H2 and H2O.3,4 By using a hydrophobic solvent, it is possible to retain moisture with high thermal energy in the H2 stream that is necessary for combustion in a gas turbine. Note that although polyethylene glycol (PEG) has good CO2/H2 solubility selectivity, which is partly due to strong CO2 interactions with −O groups of PEG, the PEG solvent itself is very hydrophilic and consequently not appropriate for precombustion carbon capture application at high temperatures, under which conditions valuable water steam in a hot/warm gas stream will be absorbed in PEG. Indeed, designing a solvent to capture CO2 from the warm or hot gas stream that contains H2O and H2 is very challenging. © XXXX American Chemical Society

The solvent must have the following properties: thermal stability, hydrophobicity, low vapor pressure, large CO2 solubility, large CO2/H2 solubility selectivity, low viscosity, no foaming tendency, safety of operations, and low price. In previous work,5,6 we have both theoretically and experimentally studied a hydrophobic solvent poly(dimethylsiloxane) (PDMS), which is very hydrophobic and even exhibits ∼8 times smaller H2O solubility compared to the H2S solubility. Unfortunately, this PDMS solvent exhibits relatively low CO2/H2 solubility selectivity (∼19).5 Consequently, other hydrophobic solvents with large CO2 solubility and large CO2/H2 solubility selectivity are still needed for precombustion CO2 capture. In one previous work,6 a new polyethylene glycol siloxane (PEGS) solvent was experimentally determined to have a high CO2/H2 solubility selectivity (∼40). Note that in that work, the PEGS solvent was given a different name with the abbreviation of HPDMS.6 This PEGS solvent is a hybrid molecule that consists of parts of a PDMS chain and parts of a PEG chain. This structure is expected to Received: July 7, 2016 Revised: August 15, 2016

A

DOI: 10.1021/acs.jpcc.6b06810 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C give the material hydrophobicity that is characteristic of PDMS combined with improved CO2 solubility and CO2/H2 solubility selectivity that is typical for PEG. Although in the previous work6 the Advanced System for Process Engineering (ASPEN) process modeling indicates that the PEGS solvent is promising for CO2 precombustion capture application, many questions remain to be addressed. For example, although PEGS contains parts of PEG and PDMS and as a result PEGS exhibits hydrophobicity and CO2/H 2 solubility selectivity between PDMS and PEG, one still needs to understand and quantify these effects. What are the gas solubility and gas diffusivity values for CO2, H2, H2S, and H2O absorption in PEGS? In any ASPEN process modeling, the solution thermophysical properties significantly affect the technoeconomic analysis. As a result, CO2 absorption effects on physical properties, such as gas diffusivity, solution viscosity, and surface tension, are important to understand. The goal of this work is to address all of the above questions using molecular simulations. To ensure that simulations are as predictive as possible, only the experimental PEGS solvent density at ambient condition was used to slightly tune the classical force field (FF) parameters for the PEGS molecule. Additionally, to verify theoretical predictions, experiments were also implemented in this work to investigate CO2 absorption effects on the solution viscosity.

Figure 1. Optimized structure for the PEGS molecule obtained from quantum AI gas-phase calculations at the B3LYP/6-311++G(d,p) level of theory. Also shown are the atomic partial charges for Si (yellow), O (red), C (silver), and H (white) atoms obtained from quantum calculations by using the CHELPG protocol16 at the same level of theory as optimization.

molecule were obtained from quantum calculations by using the CHELPG protocol16 at the B3LYP/6-311++G(d,p) level of theory. The LJ parameters were obtained from previous work5,11 except that the σ values were decreased by 1% to reproduce the experimental PEGS solvent density at ambient conditions. All FF parameters for the PEGS molecule are given in the Supporting Information. In order to calculate transport properties such as diffusivity and viscosity, typically 10−20 ns NVE molecular dynamics (MD) simulations were run by following the NPT and NVT MD simulations. More simulation details were given in previous work.9 The time step in NVE MD simulation was set to be 0.5 fs. Coordinates and pressure tensors were dumped every 1000 and 1 fs, respectively. All MD simulations in this work were performed by using the NAMD program.17 The smooth particle Mesh Ewald method18 was used to calculate the electrostatic interaction in MD simulations. 2.2. Quantum Ab Initio (AI) Calculations. To evaluate the accuracy of CO2, H2, H2O and H2S gases’ interactions with the PEGS molecule obtained from classical FF calculations, quantum AI gas-phase calculations were performed for gas− PEGS dimers by using the Gaussian 09 program.10 Geometry optimizations were performed at the B3LYP/6-311++G(d,p) level of theory followed by single-point energy calculations at the MP2/cc-pVTZ level to compute gas−PEGS interactions. To account for the basis set superposition error, we used counterpoise corrections in both geometry optimization and single-point energy calculations. The gas−PEGS interaction was calculated to be the energy required to separate two monomers in the optimized dimer structure to an infinitely large distance but keeping the monomer structures frozen at their geometries in the optimized dimer system. For each gas− PEGS dimer, typically 5−6 different initial dimer configurations were used for geometry optimization and interaction energy calculations. Only the strongest gas−PEGS interactions were reported in this work. Additionally, we also performed gas− PEGS interaction energy calculations on those optimized dimer structures at the RI-MP2 level of theory with the cc-pVDZ, ccpVTZ, and cc-pVQZ bases, and the interaction energies were extrapolated to the complete basis set limit (Supporting Information). 2.3. Gas Solubility. CO2 and H2 solubilities in PEGS solvent at 298−373 K were obtained from the Continuous Fractional Component (CFC) Monte Carlo (MC) calcula-

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

∑

kr(r − r0)2 +

bonds

+

∑ dihedrals N−1

+

∑ i=1

∑

k θ(θ − θ0)2

angles

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

∑

k ψ (ψ − ψ0)2

impropers

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

(1)

where the symbols represent their conventional meanings.7 Standard Lorentz−Berthelot combining rules were used to calculate the mixed Lennard-Jones (LJ) interaction parameters. The LJ potential cutoff distance 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,8 H2,9 H2O,8 and H2S5 were taken from previous work. For the PEGS solvent molecule (Figure 1), which contains 68 atoms and has a molecular weight of 426.767 g/mol, a flexible atomic model was used. The r0 values were obtained from quantum ab initio (AI) gas-phase geometry optimization calculations at the B3LYP/6-311+ +G(d,p) level of theory by using the Gaussian 09 program.10 The kr values were obtained from our previous work5 and the OPLS-DMEFF model11 for similar chemical bonds. For the angle bending potential, the θ0 values were obtained from quantum geometry optimization calculations. The kθ values were obtained from our previous work,5 the OPLS-DMEFF model,11 and the CHARMM par_all36_cgenff.prm FF.12 The dihedral potential parameters were obtained from the published work of others.11−15 The atomic charges for the PEGS B


Article

The Journal of Physical Chemistry C

the cubelet center and any solvent atom was less than rev = s × (σi/2 + rprobe), where σi is the LJ parameter for the ith atom of the PEGS molecule. The s is a scaling factor, and it was set to 0.8. The radius for the probe atom (rprobe) was set to 1.0 Å. More calculation details for the modified test-particle insertion method can be found in our previous work.21 Similar to the previous work, we attempted to insert5 C (CO2), O (H2O), and S (H2S) heavy atoms in the empty cubelets. The FGSF method was also used to compute electrostatic interactions in the test-particle insertion method to speed up the calculation. 2.5. Gas and Solvent Self-Diffusivity. Gas and solvent self-diffusivities were calculated by using the Einstein relation; the mean-square displacement for the center of mass of each molecule was obtained from the NVE MD simulation. More calculation details can be found in previous work.9,23 2.6. Viscosity. The viscosity could be calculated from the NVE MD simulation by using the following Einstein relation

tions8,9,19 by using in-house software. More CFC simulation details can be found in our previous work.8,9,19 For CO2 solubility calculations, simulations were performed at PCO2 = 1.25−30 bar by using a cubic solvent box; the box contained 100 PEGS molecules and had a box length of ∼43 Å. For H2 solubility calculations, simulations were performed at PH2 = 10− 400 bar, and the cubic box contained 100 PEGS molecules. The starting configurations for PEGS molecules were obtained from NPT MD simulations. In order to speed up calculations,8,9 the Fennel and Gezelter shift force (FGSF)8,20 method was used to compute the electrostatic interaction in CFC MC simulations. 2.4. Gas Henry’s Law Constant (H). In addition to CO2 and H2 full adsorption isotherms obtained from the above CFC MC calculations, the H values for CO2, H2O, and H2S gas absorption in PEGS were also calculated by using a modified test-particle insertion method,21 which is implemented in our in-house software package. The cubic simulation box contained 100 PEGS solvent molecules. When generating the snapshots for H calculation, three independent NPT MD simulations were performed by starting from different initial configurations. Simulations were carried out at 298−373 K and 1 bar. For each NPT MD simulation, about 20 000−40 000 equilibrated liquid structures corresponding to the last 10−20 ns production runs were saved for later test-particle insertion calculation. When calculating H from the saved snapshots, two million trial insertions were attempted for each snapshot. To improve the insertion efficiency, the cubelet22 and bias-removing methods21 were used. The cubelet length was set to be ∼0.8 Å. A cubelet was marked as occupied if the distance between A(t ) = ⟨(

∫t

⟨(∫ =

t0+ t

η = lim

t →∞

V k bT

∫0

0

t0+ t

Pζν(t ′) dt ′)2 ⟩t0 Pζν(t ′) dt ′ −

∫0

t0

Pζν(t ′) dt ′)2 ⟩t0 (2)

where V is the volume of the simulation box, kb is the Boltzmann constant, T is the temperature, Pζν is the symmetrized pressure tensor, and ⟨⟩t0 indicates the average obtained from multiple time origins. To further improve the statistics, the following average was used

Pζν(t ′) dt ′)2 ⟩t0

0

t0+ t

t0

Pxy(t ′) dt ′)2 ⟩t0 + ⟨(∫

t0+ t

Pxz(t ′) dt ′)2 ⟩t0 + ⟨(∫

t0+ t

t0

t0

Pyz(t ′) dt ′)2 ⟩t0

3

d log A(t ) , d log t

have computed the viscosities for water and n-decane. The simulated viscosity values agree with each other between the Einstein and the Green−Kubo methods. Additionally, our simulated viscosities are in reasonable agreement with the simulation values obtained by other groups.24,26 More details are given in the Supporting Information. These results suggest that we have correctly implemented the viscosity calculation. In the rest of this work, we only used the Einstein method to calculate the viscosity values for PEGS solvent and the CO2− PEGS solution. 2.7. Heat Capacity. The heat capacities for neat PEGS solvent at 298−373 K were calculated by using a similar procedure as that used by Maginn and co-workers.27 The molar enthalpy (Hm) is given by

a time period (t1

∞

dt ⟨Pζν(t )Pζν(0)⟩

t0+ t

∫0

≤ t ≤ t2) is identified during which the β average value is close to 1. Second, dA(t ) and η are calculated in the same time period dt t1 ≤ t ≤ t2. The viscosity could also be obtained from the Green−Kubo relation η=

∫t

1 V d = lim ⟨( t →∞ 2 k bT dt

The procedure to compute viscosity from eq 2 is the same as the self-diffusivity calculation by using a similar Einstein relation.9 First, by calculating β =

1 V d ⟨( 2 k bT dt

(3)

Compared with the Green−Kubo method, the Einstein method has a practical advantage24 in that the integration of the pressure tensor over the whole simulation time is performed when calculating η. In contrast, when the Green−Kubo method is used, one must rely on accurate integration of a rapidly varying value (⟨Pζν(t)Pζν(0)⟩) over a time interval much shorter than the entire simulation time.24 The integration in this short time interval is then fitted to a mathematical function of η − b × exp[−(t/τ)] to obtain the viscosity η.25 We have implemented both the Einstein and Green−Kubo methods in the in-house software to calculate viscosity. To verify viscosity calculations by using the in-house software, we

Hm = Umext + Umint + U kin + PVm = (Umext + PVm − RT ) + (Umint + U kin + RT ) = Hmres + Hmig

(4)

where the subscript m and superscripts ig and res indicate the molar property, the ideal gas heat capacity, and the liquid residual heat capacity, respectively. The Uext m indicates the interaction potential energy between different molecules. The Uint m indicates the potential energy within a single molecule; it C


Article

The Journal of Physical Chemistry C consists of Uint,bond and Uint,nonbond . The Uint,bond includes the m m m bond, angle bending, dihedral, and improper potential energies; the Uint,nonbond contains the 1−4, 1−5, 1−6, and so forth LJ and m electrostatic interaction energies between two atoms in the same molecule. The Ukin indicates the kinetic energy, Vm is the liquid molar volume, and R is the ideal gas constant. The constant pressure molar heat capacity (CP,m) is given by

simulation involving vapor−liquid equilibrium and the interface between them. The second step is to calculate the potential energy change for each configuration by fixing the simulation volume V but varying the interfacial surface area. Except for γsim obtained from the test-area method, a tail correction (γtail) must be added to γsim to obtain the total surface tension (γsim,total). More simulation details for γsim and (γtail) can be found in our previous work.5 For the neat PEGS solvent, about 60 000−200 000 configurations corresponding to the last 30−100 ns of NVT production runs were saved at 298−373 K. The simulation box sizes were around 75 Å × 75 Å × 165 Å, and they contained 500 PEGS molecules. For the CO2−PEGS solution at 298 K, the simulation box size was 79.084 Å × 79.084 Å × 134.810, and it contained 500 PEGS solvent molecules and 1500 CO2 molecules. About 80 000 configurations corresponding to the last 40 ns of NVT production runs were saved. Similar to our previous work,5 the ζ value was set to be 0.0001 to scale down/ up the interfacial surface area. Additionally, it was also found that the Ewald and FGSF methods give consistent γsim values for the PEGS solvent. To speed up the calculation, the FGSF method was used to calculate electrostatic interactions in the above second step.

⎛ ∂⟨H ig⟩ ⎞ ⎛ ∂⟨Hmres⟩ ⎞ ⎛ ∂⟨Hm⟩ ⎞ m ⎟ =⎜ C P,m = ⎜ ⎟ ⎟ +⎜ ⎝ ∂T ⎠ P ⎝ ∂T ⎠ P ⎝ ∂T ⎠ P res ig = C P,m + C P,m

(5)

where the symbol ⟨⟩ indicates the average. The liquid residual res = heat capacity is given by C P,m

(

∂⟨Umext + PVm⟩ ∂T

)

− R . To

P

calculate Cres P,m at T and 1 bar, three independent sets of NPT MD simulations were implemented at T − 5, T, and T + 5 and 1 bar; the resulting ⟨Uext m + PVm⟩ values versus temperatures were linearly regressed to estimate 17

Uext m

∂⟨Umext + PVm⟩ . Note that in ∂T P Uint,nonbond value was output. m

(

)

the NAMD program, the + The Uint,nonbond value for each PEGS molecule was calculated by m using the in-house software and subtracted from the NAMD output to obtain Uext m. The CigP,m = CigV,m + R was obtained from quantum AI gasphase calculations instead of from classical FF calculations; the quantum calculations are expected to be more accurate to predict CigV, m compared with the classical FF calculations. In quantum AI calculations, frequency analysis for the optimized PEGS structure was implemented to obtain the vibrational temperatures Θv,j. Along with the translational (3/2R) and rotational (3/2R for nonlinear and R for linear molecule) contributions to CigV,m, the vibrational contribution at finite temperature T was calculated to be R ∑j e Θv ,j / T

(

Θv , j / T e Θv , j / T −

3. RESULTS AND DISCUSSION 3.1. PEGS Density. The simulated PEGS solvent density values are very close to the experimental data (Figure 2), with

2

1

).

The geometry optimization and frequency calculation were both performed at the B3LYP/6-311++G(d,p) level. Note that the harmonic oscillator description is used as an approximate treatment for a low-frequency vibrational motion, such as torsional motion. To verify the above procedure for heat capacity calculation, the heat capacity values for liquid water and benzene were calculated. The calculated CigP,m values for water and benzene molecules were found to be in good agreement with the experimental data, with small differences of 0.3−3%. The CP,m values for liquid water and benzene are also comparable to the experimental data, with differences of 6−19%. More details are given in the Supporting Information. The reasonable agreement in water and benzene heat capacity values between simulations and the experimental data suggests that we have appropriately implemented the heat capacity calculation. 2.8. Surface Tension. There are two types of methods to compute the surface tension from molecular modeling, that is, the mechanical28−31 and the test-area methods.32 In this work, the test-area method, which has been implemented in our inhouse software package,5 was used to compute the surface tension values for the neat PEGS solvent and CO2−PEGS solution at 298−373 K and PCO2 = 1−30 bar. Similar to the Henry’s law constant calculation21 by using the test-particle insertion method, there are two steps in the test-area method. The first step is to generate configurations from the NVT MD

Figure 2. Simulated and experimental6 PEGS solvent density at 1 bar and different temperatures. The lines are linear fittings to the simulation (solid) and experimental values (dashed).

an absolute average difference of 0.8% at 293−343 K. Note that only the experimental density at 298 K and 1 bar was used to slightly tune the σ values for each atom of the PEGS molecule (section 2.1). All of the FF parameters were then fixed to predict PEGS densities at other temperatures and to calculate other thermodynamic and transport properties (below). The density agreement between simulation and experiment at temperatures different than 298 K indicates that the PEGS FF parameters in this work are transferrable to predict density at different temperatures. Note that slightly larger density differences between the simulation and experiment occur at higher temperatures (Figure 2). 3.2. CO2, H2, H2S, and H2O Interactions with PEGS. The optimized gas−PEGS dimer structures and the corresponding gas−PEGS interaction energies are shown in Figure 3 and Table 1. The classical FF gives very similar CO2, H2, and H2S interaction energies compared to the MP2/cc-pVTZ calculations, with small differences of 1.0−1.7 kJ/mol. For H2O, the classical FF overestimates the H2O−PEGS interaction energy by 8 kJ/mol compared to the MP2/cc-pVTZ calculation. Note D


Article


mainly due to its large electrostatic interaction (Table 1) through two hydrogen-bond interactions (Figure 3 c). In contrast, the van der Waals interaction energy predominates for the H2S−PEGS interaction, which is similar to the H2S interaction with PDMS.5 3.3. Gas Solubility. 3.3.1. CO2 and H2 Absorption Isotherms and Gas−PEGS Solution Volumes. The simulated CO2 and H2 solubilities are shown in Figure 4 and Table 2. The

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

Table 1. Interaction Energy (INT) among CO2, H2, H2O, and H2S Gases and the PEGS Molecule Obtained from Both Quantum AI and Classical FF Calculationsa FF INT (kJ/mol) system

MP2/cc-pVTZ INT (kJ/mol)

RI-MP2 INT (kJ/mol)

total

VDW

ELEC

CO2 H2 H2O H2S

−12.8 −3.0 −34.7 −24.7

−15.0 −3.3 −40.4 −29.4

−14.1 −2.0 −42.4 −23.0

−5.6 −2.0 −3.8 −14.9

−8.5 −0.0 −38.6 −8.1

Figure 4. Simulated CO2 (open squares) and H2 (open circles) absorption isotherms in PEGS solvent at 298 K. Simulation error bars are smaller than the symbols and are not shown. The simulated H2 solubilities at 10−400 bar are shown in the inset. The corresponding experimental CO2 (filled squares) and H2 (filled circles) solubility values6 at 298 K are also shown for comparison. The lines are guides for the eye.

a

In the AI calculation, the gas−PEGS dimer optimized structures (Figure 3) were 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 MP2/cc-pVTZ theory. For comparison, the interaction energy values at the complete basis set limit, which were obtained by extrapolating the RI-MP2 calculations at different basis sets, are also shown. In the classical FF calculation, full van der Waals (VDW) and electrostatic (ELEC) interactions were calculated without truncation. The optimized dimer structures and the corresponding INT energy values in classical FF calculations were obtained from 100 steps of energy minimization by using the B3LYPoptimized dimer structures (Figure 3) as the starting configurations. Note that the root-mean-square deviation values in the atomic positions between the FF-optimized and B3LYP-optimized structures were found to be very small, about 0.06−0.08 Å.

simulated CO2 solubilities at 298 K are very close to the experimental data6 (Figure 4), with differences less than 10%. For H2, the simulated gas solubilities at 298 K are also comparable to the experimental data6 (Figure 4), with differences of 30%. From the gas−PEGS solution molar volume (Vmix) in Table 2, the solution volume expansion relative to the neat PEGS solvent volume could also be calculated. It was found that CO2 absorption could appreciably increase the solution volume. For example, at PCO2 = 30 bar, the solution volume could expand by 24%. In contrast, H2 absorption has negligibly small effects on solvent volume. At 298 K and 10−400 bar, the H2−PEGS solution volume only changes by less than 0.5% compared to the neat PEGS solvent volume. 3.3.2. Henry’s Law Constants for CO2, H2, H2O, and H2S. Although full absorption isotherms (Figure 4) are important and required in detailed process modeling,6 they are more expensive to calculate compared with the Henry’s law constant obtained from the test-particle insertion method. Consequently, H values were calculated for all gases in PEGS at different temperatures. The H values for CO2, H2S, and H2O are consistent with each other among different simulation runs (Table 3), which suggest that both the number of simulation snapshots and the number of test-particle insertions for each snapshot are large enough to accurately calculate H. Note that for CO2, the average value of H for three independent testparticle insertion calculations (Table 3) was estimated to be 16.9 ± 0.3 bar at 298 K, comparable to the Henry’s law constant value of 18.7 bar, which was directly estimated from the CO2 solubility at 298 K and 1.25 bar (Table 2). The Henry’s law constant agreement between the test-particle insertion and the CFC MC methods indicates that we have correctly implemented these two methods.

that for each gas−PEGS pair, another five to six sets of classical FF and quantum AI calculations were also performed by putting the gas molecule in a different initial position relative to the PEGS molecule. For all of these calculations, a similar observation was obtained, that is, the classical FF and AI calculations give consistent CO2, H2, and H2S interactions with PEGS; for H2O, the classical FF gives stronger H2O−PEGS interaction compared to the MP2/cc-pVTZ calculation. The interaction energy results at the complete basis set limit from RI-MP2 calculations are also comparable to the classical FF calculations (Table 1 and Supporting Information). For CO2 and H2, the differences are 0.9−1.3 kJ/mol; for H2S and H2O, the differences are 6.4 and 2 kJ/mol, respectively. Additionally, it was found that CO2 interactions with each of four O atoms of the ethylene glycol segment are 3−7 kJ/mol stronger compared with CO2 interactions with each of two O atoms of the siloxane segment. Finally, both classical FF and AI calculations show that the gas−PEGS interaction decreases in the following order: H2O > H2S > CO2 > H2. The strong H2O interaction with PEGS is E


Article


Table 2. Simulated Gas Solubility (Mole Fraction xmolar and Mass Fraction xmass) and Solution Molar Volume (Vmix) Obtained from the CFC MC Methoda gas

T (K)

P (bar)

f (bar)

xmolar

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

298.2 298.2 298.2 298.2 298.2 298.2 298.2 298.2 298.2 298.2 298.2 298.2 333.2 333.2 333.2 333.2 333.2 333.2 373.2 373.2 373.2 373.2 373.2

1.25 2.5 7.004 14 20.0 30.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 10.0 20.0 40.0 100.0 400.0

1.2415 2.466 6.7385 12.9486 17.87 25.26 10.03 20.12 40.48 103.1 214.2 468.4 10.03 20.12 40.48 103.2 213.8 464.8 10.03 20.12 40.48 103.1 460.4

0.067 (3) 0.123 (4) 0.285 (4) 0.443 (3) 0.589 (3) 0.702 (5) 0.016 (1) 0.034 (1) 0.073 (1) 0.162 (2) 0.245 (4) 0.370 (3) 0.0191 (6) 0.0369 (8) 0.069 (2) 0.158 (2) 0.277 (2) 0.429 (6) 0.0250 (4) 0.0468 (8) 0.095 (2) 0.216 (2) 0.499 (4)

Vmix (cm3/mol)

xmass 0.0074 (3) 0.0143 (6) 0.0397 (8) 0.0759 (9) 0.129 (1) 0.196 (4) 78 (3) 166 (4) 371 (6) 910 (14) 1526 (29) 2759 (36) 92 (3) 181 (4) 348 (8) 884 (11) 1801 (16) 3527 (77) 121 (2) 232 (4) 499 (9) 1296 (14) 4685 (82)

× × × × × × × × × × × × × × × × ×

10−6 10−6 10−6 10−6 10−6 10−6 10−6 10−6 10−6 10−6 10−6 10−6 10−6 10−6 10−6 10−6 10−6

427.1 403.8 337.5 274.2 215.0 168.0 447.5 439.6 425.4 385.8 344.6 286.4 463.8 456.3 441.3 399.8 345.8 274.4 494.0 482.0 462.0 403.3 258.7

(12) (18) (17) (13) (11) (21) (8) (5) (3) (8) (14) (12) (5) (5) (4) (5) (6) (21) (7) (5) (3) (6) (15)

a The numbers in parentheses are uncertainties in the last digit obtained from standard block average calculations. Gas fugacity ( f) was obtained from the Peng−Robinson equation of state.

Table 3. Simulated Henry’s Law Constants for CO2, H2S, H2O, and H2 Absorption in PEGS Solvent at 298 Ka Henry’s law constant (bar) system

RUN 1

RUN 2

RUN 3

CO2 H2S H2O H2

17.3 (3) 2.45 (4) 0.018 (2) 558 (13)

16.5 (2) 2.38 (3) 0.017 (2)

16.9 (2) 2.45 (4) 0.019 (1)

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 simulations for neat PEGS, 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 the H2 pressure versus H2 solubility (mole fraction) below 40 bar (Table 2), not from the test-particle insertion method. Consequently, only one Henry’s law constant is given for H2. Uncertainties in the last digit are given in parentheses.

Figure 5. Simulated Henry’s law constant (H) versus 1/T for CO2, H2S, H2O, and H2 absorption in the PEGS solvent. Simulation error bars are smaller than the symbol size and are not shown. The lines are a linear fitting of the simulation values.

indicates that H2 solubility in PEGS increases when the temperature is increased; the increased H2 solubility is due to an increase of the solvent free molar volume at elevated temperature. For example, the PEGS molar volume and the free volume fraction at 298 K are 456.4 cm3/mol and 0.0980, respectively; at 373 K, those corresponding values are 509.8 cm3/mol and 0.1623. As a result, the PEGS free molar volume (molar volume × free volume fraction) is increased from 44.7 cm3/mol at 298 K to 82.7 cm3/mol at 373 K, by a factor of 1.85 times. The heats of absorption for H2O, H2S, and CO2 are comparable to the corresponding gas−PEGS interactions (Table 1), and their solubilities in PEGS decrease at elevated temperatures (Figure 5). Finally, it is interesting to compare the structures of the PEGS and PDMS molecules.5 The ring-like PDMS molecule is larger than the linear PEGS molecule. Additionally, the PDMS

Gas solubility in PEGS at 298 K decreases in the following order (Table 3): H2O (31000) > H2S (230) > CO2 (33) > H2 (1); the numbers in parentheses indicate the relative gas solubility compared to H2 solubility. The gas solubility decreases in the same order as the gas−PEGS interaction (Table 1). The very small water H value (0.018 bar) in PEGS implies that the PEGS solvent is more hydrophilic compared to the PDMS solvent, which exhibits a much larger water H value of 10.6 ± 0.5 bar at 298 K. By plotting log(H) vs 1/T (Figure 5), the heat of gas ∂ log(H )

absorption in PEGS was calculated as R ∂(1 / T ) . The heats of

gas absorption were found to be 1.5 ± 0.6 kJ/mol for H2, −11.7 ± 0.7 kJ/mol for CO2, −17.2 ± 0.4 kJ/mol for H2S, and −46.1 ± 0.5 kJ/mol for H2O. The positive H2 heat of absorption

F


Article


Table 4. Simulated Gas (g) and Solvent (solv) Self-Diffusivity (D) Coefficients at 298-373 K and 1 bar for CO2, H2, H2S, and H2O Absorption in the PEGS Solventa system CO2−PEGS CO2−PEGS CO2−PEGS H2−PEGS H2−PEGS H2−PEGS H2O−PEGS H2O−PEGS H2O−PEGS H2S−PEGS H2S−PEGS H2S−PEGS

T (K) 298 333 373 298 333 373 298 333 373 298 333 373

βg

Dg (m2/s) 1.75 3.26 4.9 7.8 1.98 2.21 1.41 4.6 1.49 9.6 2.65 5.82

(3) (5) (1) (3) (4) (3) (6) (1) (4) (3) (3) (8)

× × × × × × × × × × × ×

−9

10 10−9 10−9 10−9 10−8 10−8 10−10 10−10 10−9 10−10 10−9 10−9

1.06 (2) 0.97 (1) 0.90 (2) 1.04 (4) 1.08 (2) 0.95 (1) 1.05 (4) 0.95 (2) 0.94(2) 1.05 (3) 1.00 (1) 1.02 (1)

βsolv

Dsolv (m2/s) 5.1 2.04 5.06 6.2 2.10 5.46 4.5 1.64 4.8 4.1 1.94 5.82

(1) (2) (5) (1) (2) (3) (1) (2) (1) (1) (2) (6)

× × × × × × × × × × × ×

−11

10 10−10 10−10 10−11 10−10 10−10 10−11 10−10 10−10 10−11 10−10 10−10

0.93 (1) 1.00 (1) 0.98 (1) 1.02 (1) 1.01 (1) 1.007 (6) 0.92 (1) 0.94 (1) 0.96 (2) 0.87 (1) 1.00 (1) 1.03 (1)

Gas mole fractions in the solution are 0.05. The β = d(log Δr2)/d(log t) values are also shown. Uncertainties in the last digit are given in parentheses.

a

molecule contains extensive −CH3 groups, which are attached to the backbone Si atoms to form hydrophobic pockets. Out of all of the non-hydrogen heavy atoms (backbone O, Si, and branched C(−CH3) atoms), 52% (atom number based) of them are −CH3 groups. In contrast, the PEGS molecule contains a limited number of −CH3 groups, and they are attached to only one end of it (Figure 1), which causes PEGS to have fewer hydrophobic regions. Only 31% of all nonhydrogen atoms (backbone O, Si, C(−CH2), and branched C(−CH3)) in PEGS belongs to the −CH3 group. Also, the H2O−PEGS interaction is 12−20 kJ/mol stronger compared with the H2O−PDMS interaction. These factors make the PEGS solvent much more hydrophilic; water solubility in PEGS at 298 K is about 600 times larger than that in PDMS in the limit of the Henry’s law constant region. The relatively high H2O solubility in PEGS implies that H2O absorption in PEGS may affect CO2 loading, which is not addressed here. Note that the PDMS solvent has a molar volume of 1344 cm3/mol and a free volume fraction of 0.118 at 298 K,5 which leads to a free molar volume of 159 cm3/mol. The PEGS solvent has a 2.9 times smaller molar volume, 1.2 times smaller free volume fraction, and 3.6 times smaller free molar volume compared to the PDMS solvent. As a result, H2 solubility in PEGS at 298 K is 3.6 times smaller (mole-fraction-based, mol H 2 ) and 1.2 times smaller (solvent-volume-based,

based H2 solubility) but interacting strongly with CO2 (leading to large mole-fraction-based CO2 solubility) would exhibit a large CO2/H2 solubility selectivity. 3.4. Self-Diffusivity (D). For both gas and solvent diffusion in solution, most β values are close to 1 (Table 4). These findings suggest that gases and PEGS solvent molecules are in the diffusive region and simulations are long enough to reliably calculate the self-diffusivity. At 298 K, gas diffusivity in PEGS decreases in the following order: H2 (1) > CO2 (0.22) ≈ H2S (0.12) > H2O (0.018) (Table 4 and Figure 6). The numbers in parentheses denote

Figure 6. Simulated gas self-diffusivity (D) versus 1/T for CO2, H2S, H2O, and H2 gas absorption in the PEGS solvent. Simulation error bars are smaller than the symbols and are not shown. The lines are linear regressions of the simulation data.

mol H 2 + mol solvent mol H 2 ) compared solvent volume

to the H2 solubility in PDMS. The significantly smaller H2 solubility in PEGS compared to that in PDMS (mole-fraction-based) is due to the significantly smaller free molar volume for PEGS; the slightly smaller H2 solubility in PEGS (solvent-volume-based) is due to its slightly smaller free volume fraction. The CO2 solubility in PEGS at 298 K is 2.1 times smaller than that in PDMS (mole-fraction-based), which is mainly due to the significantly smaller PEGS free molar volume. On the one hand, the PEGS solvent free volume fraction is slightly smaller compared to that of PDMS. However, the CO2−PEGS interaction is ∼3.0 kJ/mol stronger compared to the CO2−PDMS interaction. These two factors result in the CO2 solubility (solvent-volume-based) in PEGS being 1.4 times larger than that in PDMS. Both the molefraction- and solvent-volume-based CO2 and H2 solubilities indicate that CO2/H2 solubility selectivity in PEGS is 1.7 times larger than that in PDMS. Our findings suggest that solvents having small free molar volume (leading to small mole-fraction-

the relative gas diffusivity over H2. Both the gas mass and the gas−solvent interaction are important to determine gas diffusivity in solvent. The lighter the gas molecule, the faster it will diffuse. The weaker the gas−solvent interaction, the longer distance the gas molecule will travel before it collides with other solvent molecules. Both of these will lead to a higher gas diffusivity. The lightest H2 molecule interacts most weakly with PEGS (Table 1). As a result, H2 exhibits the largest selfdiffusivity in the PEGS solvent. Although H2O is lighter than both H2S and CO2, its interaction with PEGS is significantly larger (Table 1). This leads to smaller H2O diffusivity in PEGS compared with CO2 and H2S diffusivities. The H2S molecule is 1.3 times lighter than CO2, but it interacts with the PEGS molecule more strongly by a factor of ∼2 (Table 1). This compromise between the gas mass and gas−PEGS interaction G


Article


Table 5. Simulated Rotation Relaxation Time τrlx (ps) for Three Principal Axes of the PEGS Solvent Molecule at 298− 373 K and 1 bar along with the Corresponding Eigenvalues (λ) for the Principal Axesa

leads to comparable CO2 and H2S diffusivity coefficients (Figure 6). By assuming a simple Arrhenius behavior, the gas selfdiffusion activation energy could be estimated from the linear fitting of log(D) vs 1/T (Figure 6). The diffusion activation energy was calculated to be 12 ± 6 kJ/mol for H2, 12.7 ± 1.2 kJ/mol for CO2, 22.3 ± 1.0 kJ/mol for H2S, and 29.1 ± 0.7 kJ/ mol for H2O. Compared with H2 diffusion (fast diffusion), the H2O diffusivity (slow diffusion) increases more sharply at elevated temperatures. Finally, it is interesting to note that gas absorption may affect the solvent diffusivity in a different way. At 298−373 K, when H2O is absorbed in PEGS, the PEGS solvent diffusivity decreases by 10−40% compared with PEGS diffusion when H2 is absorbed (Table 4). Water diffuses more slowly than H2 in PEGS, which in turn slows down the PEGS solvent diffusion. 3.5. Viscosity. Simulations overestimate the PEGS solvent viscosity by 1.7−2.5 times compared with the experimental data at 298−333 K (Figure 7). This viscosity difference is partly due

T (K)

τrlx,1 (ps)

τrlx,2 (ps)

τrlx,3 (ps)

λ1

298.2 313.2 333.2 373.2

575 313 163 61

231 120 62 24

251 131 69 26

16.52 16.58 16.54 16.63

a

λ2 (4) (4) (3) (2)

3.88 3.90 3.91 3.90

(1) (1) (1) (1)

λ3 1.57 (1) 1.59 (1) 1.603 (4) 1.614 (3)

Uncertainties in the last digit are given in parentheses.

time. As the temperature is increased, the eigenvalues (λ) only change slightly, which implies that the shape for the PEGS molecule only slightly changes. In contrast, the τrlx values significantly decrease at elevated temperatures; the PEGS molecule rotates more freely, and as a result, the PEGS solvent viscosity decreases (Figure 7). 3.6. Heat Capacity. The simulated and experimental heat capacity values for neat PEGS solvent are comparable, with differences of ∼20% (Table 6 and Figure 8). As shown in Table 6, the ideal gas heat capacity component (CigP,m) predominates, and it contributes 65−75% to the total liquid heat capacity (CP,m). It is interesting to note that CigP,m and the residual liquid heat capacity (Cres P,m) exhibit opposite temperature effects (Table 6). When the temperature is increased from 298 to 373 K, the CigP,m value increases by 104 J/(mol.K). In contrast, the Cres P,m value decreases by 80 J/(mol.K). As a result, the total heat capacity (CP,m) values do not change significantly with temperature (Table 6 and Figure 8). 3.7. Surface Tension. Similar to our previous work,5 the forward and backward calculations give consistent surface tension values; the surface tension values also agree with each other among different sets of simulations (not shown here). These findings suggest that simulations are long enough to reliably obtain the surface tension value. The simulated surface tension values are close to the experimental data (Table 7 and Figure 9), with differences of 10−20%. Additionally, the simulation and experiment data give a similar surface tension− temperature relationship (Figure 9), that is, the surface tension will decrease when the temperature is increased. 3.8. CO2 Absorption Effects on CO2−PEGS Solution Properties. 3.8.1. Diffusivity Effects. Finally, we investigated CO2 absorption effects on solution properties. It was found that CO2 diffusivity in the solution is increased by a factor of ∼2 when CO2 loading (mole fraction) is increased from 0.05 to 0.706 (Table 8). At xCO2 = 0.706, the CO2−PEGS solution volume expands by 23.5% compared to the neat PEGS solvent volume. In contrast, the volume expansion at xCO2 = 0.05 is only 0.3%. This large volume expansion at xCO2 = 0.706 will decrease the solution viscosity (below), and as a result, CO2 diffusivity will increase. Similarly, the PEGS solvent molecule diffusivity also increases by a factor of ∼5 at xCO2 = 0.706 compared to the neat PEGS solvent diffusivity. 3.8.2. Viscosity Effects. Simulations show that the CO2− PEGS solution viscosity decreases when CO2 loading is increased (Table 8 and Figure 10). For example, at xCO2 = 0.706, simulation predicts that the solution viscosity decreases by 5.6 times compared to the neat PEGS solvent viscosity. In order to confirm simulations, in this work, we have also measured the CO2−PEGS solution viscosity at 298 K and

Figure 7. Simulated neat PEGS solvent viscosities at 298−373 K and 1 bar. For comparison, the experimental viscosity data obtained in this work are also shown. Note that the experimental PEGS solvent sample contains 387 ± 101 ppm of water. The lines are guides for the eye. Simulation error bars are also shown.

to the classical FF parameters for the PEGS molecule. Note that the experimental PEGS sample contains ∼400 ppm of water, which may affect the neat PEGS solvent viscosity. Despite this viscosity difference, the simulation and experimental data give similar viscosity−temperature relationships. To further understand viscosity, we have also calculated the principal axis rotational relaxation time for the PEGS molecule by using the same procedure as that used by Borodin et al.25 and Maginn.33 Specifically, we first calculated the correlation 1

2 function, C(t ) = 2 {3[ e (0) ⃗ · e ⃗(t )] − 1} , where e(⃗ t) indicates the unit principal axis vector in x, y, and z directions. The ⟨⟩ denotes an average over all PEGS molecules, all saved configurations, and all different time origins. The rotational relaxation time was then calculated to be τrlx = ∫ ∞ 0 C(t) dt. The τrlx can be obtained either from direct numerical integration or by fitting C(t) to the Kohlrausch−Williams−Watts (KWW) exponential expression25 followed by analytic integration.34 These two methods were found to give comparable τrlx values for the PEGS solvent. In this work, only the τrlx values obtained from the second method (KWW fitting) are reported. As shown in Table 5, the largest principal axis a1, characterized by the largest eigenvalue λ1, is well separated from the other two principal axes, a2 and a3; λ1 (∼16.5) is much larger than λ2 and λ3 (∼1.5−4). The largest principal axis exhibits the slowest rotation and the longest rotation relaxation

H


Article


Table 6. Simulated Ideal Gas Heat Capacity (CigP,m), Liquid Residual Heat Capacity (Cres P,m), and the Total Liquid Heat Capacity (CP,m) Values for PEGS Solvent at 298−373 K and 1 bar along with the Experimental Total Liquid Heat Capacity (Cexp P,m) at 298 K6 Shown for Comparisona

a

T (K)

CigP,m J/(mol·K)

Cres P,m J/(mol·K)

CP,m J/(mol·K)

CP,m J/(g·K)

Cexp P,m J/(g·K)

298.2 333.2 373.2

590.08 639.01 694.09

313 (48) 291 (100) 233 (17)

903 (48) 930 (100) 927 (17)

2.1 (1) 2.2 (2) 2.17 (4)

1.77 (8) − −

Uncertainties in the last digit are given in parentheses. The − symbol indicates that the corresponding experimental data are not available.

Note that the number of first-neighboring PEGS molecules around the central PEGS molecule is 14.8 in neat PEGS solvent at ambient condition; it decreases to 11.8 in solution at xCO2 = 0.706. This decreased number of neighboring PEGS molecules will lead to decreased PEGS−PEGS interaction and decreased solution viscosity. Additionally, the rotational relaxation time (τrlx,1) corresponding to the largest PEGS principal axis also decreases when CO2 loading is increased (Figure 10). The decreased rotation relaxation time at higher CO2 loading implies that the PEGS molecule rotates more freely, and as a result, the solution viscosity decreases. In order to further understand CO2 and volume effects on viscosity, one additional simulation was performed. The model system only contains PEGS molecules at 298 K but uses the same volume as that for the CO2−PEGS solution at xCO2 = 0.706, which is 23.5% larger compared to the neat PEGS solvent volume at ambient condition. For this model system, the viscosity was found to be 3.2 ± 0.2 × 10−3 Pa.s; τrlx,1 for the PEGS molecule was calculated to be 358 ps; the number of first-neighboring PEGS molecules around the central PEGS molecule was 13.1. Compared to neat PEGS at ambient conditions, this model system has a larger volume and a smaller number of neighboring PEGS molecules, which will lead to smaller PEGS−PEGS interactions. Hence, this model system exhibits 1.6 times smaller τrlx,1 and 3.1 times smaller viscosity. When compared to the CO2−PEGS solution at xCO2 = 0.706 (Table 8), this model system does not contain CO2 and exhibits a larger number of first-neighboring PEGS molecules, which suggests that CO2 in the solution may behave as a “lubricant” to decrease PEGS−PEGS interactions. Consequently, the model system exhibits 2.8 times larger τrlx,1 and 1.8 times larger viscosity compared to the solution at xCO2 = 0.706. To conclude, CO2 absorption in PEGS expands the solution volume and decreases the number of first-neighboring PEGS molecules, which will lead to decreased PEGS−PEGS interaction. As a result, PEGS molecules will more easily rotate with a smaller τrlx,1 and the viscosity will become smaller. 3.8.3. Surface Tension Effects. At 30 bar and xCO2 = 0.706, the CO2−PEGS solution exhibits 20% smaller surface tension compared to that of the neat PEGS solvent (Table 8). This finding is similar to our previous work.5 Upon CO2 absorption, the solution volume expands, which will decrease the PEGS− PEGS interaction. As a result, it takes less energy to bring a PEGS molecule from the bulk solution phase to the interface region and the solution surface tension becomes smaller.

Figure 8. Simulated and experimental6 heat capacity values for neat PEGS solvent at 298−373 K and 1 bar. The simulation and experiment error bars are also shown.

Table 7. Simulated Liquid−Vapor Surface Tensions (in 10−3 N/m) for Neat PEGS at 298−373 K along with the Surface Tension Values (γexp) Estimated from the Experimental Data (Figure 9) by Linear Interpolation for Comparisona simulated surface tension T (K)

γsim

γtail

γsim,total

γexp

298.2 333.2 373.2

13.2 (8) 10.6 (7) 9.2 (4)

6.52 5.64 4.52

19.7 (8) 16.2 (7) 13.7 (4)

21.84 19.79 −

a The total surface tension is calculated as γsim,total = γsim + γtail. Uncertainties in the last digit are given in parentheses. The − symbol indicate that the corresponding experimental viscosity at 373 K is unavailable.

Figure 9. Simulated and experimental6 surface tension values for PEGS solvent at 298−373 K. The simulation error bars are also shown. The lines are linear fits to the simulation and experimental data.

different CO2 loadings by using a stress-controlled Anton Paar Physica MCR 302 rheometer. Although the experimental solution viscosity data are smaller compared to those from simulations (Figure 10), they exhibit a similar decrease of solution viscosity upon CO2 loading. The experimental data indicate that the solution viscosity decreases by 4.8 times at xCO2 = 0.657 compared to the neat PEGS solvent viscosity.

4. CONCLUSIONS The thermophysical properties for neat PEGS solvent were predicted from atomistic simulations. Additionally, CO2, H2, H2O, and H2S solubilities and diffusivities in PEGS along with I


Article


Table 8. Simulated CO2 and PEGS Diffusivities in the CO2−PEGS Solution and the Solution Viscosity and Surface Tension at 298 K and Different CO2 Loading (Mole Fraction, xCO2) along with the Neat PEGS (xCO2 = 0) Solvent Diffusivity, Viscosity, And Surface Tension at 298 K and 1 bar for Comparisona xCO2 0 0.05 0.123 0.444 0.706

P (bar) 1 1 2.5 14 30

DCO2 (m2/s) 1.75 1.24 2.04 3.39

(3) (1) (1) (1)

× × × ×

10−9 10−9 10−9 10−9

DPEGS (m2/s) 5.55 5.1 6.1 1.14 2.83

(3) (1) (1) (1) (3)

× × × × ×

−11

10 10−11 10−11 10−10 10−10

surface tension (10−3 N·m)

viscosity (10−3 Pa·s)

19.7 (8) − − − 15.8 (6)

9.8 (8) − 7.9 (4) 2.5 (2) 1.74 (2)

a At each CO2 loading, the corresponding CO2 pressure in the gas phase (Table 2) is also added. Uncertainties in the last digit are given in parentheses. The − symbols indicate that the simulation values were not obtained.

and CO2 gas solubilities decrease, while the H2 solubility slightly increases mainly due to the increased free molar volume for PEGS solvent at elevated temperature. Note that the H2− PEGS interaction is very weak and the H2 solubility is mainly determined by the solvent free molar volume. Gas diffusivity in PEGS decreases in an opposite order compared to gas solubility: H2 > CO2 ≈ H2S > H2O. Both the gas mass and gas−PEGS interaction are important for determining the gas diffusivity. The lightest H2 molecule interacts most weakly with PEGS, which leads to the largest H2 diffusivity. Although H2O has less mass compared to CO2 and H2S, its interaction with PEGS is much stronger. As a result, H2O exhibits the smallest diffusivity. Compared with PDMS, PEGS interacts more strongly with H2O by 12−20 kJ/mol. Additionally, PEGS has fewer branched −CH3 groups, which reduces the ability for the PEGS molecule to form hydrophobic pockets. These two factors make PEGS more hydrophilic compared to PDMS. The CO2/H2 solubility selectivity in PEGS at 298 K is 1.7 times larger than that in PDMS. This result is due to two primary reasons: the stronger CO2−PEGS interaction compared with the CO2−PDMS interaction and the smaller free molar volume and smaller free volume fraction for PEGS compared with those for PDMS. Finally, we investigated CO2 absorption effects on solution properties. At a CO2 loading of 0.706 (mole fraction), simulations show that the CO2−PEGS solution viscosity decreases by a factor of 5.6, the solution surface tension decreases by 20%, and the PEGS diffusivity increases by 5.1 times compared to the corresponding neat PEGS solvent properties. Experimental results obtained in this work also show a similar solution viscosity decrease upon CO2 absorption. At xCO2 = 0.706, the CO2−PEGS solution volume expands by 23.5% compared to neat PEGS solvent volume, which leads to fewer first-neighboring PEGS molecules around the central PEGS molecule. As a result, the PEGS−PEGS interaction decreases. Additionally, it was found that CO2 molecules could act as a lubricant to decrease the PEGS−PEGS interaction by decreasing the number of first-neighboring PEGS molecules. The above two CO2 absorption effects (the solution volume expansion and CO2 lubricant) will decrease the solution viscosity and surface tension and increase PEGS and CO2 diffusivity values at xCO2 = 0.706.

Figure 10. Simulated rotational relaxation time (squares), τrlx,1, for the largest PEGS principal axis and the solution viscosity (open triangles) at 298 K and different CO2 loadings (Table 8). The experimental solution viscosity values (filled triangles) versus CO2 loadings at 298 K obtained in this work are also shown for comparison. The lines are guides for the eye.

CO2 absorption effects on the solution properties were investigated. The simulated neat PEGS solvent density, heat capacity, surface tension, and viscosity values at 298−373 K and 1 bar are comparable to the experimental data, with differences of 0.8% for density, 20% for heat capacity, 10−20% for surface tension, and 170−250% for viscosity. The simulated CO2 and H2 solubilities in PEGS at 298 K and 1−30 bar are also close to the experimental data, with differences of 10% for CO2 and 30% for H2 solubilities. Note that all simulation results are purely predictive; only one experimental density for neat PEGS solvent at ambient conditions was used to slightly tune the LJ σ parameters for each atom of the PEGS molecule. Except for viscosity, reasonable agreement in thermophysical properties between simulations and experiments indicates that the PEGS classical FF parameters used in this work are transferrable to predict many different properties. For neat PEGS solvent, simulations show that the rotation relaxation time for the largest principal axis (τrlx,1) of the PEGS molecule correlates very well with solvent viscosity. When the temperature is increased from 298 to 373 K, the PEGS viscosity decreases by a factor of 8.3 and the τrlx,1 decreases by a factor of 9.4. The smaller τrlx,1 value at elevated temperature suggests that PEGS molecules rotate more freely, which in turn leads to smaller PEGS viscosity when the temperature is increased. For the neat PEGS solvent heat capacity, simulations indicate that the ideal gas heat capacity predominates over the residual liquid heat capacity. Interestingly, when the temperature is increased, the ideal gas heat capacity increases, while the residual liquid heat capacity decreases. Calculations show that the gas−PEGS interaction decreases in the following order: H2O > H2S > CO2 > H2. Gas solubilities in PEGS follow the same order as those for the gas−PEGS interaction. When the temperature is increased, the H2O, H2S,

■

ASSOCIATED CONTENT

* Supporting Information S

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.6b06810. Classical force field parameters for the PEGS molecule, validations of the procedure to calculate viscosity and J


Article


■

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. (10) 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. (11) Fischer, J.; Paschek, D.; Geiger, A.; Sadowski, G. Modeling of Aqueous Poly(oxyethylene) Solutions: 1. Atomistic Simulations. J. Phys. Chem. B 2008, 112, 2388−2398. (12) Vanommeslaeghe, K.; Hatcher, E.; Acharya, C.; Kundu, S.; Zhong, S.; Shim, J.; Darian, E.; Guvench, O.; Lopes, P.; Vorobyov, I.; et al. CHARMM General Force Field: A Force Field for Drug-Like Molecules Compatible with the CHARMM All-Atom Additive Biological Force Fields. J. Comput. Chem. 2009, 31, 671−690. (13) Frischknecht, A.; Curro, J. Improved United Atom Force Field for Poly(dimethylsiloxane). Macromolecules 2003, 36, 2122−2129. (14) Smith, J.; Borodin, O.; Smith, G. A Quantum Chemistry Based Force Field for Poly(dimethylsiloxane). J. Phys. Chem. B 2004, 108, 20340−20350. (15) Striolo, A.; McCabe, C.; Cummings, P. T.; Chan, E. R.; Glotzer, S. C. Aggregation of POSS Monomers in Liquid Hexane: A MolecularSimulation Study. J. Phys. Chem. B 2007, 111, 12248−12256. (16) Breneman, C. M.; Wiberg, K. B. Determining Atom-centered Monopoles from Molecular Electrostatic Potentials. The Need for High Sampling Density in Formamide Conformational Analysis. J. Comput. Chem. 1990, 11, 361−373. (17) 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. (18) 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. (19) 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. (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. (21) 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. (22) 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. (23) Shi, W.; Sorescu, D. C. Molecular Simulations of CO2 and H2 Sorption into Ionic Liquid 1-n-Hexyl-3- methylimidazolium Bis(trifluoromethylsulfonyl)amide ([hmim][Tf2N]) Confined in Carbon Nanotubes. J. Phys. Chem. B 2010, 114, 15029−15041. (24) Mondello, M.; Grest, G. Viscosity Calculations of n-Alkanes by Equilibrium Molecular Dynamics. J. Chem. Phys. 1997, 106, 9327− 9336. (25) Borodin, O.; Smith, G. D. Structure and Dynamics of N-methylN-propylpyrrolidinium bis(trifluoromethanesulfonyl)imide Ionic Liquid from Molecular Dynamics Simulations. J. Phys. Chem. B 2006, 110, 11481−11490. (26) Raabe, G.; Sadus, R. J. Molecular Dynamics Simulation of the Effect of Bond Flexibility on the Transport Properties of Water. J. Chem. Phys. 2012, 137, 104512. (27) Cadena, C.; Maginn, E. J. Molecular Simulation Study of Some Thermophysical and Transport Properties of Triazolium-based Ionic Liquids. J. Phys. Chem. B 2006, 110, 18026−18039. (28) Rowlinson, J. S.; Widom, B. Molecular Theory of Capillarity; Clarendon: Oxford, U.K., 1982. (29) Kirkwood, J.; Buff, F. The Statistical Mechanical Theory Of Surface Tension. J. Chem. Phys. 1949, 17, 338−343.

heat capacity, and the ab initio gas−PEGS interactions at the complete basis set limit obtained from the RI-MP2 calculations (PDF)

AUTHOR INFORMATION

Corresponding Author

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

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This technical effort was performed in support of the National Energy Technology Laboratory’s ongoing research in computational chemistry under RES Contract DE-FE0004000. The authors also thank Megan K. Macala for her efforts to dry the neat PEGS solvent sample. 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 AECOM. Neither the United States Government nor any agency thereof, nor any of their employees, nor AECOM, 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 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 authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof.

■

REFERENCES

(1) Sustainable Innovation Forum 2015: COP21 Paris. http://www. cop21paris.org/ (accessed May 16, 2016). (2) Clean Power Plan Home. https://www.epa.gov/cleanpowerplan (Accessed May 16, 2016). (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) Koronaios, P.; Stevenson, C.; Warman, S.; Enick, R.; Luebke, D. Thermally Stable Silicone Solvents for the Selective Absorption of CO2 from Warm Gas Streams That Also Contain H2 and H2O. Energy Fuels 2016, 30, 5901−5910. (5) Shi, W.; Siefert, N. S.; Morreale, B. D. Molecular Simulations of CO2, H2, H2O, and H2S Gas Absorption into Hydrophobic Poly(dimethylsiloxane) (PDMS) Solvent: Solubility and Surface Tension. J. Phys. Chem. C 2015, 119, 19253−19265. (6) Siefert, N.; Agarwal, S.; Shi, F.; Shi, W.; Roth, E.; Hopkinson, D.; Kusuma, V.; Thompson, R.; Luebke, D.; Nulwala, H. Hydrophobic Physical Solvents for Precombustion CO2 Capture: Experiments, Computatuonal Simulations, and Techno-economic Analysis. Int. J. Greenhouse Gas Control 2016, 49, 364−371. (7) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon: Oxford, U.K., 1987. (8) 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. (9) Shi, W.; Sorescu, D. C.; Luebke, D. R.; Keller, M. J.; Wickramanayake, S. Molecular Simulations and Experimental Studies K


Article

The Journal of Physical Chemistry C (30) Irving, J.; Kirkwood, J. The Statistical Mechanical Theory Of Transport Processes 0.4. The Equations Of Hydrodynamics. J. Chem. Phys. 1950, 18, 817−829. (31) Walton, J.; Tildesley, D.; Rowlinson, J.; Henderson, J. The Pressure Tensor at The Planar Surface of A Liquid. Mol. Phys. 1983, 48, 1357−1368. (32) Gloor, G.; Jackson, G.; Blas, F.; de Miguel, E. 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. (33) Cadena, C.; Zhao, Q.; Snurr, R. Q.; Maginn, E. J. Molecular Modeling and Experimental Studies of the Thermodynamic and Transport Properties of Pyridinium-based Ionic Liquids. J. Phys. Chem. B 2006, 110, 2821−2832. (34) Liu, H.; Maginn, E. A Molecular Dynamics Investigation of the Structural and Dynamic Properties of the Ionic Liquid 1-n-butyl-3methylimidazolium bis(trifluoromethanesulfonyl)imide. J. Chem. Phys. 2011, 135, 124507.

L


Molecular Simulations of the Thermophysical ... - ACS Publications

Recommend Documents