Thermodynamic Properties of Supercritical Mixtures of Carbon Dioxide

Article pubs.acs.org/jced

Thermodynamic Properties of Supercritical Mixtures of Carbon Dioxide and Methane: A Molecular Simulation Study Cassiano G. Aimoli,†,‡,§ Edward J. Maginn,*,§ and Charlles R. A. Abreu‡,∥ †

Santos Basin Exploration & Production Operations Unit, Petróleo Brasileiro S.A., Av. Conselheiro Nébias 175, Paquetá, Santos, SP 11015-001, Brasil ‡ School of Chemical Engineering, State University of Campinas, Av. Albert Einstein 500, Cidade Universitária Zeferino Vaz, Campinas, SP 13083-852, Brasil § Department of Chemical and Biomolecular Engineering, University of Notre Dame, 182 Fitzpatrick Hall, Notre Dame, Indiana 46556, United States ∥ School of Chemistry, Federal University of Rio de Janeiro, Av. Athos da Silveira Ramos 149, Ilha do Fundão, Rio de Janeiro, RJ 21941-909, Brasil S Supporting Information *

ABSTRACT: Volumetric and second-derivative thermodynamic properties of binary mixtures of carbon dioxide and methane were calculated under supercritical conditions at pressures up to 99.93 MPa and temperatures between (323.15 and 573.15) K using molecular dynamics simulations with the multistate Bennett acceptance ratio technique. Eleven compositions were studied, ranging from pure methane to pure carbon dioxide in 0.10 increments in carbon dioxide mole fraction. The molecular simulations utilized the TraPPE and SAFT-γ force fields to model methane and carbon dioxide. Because of the unavailability of experimental data, the results were compared against the extended range estimation and extrapolation of the GERG-2008 equation for the calculation of volume expansivity, isothermal compressibility, isobaric and isochoric heat capacities, Joule−Thomson coefficient, and speed of sound of the binary mixtures. The agreement between the simulation results and those obtained through the extrapolation of the GERG-2008 equation beyond its stated range of validity is quite good, suggesting that the force field representation of the binary mixture is accurate, especially when compared with the extrapolation of traditional cubic equations of state. Although more computationally expensive, molecular simulations relying upon only a few physically meaningful parameters are able to give similar results to those obtained with the semiempirical multiparameter GERG-2008 model. the data fitting but usually limits applications out of the original range of validity of the equations.8 In an attempt to build physically based equations of state with broader operating ranges, the statistical associating fluid theory (SAFT)9 and its derivatives10−14 combine thermodynamic perturbation theory for associating fluids with the physical parametrization of fluid interactions describing the thermodynamic behavior of complex fluids. These models are generally based on terms directly related to molecular-level interactions such as the repulsion−dispersion contribution of individual segments, the contribution of segment connectivity to form molecules, and the contribution of segment association to form molecular complexes.15

1. INTRODUCTION The measurement and modeling of thermodynamic properties of natural gas components such as methane (CH4) and carbon dioxide (CO2) are of high importance for many industrial applications such as natural gas production, processing, and transportation. They are also important for carbon capture and storage techniques1,2 and methane recovery from hydrates, which is a promising alternative fossil fuel resource.3,4 In all cases, knowledge of fluid behavior under different conditions and the availability of mathematical models to represent the thermodynamic properties accurately are needed for the engineering and performance evaluation of processing plants.5 Cubic equations of state (e.g., the Peng−Robinson6 and Soave−Redlich−Kwong7 equations) are capable of providing estimates in good agreement with experimental data in specific thermodynamic regions. Such models are normally parametrized using critical properties and other experimental data. This allows a good representation of regions considered during © XXXX American Chemical Society

Special Issue: Modeling and Simulation of Real Systems Received: February 1, 2014 Accepted: April 29, 2014

A

dx.doi.org/10.1021/je500120v | J. Chem. Eng. Data XXXX, XXX, XXX−XXX

Journal of Chemical & Engineering Data

Article

pressure ranges, molecular simulations are relatively fast and inexpensive to carry out. They may help provide a deeper understanding of fluid behavior and can be used to supplement experimental data for the development of accurate thermophysical property models. Along with the technological evolution of computational resources and algorithms, additional efforts have been made in order to improve the statistical quality of molecular simulation results and to reduce simulation times. The united-atom model is a good approximation to simulate molecular systems of specific stable chemical compounds (i.e., the carbon−hydrogen bonds found in an alkane chains32,34,35), allowing in many cases a very significant decrease in the number of interactions to be computed without compromising the quality of the physical representation. Another example is the efficiency improvement based on increasing the amount of information extracted from independent simulations,36−39 allowing the representation of equilibrium averages obtained through Monte Carlo or molecular dynamics simulations as continuous functions of the temperature and the pressure of the system.40 In this work, we demonstrate how molecular simulation techniques can be used successfully to predict volumetric and second-derivative thermodynamic properties of binary mixtures of CO2 and CH4 in regions where only experimental measurements of density are available. Two different united-atom CH4 models are studied along with three-site and single-site CO2 models. The multistate Bennett acceptance ratio (MBAR)39 technique is applied to improve the statistical quality of the simulation results. A comparison of the molecular simulation results against predictions obtained via extrapolations of the GERG-2008 and Peng−Robinson (hereafter called GERG and PR, respectively) equations of state is provided, aiming to elucidate how the physically meaningful molecular models can be used to evaluate the performance of empirical models out of their ranges of validity.

On the other hand, the broad availability of high-quality experimental data for pure compounds enables the construction of empirical mathematical models that are weakly associated with the physical description of fluid behavior but designed to predict thermodynamic properties within experimental uncertainties over a wide range of conditions. Highly accurate multiparameter equations of state that have been fitted to comprehensive collections of experimental data are available for several relevant pure compounds, including CH416 and CO2.17 Likewise, the same multiproperty fitting technique can be applied to build equations to predict mixture behavior.18,19 Despite the fact that the Groupe Europeen de Recherches Gazieres (GERG)-200420 equation of state and its expanded version GERG-200821 are able to give accurate thermodynamic property predictions of common natural gas mixtures at typical temperature and pressure ranges, the uncertainties of properties calculated outside the range of validity are unknown. Obviously the difficulties in handling highly volatile gases and the large number of individual experiments required to cover a broad range of temperatures, pressures, and mixture compositions may discourage comprehensive experimental studies in spite of the great technical importance of such information, particularly when unusual and extreme conditions are of interest. Although most industrial processes dealing with binary gaseous mixtures of CH4 and CO2 can be found under thermodynamic conditions covered by accurate multiparameter models, the interest in industrial applications that require highly compressed gases is growing rapidly. In the particular case of the so-called presalt22 oil accumulations recently found on the continental shelves off the southern coast of Brazil, the high volumes of CO2 originally found in the underground formations have to be removed from the produced natural gas and injected back into the reservoirs. The manipulation of mixtures containing mainly CO2 and CH4 at pressures that can exceed the range covered by the existing equations of state is required, as well as the ability to predict fluid behavior over a wide range of conditions. With the increase in the availability of powerful computer hardware combined with fast, highly parallelized calculation algorithms,23−26 molecular simulation is becoming an interesting option for describing the properties of fluids in regions where experimental data are unreliable, insufficient, or unavailable. Starting from the knowledge of how interactions occur at the atomic and molecular level, analytic potential energy functions (force fields) are developed to describe these interactions. Molecular simulations are then conducted to provide microscopic structural information as well as its connection to macroscopic fluid properties through the fundamentals of statistical mechanics.27,28 Besides being generally more physically meaningful, the atomistic force fields usually rely upon a much smaller number of parameters to describe thermodynamic properties of various pure compounds and mixtures29−32 compared with accurate multiparameter equations of state. This is considered an important advantage since the molecular models are often suitable to describe fluid behavior over a extended range of thermodynamic conditions, sometimes far from the original conditions used to develop the force field parameters.33 The computational demand of molecular simulations is usually higher than those of empirical models and equations of state. However, compared with acquisition of the extensive experimental data required to develop multiparameter equations of state, particularly those applicable under extreme

