Molecular Modeling Analysis of CO2 Absorption by Glymes - The

Jan 29, 2018 - The properties of diglyme + CO2 systems were analyzed through density functional theory and molecular dynamics methods with the objecti...
2 downloads 0 Views 4MB Size
Article pubs.acs.org/JPCB

Cite This: J. Phys. Chem. B 2018, 122, 1948−1957

Molecular Modeling Analysis of CO2 Absorption by Glymes Alberto Gutiérrez,† Mert Atilhan,*,‡,§ and Santiago Aparicio*,† †

Department of Chemistry, University of Burgos, 09001 Burgos, Spain Department of Chemical Engineering, Texas A&M University at Qatar, Doha, Qatar § Gas and Fuels Research Center, Texas A&M University, College Station, Texas 77843, United States ‡

S Supporting Information *

ABSTRACT: The properties of diglyme + CO2 systems were analyzed through density functional theory and molecular dynamics methods with the objective of inferring the microscopic properties of CO2 capture by glyme-based solvents and the effect of ether group regarding solvents affinity toward CO2. Calculations of diglyme + CO2 molecular clusters using density functional theory allowed accurate quantification and characterization of short-range intermolecular forces between these molecules, whereas the molecular dynamics simulation of diglyme + CO2 liquid mixtures, for different CO2 contents, were the means to infer the properties and dynamics of bulk liquid phases upon CO2 absorption. Likewise, liquid diglyme + CO2 gas interfaces were also studied using molecular dynamics methods to examine the kinetics of CO2 capture, adsorption at the gas−liquid interface, and the mechanism of interface crossing, which is of pivotal importance for the design of CO2 capturing units.



istics. Traditional amine-based CO2 capturing methods17−20 have been employed in the natural gas industry with success for many decades, but the characteristics of natural gas treatment and those of CO2 capture from flue gases are remarkably different; different pressure ranges and CO2 concentrations, which are also combined with high corrosion effect,21−23 high solvent degradation,24 and high regeneration costs,25−29 lead to a huge increase in the cost of electricity production when using amine-based carbon capture units. Hence, the use of materials such as ionic liquids,30 deep eutectic solvents,31,32 membranes,33 metal organic frameworks,34 covalent organic frameworks,35 and many other types of sorbents have been studied and combined with different technological approaches.36−39 Among the studied liquid state CO2 absorbents that work with physisorption,37,40 glymes (glycol dieters) showed a large affinity for CO2 molecules leading to high absorption performances.41 Although the enhancement of CO2 absorption by oxygenation of alkyl chains is a known effect,42 the mechanism of interaction between glymes and CO2 molecules is still not understood, especially at the microscopic level. Therefore, in order to obtain a microscopic knowledge of the mechanism of CO2 capture by glyme-based physical absorbents, a computational chemistry study was carried out in this work. Diglyme, bis(2-methoxyethyl) ether (2G), was considered as a model system of glymes, and their properties regarding CO2 absorption were studied. In the first stage, Density Functional

INTRODUCTION The problem of increasing atmospheric levels of CO21,2 and its effects on climate change3 due to excessive utilization of fossil based fuels,4 requires the development of new, sustainable, and effective carbon capture technologies in chemical process industries.5,6 Emissions from fossil-fuelled power plants are among the largest CO2 sources,7,8 and thus, developing suitable sorption technologies would contribute to a decrease and control of CO2 atmospheric levels. This can be considered as one of the pivotal challenges that the chemical sciences have faced.9 Therefore, the development of economically and technologically viable carbon capture technologies is required in the transition toward a fully decarbonized electricityproduction framework10,11 which may last for decades. The development of sustainable and economically viable CO2 capture technologies at the industrial scale for power plants must consider several engineering challenges due to the large scale of the issue.12 There is a growing technical consensus that the main challenge is the absence of suitable high-performance sorbents (either physical or chemical),13 which can perform the absorption and desorption steps favorably enough in comparison to the current state of the art method(s) both technically and economically. From the technological viewpoint, the main challenges of developing such materials for CO2 capture in fossil-fuel based power plants are based on the proper characteristics of flue gases, and low CO2 partial pressures,14 which require highly efficient sorbents to avoid the increased costs of electricity production to unfeasible levels.15,16 Therefore, many materials have been studied this past decade with respect to very different chemical and physical character© 2018 American Chemical Society

Received: October 17, 2017 Revised: December 22, 2017 Published: January 29, 2018 1948

DOI: 10.1021/acs.jpcb.7b10276 J. Phys. Chem. B 2018, 122, 1948−1957

Article

The Journal of Physical Chemistry B Theory (DFT) calculations on 2G + CO2 molecular clusters were carried out to obtain an accurate characterization of glymes−CO2 intermolecular interactions. However, the reported DFT calculations in the gas phase lack the data to infer additional effects in bulk liquid phases, such as larger molecular clusters, development of solvation shells, or molecular packing effects. Thus, molecular dynamics (MD) simulations were carried out. The study of liquid phase properties via MD was done considering 2G + CO2 liquid mixtures for different mixture compositions, that is, increasing CO2 mole fraction, according to experimental compositions reported by Kodama et al.41 The analysis of intermolecular forces, molecular arrangements and distributions, and dynamic properties allow a characterization of the effect of CO2 on 2G properties and thus the CO2 capture mechanism as well. Likewise, to infer the kinetics of CO2 capture at the microscopic level, MD simulations on liquid 2G in contact with CO2 gas phases were carried out, inferring the properties at liquid−gas interface (i.e., CO2 adsorption), diffusion across the interface, and mass transfer form gas to liquid phases. The reported computational study herein focuses on the microscopic characterization of glymes−CO2 systems, analyzing the effect of ether-functionalization on CO2 capture, and contributing to the development of suitable solvents for effective CO2 capture.

Table 1. Systems Considered for NPT Molecular Dynamics Simulations of x CO2 + (1 − x) 2Ga x

N (2G)

0 0.19 0.31 0.44 0.49 0.59 0.68 0.77 0.86 1

1000 1000 1000 1000 1000 1000 1000 1000 1000 −

N (CO2)

