Molecular Simulations of CO2 and H2 Solubility, CO2 Diffusivity, and

Article Cite This: J. Chem. Eng. Data XXXX, XXX, XXX−XXX

pubs.acs.org/jced

Molecular Simulations of CO2 and H2 Solubility, CO2 Diffusivity, and Solvent Viscosity at 298 K for 27 Commercially Available Physical Solvents Wei Shi,*,†,‡ Robert L. Thompson,†,‡ Megan K. Macala,†,‡ Kevin Resnik,†,‡ Janice A. Steckel,† Nicholas S. Siefert,† 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

‡

J. Chem. Eng. Data Downloaded from pubs.acs.org by LUND UNIV on 03/05/19. For personal use only.

S Supporting Information *

ABSTRACT: CO2 and H2 solubilities, CO2/H2 solubility selectivities, CO2 diffusivities, and solvent viscosities in 27 commercially available physical solvents at 298 K were calculated from molecular simulations using the CHARMM36 all-atom force field for most solvents, and the simulation results were compared with available experimental data. The van der Waals radius parameters for solvents were slightly tuned to reproduce the experimental solvent density. The simulated CO2 solubilities are comparable with the experimental data, with an average absolute difference of 28%. For the homologous compounds containing the −(OCH2CH2)− repeat unit, both simulated and experimental data show that CO2 solubility decreases when the number of repeat units is increased; CO2 solubilities in these homologous compounds exhibit almost a perfect positive linear correlation with the solvent free-volume fractions. The simulated H2 solubilities and CO2/H2 solubility selectivities are also comparable with the experimental data, with differences of 22% and 17%, respectively. The H2 solubilities in all solvents studied in this work correlate very well with the solvent freevolume fractions, exhibiting a positive linear correlation coefficient of 0.84. Additionally, simulations show that CO2 solubility decreases when the temperature is increased. In contrast, H2 solubility increases at elevated temperature, which is partly due to the increased solvent free-volume fraction at elevated temperature. Finally, although the viscosity difference tends to be large (30%−246%) between simulation and experiment, both simulated and experimental data exhibit a similar solvent viscosity trend. Furthermore, simulations show that CO2 diffusivities in solvents are very strongly correlated with the solvent viscosities ± 0.03 . and the relationship between them is given by DCO2 = (2.6 ± 0.3) × 10−9/η0.59 solvent

1. INTRODUCTION Capturing CO2 from coal power plants has environmental benefits, and potential economic benefits in locations when CO2 is needed for industrial applications. Postcombustion CO2 capture is typically done using chemical solvents because of the low CO2 partial pressure (∼0.1 bar) in a nitrogen-rich flue gas. At integrated gasification combined cycle (IGCC) power plants, CO2 has significantly higher partial pressure (∼25 bar) in a hydrogen-rich syngas, and CO2 can be captured by using physical solvents in the IGCC power plant. Recently, we have developed an approach that integrates database, molecular modeling, and machine learning methods to screen physical solvents for gas separation and CO 2 precombustion capture applications.1,2 In the integrated approach, molecular simulation was used to calculate important physical properties, such as CO2 and H2 solubilities and selectivities in solvents, CO2 diffusivities, solvent viscosities, and solvent surface tensions. Furthermore, molecular simulations © XXXX American Chemical Society

will be applied to predict physical properties for new physical solvents or for properties which are unavailable from experiments. These properties will be used in Advanced System for Process Engineering (ASPEN) for techno-economic analysis. Hence, it is very important to evaluate the accuracy for the above simulated physical properties. In this work, we summarized the molecular simulation results for four crucial properties of 27 commercially available physical solvents, that is, CO2 solubility, CO2/H2 solubility selectivity, CO2 diffusivity, and solvent viscosity. These four properties will affect the levelized cost of precombustion carbon capture at an IGCC power plant. A process flow diagram of a 3-flash Special Issue: FOMMS Received: December 19, 2018 Accepted: February 19, 2019

A

DOI: 10.1021/acs.jced.8b01228 J. Chem. Eng. Data XXXX, XXX, XXX−XXX

Journal of Chemical & Engineering Data

Article

Table 1. Chemical Name, Chemical Formula, CAS Registry Number (CASRN), and the Schematic Molecular Structure for 27 Solvents Studied in This Work

B



Article

regeneration process (Figure 5 in our previous work3) shows that a solvent with high CO2 uptake between 2 and 25 bar can reduce the solvent flow rate and lower the pumping cost. A solvent with high CO2/H2 solubility selectivity would result in the high pressure flash tank being operated with minimal need for recycle, and a solvent with high CO2 diffusivity and low viscosity would lead to a small and less costly absorber. Note that for fair comparison between different solvents, we will calculate the four crucial properties at the same temperature of 298 K. The best solvents will be selected and the ASPEN process modeling will be used to optimize the operating conditions (such as solvent flow rate and the absorption temperature) for each solvent. For most solvent molecules, the CHARMM36 all-atom force field was used in this work.4 Only the van der Waals radius parameters for solvent molecules were tuned to match the solvent densities at ambient conditions. The force field parameters were then fixed to predict all other properties, such as CO2 and H2 solubility, CO2 diffusivity, and solvent viscosity. The simulation results were compared with available experimental data to evaluate the simulation accuracy. Additionally, it was found that H2 solubilities in all solvents studied in this work exhibit a strong positive linear correlation with solvent free-volume fractions; CO2 solubilities in the homologous compounds almost exhibit a perfect positive linear correlation with the solvent free-volume fractions, and CO2 diffusivity is very strongly correlated with the solvent viscosity.