2. METHODOLOGY Initially, the densities of binary mixtures of CO2 and CH4 were calculated through molecular simulations over a wide range of compositions and several conditions of pressure and temperature in the supercritical region, using two distinct models to represent the atomic and molecular interactions in the systems. The GERG multiparameter equation and the PR equation of state were also used to provide estimates of density in the same region of thermodynamic conditions. All of the results were compared against available experimental data.41 Afterward, additional volumetric and second-derivative properties were calculated through molecular simulations using isobaric−isothermal volume and energy fluctuations. Because of the lack of experimental data describing these properties under the conditions studied, extrapolations of the GERG equation beyond its generally accepted range of validity were tested against the molecular simulations results in order to provide information regarding the level of agreement between the two estimation approaches. Lastly, a comparison between the computational time requirements found for the two molecular models employed in the simulations was made. The MBAR technique is introduced as an efficient method to provide continuous estimation of thermodynamic properties based on results of independent simulations at multiple thermodynamic states. 2.1. Interaction Potentials. All of the molecules involved in this study were either considered to be rigid (in the case of B



Article

Table 2. Functional Forms of the Peng−Robinson, SAFT-VR Mie, and GERG-2008 Equations of State

Table 1. Force Field Parameters for Each CO2 and CH4 Model CO2 model type (ε/kB)/K σ/Å q/e lC−O/Å θ/deg λm λn

CH4

equation

TraPPE32

SAFT-γ43

TraPPE34

SAFT-γ44

three-site 27.00 (C−C) 79.00 (O−O) 2.80 (C−C) 3.05 (O−O) +0.70 (C) −0.35 (O) 1.16 180.0 12.00 6.00

single-site 361.69

single-site 148.00

single-site 153.36

3.7410

3.7300

3.7412

−

−

−

RT a − (V − b) [V (V + b) + b(V − b)]

Peng−Robinsona

p=

SAFT-VR Mieb

A Aideal Amono Achain Aassoc = + + + NkBT NkBT NkBT NkBT NkBT

GERG-2008c

α = α ideal +

NG

∑ xiα0ri + Δα r i=1

a

− − 23.00 6.66

− − 12.00 6.00

Explicit in pressure, defined as the sum of the repulsion and attraction pressures. bExplicit in the Helmholtz free energy, defined as the sum of free energy contributions due to ideal gas behavior, segment−segment intermolecular interactions, formation of chain molecules, and molecular association, respectively. cExplicit in the reduced Helmholtz free energy α = A/NkBT, defined as the sum of an ideal gas part, the contribution of the pure substances, and the so-called departure function to account for mixtures.

− − 12.65 6.00

the three-site CO2 model) or described as single particles (for both united-atom CH4 and CO2 models), and therefore, no intramolecular energy was taken into account; intermolecular interactions were calculated only with dispersion−repulsion and electrostatic contributions. The most common model for pairwise dispersion−repulsion interactions is the Lennard-Jones potential, which follows the same functional form as the Mie potential42 but has fixed values for the exponents in the dispersion and repulsion terms. The Mie pairwise-additive potential is given by UMie

common functional form adopted

⎡⎛ ⎞ λ m ⎛ ⎞ λ n ⎤ σ σ = Urep + Udisp = Cε⎢⎢⎜⎜ ⎟⎟ − ⎜⎜ ⎟⎟ ⎥⎥ r r ⎝ ij ⎠ ⎦ ⎣⎝ ij ⎠

Table 3. Average Absolute Relative Deviations (AARDs) between Densities Calculated Using Molecular Simulations or Equations of State and Experimental Data AARD/%

(1)

where rij represents the distance between a pair of atoms i and j, ε is the potential well depth, σ is related to the atomic diameters, and λm and λn are the repulsion and dispersion exponents that characterize the potential energy behavior. The constant C is a function of the given exponents, defined as λn / λm − λn λm ⎛ λm ⎞ C= ⎜ ⎟ λm − λn ⎝ λn ⎠

GERG

TraPPE

SAFT-γ

PR

0.31

0.98

0.93

1.58

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

0.03 0.86 0.69 0.44 0.36 0.23 0.16 0.21 0.31 0.25 0.03

0.07 0.67 0.67 0.61 0.88 0.87 0.95 1.01 1.24 1.32 2.42

0.55 1.28 1.08 0.91 0.97 0.92 0.92 0.96 0.89 0.82 1.06

0.79 1.45 1.43 1.30 1.40 1.51 1.53 1.77 1.79 1.96 2.49

19.94 29.94 39.94 59.93 79.93 99.93

0.52 0.33 0.34 0.24 0.13 0.18

1.52 1.25 1.02 0.66 0.32 0.68

1.35 1.01 0.87 0.62 0.48 1.01

2.06 1.97 1.46 0.89 1.55 1.64

323.15 373.15 473.15 573.15

0.23 0.24 0.44 0.43

1.49 0.77 0.80 1.03

1.30 0.84 0.77 1.00

1.93 1.77 1.57 1.13

global xCO2

p/MPa

(2)

When λm = 12 and λn = 6, eq 1 reduces to the familiar LennardJones 12−6 potential. The effects of Coulombic interactions were taken into account through the incorporation of partially charged sites into the potential model according to the expression qiqj Uelec = 4π ϵ0rij (3)

T/K

where qi and qj are the partial charges of the atoms and ϵ0 is the vacuum permittivity. In this work, two different sets of force field parameters, originally developed to model vapor−liquid equilibrium, were employed to represent the molecular interactions in the binary mixtures of CO2 and CH4. First, the well-established transferable potentials for phase equilibria (TraPPE) models of CO232 and CH434 were used. TraPPE represents CO2 with three Lennard-Jones sites united by two collinear chemical bonds of fixed length and with partial charges on each atom center, while CH4 is modeled as a single, uncharged LennardJones site. Standard Lorentz−Berthelot mixing rules provide cross-interaction parameters between unlike atoms. Second, two SAFT-γ single-site models for CO243 and CH444 were

utilized, in which pairwise interactions are calculated using eq 1 with specific parameters and exponents for each molecule.43−45 The SAFT-γ models were originally developed via a top-down concept for the development of accurate coarse-grained intermolecular potentials using the SAFT-VR Mie equation of state to fit experimental data and generate the parameters of interaction. The unlike size parameters are obtained using the arithmetic-mean combining rule for the molecular diameter: σii + σjj σij = (4) 2 C



Article

Figure 1. Density of CO2 + CH4 binary mixtures as a function of pressure and CO2 mole fraction at (a) 323.15 K, (b) 373.15 K, (c) 473.15 K, and (d) 573.15 K. The colored surfaces correspond to smoothed estimations from the GERG equation of state. The solid black triangles correspond to the experimental data available for binary mixtures and pure compounds. The open symbols represent molecular simulation results obtained using the TraPPE (black squares) and SAFT-γ (blue circles) force fields.

expression is consistent with unlike Mie attractive energies εij obtained using the geometric-mean rule for the van der Waals attractive constants αij,λn,ij = (αii,λn,iiαjj,λn,jj)1/2 of a hard-core inter-