p/bar

T/K

235 449 786 961 1439 2125 3348 6143 1000

1 12 21 31 35 45 55 64 71 50

313 313 313 313 313 313 313 313 313 313

a Notation: x, mole fraction; p, pressure; T, temperature; N, number of molecules. Composition, pressure, and temperature of the mixtures were selected according to experimental CO2 solubility in 2G as reported by Kodama et al.41

production (10 ns) runs. The MD simulations were carried out in a NPT ensemble at 313 K, and pressures are reported in Table 1. The Nose−Hoover method was used for pressure and temperature control, the Ewald method58 (15 Å for cutoff radius) was used for treating Coulombic interactions, Lorentz− Berthelot mixing rules were used for Lennard-Jones cross terms, and the Tuckerman−Berne double time step algorithm59 (1 and 0.1 fs for long and short time steps) was used for solving equations of motion. Periodic boundary conditions were applied. Regarding simulations of the 2G−CO2 systems, a 2G liquid phase was put in contact (i.e., developing an interface) with a gas phase containing 2000 CO2 molecules, with both phases being previously equilibrated at 313 K. Simulations as a function of time (10 ns long) were carried out in the NVT ensemble, with the same methods as for the mentioned simulations for the 2G + CO2 mixtures.



METHODS The modeling of 2G-CO2 intermolecular interactions was performed with the use of molecular clusters with one 2G molecule and n (=1, 2, or 3) CO2 molecules, which were optimized using the B3LYP43−45 functional, with Grime’s method (DFT-D3) for treating dispersion interactions,46 and 6311++G(d,p) basis set. For the calculation of interaction energies, ΔE, a counterpoise procedure was applied to correct the Basis Set Superposition Error (BSSE).47 All DFT calculations were carried out by using the ORCA program.48 DFT calculations were carried out only at the B3LYP-D3/6311++G(d,p) theoretical level, and no comparison was carried out with higher level approaches such as MP2. Previous studies have shown that differences for ΔE between B3LYP with dispersion corrections and MP2 methods are lower than 0.5 kcal mol−1,49 and thus, methods reported in this work are a reasonably accurate approach to ΔE50−52 for the studied systems. Intermolecular interactions were also analyzed based on Bader’s theory (Atoms in a Molecule, AIM) using the Multiwfn code.53 MD simulations were carried out using the MDynaMix v.5.2 program.54 The force field parametrization for 2G was previously reported by our group,55 whereas that for CO2 was obtained from the literature.56 Two different types of systems were considered for studying the behavior of 2G + CO2 by using MD: (i) liquid mixtures of 2G + CO2 with different compositions, to study the mechanism of interaction, intermolecular forces, and other effects controlling 2G affinity for CO2 molecules; (ii) liquid 2G in contact with a CO2 gas phase, to study the kinetics of CO2 absorption, adsorption of CO2 molecules at 2G liquid interface, and other interfacial phenomena controlling CO2 capture by 2G. 2G + CO2 mixtures were modeled using the systems reported in Table 1, considering a fixed number of 2G molecules and different contents of CO2 to mimic experimental absorption data by Kodama et al.41 Simulation boxes were built using the Packmol57 program. MD simulations were developed in two stages: (i) initial equilibration (10 ns); followed by (ii)



RESULTS AND DISCUSSION The characterization of 2G-CO2 intermolecular forces carried out using DFT was done for 2G + n (=1,2,3) CO2 molecular clusters, at different interaction sites. The presence of ether group is considered to enhance the affinity toward CO2 molecules, and thus several configurations were considered with CO2 molecules placed close to oxygen atoms in 2G, Figure 1. Two different types of oxygen atoms may be considered in glymes: those placed at the end of the chains (O1), between a methylene and a terminal methyl group, and those in the chain between two methylene groups (O2). Therefore, the interaction of CO2 molecules with these two types of oxygen atoms was studied. In the case of a single CO2 molecules interacting with 2G, results that are presented in Figure 1 show that larger binding energies for interactions through O1 than through O2 sites, although the low values for both sites are indicative of van der Waals interactions of moderate strength. Interaction through the O2 site is characterized by CO2 molecules adopting a parallel arrangement to 2G chains with a carbon atom on top of an oxygen atom (at 3.07 Å) and oxygen atoms on top of the two-methylene carbons neighboring O2 site. In the case of the O1 site, the distance between the carbon atom in CO2 and the O1 site is 2.75 Å, which is lower than that for the interactions through the O2 site, which would justify the larger binding energies through O1 than through O2. Likewise, the interaction through the O1 site is characterized by the interaction between the hydrogen atom 1949

DOI: 10.1021/acs.jpcb.7b10276 J. Phys. Chem. B 2018, 122, 1948−1957

Article

The Journal of Physical Chemistry B

Figure 1. Structures of 2G + n CO2 (n = 1, 2, 3) clusters, for several possible positions (Pi) optimized in the gas phase at B3LYP-D3/6-311++g(d,p) theoretical level. Counterpoise corrected binding energies, ΔE, are reported inside each panel.

in the terminal methyl group and oxygen atom in CO2 which is stronger than the interaction with the methylene hydrogens when CO2 interacts through the O2 site. This is confirmed by the bond paths reported in the AIM analysis (Figure S1, Supporting Information). Additional details of AIM analysis for intermolecular interactions are reported in the Supporting Information. Larger clusters were considered to infer the effect of several CO2 molecules competing and interacting with a single 2G molecule, Figure 1. The interaction of two CO2 molecules with one 2G molecule leads to binding energies as the sum of those corresponding to single interactions, but when one of the molecules is out of the main interaction sites (O1 and O2) binding energy decreases remarkably, as reported in Figure 1. Regarding the larger number of CO2 molecules, results in Figure 1 confirm additive binding energies and no perturbation on interactions between neighbor sites. Therefore, each 2G molecule can fit one CO2 molecule per O1 or O2 site, confirming the sorption capacity experimentally reported by Kodama et al.41 The characterization of bulk 2G liquid phases upon CO2 absorption is also analyzed in this work. The absorption of CO2 may be accompanied by solvent swelling, that is, volume expansion upon gas sorption, %Vexp, which was quantified according to a method proposed by Gallager et al.60 and reported in Figure 2. The large %Vexp reported in Figure 2 should be analyzed considering the high CO2 concentrations in 2G, which were also confirmed by Kodama et al.41 This expansion shows a disruption of the 2G structure when gas is captured. CO2 concentrations lower than 0.6 yielded a percentage volume expansion of 40%, which could be considered as reasonable with respect to possible industrial scale applications of the solvent. The 2G + CO2 mixtures showed nonideal behavior of the evolution of density with composition, as shown in Figure 3. When the density profile is examined, there is a region at which a sudden increase take place between the 0 and 0.2 CO2