work contain the N atom, such as 2-benzylpyridine, 2,2methylene bispyridine, and n-methyl-2-pyrrolidone, both our quantum mechanical calculations and the experimental data show that they do not chemically react with CO2; they still physically absorb CO2. 2.2. Classical Force Field (FF). The classical FF potential used to simulate CO2, H2, the 27 solvents, and their interactions is given by =(r) =

∑

kr(r − r0)2 +

bonds

+

∑

kθ(θ − θ0)2

angles

∑

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

dihedrals

+

∑

k ψ (ψ − ψ0)2

ÄÅ É l ÅÅi y12 i y6ÑÑÑ q q | o o o o Å σ σ o o j z j zÑ ∑ omoo4ϵijÅÅÅÅÅjjjjj ij zzzzz − jjjjj ij zzzzz ÑÑÑÑÑ + i j o}oo r Ñ rij o o ÅÅÅk rij { o j=i+1 o k ij { ÑÑÑÖ n ÅÇ ~

impropers N−1

+

∑ i=1

N

(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 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 H2O,8 and H29 were obtained from our previous work. For solvent molecules, the r0 and θ0 values were obtained from quantum mechanical (QM) gas-phase geometry optimization calculations at the B3LYP/6311++G(d,p) level of theory by using the Gaussian 09 program.10 Within a single optimized molecule, if the same type of bond has different r0 values, the average for r0 will be used for that bond. Similarly, the average for θ0 will be used for an angle if different θ0 values occur for the same type of angle. The parameters for kr, kθ, dihedral potential energy, and LJ interaction were typically obtained from the CHARMM par−all36−cgenff.prm file.4 If the parameters for kr, kθ, and dihedral were unavailable in the par−all36−cgenff.prm, we manually assigned these values by choosing parameters for similar bonds, angles, or dihedrals. Following geometry optimization calculations, atomic charges for the solvent molecules were obtained from QM calculations using the CHELPG protocol11 at the B3LYP/6-311++G(d,p) level of theory. The atomic charges for the topologically equivalent atoms12 were set to be equal by averaging the charges. Note that some FF parameters for tributyl phosphate and 2,5,8,11,14pentaoxapentadecane compounds were obtained from the previous work.13,14 Additionally, the homologous compounds, such as 1-methoxy-2-(2-methoxyethoxy)ethane, 1,2-bis(2methoxyethoxy)ethane, 2,5,8,11,14-pentaoxapentadecane, and 2,5,8,11,14,17- hexaoxaoctadecane have the same FF parameters except the atomic charges. In-house codes were developed to automatically set up all the above FF parameters except that atom types were manually assigned. Additionally, the psf structure and crd coordinate files, which are needed for any molecular dynamics (MD) simulation, were also generated using the in-house codes. For example, by starting from the charge calculation output file obtained from QM calculation on a single solvent molecule, the in-house codes

2. SIMULATION DETAILS 2.1. Chemical System Information. The chemical name, chemical formula, CAS registry number, and the schematic molecular structure for all 27 solvents studied in this work are given in Table 1. The methyl 4-methylbenzenesulfonate compound has a melting point between 298 and 301 K.5 This compound may be a solid at ambient conditions (298 K and 1 bar). Simulations for this compound were conducted by using irregular amorphous configurations, which implies that this compound is in the supercooled liquid state at ambient conditions. The above 27 physical solvents were not chosen randomly. For example, some solvents were obtained from our in-house database calculations,1,2 and they were found to interact strongly with CO2, such as 2-benzylpyridine. Methanol, propylene carbonate, and n-methyl-2-pyrrolidone have been studied to physically absorb CO2,6 and as a result they were included in this study. Similarly, the Selexol solvent, a mixture of dimethyl ethers of polyethylene glycol, has also been used for CO2 absorption.6 Consequently, the homologous compounds 1-methoxy-2-(2methoxyethoxy)ethane, 1,2-bis(2-methoxyethoxy)ethane, 2,5,8,11,14-pentaoxapentadecane, and 2,5,8,11,14,17-hexaoxaoctadecane, which resemble the Selexol solvent, were also included in this work. Water and decane were included for comparison. Note that many solvents studied in this work are not commonly used for CO2 capture. For example, water and decane have small CO2 solubilities and they are not good for CO2 capture. However, they are still included in this work because we want to compare the simulated CO2 solubility with the experimental data for these two solvents to evaluate the accuracy of molecular simulations to predict CO2 solubility. In addition, some other solvents, such as sulfolane, which is also a commonly used physical solvent for CO2 capture, were not included in the current work. They will be included in our future work. Note that even though several compounds studied in this C



Article

can generate the crd file, the atom, bond, angle, dihedral, aromatic improper dihedral, and the atomic charges contained in the psf file. These in-house codes have been validated for many test molecules; the in-house codes were found to give the same FF parameters as would be obtained manually, and the inhouse codes also generate the same psf and crd files as the VEGA ZZ software.15 For each solvent, the atomic radii (σ in eq 1) was slightly tuned such that the simulated density reproduced the experimental solvent density at 298 K and 1 bar. In-house codes were developed to tune σ values automatically instead of manually, and the codes are provided in the sigma_tune_example directory (Supporting Information). Additionally, for each solvent, five files are also given in the run_sim_file directory (Supporting Information). The five files are the FF parameter file, the psf and crd files for one solvent molecule, and psf and PDB files for the system, which contains 100−1024 solvent molecules in a cubic box with box sizes of 30−47 Å. The simulation box sizes are larger than twice the LJ cutoff value (12.0 Å). To be more rigorous, one should systematically study the system size effects on the properties studied in this work, such as gas solubility, solvent viscosity, and gas diffusivity. For example, it has been shown that the diffusivity obtained by using a larger simulation box is larger than that obtained by using a smaller simulation box.16,17 However, the studies of system size effects are beyond the scope of this paper and they were not performed. We will study these system size effects in another paper. Isothermal−isobaric (NPT) MD simulations were used to calculate solvent density and obtain trajectories from which gas Henry’s law constant (see below) calculations were performed. NPT MD simulations were performed for 40 ns, and 40 000 simulation snap shots were saved every 1 ps. The simulation time step was set to be 1 fs. For each solvent, the system contains 100−1024 solvent molecules. Langevin dynamics with a Nosé− Hoover barostat/thermostat were applied for temperature and pressure controls. The barostat oscillation time and damping factors were set to be 1 ps, and the thermostat damping factor was set to be 5 ps−1. A typical configuration file (NPT.conf) containing the NPT simulation details may be found in the sigma_tune_example directory (Supporting Information). To calculate transport properties such as gas diffusivity and solvent viscosity, typically 20 ns microcanonical (NVE) MD simulations were performed subsequent to 20 ns NPT and 20 ns canonical (NVT) MD simulations. The simulation time step in NVE simulation was set to be 0.5 fs, and 1 fs in both NPT and NVT simulations. The volume V in NVE and NVT was set to be the average volume obtained from NPT MD simulations. In each NVE MD simulation, 2000 minimization steps were initially performed, followed by heating the system with a rate of 0.01 K per step from 0 to 298 K; the equilibration and final production were typically set to be 2 and 20 ns, respectively. In the final production runs, coordinates and pressure tensors were dumped every 1000 and 1 fs, respectively, for later gas diffusivity and solvent viscosity calculations (see below). For each solvent, the simulation system contains five CO2 molecules and 100−1024 solvent molecules for CO2 diffusivity and solvent viscosity calculations. For all the above MD simulations, the smooth particle mesh Ewald method18 was used to calculate the electrostatic interaction. The MD simulations were performed using the NAMD program.19 Further simulation details can be found in our previous work.9

2.3. Henry’s Law Constant and Gas Solubility. The gas Henry’s law constant is defined as Hgas =

lim

Pgas

Pgas , xgas → 0 xgas ngas

denoted as a mole fraction, that is, xgas =

ngas + nsolv

. If xgas is (n is the

number of moles), H will be given by the following equation20 H=

ρmolar,solv k b × T exp( −βμex )

where ρmolar,solv =

(2)

nsolv , k is the Boltzmann constant, β = 1/(kb × Vsolv b ex

T), T is the temperature. The μ is given by i exp( −βμex ) = jjj k

∫ exp(−βUsolv,gas) × r12 sin θ1

y ··· rK2 − 1 sin θK − 1 drsolv dr1 dθ1 dϕ1 ··· drK − 1 dθK − 1 dϕK − 1zzz { /

( ∫ exp(−βU

solv )

)

drsolv × Cintra,gas

(3)

where K is the number of atoms for one gas molecule, for example, K is 3 for CO2. The rsolv denotes the x, y, and z coordinates for all the solvent molecules in the system. The Usolv,gas is the total potential energy for the system containing nsolv solvent molecules and one gas molecule. The Usolv is the potential energy for the solvent system containing nsolv solvent molecules. The Cintra is given by Cintra =

∫ exp(−βUgas,intra) × r12 sin θ1 ··· rK2 − 1 sin θK − 1 dr1 dθ1 dϕ1 ··· drK − 1 dθK − 1 dϕK − 1 (4)

where Ugas,intra denotes the potential energy for a single gas molecule, which includes the bond, angle, dihedral, improper energies, and the interaction between different atoms of the same gas molecule separated by more than three consecutive bonds. For CO2, Ugas,intra only consists of two bonds and one angle energy terms. More details were given in our previous work.20 mgas If gas mass fraction is used, that is, xgas = m + m (m is the gas

solv

mass), the Henry’s law constant is given by Hm =

ρmass,solv k b × T Mgasexp( −βμex )

where ρmass,solv =

msolv , Vsolv

(5)

Mgas is the molecular mass for the gas

molecule. If volume-based gas solubility is used, that is, xgas =

ngas Vsolv

, the

Henry’s law constant is given by Hv =

kb × T exp( −βμex )

(6)

It can be concluded that the gas Henry’s law constant depends on solvent density (eqs 2 and 5) if either gas mole fraction or gas mass fraction is used. As a result, it is inaccurate to compare gas solubilities in different solvents using either H or Hm, both of which have units of bar. In contrast, when the volume-based gas solubility is used, the Henry’s law constant Hv only depends on μex (eq 6), which in turn only depends on gas interaction with D



Article

solvent molecules. Consequently, it is appropriate to use Hv for accurate comparison of gas solubilities in different solvents. Hv has units of bar × cm3/mol, or other similar units, such as 106Pa(MPa) × L/mol. In this work, Hv was computed and gas solubility was calculated to be 1/Hv, which indicates how many moles of gas are absorbed in one cm3 solvent at one bar of the gas pressure. To facilitate gas solubility comparison between simulation and experiment, the experimental Henry’s law constant H (mole fraction-based) was converted to H v (volume-based) by multiplying the experimental H and the experimental solvent molar volume. When calculating the Henry’s law constant from the snapshots obtained from NPT MD simulations by using the modified Widom test particle insertion method,20,21 2 × 106 trial insertions were attempted for each snapshot. To improve the insertion efficiency, the cubelet21 and bias-removing methods20 were used. Further calculation details were reported in our previous work.20 This Henry’s law constant calculation method has been implemented in our in-house codes, which14,20 have been shown to give the same gas Henry’s law constant as obtained from the continuous fractional component (CFC) Monte Carlo (MC) method. Additionally, in this work, CO2 solubilities in 2-benzylpyridine were calculated at 298 K by using the CFC MC method8,22 and they were found to be 0.022 ± 0.001 at 1.25 bar and 0.045 ± 0.002 at 2.5 bar, respectively. By linear fitting of CO2 solubilities to CO2 pressures (including the origin), CO2 Henry’s law constant was estimated to be 56 ± 1 bar and it was found to be consistent with the value of 58.2 ± 0.4 bar obtained from the test particle insertion method, with a small difference of 4%. Similarly, H2 solubilities in 2-benzylpyridine at 298 K and 10−40 bar were also calculated from the CFC MC method and they were found to be 0.0036 ± 0.0001 at 10 bar, 0.0070 ± 0.0002 at 20 bar, and 0.0152 ± 0.0005 at 40 bar, respectively. By linear fitting, the H2 Henry’s law constant was estimated to be 2660 ± 76 bar, and it was found to be consistent with the value of 2821 ± 8 bar obtained from the test particle insertion method, with a difference of 6%. Although the modified Widom test particle insertion and the CFC MC methods give consistent gas Henry’s law constant, we chose to use the modified Widom test particle insertion method to calculate gas Henry’s law constant. In practice, it was convenient to split the snapshots saved during MD simulations among different computing nodes to calculate the gas Henry’s law constant in a parallel way when using the modified Widom test particle insertion method. As a result, it is much faster to obtain the gas Henry’s law constant by using the modified Widom test particle insertion method. Although we have found good agreement in Henry’s law constant between the modified Widom test particle insertion and the CFC MC methods, and H2 solubility agreement in the [C6mim][Tf2N] ionic liquid between the CFC MC9 and the configurational-biased Gibbs ensemble methods,23 it is interesting to conduct more extensive gas solubility simulation comparisons between different groups, which will be performed in our other papers. In the test particle insertion method, the cubelet length was set to be ∼0.8 Å. A cubelet was marked as occupied if the distance between the cubelet center and any solvent atom is less than rev = s × (σi/2 + rprobe), where σi is LJ parameter for the ith atom of the solvent molecule. The variable s is a scaling factor and it was set to be 0.8. The radius for the probe atom (rprobe) was set to be 1.0 Å for CO2 insertion. For H2 insertion, rprobe was set to be 0.9 Å. The Fennel and Gezelter shift force (FGSF)8,24 method was used to compute the electrostatic interaction

between gas and solvent molecules to speed up calculations when calculating the Henry’s law constant. It has been shown that the FGSF method gives CO2−ionic liquid electrostatic interaction energies very close to the values obtained from the standard Ewald method.20 To correlate gas solubility with solvent free-volume fraction, the solvent free-volume fraction was also calculated by using the in-house code.20 The freevolume fraction was calculated to be the number of unoccupied cubelets divided by the total (occupied + unoccupied) number of cubelets. 2.4. CO2 Diffusivity. CO2 self-diffusivity was calculated by using the Einstein relation. The mean-squared displacement for the center of mass of each CO2 molecule was obtained from the NVE MD simulation. Further calculation details can be found in our previous work.9 2.5. Solvent Viscosity. Solvent viscosity was calculated from the NVE MD simulation by using the Einstein relation. Pressure tensors saved in NVE MD simulations were used to calculate viscosity. Similar to the CO2 diffusivity calculation, a βviscosity value was first calculated to identify the time period in which the βviscosity average is close to 1. Viscosity was then calculated in the same time period. Further calculation details were reported in our previous work.14,25 2.6. Simulated Uncertainty Quantification. In this work, we used the standard block average method to estimate the simulated uncertainty.26 Ten blocks were used to calculate the standard deviation of the block averages (s(xb)).26 The simulated error bars and uncertainties reported in this work were set to be s(xb), which is larger than the standard deviation of the mean for block averages (s(x̅b))26 by a factor of 10.0 = 3.162. By using the t distribution,27,28 the s(xb) value corresponds to a confidence level of 98−99%, that is, we are 98− 99% certain that the simulated average properties deviate from the means of the block averages by less than s(xb). Note that in this work we have not performed the rigorous autocorrelation analysis among the block averages. We only investigated the correlation between adjacent block averages for several systems. The block averages were found to fluctuate about mean values. For density and gas solubility, the adjacent block averages were found to exhibit weak linear correlation with correlation coefficients between 0.18−0.4. For gas diffusivity and solvent viscosity, the adjacent block averages exhibit moderate linear correlation with correlation coefficients of about 0.5. These calculations suggest that the block averages should not exhibit strong autocorrelation.

3. RESULTS AND DISCUSSION 3.1. CO2 Solubility. The simulated CO2 solubilities in 27 solvents at 298 K are given in Figure 1 and the corresponding numerical data are shown in Table S1 (Supporting Information). The simulated CO2 solubilities are between 0.152 and 2.37 mol/MPa·L. For 12 solvents for which experimental data are available, a comparison between simulated and measured CO2 solubility is shown (Figure 2 and Table S1). The average absolute CO2 solubility difference between simulation and experiment is 28%, with minimum and maximum differences of 0.6% and 50%, respectively. Considering that there is a variation among reported experimental CO2 solubilities, such as 15% for CO2 in methanol,29 the CO2 solubility difference between simulation and experiment is acceptable. This reasonable CO2 solubility comparison is partly due to a favorable CO2−solvent interaction E



Article

CO2 solubility in water (0.152 mol/(MPa·L)) is only for physical absorption; it is about half of the experimental value of 0.341 mol/(MPa·L) (Supporting Information, Table S1), which consists of both physical and chemical absorption. To evaluate the solvent FF parameter effects on CO2 solubility, we have used two different sets of FF parameters for each of the 2,5,8,11,14-pentaoxapentadecane, 2,4,6,8-tetraoxanonane, and 2-(2-methoxyethoxy)ethanol 1,1′,1″-phosphate solvents. The only difference between these two sets of FF parameters are their LJ parameters. For 2,5,8,11,14-pentaoxapentadecane, the two FF sets with both different σ and different ϵ values give solvent densities very close to each other and the simulated densities are in agreement with the experimental data. However, the second set gives 16% lower CO2 solubility than the first set does. For 2,4,6,8-tetraoxanonane, the second FF set has the same ϵ but 4.5% higher σ values compared with the first set and it gives 6% lower solvent density. The second FF set also gives 19% lower CO2 solubility compared with the first set. For 2-(2-methoxyethoxy)ethanol 1,1′,1″-phosphate, the second FF set has the same ϵ values but different σ values compared with the first set. The second FF set has larger σ values for some atoms but smaller σ values for other atoms compared with the first set. Both FF sets give solvent densities in good agreement with the experimental data, while the second set gives 20% larger CO2 solubility compared with the first set. All these findings suggest that the LJ parameters for solvent molecules will affect CO2 solubility by less than 20%. The 1-methoxy-2-(2-methoxyethoxy)ethane, 1,2-bis(2methoxyethoxy)ethane, 2,5,8,11,14-pentaoxapentadecane, and 2,5,8,11,14,17-hexaoxaoctadecane are homologous (Table 1). Both simulation and experiment (Figures 1 and 2, Table S1) show that when the number of −(OCH2CH2)− repeat units increases, the CO2 solubility decreases. Note that the same FF parameters were used for all those four solvents except the atomic charges. When the number of repeat units increases, the solvent density increases (Table S1) and the solvent free-volume fraction decreases (Figure 3), which in turn leads to a smaller CO2 solubility (Figure 3). For these homologous solvents, CO2 solubilities in them exhibit almost perfect positive linear correlation with the solvent free-volume fractions (Figure 3), with a correlation coefficient of 0.999. This perfect correlation is partly due to the similar CO2 interactions with the homologous

Figure 1. Simulated CO2 solubility values in 27 solvents (Table 1) at 298 K obtained from Henry’s law constant calculations. Also shown are error bars obtained from standard block average calculations.

Figure 2. Comparison of CO2 solubility values in 12 solvents between simulations and the experimental data. The experimental CO2 solubility data for CASRN 124-18-5, 126-73-8, 67-56-1, 108-32-7, 107-12-0, 10999-9 were obtained from the book compilation.29 For CASRN 143-248, 112-49-2, 111-96-6, 872-50-4, the experimental CO2 solubilities were obtained from the NIST collection.30 For CAS 101-82-6 and 1132-37-2 solvents, the experimental CO2 solubilities were estimated by linearly fitting CO2 solubilities obtained in this work versus pressures below 5 bar (including the origin). The simulation error bars obtained from standard block average calculations are also shown.

comparison between classical FF and QM calculations. By using the same procedure as described in our previous work,14 the classical FF was found to give reasonable CO2−solvent interactions compared with the QM calculations. For all solvents studied in this work, the absolute difference for CO2−solvent interaction between classical FF and QM calculation was found to be 0.16−3.59 kJ/mol, with an average of 1.4 kJ/mol. The three homologous solvents (2,5,8,11,14pentaoxapentadecane, 1,2-bis(2-methoxyethoxy)ethane, 1-methoxy-2-(2-methoxyethoxy)ethane) exhibit the largest CO2 solubility differences between simulation and experiment. This finding suggests that a more accurate FF is needed for this type of solvent to better predict CO2 solubility. In spite of this CO2 solubility difference between simulation and experiment, simulations give a similar CO2 solubility trend in different solvents as the experimental data (Figure 2), especially for the homologous compounds. Finally, it was noted that the simulated

Figure 3. Simulated CO2 solubility data at 298 K in homologous solvents 1-methoxy-2-(2-methoxyethoxy)ethane, 1,2-bis(2methoxyethoxy)ethane, 2,5,8,11,14-pentaoxapentadecane, and 2,5,8,11,14,17-hexaoxaoctadecane versus the simulated solvent freevolume fraction. The chemical formulas are also included in the figure. The error bars obtained from standard block average calculations are also shown. The line is the linear fitting of the simulation data with a linear correlation coefficient of 0.999. F



Article

compounds. As a result, CO2 solubilities in the homologous compounds are mainly determined by the solvent free-volume fractions. In contrast, CO2 interactions are different between different types of solvents, which leads to the CO2 solubility dependence on both CO2−solvent interaction and the solvent free-volume fraction. Consequently, CO2 solubilities in all 27 solvents studied in this work only show a very weak linear correlation with the solvent free-volume fractions (not shown here), with a linear correlation coefficient of 0.31. CO2 solubilities in common solvents, namely, methanol, propionitrile, and tetrahydrofuran, were also calculated at 233 and 253 K (Figure 4). When the temperature is decreased, CO2

solubility is also comparable with the available experimental data (Figure 6). The absolute H2 solubility difference between

Figure 6. Comparison of H2 solubility values in seven solvents at 298 K between simulation and experimental data.30 The simulation error bars obtained from standard block average calculations are also shown.

simulation and experiment is 10−41%, with an average of 22%. The simulation and experimental H2 solubilities in different solvents exhibit similar trends. Note that experiments also show H2 solubility differences in the same solvent, for example, H2 solubility in decane exhibits a difference of 20% between different experimental data,30 which is larger than the solubility difference between simulation and experiment. The simulated H2 solubilities in all solvents exhibit a strong positive linear correlation with the solvent free-volume fractions (Figure 7), with a linear correlation coefficient of 0.84. This is

Figure 4. Simulated CO2 (open) and H2 (filled) solubility values in methanol (circles), propionitrile (squares), and tetrahydrofuran (triangles) solvents at 233 K, 253 K, and 298 K. The lines are linear regressions of the simulation data. The error bars are smaller than the symbols and are not shown.

solubilities in these solvents increase. From the slopes of the lines (Figure 4), CO2 heats of absorption in these three solvents were estimated to be −15 kJ/mol. Although these three solvents exhibit similar CO2 heats of absorption, both propionitrile and tetrahydrofuran exhibit a larger CO2 solubility than methanol does. The tetrahydrofuran solvent exhibits about two times larger CO2 solubility compared with methanol both at 298 and 253 K. 3.2. H2 Solubility. The simulated H2 solubilities in 27 solvents at 298 K are given in Figure 5 and Table S1. The H2 solubilities are 0.0038−0.0534 mol/(MPa·L). The simulated H2 Figure 7. Simulated H2 solubility values in 27 solvents versus the corresponding solvent free-volume fractions at 298 K. The line is the linear fitting of the simulation data. The simulation error bars are smaller than the symbols and are not shown.

expected since H2 interacts very weakly with all solvent molecules.2 As a result, H2 solubility in a solvent is primarily determined by the solvent free-volume fraction. When temperature is increased, H2 solubilities in methanol, propionitrile, and tetrahydrofuran (Figure 4) increase. From the slopes of the lines (Figure 4), H2 heats of absorption in these three solvents were estimated to be 3−5 kJ/mol. When the temperature is increased, simulations show that the solvent freevolume fraction also increases, which in turn leads to increased H2 solubility in solvent. The simulated mole fraction-based H2 solubilities in these three solvents also increase at elevated temperatures (not shown). Similarly, the experimental H2 solubilities (mole fraction-based) in methanol, n-methyl-2pyrrolidone, propylene carbonate, decane, and tetrahydrofuran also increase when the temperature increases.30

Figure 5. Simulated H2 solubility values in 27 solvents (Table 1) at 298 K obtained from Henry’s law constant calculations. Also shown are error bars obtained from standard block average calculations. G



Article

3.3. CO2/H2 Solubility Selectivity. The CO2/H2 solubility selectivity was estimated from the pure CO2 and H2 solubility values. The selectivity values are between 16.7−263 (Figure 8),

Figure 8. CO2/H2 solubility selectivity data in 27 solvents at 298 K estimated from their respective pure gas solubilities (Table S1). Also shown are error bars obtained from error propagation calculations based on the error bars for pure gas solubilities.

Figure 10. Simulated CO2 diffusivities in 25 solvents at 298 K. Also shown are error bars obtained from standard block average calculations.

larger than 1, which suggests that CO2 exhibits larger gas solubility than H2 in all solvents studied in this work. This result is due to the stronger CO2−solvent interaction compared with H2-solvent interaction.2 The simulated CO2/H2 selectivity values are comparable with the experimental data (Figure 9),

with CO2 mole fractions of 0.5−5%. The CO2 diffusivities are between 10−12 and 10−8 m2/s, and most solvents exhibit CO2 diffusivities larger than 10−9 m2/s. For homologous compounds (1-methoxy-2-(2-methoxyethoxy)ethane, 1,2-bis(2methoxyethoxy)ethane, 2,5,8,11,14-pentaoxapentadecane, 2,5,8,11,14,17-hexaoxaoctadecane), when the number of repeat units is increased, CO2 diffusivity in the solvent decreases, which is due to decreased solvent free-volume fraction (Figure 3) and increased solvent viscosity (see below). 3.5. Solvent Viscosity. Simulated solvent viscosities are shown in Figure 11 and Table S2. The solvent viscosities are

Figure 9. Comparison of CO2/H2 solubility selectivities in six solvents at 298 K between simulation and experimental data. The CO2/H2 selectivity was estimated from their respective pure gas solubilities (Table S1). Also shown are error bars obtained from error propagation calculations based on the error bars for pure gas solubilities.

with an average absolute difference of 17%. Both simulation and experiment show that n-methyl-2-pyrrolidone and tetrahydrofuran solvents exhibit both larger CO2 solubility (Figure 2) and larger CO2/H2 solubility selectivity (Figure 9) than methanol. Additionally, simulations show that 2,4,6,8-tetraoxanonane and propionitrile exhibit both larger CO2 solubility (Figure 1) and larger CO2/H2 solubility selectivity (Figure 8) than methanol. These theoretical predictions need to be validated experimentally. 3.4. CO2 Diffusivity. Simulated CO2 diffusivities in solvents are shown in Figure 10 and Table S2. Note that the CO2 diffusivities correspond to low CO2 concentrations in solvents,

Figure 11. Simulated viscosity values for 25 solvents at 298 K. Also shown are error bars obtained from standard block average calculations. H



Article

between 10−4 and 102 Pa·s at 298 K. When the solvent viscosity is smaller than 10−3 Pa·s, the simulated and experimental solvent viscosity values31 are comparable to each other (Figure 12 and

experimental data, this relationship will be very useful to predict CO2 diffusivity from the solvent viscosity, which can easily be determined from the experiment. Note that it is difficult to directly measure CO2 diffusivity in a solvent.

4. CONCLUSIONS CO2 and H2 solubility, CO2/H2 solubility selectivity, CO2 diffusivity, and solvent viscosity in 27 commercially available physical solvents at 298 K were calculated from molecular simulation by using the CHARMM36 all-atom force field for most solvent molecules. All the simulation results are comparable to the available experimental data with the exception of highly viscous solvents. All of the aforementioned simulation results are predictive, since only the van der Waals radius parameters for all solvent molecules were slightly tuned to reproduce the experimental solvent density and the FF parameters were then fixed to calculate other properties. The average absolute CO2 solubility difference between simulation and experiment is 28%. Additionally, simulations exhibited similar CO2 solubility trends in different solvents as the experimental data. The reasonable CO2 solubility agreement between simulation and experiment is partly due to the reasonably accurate CO2−solvent interaction obtained from the classical FF compared with the QM calculation. The CO2 solubilities in the homologous compounds were found to exhibit nearly a perfect positive linear correlation with the solvent freevolume fractions, which is partly due to the similar CO2 interactions with the homologous compounds. The simulated H2 solubility is also comparable with the experimental data and exhibits a trend similar to that of the experiment. The average absolute H2 solubility difference between simulation and experiment is 22%. The simulated H2 solubilities in all solvents exhibit a strong positive linear correlation with the solvent free-volume fractions (correlation coefficient of 0.84). This is expected since H2 interacts very weakly with solvent molecules. As a result, H2 solubility in a solvent is mainly determined by the solvent free-volume fraction. The simulated CO2/H2 solubility selectivity is also comparable with the experimental data, with an average absolute difference of 17%. In addition, simulations predict that nmethyl-2-pyrrolidone, tetrahydrofuran, 2,4,6,8-tetraoxanonane, and propionitrile exhibit both larger CO2 solubility and larger CO2/H2 solubility selectivity than methanol. Some predictions have been confirmed experimentally while other predictions (2,4,6,8-tetraoxanonane and propionitrile) remain to be validated. Simulated solvent viscosity differ from the experimental values for viscous solvents. When the solvent viscosity is less than 10−3 Pa·s, the viscosity difference between simulation and experiment is only 30%. However, when the solvent viscosity is larger than 10−3Pa·s, the average absolute viscosity difference is large (236%). In spite of this difference, the simulation gives a viscosity trend similar to that of the experiment. For example, for the homologous compounds containing the −(OCH2CH2)− repeat unit, when the repeat number is increased, both simulation and experiment show that the solvent viscosity increases. Finally, the simulated CO2 diffusivity was found to correlate very well with the simulated solvent viscosity, and their ± 0.03 . relationship is given by DCO2 = (2.6 ± 0.3) × 10−9/η0.59 solvent

Figure 12. Comparison of viscosity values for 13 solvents at 298 K between simulation and experimental data.31 The simulation error bars obtained from standard block average calculations are also shown.

Table S2), with an average absolute difference of 30%. When the solvent viscosity is larger than 10−3 Pa·s, the viscosity difference between simulation and experimental data is large, with an average absolute difference of 236%. Despite this large difference, the simulation gives a viscosity trend similar to that of the experimental data. For example, for the homologous compounds (1-methoxy-2-(2-methoxyethoxy)ethane, 1,2-bis(2-methoxyethoxy)ethane, 2,5,8,11,14-pentaoxapentadecane, 2,5,8,11,14,17-hexaoxaoctadecane), when the solvent −(OCH2CH2)− repeat unit is increased, both simulation and experiment show that the solvent viscosity increases. Because of the similar viscosity trend between simulation and the experimental data, it is reasonable to use molecular modeling to predict the solvent viscosity trend even before new solvents are synthesized. CO2 diffusivity (D) correlates very well with the solvent viscosity (η) (Figure 13), with a linear correlation coefficient of −0.97, which suggests that log D exhibits a very strong negative linear correlation with log η. It would be interesting to compare this simulated relationship between D and η with experimental data. If this simulated relationship holds true for the

Figure 13. Simulated CO2 diffusivity (D) values in 25 solvents versus the corresponding simulated solvent viscosities (η) at 298 K. The line indicates a linear fitting of log D to log η with a linear negative correlation coefficient of −0.97. The fitted equation is also given in the figure. I



■

Article

Computatuonal Simulations, and Techno-economic Analysis. Int. J. Greenhouse Gas Control 2016, 49, 364−371. (4) Vanommeslaeghe, K.; Hatcher, E.; Acharya, C.; Kundu, S.; Zhong, S.; Shim, J.; Darian, E.; Guvench, O.; Lopes, P.; Vorobyov, I.; MacKerell, A. D., Jr. 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. (5) Chemical Book. https://www.chemicalbook.com/ (accessed 2018, November, 02). (6) Burr, B.; Lyddon, L. Which Physical Solvent is Best for Acid Gas Removal? Hydrocarb. Process. 2009, 88, 43−50. (7) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon: Oxford, 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 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.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, O.; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian 09, revision A.1. Gaussian, Inc.: Wallingford, CT, 2009. (11) 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. (12) Antechamber, An Accessory Software Package For Molecular Mechanical Calculations. http://ambermd.org/antechamber/ antechamber.pdf (accessed 2018, November, 29. (13) Cui, S.; de Almeida, V. F.; Hay, B. P.; Ye, X.; Khomami, B. Molecular Dynamics Simulation of Tri-n-butyl-Phosphate Liquid: A Force Field Comparative Study. J. Phys. Chem. B 2012, 116, 305−313. (14) Shi, W.; Siefert, N. S.; Baled, H. O.; Steckel, J. A.; Hopkinson, D. P. Molecular Simulations of the Thermophysical Properties of Polyethylene Glycol Siloxane (PEGS) Solvent for Precombustion CO2 Capture. J. Phys. Chem. C 2016, 120, 20158−20169. (15) VEGA ZZ Manual. http://www.msg.ucsf.edu/local/programs/ Vega/pages/gl_index.htm (accessed 2018November, 29). (16) Jamali, S. H.; Wolff, L.; Becker, T. M.; Bardow, A.; Vlugt, T. J. H.; Moultos, O. A. Finite-Size Effects of Binary Mutual Diffusion Coefficients from Molecular Dynamics. J. Chem. Theory Comput. 2018, 14, 2667−2677. (17) Jamali, S. H.; Hartkamp, R.; Bardas, C.; Sohl, J.; Vlugt, T. J. H.; Moultos, O. A. Shear Viscosity Computed from the Finite-Size Effects of Self-Diffusivity in Equilibrium Molecular Dynamics. J. Chem. Theory Comput. 2018, 14, 5959−5968. (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) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kale, L.; Schulten, K. Scalable

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jced.8b01228. Summary of the numerical values for CO2 and H2 solubilities, and solvent density obtained both from simulation and experiment; summary of the numerical values for simulated CO2 diffusivity, and solvent viscosity both from simulation and experiment (PDF) run_sim_file directory containing psf and crd files for a single solvent molecule, the force field parameter file for that solvent, and psf and PDB files for the system (ZIP) sigma_tune_example directory containing the codes used to automatically tune σ values (ZIP)

■

AUTHOR INFORMATION

Corresponding Author

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

Wei Shi: 0000-0002-7295-9443 Robert L. Thompson: 0000-0003-3266-6600 Funding

This technical effort was performed in support of the National Energy Technology Laboratory’s ongoing research in carbon capture under the 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 AECOM. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors thank for Jeffrey Culp for measuring CO2 solubilities in 2-benzylpyridine and 2,2-methylene bispyridine by using the gravimetric method,32 and Surya Tiwari for making a long-table template by using the LYX software. 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) The Development of Machine Learning, Group Contribution and Molecular Modeling Approach to Screen Physical Solvents for Gas Separation. https://aiche.confex.com/aiche/2018/meetingapp.cgi/ Paper/519751 (accessed 2018, December, 2). (2) Shi, W., 2018, unpublished. (3) 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, J



Article

Molecular Dynamics with NAMD. J. Comput. Chem. 2005, 26, 1781− 1802. (20) 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. (21) 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. (22) 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. (23) Singh, R.; Marin-Rimoldi, E.; Maginn, E. J. A Monte Carlo Simulation Study To Predict the Solubility of Carbon Dioxide, Hydrogen, and Their Mixture in the Ionic Liquids 1-Alkyl-3methylimidazolium bis(trifluoromethanesulfonyl)amide ([C(n)mim(+)][Tf2N(−)], n = 4, 6). Ind. Eng. Chem. Res. 2015, 54, 4385−4395. (24) 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. (25) Shi, W.; Widger, L. R.; Sarma, M.; Lippert, C. A.; Alman, D. E.; Liu, K. Molecular Modeling of the Physical Properties for Aqueous Amine Solution Containing a CO2 Hydration Catalyst. Ind. Eng. Chem. Res. 2017, 56, 11644−11651. (26) Grossfield, A.; Patrone, P.; Roe, D.; Schultz, A.; Siderius, D.; Zuckerman, D. Best Practices for Quanti1cation of Uncertainty and Sampling Quality in Molecular Simulations [Article v1.0]. Living J. Comp. Mol. Sci. 2019, 1, 5067. (27) Confidence Interval on the Mean. http://onlinestatbook.com/ 2/estimation/mean.html, (accessed 2019, February, 10). (28) Distribution Critical Values. https://www.stat.tamu.edu/ ~lzhou/stat302/T-Table.pdf (accessed 2019February, 10). (29) Hansen, C. M. Hansen Solubility Parameters: A User’s Handbook, 2nd ed.; CRC: Boca Raton, FL, 2007; pp 179−182. (30) IUPAC Solubility Data Series (Online). https://srdata.nist.gov/ solubility/IUPAC/iupac.aspx (accessed 2018November, 02). (31) Diky, V.; Chirico, R. D.; Frenkel, M.; Bazyleva, A.; Magee, J. W.; Paulechka, E.; Kazakov, A.; Lemmon, E. W.; Muzny, C. D.; Smolyanitsky, A. Y.; Townsend, S.; Kroenlein, K. NIST Standard Reference Database 103a, ThermoData Engine (TDE), version 10.0a (Pure Compounds); Standard Reference Data Program, National Institute of Standards and Technology: Gaithersburg, MD, 20899. (32) Thompson, R. L.; Culp, J.; Tiwari, S. P.; Basha, O.; Damodaran, K.; Shi, W.; Siefert, N.; Hopkinson, D. Effect of Molecular Structure on Hydrophobic Carbon Capture Solvents Consisting of Hybrid Poly Ethylene Glycol/Poly Dimethyl Siloxane Units. 2018, unpublished.

K


Molecular Simulations of CO2 and H2 Solubility, CO2 Diffusivity, and

Recommend Documents