By application of the standard Lorentz−Berthelot rules to combine the parameters of common cubic equations of state (i.e., the geometric-mean rule for the cross-energy a and the arithmetic-mean rule for the cross-covolume parameter b) and recognition of the proportionality relationships between microscopic and macroscopic properties (εσ3 ∝ a and σ3 ∝ b), the parameters for unlike-pair energy interactions are calculated using the following equation:46 εij = (1 − k ij)

σii 3σjj3 σij3

action model such as the Sutherland potential, where αij,λn,ij = λn,ij 2 2πεij∫ ∞ σij (σij/r) r dr. The combined exponents are used in eq 2 to give the combined constant C, and the unlike-pair interactions are finally calculated using eq 1. The force field parameters for both TraPPE and SAFT-γ are summarized in Table 1. Whereas the two single-site CH4 models studied are relatively similar, significant differences between the TraPPE and SAFT-γ CO2 force fields can be noted. In an effort to account for all of the molecular interactions in a single-site model, the values of ε and λm adopted by the SAFT-γ model are higher, since this model approximates molecules as single-site uncharged particles and the Mie parameters are therefore fitted to account also for electrostatic interactions and intramolecular energy modes that might be relevant. Although leading to unrealistically high values of ε and λm, the single-site approach can generally allow for a significant decrease in computer demand and simulation duration compared with all-atom models that account for the presence of electric charges explicitly. 2.2. Thermodynamic Properties. The density (ρ) of each pure compound and each mixture was computed in a

εiiεjj (5)

where kij is a correction term to address deviations that arise from the difference in the chemical natures of the two compounds considered. In this study, the value kij = 0.02 was adopted, according to preliminary high-pressure vapor−liquid equilibrium data regression carried out by the researchers who developed the SAFT-γ force field.43−45 The parameters λm and λn for pure compounds can also be combined to give crossinteraction exponents using the expression (λk,ij − 3) =

(λk,ii − 3)(λk,jj − 3)

(6)

where the subscript k is either m (repulsion) or n (dispersion). As pointed out recently by Lafitte and co-workers,44 this D



Article

Figure 2. Volume expansivity of CO2 + CH4 binary mixtures as a function of pressure and CO2 mole fraction at (a) 323.15 K, (b) 373.15 K, (c) 473.15 K, and (d) 573.15 K. Legend as in Figure 1.

In order to allow direct comparison with experimental data, the ideal contribution to the isobaric heat capacity was added to the reid sidual part to give the overall isobaric heat capacity: Cp = Cres p + Cp , where Cidp is defined as the sum of the ideal isobaric heat capacities of CO2 and CH4 obtained individually from experimental correlations,16,17 weighted by the numbers of molecules of each compound used to define the mixture composition in the simulations. The value of Cres p was obtained through isobaric−isothermal ensemble fluctuations according to the following expression:

straightforward manner through the ensemble average of the system volume V according to ρ=

NCH4 MWCH4 + NCO2 MWCO2 NA⟨V ⟩

(7)

where Nj is the number of molecules of compound j used in the simulation, MWj is the corresponding molecular weight, NA is the Avogadro number, and the angle brackets indicate an ensemble average. Thermodynamic properties defined as temperature or pressure derivatives, such as the volume expansivity (αp) and isothermal compressibility (κT), were evaluated by isobaric−isothermal ensemble fluctuations using the following expressions43,47−49 with proper unit conversion: αp =

κT =

⟨VH conf ⟩ − ⟨V ⟩⟨H conf ⟩ ⟨V ⟩kBT 2

(8)

⟨V 2⟩ − ⟨V ⟩2 ⟨V ⟩kBT

(9)

Cpres =

⟨U conf H conf ⟩ − ⟨U conf ⟩⟨H conf ⟩ kBT 2 +p

⟨VH conf ⟩ − ⟨V ⟩⟨H conf ⟩ − (NCH4 + NCO2)kB kBT 2 (10)

In addition, second-derivative properties, including isochoric heat capacity (CV), Joule−Thomson coefficient (μJT), and speed of sound (csound), were obtained using standard thermodynamic relationships:

in which kB is the Boltzmann constant, p and T correspond to the specified thermodynamic pressure and temperature, respectively, and Hconf = Uconf + pV, where Uconf is the potential energy of the system due to intermolecular interactions (equal to the total potential energy in this case, since no intramolecular terms are considered).

CV = Cp −

E

T ⟨V ⟩αp 2 κT

(11)

⎛ ∂T ⎞ ⟨V ⟩(Tαp − 1) μJT = ⎜ ⎟ = Cp ⎝ ∂p ⎠ H

(12)



Article

Figure 3. Isothermal compressibility of CO2 + CH4 binary mixtures as a function of pressure and CO2 mole fraction at (a) 323.15 K, (b) 373.15 K, (c) 473.15 K, and (d) 573.15 K. Legend as in Figure 1.

csound =

Cp CV κTρ

used, followed by a production run of 5 ns during which the properties of interest were recorded every 1 ps (resulting in 5001 samples for each simulation). The production simulations were divided into five equal intervals, and the standard deviations of the interval averages were used to estimate the uncertainties. Parameters such as integration time step and simulation length were chosen by conducting preliminary simulations and ensuring the reproducibility and reliability of the results for all of the force fields over the temperature, pressure, and composition ranges studied. Standard long-range corrections to the energy and pressure virial were also included, and conventional periodic boundary conditions were applied.47 The simulation length was kept the same in all cases to ensure reliable comparison between the TraPPE and SAFT-γ models. The potential cutoff was also standardized and fixed to a value of 4σ (considering the higher value for each system), and the number of molecules in each of the simulations was chosen to ensure that the box length was at least 8σ on each side. Although many force fields specify the cutoff (so technically this parameter is part of the actual force field), we decided to use a standard cutoff for all of the force fields to facilitate comparison. It should be noted that a cutoff of 4σ does not differ much from (and is never smaller than) the values used in the original force field publications. The thermodynamic states to be simulated were defined on the basis of the availability of experimental density measurements41

(13)

2.3. Simulation Details. All of the molecular simulations were performed using a time step of 1.0 fs, and the equations of motion were integrated with the standard velocity-Verlet algorithm50 implemented in the LAMMPS package.25 In the equilibration and production simulations of single-site models, temperature and pressure were specified with the Nosé−Hoover thermostat and barostat, respectively, according to the formulation given by Shinoda et al.,51 whereas the algorithm used in the simulation of the rigid three-site CO2 model was the one described by Kamberaj et al.52 The molecular configurations used to initialize the simulations were generated using the Packmol package.53 A particle-particle-particle mesh (PPPM) solver was used with a desired relative error in forces of 1 × 10−4 to calculate the long-range Coulombic interactions. In contrast to the Ewald method,47 which evaluates the Fourier series directly, the PPPM method takes advantage of computationally efficient fast Fourier transform algorithms. The short-range contribution to the nonbonded interactions (including dispersion, repulsion, and the short-range component of electrostatic interactions) is calculated by explicit particle−particle summation with (minimum-image) spherical truncation.54 The simulations were carried out using a total of 500 molecules placed in a cubic box. A 5 ns equilibration period was F



Article

Figure 4. Isobaric heat capacity of CO2 + CH4 binary mixtures as a function of pressure and CO2 mole fraction at (a) 323.15 K, (b) 373.15 K, (c) 473.15 K, and (d) 573.15 K. Legend as in Figure 1.