Figure 2. Percentage of volume expansion, %Vexp, upon CO2 absorption for x CO2 + (1 − x) 2G mixtures at 313 K and pressures reported in Table 1.

Figure 3. Density, ρ, and dielectric constant, ε, for x CO2 + (1 − x) 2G mixtures at 313 K and the pressures reported in Table 1.

1950

DOI: 10.1021/acs.jpcb.7b10276 J. Phys. Chem. B 2018, 122, 1948−1957

Article

The Journal of Physical Chemistry B fractions. For CO2 concentrations between 0.2 and 0.6, there is a transitional density region observed, and there is a linear evolution in the density for mole fractions higher than 0.6. This is in agreement with the behavior of %Vexp; therefore, three regions in the behavior of 2G + CO2 mixtures are inferred, which is confirmed by the evolution of dielectric constant, ε. The force field parametrizations that are used are able to predict ε accurately for both pure 2G and CO2; 6.886 for 2G at 318.15 K was reported experimentally by Lago et al.,61 and the value was 6.1 as a result of MD simulations at 313 K (Figure 3). Yet, experimental ε values were reported close to 1 at 313.15 K for pure CO2,62 which is also in agreement with values for CO2rich mixtures as presented in Figure 3. One of the main origins of CO2 capturing by liquid sorbents is based on the rearrangement of available empty molecular space or cavities to fit the absorbed molecules. The spontaneous rearrangement of fluids cavities has been considered as a pivotal factor for CO2 capture in fluids such as ionic liquids,63,64 but there are not enough scientific studies that handle 2G at molecular level for the same purpose. To analyze the distribution of cavities in pure 2G and in 2G + CO2 mixtures, the method studied by Margulis65 was considered; and to do that, a network of random points was distributed in simulation boxes and the distance between each point and the edge of the nearest atom (defined by its van der Waals radius) is considered as the size of the cavity defined as spheres. Probabilities distributions of cavities are reported in Figure 4

Figure 5. Center-of-mass radial distribution functions, g(r), for x CO2 + (1 − x) 2G mixtures at 313 K and pressures reported in Table 1. Values inside the panels show the position of maxima and minima (bold) for the most relevant peaks.

Figure 5b also show a trend to self-associate as the presence of at least two well-defined RDFs peaks indicate. RDFs for 2G− CO2 pairs, Figure 5c, show complex behavior, five RDF peaks are obtained for distances lower than 20 Å, with their maxima remaining constant upon the increase of CO2 concentration. Nevertheless, the intensity of RDFs for 2G−CO2 pairs is remarkably lower than those for 2G−2G or CO2−CO2 pairs, which confirms the large trend of both 2G and CO2 to selfassociate. The behavior of RDFs shows the existence of microdomains of 2G and CO2 molecules, with some regions allowing the development of 2G−CO2 interactions but with the trend to self-associate being the pivotal role in the structuring of these mixed fluids. The formation of microdomains have also been reported for other types of CO2 absorbents, such as ionic liquids,66,67 although in this case these microdomains originated from the presence of large apolar chains with self-aggregation trends leading to apolar microdomains,68 whereas in the case of 2G, microdomains are produced even with the presence of a polar ether group because of the trend of 2G molecules to develop parallel arrangements for improving van der Waals interactions. The number of molecules in the solvation shells obtained from RDFs is reported in Figure 6; regarding CO2 molecules around 2G (Figure 6a,b), a nonlinear evolution with a mixture composition is inferred, with three well-defined regions corresponding to those for macroscopic properties in Figures 2 and 4. The number of CO2 molecules in the first solvation shell (with radius of 4.2 Å, Figure 5c) is almost constant for 2G-rich mixtures (x(CO2) < 0.4), but its value of ∼1.8 CO2 molecules is low in contrast with the large value of CO2 molecules in the second shell (Figure 6c). This is in agreement with the picture of a low number of CO2 molecules

Figure 4. Probability distribution plots of spherical cavities radius, R, for x CO2 + (1 − x) 2G mixtures at 313 K and pressures reported in Table 1.

for pure 2G and for mixtures with different CO2 mole fractions. These results show that a rearrangement of 2G spherical cavities is produced upon CO2 absorption, with a decrease of cavities sizes, in agreement with the increase of density reported in Figure 3. The molecular arrangements at the molecular level were analyzed using radial distribution functions (RDFs), Figure 5. Pure 2G is a largely structured fluid and its center-of-mass RDF reported in Figure 5a shows a first and narrow solvation peak (with maximum at 4.6 Å) which is followed by three additional peaks that are wider and less intense than the first one and separated roughly 4 Å. This structure is maintained upon CO2 absorption, with a slight decrease in the intensity of RDFs peaks because of dilution of 2G molecules, and thus shows the strong trend of 2G molecules to self-associate. Regarding the behavior of CO2 molecules when absorbed in 2G, results in 1951

DOI: 10.1021/acs.jpcb.7b10276 J. Phys. Chem. B 2018, 122, 1948−1957

Article

The Journal of Physical Chemistry B

Figure 6. Number of (a) CO2 molecules around 2G and (c,d) CO2 molecules around CO2 molecules for x CO2 + (1 − x) 2G mixtures at 313 K and pressures reported in Table 1. Values are obtained from radial distribution functions reported in Figure 5 for the corresponding first and second solvation shells, as defined by the first and second minimum in each function.