for binary mixtures with CO2 mole fractions between 0.10 and 0.90 at temperatures of (323.15, 373.15, 473.15, and 573.15) K and pressures of (19.94, 29.94, 39.94, 59.93, 79.93, and 99.93) MPa. Additional simulations were also performed at the intermediate pressures of (25, 35, 45, 50, 55, 65, 70, and 90) MPa, and calculations on pure CO2 and pure CH4 were included for the sake of comparison. Thus, a total of 1232 independent simulations were performed using both the TraPPE and SAFT-γ force fields as part of this study. Thermodynamic property calculations using the GERG and PR equations of state were also performed for the aforementioned conditions. Table 2 shows these equations’ functional forms along with the SAFT-VR Mie equation for comparison. The equation-of-state estimates of thermodynamic properties were obtained through the implementation of the original equations in the commercial simulator REFPROP 9.0. The mixture parameters a and b used for the PR equation of state were calculated as follows: a=

∑ ∑ xixj(1 − δij) i

b=

j

∑ xibi i

The molecular simulation results for density as a function of temperature, pressure, and mixture composition and the results obtained for the same property using the GERG and the PR equations of state were compared against the experimental data available. For all of the other thermodynamic properties, because of the lack of experimental data, the GERG results were arbitrarily chosen as the reference values against which the molecular simulation results were compared. In all cases, the average absolute relative deviation (AARD) was calculated via the expression AARD =

1 Ns

Ns

∑ i=1

Airef − Aicalc Airef

·100 % (16)

where Aref is either the experimental value for density or GERG calculation for all the other properties and Acalc is the calculated thermodynamic property obtained by molecular simulations. The summation runs over all of the states to be considered, which can include all of the states at once or may include only states found to be under some specific constant condition of temperature, pressure, or concentration in order to provide information regarding the deviation as a function of the thermodynamic state. The molecular dynamics results of observables such as volume and potential energy were used to provide optimal estimates of thermodynamic expectations using the MBAR estimator.39 For each fluid composition, the MBAR analyses

ai aj (14)

(15)

where the a and b values of pure compounds are defined elsewhere6 as functions of the critical properties and acentric factor and the binary interaction coefficient δij was set to 0.115.55 G



Article

Figure 5. Isochoric heat capacity of CO2 + CH4 binary mixtures as a function of pressure and CO2 mole fraction at (a) 323.15 K, (b) 373.15 K, (c) 473.15 K. and (d) 573.15 K. Legend as in Figure 1.

3.1. Density Validation. Molecular simulations using the TraPPE and SAFT-γ force fields along with equation-of-state calculations using the GERG and PR models were performed in order to provide discrete density results under supercritical conditions for pure CO2, pure CH4, and their mixtures. The thermodynamic states studied were chosen following values previously reported by Seitz and co-workers41 in their experimental measurements of mixture densities in order to allow direct comparison between the simulated and experimental data. Table 3 shows the AARDs between the molecular simulation and equation-of-state results and the experimental data. The first column lists the variables held constant during the averaging. For instance, the AARDs shown in the first line are the ones calculated for each force field or equation of state over the whole range of temperatures, pressures, and compositions considered in this study, while those in the second line were calculated over all of the simulated temperatures and pressures for pure CH4 only. All of the subsequent lines show similar values calculated at different compositions, pressures, and temperatures, respectively. It can be noted that the GERG equation of state provides the most accurate results under all conditions studied, including temperatures above 450 K and pressures above 35 MPa, which are above the reported range of validity of the GERG model. This is to be expected, since the experimental data considered here were the same as those used to fit the parameters of the GERG equation. Taking advantage of very accurate multiparameter

were performed by combining the simulation results obtained at all different temperatures and pressures, seeking to improve the statistical quality of the individual thermodynamic property estimates. The MBAR estimator is a direct extension of Bennett’s acceptance ratio (BAR) estimator,36 as it allows for assessment of data from different simulated state points to improve the results in a statistically optimal way and predict free energies of states were no simulation was performed. The MBAR method is also equivalent to the weighted histogram analysis method (WHAM)38 in the limit that the histogram bin widths approach zero, but it does not invoke histograms and thus provides a direct estimation of uncertainties with a modest computational expense.39 The uncorrelated data sets required for sensible MBAR estimation were ensured by subsampling the original molecular dynamics results according to the methodology proposed by Shirts and Chodera.39 All of the MBAR calculations were carried out using the pyMBAR package.56

3. RESULTS AND DISCUSSION The results of thermodynamic property calculations using molecular simulations and equations of state are shown in this section. Because of the large amount of results obtained, they are presented graphically and the source data are provided in the Supporting Information, including all of the estimated uncertainties. H



Article

Figure 6. Joule−Thomson coefficient of CO2 + CH4 binary mixtures as a function of pressure and CO2 mole fraction at (a) 323.15 K, (b) 373.15 K, (c) 473.15 K, and (d) 573.15 K. Legend as in Figure 1.

Figure 1 shows the density results for pure CO2, pure CH4, and their binary mixtures as functions of pressure and CO2 mole fraction at temperatures of (323.15, 373.15, 473.15, and 573.15) K. The results for the GERG equation and both molecular models studied are shown, whereas the PR results have been omitted for simplicity. Although generally smaller than the symbol sizes, the error bars on molecular simulation calculations are also included. The density graphs shown in Figure 1 include the values considered in the deviation calculations presented in Table 3. The uncertainties of the GERG equation for density calculations are reported to be 0.1 % inside its normal range of validity (going as high as 0.3 % in some states),21 which are consistent with our findings. Outside the model’s normal range of validity, in the region called the extended range of validity and beyond, the uncertainties of the GERG equation are roughly estimated to be between (0.5 and 1.0) %, which are also within most of the deviation values estimated for the molecular simulation results. 3.2. Thermodynamic Properties. Inspired by the excellent estimation of densities obtained with the GERG model and also because of the lack of experimental data available for other thermodynamic properties of CO2 + CH4 binary mixtures over the ranges of pressures and temperatures studied, the molecular simulation results for both the TraPPE and the SAFT-γ force fields were compared to the GERG results for volume expansivity (αp), isothermal compressibility (κT), isobaric heat capacity (Cp),

representations for both pure CO2 and CH4, this model provides deviations from experimental data that are within the uncertainties typically observed for measuring techniques of pure compounds, as pointed out by the developing group.20,21 For mixtures, the deviations found suggest an increase in the model accuracy as the CO2 mole fraction and the system pressure increase, and lower temperatures also seem to allow slightly better estimations. Generally, the molecular simulations utilizing the TraPPE and SAFT-γ force fields provided similar results that are consistently more accurate than those for the PR equation of state. It is worth noting that despite the simplicity of the coarsegrained representation of CO2 molecules employed in the SAFT-γ model, it is still able to provide results similar to those obtained using the more complex and physically meaningful TraPPE representation of CO2. Although generally similar, the force field deviations differ with respect to how the values are distributed over the conditions studied. Whereas the deviations appears to rise as the fluid becomes more concentrated in CO2 using the TraPPE representation, the composition does not seem to play a major role in the density deviations obtained with the SAFT-γ model. The gradual increase in the CH4 mole fraction helps to overcome the deviations observed for pure CO2 using the TraPPE model, though the same behavior is not easily observed with the SAFT-γ model. I



Article

Figure 7. Speed of sound for CO2 + CH4 binary mixtures as a function of pressure and CO2 mole fraction at (a) 323.15 K, (b) 373.15 K, (c) 473.15 K, and (d) 573.15 K. Legend as in Figure 1.

were also investigated for pure CO2, pure CH4, and their binary mixtures, and the results are shown in Figures 5 to 7. Substantial differences between the isochoric heat capacities estimated using the molecular models studied can be seen in Figure 5, where the TraPPE representation provides values consistently higher and closer to the GERG results than those obtained using the SAFT-γ force field. The discrepancies between the molecular models are more apparent at higher pressures and higher CO2 mole fractions, while increasing the temperature mitigates the differences. The differences between the two force fields observed for the isochoric heat capacity are related to their distinct molecular representations. In this sense, the CO2 molecular representation appears to be of major importance, not only at higher CO2 mole fractions, as expected, but also at lower CO2 mole fractions and higher pressures. The three-site CO2 TraPPE model has a rotational degree of freedom that increases its heat capacity relative to the SAFT-γ single-particle model, where all of the kinetic energy of the system is translational. This difference becomes more evident under conditions that have lower contributions of the translational kinetics to the total molecular movement, such as at lower temperatures and denser (higher pressure) conditions. Although being subject to error propagation, the results obtained for the speed of sound using both molecular models studied are in excellent agreement with those from the GERG equation (Figure 7). Despite the fact that significant differences between the Joule−Thomson coefficients calculated using

isochoric heat capacity (CV), Joule−Thomson coefficient (μJT), and speed of sound (csound). This strategy can be understood as an attempt to fill in the gaps in the availability of experimental data through the use of molecular simulation results based on two different molecular interaction models with varying levels of adherence to the real fluid structures. Moreover, the physically based calculation of thermodynamic properties provided by molecular simulations can be very useful in assessing the performance of empirical models such as GERG, especially under fluid conditions not included in the original regression of the model parameters. Figures 2 to 4 show the results for the calculated volume expansivities, isothermal compressibilities, and isobaric heat capacities, respectively. Reference values of the pure-compound properties from the National Institute of Standards and Technology (NIST)57 are also shown for comparison. It can be observed that both the TraPPE and SAFT-γ results approach the GERG results as the temperature of the simulations increases, while larger deviations appear to be concentrated in the low pressure/high CO2 mole fraction region of the figure. The molecular simulation uncertainties for results in the low pressure/high CO2 mole fraction region are also generally higher as the fluid conditions evolve toward the pure CO2 region. Overall, however, the simulation results match the values obtained from the GERG equation to a remarkable degree. Second-derivative properties, namely, the isochoric heat capacity, the Joule−Thomson coefficient, and the speed of sound, J



Article

Table 4. AARDs between Thermodynamic Properties Calculated through Molecular Simulations and Results Obtained from the GERG Equation of State AARD/% TraPPE

SAFT-γ

ρ

αp

κT

Cp

CV

μJT

csound

ρ

α

κT

Cp

Cv

μJT

csound

0.82

2.32

2.76

1.76

0.75

35.53

2.28

0.66

1.84

3.13

1.07

2.98

20.88

3.27

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

0.08 0.35 0.46 0.54 0.56 0.60 0.70 0.84 1.09 1.35 2.47

0.68 1.06 1.27 1.33 1.63 1.95 2.52 2.96 4.70 4.48 2.96

0.68 1.52 1.83 2.20 2.33 2.50 2.90 3.27 5.47 5.85 1.81

0.42 0.38 0.50 0.78 1.30 1.84 2.36 2.62 3.27 3.51 2.41

0.64 0.68 0.68 0.63 0.53 0.49 0.47 0.58 0.63 1.03 1.85

10.40 32.99 53.57 25.35 58.83 28.74 18.69 21.04 62.03 57.81 21.36

0.34 0.77 1.09 1.42 1.74 2.19 2.53 3.06 3.66 4.63 3.62

0.59 0.50 0.50 0.55 0.59 0.65 0.71 0.76 0.79 0.81 0.86

0.67 0.78 1.39 1.43 1.79 2.10 1.88 1.99 2.00 2.74 3.42

0.97 1.51 1.81 2.02 2.30 2.63 3.07 3.58 4.44 5.54 6.58

0.24 0.39 0.78 0.92 1.08 1.14 0.96 0.99 1.09 1.63 2.49

0.65 0.83 1.06 1.35 1.75 2.24 2.90 3.72 4.76 6.01 7.52

8.08 9.75 27.14 19.13 33.59 25.67 13.43 17.05 23.27 27.39 25.20

1.04 1.52 1.82 2.24 2.61 2.95 3.43 3.97 4.62 5.43 6.40

19.94 25.00 29.94 35.00 39.94 45.00 50.00 55.00 59.93 65.00 70.00 79.93 90.00 99.93

1.31 1.30 1.20 1.07 0.96 0.84 0.77 0.70 0.64 0.60 0.58 0.51 0.48 0.53

2.40 3.41 2.15 2.45 1.93 2.35 2.03 1.91 1.82 2.30 2.00 1.99 1.94 3.84

3.31 3.28 3.25 3.05 2.80 2.34 2.60 2.75 2.21 2.31 2.18 2.34 2.85 3.36

1.84 2.15 1.76 1.76 1.56 1.94 1.75 1.63 1.76 1.79 1.79 1.58 1.28 2.10

0.80 0.77 0.76 0.74 0.77 0.77 0.77 0.75 0.75 0.74 0.72 0.75 0.72 0.64

6.45 7.38 6.79 7.58 8.67 13.94 55.91 33.60 135.28 84.98 42.23 61.63 19.13 13.80

2.22 2.16 2.20 2.40 2.29 2.11 2.19 2.35 2.26 2.25 2.26 2.27 2.43 2.49

0.98 0.94 0.84 0.75 0.64 0.57 0.51 0.47 0.45 0.45 0.48 0.61 0.75 0.88

2.84 2.44 2.29 1.98 1.99 1.62 1.63 1.65 1.56 1.35 1.39 1.68 1.64 1.63

2.28 2.03 2.40 2.81 2.94 3.03 3.10 3.18 3.40 3.53 3.55 3.66 3.85 4.07

1.44 1.27 1.13 0.96 1.04 0.91 0.94 1.00 0.98 0.92 0.97 1.10 1.10 1.15

1.91 2.17 2.40 2.60 2.76 2.87 3.00 3.12 3.21 3.31 3.39 3.53 3.64 3.79

4.94 4.85 5.39 6.92 8.49 11.91 36.10 28.77 53.04 41.23 32.73 31.83 18.67 7.49

1.48 1.91 2.31 2.70 2.93 3.18 3.34 3.49 3.69 3.85 3.96 4.10 4.37 4.53

323.15 373.15 473.15 573.15

1.11 0.76 0.64 0.77

2.74 2.60 1.99 1.96

4.73 2.79 1.71 1.82

2.99 2.28 1.06 0.71

1.12 1.03 0.59 0.24

21.09 43.03 44.63 33.37

3.91 2.58 1.41 1.21

0.86 0.64 0.45 0.70

2.15 2.10 1.68 1.40

6.32 3.65 1.49 1.07

1.92 1.29 0.58 0.48

5.38 3.78 1.85 0.91

15.08 17.24 22.42 28.79

5.91 3.81 1.96 1.41

global xCO2

p/MPa

T/K

molecular simulations and the GERG equation can be observed in Figure 6, it is still possible to note a good agreement among all of the methods in the description of the Joule−Thomson inversion curves under all conditions, which also agrees well to the findings of Vrabec and co-workers.58 Table 4 shows a thorough comparison of the observed AARDs for all of the property estimation methods studied. The GERG equation was taken as the reference method, and the AARDs of each molecular model against the GERG estimates are shown for all of the properties, following the same structure previously adopted in Table 3. The differences between the density values calculated through molecular simulations and the GERG equation are comparable to the deviations observed between the experimental values and the GERG equation in its extended range of validity. In regard to the isochoric heat capacity, the TraPPE model provides results in better agreement with those calculated using the GERG equation over the full range of compositions. In this case, the deviation values are within the uncertainties estimated for experimental measurements of binary and multicomponent mixtures.21 For the same property, the SAFT-γ model shows a decrease in accuracy as the CO2