Figure 7. Spatial distribution functions of CO2 center-of-mass around 2G molecules in x CO2 + (1 − x) 2G mixtures at 313 K and pressures reported in Table 1.

Figure 8. Snapshot of structure after 10 ns of simulation for x CO2 + (1 − x) 2G mixtures at 313 K and pressures reported in Table 1. CO2 molecules are reported as van der Waals spheres and 2G molecules as green lines. Plots in the bottom row show the same results as in the top row but CO2 molecules are omitted to improve visibility of microheterogeneities.

The distribution of CO2 molecules is characterized spatially by the presence around 2G-oxygen atoms, as shown as spatial distributions functions (SDFs) in Figure 7 for the full composition range. Likewise, SDFs around the terminal methyl groups show the presence of CO2 molecules around these sites,

directly interacting with 2G molecules and a large number of CO2 molecules in neighbor CO2-rich clusters self-associated by quadrupolar interactions but not directly interacting with 2G molecules. This behavior is confirmed by results for solvation shells of CO2−CO2 type reported in Figure 6c,d. 1952

DOI: 10.1021/acs.jpcb.7b10276 J. Phys. Chem. B 2018, 122, 1948−1957

Article

The Journal of Physical Chemistry B

CO2 interactions when compared with those of CO2−2G. DFT results in Figure 1 reported interaction energies for CO2−2G clusters, and thus, for comparison purposes the interaction energy for the CO2−CO2 dimer was calculated at the same theoretical level. The properties of the CO2−CO2 dimers have been widely studied in the literature using many different theoretical approaches,70 but all of them have shown that the slipped parallel configuration is preferred over the T-shaped or other orientations for CO2−CO2 dimers, thus only the slipped parallel dimer was considered in this work. The ΔE for the CO2−CO2 slipped parallel dimer at B3LYP-D3/6-311++g(d,p) level is −7.0 kJ mol−1, which is less than the half of the ΔE for the CO2−2G clusters reported in Figure 1, and thus, this is not the reason for the larger residence times for CO2−CO2 reported in Figure 9. These larger residence times can be attributed to larger CO2−CO2 clusters in comparison with CO2−2G, which decrease the molecular mobility of CO2 molecules in the solvation sphere of central CO2 molecules. Regarding the macroscopic dynamics, self-diffusion coefficients (D) were calculated from mean square displacements (msd) using Einstein’s equation for fully diffusive regimes, Figure 10. The β exponent of msd as a function of simulation time was calculated to ensure that a fully diffusive (Fickian) regime was reached in the time intervals from which D coefficients were calculated; this method has been used in the literature for confirming diffusive regimes in highly viscous fluids such as ionic liquids.71−73 The obtained β = 1 values for all the systems reported in Figure 10 confirm that proper equilibrations were carried out, and D values are obtained under a diffusive regime. This is in agreement with the low viscosity of 2G (0.777 mPa s at 303.15 K16), which leads to short simulation times for reaching a fully diffusive regime. The reported results show larger D values for CO2 than for 2G, although the evolution of D with mixture composition follows the same trend for both compounds. Three regions are inferred for D as a function of CO2 mole fraction: (i) first, D values decrease (i.e., molecular mobility is reduced with increasing CO2 concentration for 2Grich mixtures), which corresponds to mixtures with large 2G domains with remarkable microheterogeneities as reported in Figure 8 and it hinders the molecular mobilities for both types of molecules, (ii) in the second composition region, D values remains almost constant, corresponding to the aggregation stages of small CO2-rich domains into a larger one but still with the presence of 2G-self-aggregated domains (Figure 8), and (iii) a third composition region, corresponding to fluids dominated by CO2, in which self-aggregation of 2G molecules is less remarkable and fluids are dominated by CO2-like domains, thus remarkably increasing the molecular mobility as showed in Figure 10. The evolution of intermolecular interaction energies, Einter, with mixtures composition is reported in Figure 11. First, it should be remarked that all the reported Einter follow nonlinear trends with mixtures composition, in agreement with the evolution of mixtures’ structuring according to the three composition regions reported in previous paragraphs. Considering that Einter are of van der Waals type, as confirmed by DFT in Figure S1 (Supporting Information) and Table 1, those Einter involving a larger number of atoms are larger than those with a lower number of atoms; that is, Einter for CO2−CO2 interactions are lower than the remaining ones. The fluids’ structuring is characterized by 2G-microdomains in 2G-rich mixtures and its evolution toward smaller 2G clusters diluted in CO2-rich mixtures is confirmed by the decrease (in absolute value) of

but at larger distances than CO2 molecules placed around the 2G-oxygen atoms, which leads to weaker interactions. This is confirmed by DFT calculations of a single CO2 molecule interacting directly with a terminal methyl group in 2G leading to ΔE = −6.2 kJ mol−1, which is remarkably lower than the values reported in Figure 1. This arrangement is maintained in the full composition range without remarkable changes, which shows a small number of CO2 molecules are interacting closely with 2G ether groups and most of the CO2 molecules selfassociated. Additional information on the structure of 2G−CO2 mixtures may be inferred from the simulation snapshots that are reported in Figure 8. These snapshots show the formation of microaggregates of 2G molecules even at low CO2 concentrations (see panel corresponding to x(CO2) = 0.19 in Figure 8). These aggregates are especially well-defined for 2Grich mixtures, but are present even for high CO2 concentrations. Likewise, these 2G aggregates are characterized by parallel arrangements of 2G molecules adopting pillar-like clustered structures. Moreover, the heterogeneous distribution of CO2 molecules is also inferred from Figure 8. CO2-rich domains are present for low CO2 concentrations, which are aggregated upon an increase of CO2 concentration leading to CO2-dominated fluid. This confirms the three regimes inferred for 2G-CO2 mixtures as a function of mixture composition: (i) first region, low CO2 concentrations, self-aggregated CO2 domains dispersed in 2G microdomains, (ii) second region, CO2-domains start to self-aggregate and 2G-domains are dispersed in them, and (iii) third region CO2-dominated fluids with small 2G-clusters, less organized than for 2G-rich mixtures. The dynamic properties of the studied mixtures were characterized at the local level around each type of molecule by the so-called residence time, tres, which were calculated as previously reported.69 tres defines the average time that a molecule remains inside a sphere of defined radius around another central molecule. The radius of spheres used to define t res were obtained from RDFs reported in Figure 5, corresponding to the first minima used to define first solvation spheres. These results show slower dynamics of CO2 molecules for CO2−CO2 interactions than for those of CO2−2G, Figure 9, although for both cases tres values increase with increasing CO2 concentration in a nonlinear trend. The slower dynamics of CO2 in CO2 clusters could lead to the idea of stronger CO2−