concentration increases, and the agreement with the GERG results at CO2 mole fractions above 0.5 is only moderate. The deviations in the speed of sound calculations are generally higher than those reported for gas-phase experimental measurements [ranging from (0.02 to 0.1) % at temperatures between (250 and 350) K and pressures up to 12 MPa21]. Increasing the temperature appears to allow better agreement between the molecular simulations and the GERG results, whereas mixtures with higher CO2 mole fractions show an increase in the deviations. In addition, for the SAFT-γ model, the relation between the deviations observed and the pressure is almost a direct proportion, while for the TraPPE model the deviation is not significantly affected by differences in pressure. In both cases, the monotonic behavior of the deviations considering the full range of conditions studied suggests that the semiempirical GERG approach and the physically meaningful representation employed in the molecular simulations maintain approximately the same level of agreement even in the high-pressure and -temperature regions outside the range of validity of the GERG equation. The larger differences in the Joule−Thomson coefficients calculated using molecular simulation results and the GERG equation can be attributed to the high sensitivity of individual K



Article

values of this property calculated through molecular simulation and eq 12, where small variations in properties such as the volume expansivity can lead to very large uncertainties in individual Joule−Thomson coefficient values. This is a well-known difficulty about which other authors have commented.59−63 Although not quantitatively precise for discrete values, the good representation of the overall qualitative results, including the Joule−Thomson inversion curves, is remarkable and can be used to stress the good adherence between the GERG semiempirical model and physically meaningful molecular models. 3.3. Computational Cost. Besides being generally more physically meaningful, the atomistic approach usually relies upon a smaller number of parameters to describe the thermodynamic properties of various pure compounds and mixtures, even when compared with simple equations of state. Even though the continuous advances in computational resources and methods have been trimming down simulation times, the usually long duration of molecular simulations, especially compared with the almost on-the-fly calculations of properties provided by equations of state and multiparameter models, suggests that molecular simulations will not supplant these equations when many rapid evaluations are required, such as in the use of a process simulator. However, compared with the time and expense required to measure thermodynamic properties such as those reported here, simulations are significantly faster and less costly. This suggests that once verified against a limited amount of experimental data, simulations can be used to augment the data to enable the development of better property models. In addition, the time and cost associated with simulations will only continue to decrease with time. The coarse-grained representation of either simple or complex fluids appears as a promising and valuable option to help bring molecular simulations down to the level of a routine calculation of use to industry. For instance, the property calculations presented in this work show that similar (and sometimes even more accurate) results can be achieved by the utilization of a simplified representation of molecules that enables a drastic decrease in computational resources and simulation times. Table 5 shows a comparison of the average simulation times required to perform the thermodynamic property calculations previously shown for the TraPPE and the SAFT-γ models. The global value corresponds to the average over all of the simulations performed. The subsequent lines show the averages taken under specific conditions of constant composition, pressure, or temperature, following a structure similar to that presented in Tables 3 and 4. For systems consisting of pure CH4, both TraPPE and SAFT-γ are expected to give similarly low simulation times, as only interactions of uncharged single particles are considered. However, the implementation of the Lennard-Jones potential in LAMMPS takes advantage of the factor of 2 between the repulsive and dispersive exponents to provide very efficient calculations. To account for real numbers that can arise from the parametrization of the exponents, the Mie implementation relies upon less-optimized mathematical operations, leading to longer simulation times, as can be noted in Table 5. The superiority of the Lennard-Jones potential found in the pure CH4 simulations immediately vanishes as the CO2 mole fraction increases to values as low as 0.1. In the TraPPE representation, the presence of CO2 molecules adds three new interaction sites for each molecule and also requires the expensive calculation of electrostatic interactions, while the SAFT-γ model describes the CO2 in a very similar way to CH4. Hence, the simulation times become almost insensitive to the mixture composition.

Table 5. Comparison of Average Simulation Times for the TraPPE and SAFT-γ Representations of CO2 and CH4 Pure Compounds and Binary Mixturesa simulation time (h) TraPPE

SAFT-γ

ratio

2.36

0.91

2.59

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

0.39 1.11 1.42 1.74 2.08 2.46 2.79 3.11 3.51 3.93 3.45

0.84 0.84 0.85 0.86 0.88 0.90 0.92 0.94 0.97 0.99 1.04

0.47 1.33 1.68 2.02 2.36 2.74 3.05 3.31 3.61 3.95 3.32

19.94 25.00 29.94 35.00 39.94 45.00 50.00 55.00 59.93 65.00 70.00 79.93 90.00 99.93

1.95 2.07 2.12 2.20 2.27 2.32 2.36 2.42 2.49 2.48 2.51 2.57 2.65 2.67

0.61 0.68 0.75 0.80 0.84 0.88 0.92 0.95 0.97 1.00 1.03 1.07 1.11 1.14

3.21 3.03 2.84 2.74 2.70 2.63 2.57 2.55 2.56 2.48 2.44 2.40 2.38 2.34

323.15 373.15 473.15 573.15

2.73 2.50 2.19 2.03

1.14 1.00 0.81 0.70

2.40 2.50 2.70 2.92

global xCO2

p/MPa

T/K

a

The simulations were performed on a dual eight-core 2.3 GHz AMD Opteron (24 GB RAM) cluster.

The system density is also important in the simulation time evaluation, as it determines the number of neighbors of each molecule and therefore is related to the computational effort required to calculate the molecular interactions at each time step. As shown in Table 5, upon an increase in the pressure or, equivalently, a decrease in the temperature of the molecular system, the density increase leads to longer simulation times for both the TraPPE and SAFT-γ force fields. It should be noted that every simulation reported here was independent, so in principle all of the calculations could be carried out at the same time if enough parallel computing resources are available.

4. CONCLUSIONS Extensive molecular dynamics simulations were carried out to compute thermodynamic properties of CO2, CH4, and their mixtures over a wide range of temperatures and pressures. Carbon dioxide was modeled using two different force fields, the popular three-site TraPPE model as well as a single-site Mie potential parametrized using the SAFT-γ approach. Methane was modeled using the TraPPE single-site Lennard-Jones model and a single-site SAFT-γ Mie potential. The simulations L



Article

and the mixing rule parameters with CO2 prior to publication. We are also grateful for the computational resources provided by the Center for Research Computing (CRC) at the University of Notre Dame.

were compared against experimental density data as well as results from the multiparameter GERG-2008 equation and the Peng−Robinson equation of state. The simulations and the GERG equation all modeled the pure-compound and mixture densities with high accuracy, while the PR equation of state performed the worst, especially at higher pressures. Other properties, including volume expansivity (αp), isothermal compressibility (κT), isobaric heat capacity (Cp), isochoric heat capacity (CV), Joule−Thomson coefficient (μJT), and speed of sound (csound) were computed using the two different force fields. Because no experimental data for these properties for mixtures were found in the literature, comparison was made with results obtained using the GERG equation. With the exception of the Joule−Thomson coefficient, both force fields yielded overall property estimates within (1 to 3) % of those obtained from the GERG equation. The deviation for the Joule−Thomson coefficient, a property notoriously difficult to compute because of error propagation issues, was on the order of (20 to 35) %. The results suggest that one can confidently use the GERG equation or either the TraPPE or SAFT-γ molecular model to predict thermodynamic properties of CO2 and CH4 mixtures over the range of conditions considered here, including the regions outside the normal range of validity of the GERG equation of state. The single-site representation of CO2 in the SAFT-γ force field was able to give results similar to those for the charged three-site CO2 model used in the TraPPE force field while allowing a very significant reduction in simulation time. The reduction in computational demand is an imminent need that must be met in order to allow molecular simulations to become a feasible and reliable tool available for daily industrial purposes. The multistate Bennett acceptance ratio (MBAR) method was used to provide a statistically optimal estimation of equilibrium observables, improving the accuracy and reducing the calculation uncertainties. The application of the MBAR technique in thermodynamic property calculations offers a number of benefits ranging from the reduction of the number of simulations required to the calculation of properties over a quasicontinuous range of temperatures, thus saving resources and providing fast results. In this sense, the MBAR technique can be viewed as a numerical equation of state that relies upon a predetermined set of individual simulations in the surroundings of some thermodynamic state of interest in a given ensemble to provide thermodynamic property estimations over a welldefined range of conditions and uncertainties known a priori.

■

■

(1) Oldenburg, C. M.; Pruess, K.; Benson, S. M. Process Modeling of CO2 Injection into Natural Gas Reservoirs for Carbon Sequestration and Enhanced Gas Recovery. Energy Fuels 2001, 15, 293−298. (2) Class, H.; Ebigbo, A.; Helmig, R.; Dahle, H. K.; Nordbotten, J. M.; Celia, M. A.; Audigane, P.; Darcis, M.; Ennis-King, J.; Fan, Y.; Flemisch, B.; Gasda, S. E.; Jin, M.; Krug, S.; Labregere, D.; Naderi Beni, A.; Pawar, R. J.; Sbai, A.; Thomas, S. G.; Trenty, L.; Wei, L. A Benchmark Study on Problems Related to CO2 Storage in Geologic Formations. Comput. Geosci. 2009, 13, 409−434. (3) Kvenvolden, K. A. Gas HydratesGeological Perspective and Global Change. Rev. Geophys. 1993, 31, 173−187. (4) Deusner, C.; Bigalke, N.; Kossel, E.; Haeckel, M. Methane Production from Gas Hydrate Deposits through Injection of Supercritical CO2. Energies 2012, 5, 2112−2140. (5) Ungerer, P.; Tavitian, B.; Boutin, A. Applications of Molecular Simulation in the Oil and Gas Industry: Monte Carlo Methods; Editions Technip: Paris, 2005. (6) Peng, D. Y.; Robinson, D. B. A New Two-Constant Equation of State. Ind. Eng. Chem. Fundam. 1976, 15, 59−64. (7) Soave, G. Equilibrium Constants from a Modified Redlich− Kwong Equation of State. Chem. Eng. Sci. 1972, 27, 1197−1203. (8) Vera, J. H.; Huron, M. J.; Vidal, J. On the Flexibility and Limitations of Cubic Equations of State. Chem. Eng. Commun. 1984, 26, 311−318. (9) Chapman, W. G.; Gubbins, K. E.; Jackson, G.; Radosz, M. New Reference Equation of State for Associating Liquids. Ind. Eng. Chem. Res. 1990, 29, 1709−1721. (10) Huang, S. H.; Radosz, M. Equation of State for Small, Large, Polydisperse, and Associating Molecules. Ind. Eng. Chem. Res. 1990, 29, 2284−2294. (11) Huang, S. H.; Radosz, M. Equation of State for Small, Large, Polydisperse, and Associating Molecules: Extension to Fluid Mixtures. Ind. Eng. Chem. Res. 1991, 30, 1994−2005. (12) Gil-Villegas, A.; Galindo, A.; Whitehead, P. J.; Mills, S. J.; Jackson, G.; Burgess, A. N. Statistical Associating Fluid Theory for Chain Molecules with Attractive Potentials of Variable Range. J. Chem. Phys. 1997, 106, 4168−4186. (13) Jog, P. K.; Garcia-Cuellar, A.; Chapman, W. G. Extensions and Applications of the SAFT Equation of State to Solvents, Monomers, and Polymers. Fluid Phase Equilib. 1999, 158, 321−326. (14) Gross, J.; Sadowski, G. Perturbed-Chain SAFT: An Equation of State Based on a Perturbation Theory for Chain Molecules. Ind. Eng. Chem. Res. 2001, 40, 1244−1260. (15) Müller, E. A.; Gubbins, K. E. Molecular-Based Equations of State for Associating Fluids: A Review of SAFT and Related Approaches. Ind. Eng. Chem. Res. 2001, 40, 2193−2211. (16) Setzmann, U.; Wagner, W. A New Equation of State and Tables of Thermodynamic Properties for Methane Covering the Range from the Melting Line to 625 K at Pressures up to 1000 MPa. J. Phys. Chem. Ref. Data 1991, 20, 1061−1155. (17) Span, R.; Wagner, W. A New Equation of State for Carbon Dioxide Covering the Fluid Region from the Triple-Point Temperature to 1100 K at Pressures up to 800 MPa. J. Phys. Chem. Ref. Data 1996, 25, 1509−1596. (18) Lemmon, E. W.; Jacobsen, R. T. A Generalized Model for the Thermodynamic Properties of Mixtures. Int. J. Thermophys. 1999, 20, 825−835. (19) Span, R.; Wagner, W.; Lemmon, E. W.; Jacobsen, R. T. Multiparameter Equations of StateRecent Trends and Future Challenges. Fluid Phase Equilib. 2001, 183, 1−20. (20) Kunz, O.; Klimeck, R.; Wagner, W.; Jaeschke, M. The GERG2004 Wide Range Equation of State for Natural Gases and Other

ASSOCIATED CONTENT

* Supporting Information S

Source data used to generate the figures. This material is available free of charge via the Internet at http://pubs.acs.org.

■

REFERENCES

AUTHOR INFORMATION

Corresponding Author

*Phone: +1 574-631-5687. Fax: +1 574-631-8366. E-mail: ed@ nd.edu. Funding

This work was supported by Petrobras. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We thank Dr. Erich Müller (Imperial College London) for providing us with the CH4 SAFT-γ Mie force field parameters M



Article