Figure 9. Residence time, tres, for the reported molecular pairs, in the first solvation sphere, defined from radial distribution functions reported in Figure 6 for x CO2 + (1 − x) 2G mixtures at 313 K and pressures reported in Table 1. 1953

DOI: 10.1021/acs.jpcb.7b10276 J. Phys. Chem. B 2018, 122, 1948−1957

Article

The Journal of Physical Chemistry B

Figure 10. Mean square displacements, msd, and self-diffusion coefficients (calculated from msd using Einsten’s equation), D, for each type of molecule in x CO2 + (1 − x) 2G mixtures at 313 K and pressures reported in Table 1.

Figure 11. Intermolecular interaction energies, Einter (sum of LennardJones and Coulombic contributions), for the reported molecular pairs in x CO2 + (1 − x) 2G mixtures at 313 K and pressures reported in Table 1.

Figure 12. Number density, ρ, profiles for CO2 molecules at 2G interface as a function of simulation time, t. z stands for the coordinate perpendicular to the 2G−CO2 interface, and zGDS stands for the coordinate of Gibbs dividing surface. Values are calculated at 313 K.

Einter for 2G−2G pairs and the increase (in absolute value) for 2G−CO2 interactions (Figure 11). The large increase of 2G− CO2 Einter on going to CO2-rich mixtures may be justified considering that the decrease of 2G self-association, that is, the decrease in the size of 2G clusters and the vanishing of the large pillar-like microheterogeneities (Figure 8) allows a better interaction with CO2 molecules at the cost of decreasing those of the 2G−2G. This behavior is confirmed by the relationship between the evolution of Einter and D coefficients data with mixtures composition, Figures 10 and 11. For CO2-rich mixtures CO2−CO2 interactions prevail, which are remarkably weaker than those of CO2−2G, which combined with the vanishing of 2G−2G interactions leads to the sudden increase in D coefficients for CO2-rich mixtures. Nevertheless, the large Einter obtained for most of the composition range should be noted as it should lead to largely nonideal mixtures. The kinetics of CO2 absorption by 2G was studied by the MD simulation and analysis of interfaces between liquid 2G and CO2 in the gas phase, and their evolution with simulation time. The large affinity of CO2 molecules for 2G is confirmed by the development of an adsorbed and dense layer of CO2 molecules on top of liquid 2G for short simulation times. Results that are presented in Figure 12 show that just after 0.2

ns of simulation CO2, a high number of CO2 molecules are adsorbed at the interface and some of them have crossed the Gibbs dividing surface starting to penetrate inside the 2G liquid layer in the subsurface regions. The density at the interface layer increases with simulation time and also the number of subsurface CO2 molecules, but a concentration limit at the interface is reached at around 5 ns, Figure 12. The time evolution of Einter reported in Figure 13 confirms the three stages of the CO2 capture mechanism: absorption that is followed by diffusion to subsurface regions and then migration to bulk liquid 2G phase. The decrease (in absolute value) in Einter for 2G−2G interactions shows structural changes in the 2G liquid at the interface to accommodate adsorbed CO2 molecules, but this rearrangement after diffusion of CO2 molecules to subsurface regions is not large because when the stationary regime is reached 2G−2G Einter are almost the same as before the CO2 adsorption. The formation of microheterogeneities in the liquid 2G reported in Figure 8 for 2G + CO2 mixtures is confirmed from the simulations of interfaces. Results in Figure 14 show how 2G self-association, pillar-like domains with 2G in parallel arrangements, is inmediately formed when CO2 molecules diffuse from the subsurface regions toward the bulk liquid phase, thus being a 1954

DOI: 10.1021/acs.jpcb.7b10276 J. Phys. Chem. B 2018, 122, 1948−1957

Article

The Journal of Physical Chemistry B

interface and by considering the density profiles: low residence when adsorbed, then very quick diffusion to the subsurface layer, and then migration top bulk liquid phase. Therefore, largely heterogeneous mixed fluids are formed with a microstructure dependent on mixture composition and a trend to self-aggregate for both types of molecules combined with large affinity by 2G toward CO2 molecules.



ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.7b10276. Results of AIM analysis of intermolecular interactions (PDF)



Figure 13. Time evolution of intermolecular interaction energies, Einter (sum of Lennard-Jones and Coulombic contributions), for the reported molecular pairs for 2G in contact with CO2 interface at 313 K.

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected]. *E-mail: [email protected]. ORCID

Santiago Aparicio: 0000-0001-9996-2426 Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS This work was funded by Junta de Castilla y León (Spain, project BU324U14). We also acknowledge The Foundation of Supercomputing Center of Castile and León (FCSCL, Spain) for providing supercomputing facilities. The statements made herein are solely the responsibility of the authors.

Figure 14. Snapshot of structure after 10 ns of simulation for (a) pure 2G and (b) 2G with CO2 interface, both at 313 K. CO2 molecules are reported as van der Waals spheres and 2G molecules as green lines. Arrows indicate some relevant microheterogeneities with 2G molecules adopting large parallel arrangements.



REFERENCES