Mixtures; GERG Technical Monograph 15; VDI Verlag: Düsseldorf, Germany, 2007. (21) Kunz, O.; Wagner, W. The GERG-2008 Wide-Range Equation of State for Natural Gases and Other Mixtures: An Expansion of GERG-2004. J. Chem. Eng. Data 2012, 57, 3032−3091. (22) Melo, C. L.; Thedy, E. A.; Rocha, P. S.; Almeida, A. S.; Musse, A. P. The Challenges on the CCGS Monitoring in the Development of Santos Basin Pre-salt Cluster. Energy Procedia 2011, 4, 3394−3398. (23) Anderson, J. A.; Jankowski, E.; Grubb, T. L.; Engel, M.; Glotzer, S. C. Massively Parallel Monte Carlo for Many-Particle Simulations on GPUs. J. Comput. Phys. 2013, 254, 27−38. (24) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435−447. (25) Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys. 1995, 117, 1−19. (26) Fincham, D. Parallel Computers and Molecular Simulation. Mol. Simul. 1987, 1, 1−45. (27) Frenkel, D.; Smit, B. Understanding Molecular Simulation; Academic Press: New York, 2002. (28) McQuarrie, D. A. Statistical Mechanics; University Science Books: Sausalito, CA, 2000. (29) Jorgensen, W. L.; Madura, J. D.; Swenson, C. J. Optimized Intermolecular Potential Functions for Liquid Hydrocarbons. J. Am. Chem. Soc. 1984, 106, 6638−6646. (30) MacKerell, A. D.; Banavali, N.; Foloppe, N. Development and Current Status of the CHARMM Force Field for Nucleic Acids. Biopolymers 2001, 56, 257−265. (31) Yang, J.; Ren, Y.; Tian, A.; Sun, H. COMPASS Force Field for 14 Inorganic Molecules, He, Ne, Ar, Kr, Xe, H2, O2, N2, NO, CO, CO2, NO2, CS2, and SO2, in Liquid Phases. J. Phys. Chem. B 2000, 104, 4951−4957. (32) Potoff, J. J.; Siepmann, J. I. Vapor−Liquid Equilibria of Mixtures Containing Alkanes, Carbon Dioxide, and Nitrogen. AIChE J. 2001, 47, 1676−1682. (33) Panagiotopoulos, A. Z. Force-Field Development for Simulations of Condensed Phases. AIChE Symp. Ser. 2001, 97, 61−70. (34) Martin, M. G.; Siepmann, J. I. Transferable Potentials for Phase Equilibria. 1. United-Atom Description of n-Alkanes. J. Phys. Chem. B 1998, 102, 2569−2577. (35) Maerzke, K. A.; Siepmann, J. I. Transferable Potentials for Phase EquilibriaCoarse-Grain Description for Linear Alkanes. J. Phys. Chem. B 2011, 115, 3452−3465. (36) Bennett, C. H. Efficient Estimation of Free Energy Differences from Monte Carlo Data. J. Comput. Chem. 1976, 22, 245−268. (37) Ferrenberg, A. M.; Swendsen, R. H. Optimized Monte Carlo Data Analysis. Phys. Rev. Lett. 1989, 63, 1195−1198. (38) Kumar, S.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A.; Rosenberg, J. M. The Weighted Histogram Analysis Method for FreeEnergy Calculations on Biomolecules. I. The Method. J. Comput. Chem. 1992, 13, 1011−1021. (39) Shirts, M. R.; Chodera, J. D. Statistically Optimal Analysis of Samples from Multiple Equilibrium States. J. Chem. Phys. 2008, 129, No. 124105. (40) Paliwal, H.; Shirts, M. R. Multistate Reweighting and Configuration Mapping Together Accelerate the Efficiency of Thermodynamic Calculations as a Function of Molecular Geometry by Orders of Magnitude. J. Chem. Phys. 2013, 138, No. 154108. (41) Seitz, J. C.; Blencoe, J. G.; Bodnar, R. J. Volumetric Properties for {(1−x)CO2 + xCH4}, {(1−x)CO2 + xN2}, and {(1−x)CH4 + xN2} at the Pressures (9.94, 19.94, 29.94, 39.94, 59.93, 79.93, and 99.93) MPa and Temperatures (323.15, 373.15, 473.15, and 573.15) K. J. Chem. Thermodyn. 1996, 28, 521−538. (42) Mie, G. Zur kinetischen Theorie der einatomigen. Ann. Phys. 1903, 11, 657−697. (43) Avendaño, C.; Lafitte, T.; Galindo, A.; Adjiman, C. S.; Jackson, G.; Müller, E. A. SAFT-γ Force Field for the Simulation of Molecular Fluids. 1. A Single-Site Coarse Grained Model of Carbon Dioxide. J. Phys. Chem. B 2011, 115, 11154−11169.

(44) Lafitte, T.; Apostolakou, A.; Avendaño, C.; Galindo, A.; Adjiman, C. S.; Müller, E. A.; Jackson, G. Accurate Statistical Associating Fluid Theory for Chain Molecules Formed from Mie Segments. J. Chem. Phys. 2013, 139, No. 154504. (45) Avendaño, C.; Lafitte, T.; Adjiman, C. S.; Galindo, A.; Müller, E. A.; Jackson, G. SAFT-γ Force Field for the Simulation of Molecular Fluids: 2. Coarse-Grained Models of Greenhouse Gases, Refrigerants, and Long Alkanes. J. Phys. Chem. B 2013, 117, 2717−2733. (46) Kontogeorgis, G. M.; Folas, G. K. Thermodynamic Models for Industrial Applications: From Classical and Advanced Mixing Rules to Association Theories; John Wiley & Sons: New York, 2010. (47) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, U.K., 1987. (48) Lagache, M.; Ungerer, P.; Boutin, A.; Fuchs, A. H. Prediction of Thermodynamic Derivative Properties of Fluids by Monte Carlo Simulation. Phys. Chem. Chem. Phys. 2001, 3, 4333−4339. (49) Colina, C. M.; Olivera-Fuentes, C. G.; Siperstein, F. R.; Lísal, M.; Gubbins, K. E. Thermal Properties of Supercritical Carbon Dioxide by Monte Carlo Simulations. Mol. Simul. 2003, 29, 405−412. (50) Verlet, L. Computer “Experiments” on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones Molecules. Phys. Rev. 1967, 159, 98−103. (51) Shinoda, W.; Shiga, M.; Mikami, M. Rapid Estimation of Elastic Constants by Molecular Dynamics Simulation under Constant Stress. Phys. Rev. B 2004, 69, No. 134103. (52) Kamberaj, H.; Low, R. J.; Neal, M. P. Time Reversible and Symplectic Integrators for Molecular Dynamics Simulations of Rigid Molecules. J. Chem. Phys. 2005, 122, No. 224114. (53) 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. (54) Gargallo, R.; Hünenberger, P. H.; Avilés, F. X.; Oliva, B. Molecular Dynamics Simulation of Highly Charged Proteins: Comparison of the Particle-Particle Particle-Mesh and Reaction Field Methods for the Calculation of Electrostatic Interactions. Protein Sci. 2003, 12, 2161−2172. (55) Nishiumi, H.; Arai, T.; Takeuchi, K. Generalization of the Binary Interaction Parameter of the Peng−Robinson Equation of State by Component Family. Fluid Phase Equilib. 1988, 42, 43−62. (56) Shirts, M. R.; Chodera, J. D. A Python Implementation of the Multistate Bennett Acceptance Ratio (MBAR). https://simtk.org/ home/pymbar (accessed Feb 22, 2013). (57) NIST Chemistry WebBook; Linstrom, P. J., Mallard, W. G., Eds.; NIST Standard Reference Database Number 69; National Institute of Standards and Technology: Gaithersburg, MD, 2005; http://webbook. nist.gov (accessed Sept 17, 2013). (58) Vrabec, J.; Kumar, A.; Hasse, H. Joule−Thomson Inversion Curves of Mixtures by Molecular Simulation in Comparison to Advanced Equations of State: Natural Gas as an Example. Fluid Phase Equilib. 2007, 258, 34−40. (59) Colina, C. M.; Müller, E. A. Molecular Simulation of Joule− Thomson Inversion Curves. Int. J. Thermophys. 1999, 20, 229−235. (60) Chacín, A.; Vázquez, J. M.; Müller, E. A. Molecular Simulation of the Joule−Thomson Inversion Curve of Carbon Dioxide. Fluid Phase Equilib. 1999, 165, 147−155. (61) Colina, C. M.; Lísal, M.; Siperstein, F. R.; Gubbins, K. E. Accurate CO2 Joule−Thomson Inversion Curve by Molecular Simulations. Fluid Phase Equilib. 2002, 202, 253−262. (62) Kioupis, L. I.; Maginn, E. J. Pressure−Enthalpy Driven Molecular Dynamics for Thermodynamic Property Calculation: I. Methodology. Fluid Phase Equilib. 2002, 200, 75−92. (63) Kioupis, L. I.; Arya, G.; Maginn, E. J. Pressure−Enthalpy Driven Molecular Dynamics for Thermodynamic Property Calculation II: Applications. Fluid Phase Equilib. 2002, 200, 93−110.

N


Thermodynamic Properties of Supercritical Mixtures of Carbon Dioxide

Recommend Documents