(1) Monastersky, R. Global Carbon Dioxide Levels Near Worrisome Milestone. Nature 2013, 497, 13−14. (2) Jackson, R. B.; Canadell, J. G.; Le Quéré, C.; Andrew, R. M.; Korsbakken, J. I.; Peters, G. P.; Nakicenovic, N. Reaching Peak Emissions. Nat. Clim. Change 2016, 6, 7−10. (3) Peters, G. P.; Andrew, R. M.; Boden, T.; Canadell, J. G.; Ciais, P.; LeQuéré, C.; Marland, G.; Raupach, M. R.; Wilson, C. The Challenge to Keep Global Warming Below 2 °C. Nat. Clim. Change 2013, 3, 4−6. (4) Climate Change 2014: Synthesis Report; Contribution of Working Groups I, II, and III to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change; Pachauri, R. K., Meyer, L.A., Eds.; IPCC: Geneva, Switzerland, 2014; 151 pp. (5) Koytsoumpa, E. I.; Bergins, C.; Kakaras, E. The CO2 Economy: Review of CO2 Capture and Reuse Technologies. J. Supercrit. Fluids 2017, 132, 3−16. (6) Senftle, T. P.; Carter, E. A. The Holy Grail: Chemistry Enabling and Economically Viable CO2 Capture, Utilization, and Storage Strategy. Acc. Chem. Res. 2017, 50, 472−475. (7) Davis, S. J.; Socolow, R. H. Commitment accounting of CO2 emissions. Environ. Res. Lett. 2014, 9, 084018. (8) Echevarria-Huaman, R. N.; Jun, T. X. Energy Related CO2 Emissions and the Progress on CCS Projects: A Review. Renewable Sustainable Energy Rev. 2014, 31, 368−385. (9) Sholl, D. S.; Lively, R. P. Seven Chemical Separations to Change the World. Nature 2016, 532, 435−437. (10) Chu, S. Carbon Capture and Sequestration. Science 2009, 325, 1599. (11) Chu, S.; Cui, Y.; Liu, N. The Path Towards Sustainable Energy. Nat. Mater. 2017, 16, 16−22. (12) Cebrucean, D.; Cebrucean, V.; Ionel, I. CO2 Capture and Storage from Fossil Fuel Power Plants. Energy Procedia 2014, 63, 18− 26.

response of 2G molecules to the presence of absorbed CO2 molecules to maintain the large 2G−2G Einter reported in Figure 14, and confirming that 2G is able to absorb large amounts of CO2 but leading to highly heterogeneous mixed fluids.



CONCLUSIONS The absorption of CO2 by liquid 2G was studied initially by quantum chemistry methods, showing that CO2 molecules can develop effective van der Waals interactions with ether groups in glymes with stronger interactions via terminal oxygen atoms. Nevertheless, molecular dynamics simulations have shown that 2G + CO2 mixtures are complex fluids with structuring being largely dependent on the mixture composition. The 2G-rich fluids are characterized by large microheterogeneities with 2G molecules developing domains with parallel arrangements. As the composition increases in CO2, CO2-rich domains start to self-aggregate dispersing 2G-clusters, which leads to small 2G clusters dispersed in CO2-dominated fluids. The intermolecular interaction mechanism is characterized by CO2 molecules interacting with ether groups, and it was observed that their trend to self-associate hinders larger clustering around 2G. The absorption is also characterized by a large volume expansion and by the rearrangement of the available empty volume, leading to a decrease of the size of available cavities with increasing CO2 concentration. The kinetics of capture is characterized by very fast adsorption of CO2 molecules at the 1955

DOI: 10.1021/acs.jpcb.7b10276 J. Phys. Chem. B 2018, 122, 1948−1957

Article

The Journal of Physical Chemistry B (13) Bhattacharyya, D.; Miller, D. C. Post-Combustion CO2 Capture Technologies − A Review of Processes for Solvent − Based and Sorbent − Based CO2 Capture. Curr. Opin. Chem. Eng. 2017, 17, 78− 92. (14) Heldebrant, D.; Koech, P. K.; Glezakou, V. A.; Rousseau, R.; Malhotra, D.; Cantu, D. C. Water-Lean Solvents for Post-Combustion CO2 Capture: Fundamentals, Uncertainties, Opportunities, and Outlook. Chem. Rev. 2017, 117, 9594−9624. (15) Rao, A. B.; Rubin, E. S. A Technical, Economic, and Environmental Assesment of Amine-Based CO2 Capture Technology for Power Plant Greenhouse Gas Control. Environ. Sci. Technol. 2002, 36, 4467−4475. (16) Kanniche, M.; Le Moullec, Y.; Authier, O.; Hagi, H.; Bontemps, D.; Neveux, T.; Louis-Louisy, M. Up-to-date CO2 Capture in Thermal Power Plants. Energy Procedia 2017, 114, 95−103. (17) Rao, A. B.; Rubin, E. S. A Technical, Economic, and Environmental Assessment of Amine-Based CO2 Capture Technology for Power Plant Greenhouse Gas Control. Environ. Sci. Technol. 2002, 36, 4467−4475. (18) Rochelle, G. T. Amine Scrubbing for CO2 Capture. Science 2009, 325, 1652−1654. (19) Puxty, G.; Rowland, R.; Allport, A.; Yang, Q.; Bown, Q.; Burns, R.; Maeder, M.; Attalla, M. Carbon Dioxide Postcombustion Capture: A Novel Screening Study of the Carbon Dioxide Absorption Performance of 76 Amines. Environ. Sci. Technol. 2009, 43, 6427− 6433. (20) Chowdhury, F. A.; Yamada, H.; Higashii, T.; Goto, K.; Onoda, M. CO2 Capture by Tertiary Amine Absorbents: A Performance Comparison Study. Ind. Eng. Chem. Res. 2013, 52, 8323−8331. (21) Kittel, J.; Idem, R.; Gelowitz, D.; Tontiwachwuthikul, P.; Parrain, G.; Bonneau, A. Corrosion in MEA Units for CO2 Capture: Pilot Plant Studies. Energy Procedia 2009, 1, 791−797. (22) Kittel, J.; Fleury, E.; Vuillemin, B.; Gonzalez, S.; Ropital, F.; Oltra, R. Corrosion in Alkanolamine Used for Acid Gas Removal: From Natural Gas Processing to CO2 Capture. Mater. Corros. 2012, 63, 223−230. (23) Rafat, A.; Atilhan, M.; Kahraman, R. Corrosion Behavior of Carbon Steel in CO 2 Saturated Amine and Imidazolium-, Ammonium-, and Phosphonium-Based Ionic Liquid Solutions. Ind. Eng. Chem. Res. 2016, 55, 446−454. (24) Martin, S.; Lepaumier, H.; Picq, D.; Kittel, J.; de Bruin, T.; Faraj, A.; Carrett, P. L. New Amines for CO2. Capture. IV. Degradation, Corrosion, and Quantitative Structure Property Relationship Model. Ind. Eng. Chem. Res. 2012, 51, 6283−6289. (25) Rao, A.; Rubin, E. S. Identifying Cost-Effective CO2 Control Levels for Amine-Based CO2 Capture Systems. Ind. Eng. Chem. Res. 2006, 45, 2421−2429. (26) Husebye, J.; Brunsvold, A. L.; Roussanaly, S.; Zhang, X. Techno Economic Evaluation of Amine based CO2 Capture: Impact of CO2 Concentration and Steam Supply. Energy Procedia 2012, 23, 381−390. (27) El Nasr, A. S.; Nelson, T.; Abu-Zahra, M. R. M. Technoeconomic Evaluation Methodology and Preliminary Comparison of an Amine-Based and Advanced Solid Sorbent-Based CO2 Capture Process for NGCC Power Plants. Energy Procedia 2013, 37, 2432− 2442. (28) Vaccarelli, M.; Carapellucci, R.; Giordano, L. Energy and Economic Analysis of the CO2 Capture from Flue Gas of Combined Cycle Power Plants. Energy Procedia 2014, 45, 1165−1174. (29) Manzolini, G.; Sanchez-Fernandez, E.; Rezvani, S.; Macchi, E.; Goether, E. L. V.; Vlugt, T. J. H. Economic assessment of novel amine based CO2 capture technologies integrated in power plants based on European Benchmarking Task Force methodology. Appl. Energy 2015, 138, 546−558. (30) Zeng, S.; Zhang, X.; Bai, L.; Zhang, X.; Wang, H.; Wang, J.; Bao, D.; Li, M.; Zhang, S.; Liu, X. Ionic-Liquid-Based CO2 Capture Systems: Structure, Interaction and Process. Chem. Rev. 2017, 117, 9625−9673.

(31) García, G.; Aparicio, S.; Ullah, R.; Atilhan, M. Deep Eutectic Solvents: Physicochemical Properties and Gas Separation Applications. Energy Fuels 2015, 29, 2616−2644. (32) Sarmad, S.; Mikkola, J. P.; Ji, X. Carbon Dioxide Capture with Ionic Liquids and Deep Eutectic Solvents: A New Generation of Sorbents. ChemSusChem 2017, 10, 324−352. (33) Rafiq, S.; Deng, L.; Hägg, M. B. Role of Facilitated Transport Membranes and Composite Membranes for Efficient CO2 Capture − A Review. ChemBioEng Rev. 2016, 3, 68−85. (34) Yu, J.; Xie, L. H.; Li, J. R.; Ma, Y.; Seminario, J. M.; Balbuena, P. B. CO2 Capture and Separations Using MOFs: Computational and Experimental Studies. Chem. Rev. 2017, 117, 9674−9754. (35) Zeng, Y.; Zou, R.; Zhao, Y. Covalent Organic Frameworks for CO2 Capture. Adv. Mater. 2016, 28, 2855−2873. (36) Lee, S. Y.; Park, S. J. A Review on Solid Adsorbents for Carbon Dioxide Capture. J. Ind. Eng. Chem. 2015, 23, 1−11. (37) Ben-Mansour, R.; Habib, M. A.; Bamidele, O. E.; Basha, M.; Qasem, N. A. A.; Peedikakkal, Q. A.; Laoui, T.; Ali, M. Carbon Capture by Physical Adsorption: Materials, Experimental Investigations and Numerical Modeling and Simulations − A Review. Appl. Energy 2016, 161, 225−255. (38) Creamer, A. E.; Gao, B. Carbon-Based Adsorbents for Postcombustion CO2 Capture: A Critical Review. Environ. Sci. Technol. 2016, 50, 7276−7289. (39) Leung, D. Y. C.; Caramanna, G.; Maroto-Vader, M. M. An Overview of Current Status of Carbon Dioxide Capture and Storage Technologies. Renewable Sustainable Energy Rev. 2014, 39, 426−443. (40) De Riva, J.; Suarez-Reyes, J.; Moreno, D.; Díaz, I.; Ferro, V.; Palomar, J. Ionic Liquids for Post-combustion CO2 Capture by Physical Absorption: Thermodynamic, Kinetic and Process Analysis. Int. J. Greenhouse Gas Control 2017, 61, 61−70. (41) Kodama, D.; Kanakubo, M.; Kokubo, M.; Hashimoto, S.; Nanjo, H.; Kato, M. Density, Viscosity, and Solubility of Carbon Dioxide in Glymes. Fluid Phase Equilib. 2011, 302, 103−108. (42) Van Ginderen, P.; Herrebout, W. A.; van der Veken, B. J. van der Waals Complex of Dimethyl Ether with Carbon Dioxide. J. Phys. Chem. A 2003, 107, 5391−5396. (43) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti Correlation-Energy Formula into a Functional of the Electron Density. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37, 785−789. (44) Becke, A. D. Density-Functional Exchange-Energy Approximation with Correct Asymptotic Behavior. Phys. Rev. A: At., Mol., Opt. Phys. 1988, 38, 3098−3100. (45) Becke, A. D. Density-Functional Thermochemistry. III. The Role of Exact Exchange. J. Chem. Phys. 1993, 98, 5648−5652. (46) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A. Consistent and Accurate Ab Initio Parametrization of Density Functional Dispersion Correction (DFT-D) for the 94 Elements H-Pu. J. Chem. Phys. 2010, 132, 154104. (47) Simon, S.; Duran, M.; Dannenberg, J. How Does Basis Set Superposition Error Change the Potential Surfaces for HydrogenBonded Dimers? J. Chem. Phys. 1996, 105, 11024. (48) Neese, F. The ORCA Program System. WIREs Comput. Mol. Sci. 2012, 2, 73−78. (49) García, G.; Atilhan, M.; Aparicio, S. Assessment of DFT Methods for Studying Acid Gas Capture by Ionic Liquids. Phys. Chem. Chem. Phys. 2015, 17, 26875−26891. (50) Damas, G. B.; Dias, A. B. A.; Costa, L. T. A Quantum Chemistry Study for Ionic Liquids Applied to Gas Capture and Separation. J. Phys. Chem. B 2014, 118, 9046−9064. (51) Day, C.; Yang, Y.; Fisher, A.; Liu, Z.; Cheng, S. Interaction of CO2 with Metal Cluster-Functionalized Ionic Liquids. J. CO2 Util. 2016, 16, 257−263. (52) Seyedhosseini, B.; Izadyar, M.; Housaindokht, M. R. A Computational Exploration of H2S and CO2 Capture by Ionic Liquids Based on α-Amino Acid Anion and N7,N9-Dimethyladeninium Cation. J. Phys. Chem. A 2017, 121, 4352−4362. (53) Lu, T.; Chen, F. Multiwfn: A Multifunctional Wavefunction Analyzer. J. Comput. Chem. 2012, 33, 580−592. 1956

DOI: 10.1021/acs.jpcb.7b10276 J. Phys. Chem. B 2018, 122, 1948−1957

Article

The Journal of Physical Chemistry B (54) Lyubartsev, A. P.; Laaksonen, A. MDynaMix - A Scalable Portable Parallel MD Simulation Package for Arbitrary Molecular Mixtures. Comput. Phys. Commun. 2000, 128, 565−589. (55) Alcalde, R.; Gutiérrez, A.; Atilhan, M.; Trenzado, J. L.; Aparicio, S. Insights into Glycol Ether−Alkanol Mixtures from a Combined Experimental and Theoretical Approach. J. Phys. Chem. B 2017, 121, 5601−5612. (56) Shi, W.; Maginn, E. 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. (57) Martínez, L.; Andrade, R.; Birgin, E. G.; Martínez, J. M. Packmol: A Package for Building Initial Configurations for Molecular Dynamics Simulations. J. Comput. Chem. 2009, 30, 2157−2164. (58) Essmann, U. L.; Perera, M. L.; Berkowitz, T.; Darden, H.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103, 8577−8593. (59) Tuckerman, M.; Berne, B. J.; Martyna, G. J. Reversible Multiple Time Scale Molecular Dynamics. J. Chem. Phys. 1992, 97, 1990−2001. (60) Gallagher, P. M.; Coffey, M. P.; Krukonis, V. J.; Klasutis, N. ACS Symposium Series 406; In Supercritical Fluid Science and Technology; Johnson, K. P., Penninger, J. M. L., Eds.; American Chemical Society: Washington, DC, 1989; pp 334−354. (61) Lago, A.; Rivas, M. A.; Legido, J.; Iglesias, T. P. Study of Static Permittivity and Density of the Systems {(n-Nonane + Monoglyme or Diglyme)} at Various Temperatures. J. Chem. Thermodyn. 2009, 41, 257−264. (62) Moriyoshi, T.; Kita, T.; UOsaki, Y. Static Relative Permittivity of Carbon Dioxide and Nitrous Oxide up to 30 MPa. Ber. Bunsenges. Phys. Chem. 1993, 97, 589−596. (63) Huang, X.; Margulis, C. J.; Li, Y.; Berne, B. J. Why Is the Partial Molar Volume of CO2 So Small When Dissolved in a Room Temperature Ionic Liquid? Structure and Dynamics of CO2 Dissolved in [Bmim+] [PF6−]. J. Am. Chem. Soc. 2005, 127, 17842−17851. (64) Klähn, M.; Seduraman, A. What Determines CO2 Solubility in Ionic Liquids? A Molecular Simulation Study. J. Phys. Chem. B 2015, 119, 10066−10078. (65) Margulis, C. J. Computational Study of Imidazolium-based Ionic Solvents with Alkyl Substituents of Different Lengths. Mol. Phys. 2004, 102, 829−838. (66) Sheridan, Q. R.; Oh, S.; Morales-Collazo, O.; Castner, E. W.; Brennecke, J. F.; Maginn, E. J. Liquid Structure of CO2−Reactive Aprotic Heterocyclic Anion Ionic Liquids from X-ray Scattering and Molecular Dynamics. J. Phys. Chem. B 2016, 120, 11951−11960. (67) Pereiro, A. B.; Pastoriza-Gallego, M. J.; Shimizu, K.; Marrucho, I. M.; Lopes, J. N. C.; Piñeiro, M. M.; Rebelo, L.P. N. On the Formation of a Third, Nanostructured Domain in Ionic Liquids. J. Phys. Chem. B 2013, 117, 10826−10833. (68) Canongia-Lopes, J. N. A.; Pádua, A. A. H. Nanostructural Organization in Ionic Liquids. J. Phys. Chem. B 2006, 110, 3330−3335. (69) Atilhan, M.; Anaya, B.; Ullah, R.; Costa, L. T.; Aparicio, S. Double Salt Ionic Liquids Based on Ammonium Cations and Their Application for CO2 Capture. J. Phys. Chem. C 2016, 120, 17829− 17844. (70) Kalugina, Y. N.; Buryak, I. A.; Ajili, Y.; Vigasin, A. A.; Jaidane, N. E.; Hochlaf, M. Explicit Correlation Treatment of the Potential Energy Surface of CO2 Dimer. J. Chem. Phys. 2014, 140, 234310. (71) Del Pópolo, M. G.; Voth, G. A. On the Structure and Dynamics of Ionic Liquids. J. Phys. Chem. B 2004, 108, 1744−1752. (72) 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. (73) Franca, J. M. P.; Nieto de Castro, C. A.; Pádua, A. H. Molecular Interactions and Thermal Transport in Ionic Liquids with Carbon Nanomaterials. Phys. Chem. Chem. Phys. 2017, 19, 17075−17087.

1957

DOI: 10.1021/acs.jpcb.7b10276 J. Phys. Chem. B 2018, 122, 1948−1957