CH4 Mixtures from

Article pubs.acs.org/jced

Thermodynamic Properties of Supercritical CO2/CH4 Mixtures from the Virial Equation of State Shu Yang, Andrew J. Schultz, and David A. Kofke* Department of Chemical and Biological Engineering, University at Buffalo, The State University of New York, Buffalo, New York 14260-4200, United States S Supporting Information *

ABSTRACT: Mixture virial coefficients to seventh order are presented for the system CO2/CH4 at four supercritical temperatures: 323.15, 373.15, 473.15, and 573.15 K. Values are evaluated via the Mayer sampling Monte Carlo method using a three-site TraPPE model for CO2 and a one-site model for CH4. The coefficients are used to compute seven thermodynamic properties (viz., compressibility factor, isothermal compressibility, volume expansivity, isochoric and isobaric heat capacities, Joule−Thomson coefficient, and speed of sound) as a function of mole fraction and density for these temperatures. Comparison is made with corresponding data in the literature as obtained by molecular dynamics simulation, covering densities up to about twice the critical density. Key conclusions are as follows, noting that some exceptions are observed in each case: (a) The virial equation of state (VEOS) to fourth or fifth order describes all properties to within the simulation uncertainty for densities up to at least the critical density, and the addition of terms up to seventh order extends this range considerably. (b) The accuracy of the VEOS is severely diminished for conditions approaching the critical point (the present work extends down to a reduced temperature of 1.06 for CO2), and the study of the pure component behavior suggests the critical singularity blocks convergence for conditions at considerably higher temperatures, albeit at correspondingly higher pressures. (c) Comparison of the VEOS at different orders provides a reliable guide to its accuracy at a given order, so the VEOS can provide a self-assessment of its accuracy when independent data for comparison are unavailable. (d) The VEOS provides a good description of the Joule−Thomson coefficient, including the inversion point in particular. The third-order series is needed to obtain behavior that is qualitatively correct, and the addition of higher-order terms steadily improves the accuracy quantitatively. (e) Under conditions where the seventh-order series is converged, properties can be computed to a given precision with VEOS using much less computational effort in comparison to molecular simulation.

■

INTRODUCTION In engineering applications, it is common practice to represent thermophysical properties via explicit equations, which can quickly yield property values for arbitrary thermodynamic states. Often such models are based on a common form that is parametrized against experimental data taken from the laboratory. The ability of such models to describe data under conditions between or beyond the experimental data depends on the amount of data used for the fit and the correctness or rigor of the form of the model. Most technological applications involve mixtures, so it is here that the need for models is greatest. Mixtures are particularly challenging, however, because (1) each added component introduces a new degree of freedom and therefore an exponential increase in the amount of data needed to fully characterize behavior and (2) the dependence of properties on mixture composition is nontrivial and difficult to formulate correctly in a mathematical model. Improvements in molecular modeling and computer-based data-generation methods are leading to the emergence of computational methods as a surrogate for laboratory-based data. Molecular simulation is playing a key role in this development, © XXXX American Chemical Society

which can be used to generate thermophysical property data with accuracy that is limited only by the quality of the underlying molecular models on which the simulation is based. Of course, the formulation of a good-quality molecular model is a key obstacle in this process, but steady progress is being made through the development of force fields that are fit to experimental data and/or the results of accurate ab initio computational chemistry calculations. Although this pushes the fitting problem down to the molecular scale, this is nevertheless viewed as progress because models formulated at this level should be more broadly transferable across properties and state conditions. Thus, the availability of accurate computationally based properties helps to address the need for data, but still it does not otherwise advance the problem of formulating accurate models for mixture properties. Special Issue: Proceedings of PPEPPD 2016 Received: August 5, 2016 Accepted: October 7, 2016

A

DOI: 10.1021/acs.jced.6b00702 J. Chem. Eng. Data XXXX, XXX, XXX−XXX

Journal of Chemical & Engineering Data

Article

The virial equation of state (VEOS)1 holds a position intermediate between molecular simulation on the one hand and parametrized thermophysical property models on the other. The VEOS expresses the pressure p as a power series in molar density ρ and may be written in terms of the compressibility factor Z p Z≡ = 1 + B2 ρ + B3ρ2 + ... = ∑ Bnρn − 1 ρRT (1a) n=1

This picture is changing because of advances in methods for computing virial coefficients.9−12 These developments, coupled with the increasing prospect for molecular models to provide quantitative predictions, make it appropriate to examine more closely the capabilities of the VEOS for describing thermophysical properties, particularly for mixtures. Studies have appeared recently13−16 that consider this issue, but their attention is largely restricted to low-order coefficients. A notable exception is found in the work of Jäger,17 who computed a broad set of thermodynamic properties for argon based on the seventh-order VEOS; this contribution is particularly noteworthy because all coefficients used by it were calculated from ab initio-based potentials,18 including multibody and nuclear quantum effects. Moreover, the computed properties were found to agree well with data for argon as obtained experimentally. In the present work, we also examine the performance of the VEOS when applied at higher orders, also considering coefficients up to seventh order in density. Specifically, we examine the CO2(1)/CH4(2) mixture (CAS registry numbers 124-38-9 and 74-82-8, respectively). These gases find many uses in industrial chemical, petroleum, and food processes, and they play a key role in many problems relating to energy and the environment. Most important to the present context, Aimoli et al.19 recently reported results from molecular dynamics (MD) simulations of a realistic model (TraPPE-UA) of CO2/CH4 mixtures, providing values for a comprehensive set of thermodynamic properties over a wide range of states and mixture compositions in the supercritical region. These data are ideally suited as a benchmark for gauging the performance of the VEOS, which can be developed for the same molecular models used for the MD study. In separate work,20 this group performed a comparative study for the CO2/CH4 system using a variety of molecular models and concluded that the ones providing the focus of their original study, which are the ones we examine here, are among the best in describing experimental data. The remainder of this article is organized as follows. In the next section, we describe the molecular models and methods used for the calculation of the virial coefficients. Then we write formulas derived from the VEOS for the full set of properties reported by Aimoli et al.19 Following that, we report our results and provide an analysis of the VEOS based upon them, giving attention as well to the efficiency of the VEOS relative to molecular simulation as a means to compute properties from a molecular model. The final section provides conclusions and guidelines for further study.

where R is the gas constant, T is the temperature, and Bn is the nth-order virial coefficient, defining B1 ≡ 1. We sometimes use VEOSn to indicate the application of the virial series including coefficients up to order n. The virial coefficients are functions of temperature and, for mixtures, the mole fractions yi. The composition dependence is expressed in a simple analytical form, and in particular for a binary mixture it is given as B2 = B20 y12 + 2B11y1y2 + B02 y22 B3 = B30 y13 + 3B21y12 y2 + 3B12 y1y22 + B03y23

or for arbitrary n, n

Bn =

⎛n⎞

∑ ⎜⎝ k ⎟⎠B(n− k)k y1n− k y2k k=0

(1b)

The Bij coefficients are functions only of temperature and characterize the interactions of i species-1 molecules and j species-2 molecules. As with other property models, values for these coefficients can be obtained by the appropriate analysis of experimental or simulation-based PVTy data, but this becomes extremely difficult to accomplish well for n ≥ 4. However, unlike other thermodynamic models, the parameters of the VEOS (i.e., the virial coefficients) can alternatively be computed directly from a given molecular model. Indeed, the VEOS is unique among engineering equations of state in its ability to determine its parameters rigorously from molecular considerations. In this respect, the VEOS is very much like molecular simulation, which provides accurate thermodynamic properties for a given molecular model. But unlike molecular simulation, the VEOS is an equation of state, and as such it may be manipulated to yield other thermodynamic properties of interest. Moreover, the dependence of the VEOS on mixture composition (yi) is rigorous; it is exactly the correct form for a given order in density. Despite these appealing features and despite its long history, the VEOS has seen relatively little application. Part of the reason for this is the necessarily limited range of validity of the VEOS; it is not applicable to condensed phases because any series-based approach has difficulty traversing the mathematical singularities relating to condensation (although recently there have been interesting developments in connection to this2−8). The real crux of the problem, though, is the difficulty of actually computing virial coefficients from a molecular model using the formal expressions that relate them. Throughout much of its history, virial coefficients have been computed only up to second or maybe third order, and coefficients computed this way for realistic molecular models (not to mention mixtures thereof) have been particularly rare. Without values for the coefficients, it is not possible to study the convergence of the VEOS so as to understand the conditions under which it can be applied. Consequently, lacking evidence-based knowledge of its strengths and limitations (and lacking coefficient values), thermodynamic modelers have been reluctant to consider the VEOS for their applications.

■

FORMALISM AND METHODS Properties. Aimoli et al.19 measured a large set of thermodynamic properties in their MD simulations. An expression for each property can be derived from the VEOS, eq 1, and in this section, we summarize the relevant formulas. These may be derived from standard thermodynamic manipulations on the VEOS or, in the case of eqs 5−7, are available from many thermodynamics texts. Isothermal Compressibility, κT. κT ≡ −

1 ⎛ ∂V ⎞ ⎜ ⎟ V ⎝ ∂p ⎠

T

= ⎡⎣RT ∑n = 1 nBnρn ⎤⎦−1 = ⎣⎡RT (ρ + 2B2 ρ2 + 3B3ρ3 + ...)⎦⎤−1

(2)

Volume Expansivity, αp. B


Journal of Chemical & Engineering Data αp ≡

Article

1 ⎛⎜ ∂V ⎞⎟ V ⎝ ∂T ⎠ p

=

∑n = 1 (Bn + TBn′)ρn − 1 T ∑n = 1 nBnρn − 1

=

1 + (Bn + TB2′ )ρ + (B3 + TB3′)ρ2 + ... T (1 + 2B2 ρ + 3B3ρ2 + ...)

expansion results in a temperature increase versus a decrease. The behavior of μJT and the inversion curve in particular have received some attention in the context of the VEOS.23−26 Speed of Sound, u. The speed of sound relates to the adiabatic compressibility, which in turn is given in terms of the quantities above u2 =

(3)

where B′n ≡ dBn/dT, which is computed from the temperature derivatives Bij′ (T) via the differentiation of eq 1b. Molar Isochoric Heat Capacity, cV. 1 (2TBn′ n−1

cV = cVid(T ) − ∑n = 2

1 (2TB3′ + T 2B3′′)ρ2 − ... 2

(4a)

where B″n ≡ d Bn/dT , which is again given in terms of B″ij (T) via eq 1b. Also, cidV is the ideal gas contribution to the heat capacity, which of course is obtained independently of the VEOS. For this contribution, we use the same correlations employed by Aimoli et al.19 Specifically, for CO2 the contribution is given by the correlation of Span and Wagner21 2

cVid,CO2(T ) R

2

= 2.5 +

23.1996e6.1119/ Tr 18.9205e6.77708/ Tr + 6.77708/ T 6.1119/ Tr 2 2 r − 1)2 T 2 (e (e − 1) Tr r

+

61.1048e 27.0879/ Tr 133.396e11.3238/ Tr + 27.0879/ T 11.3238/ Tr 2 2 r − 1)2 T 2 (e (e − 1) Tr r

+

19.8086e3.15163/ Tr (e3.15163/ Tr − 1)2 Tr 2

(4b)

where Tr = T/Tc, with Tc = 304.25 K. For CH4, we use the formula from Setzmann and Wagner:22 cVid,CH4(T ) R +

= 3.0016 +

1.79781 × 107e1957/ T 5.28938 × 107e3895/ T + 1957/ T 2 2 − 1) T (e (e3895/ T − 1)2 T 2

5.39369 × 107e5705/ T 3.20984 × 108e15080/ T + 5705/ T 2 2 − 1) T (e (e15080/ T − 1)2 T 2

3547.77e648/ T + 648/ T − 1)2 T 2 (e

Bijfit (T ) = aij ,0 +

(5)

Joule−Thomson Coefficient, μJT. The Joule−Thomson coefficient is key to the description of throttling processes, such as those at the core of refrigeration systems. There is particular interest in this property for mixtures, as part of ongoing efforts to formulate new refrigerants by blending different substances. The Joule−Thomson coefficient can be evaluated in terms of the properties given above, via the formula ⎛ ∂T ⎞ Tαp − 1 μJT ≡ ⎜ ⎟ = ρcp ⎝ ∂p ⎠ H

aij ,1 T

+

aij ,2 T

2

+

aij ,3 T3

(8)

As discussed below, we evaluate each Bij(T) coefficient for four temperatures. These data are computed using stochastic methods, hence they exhibit some uncertainty. However, the four-parameter model for the temperature dependence (eq 8) will describe them exactly and thus overfits the data. We examined a three-parameter quadratic polynomial as well and found that it is not adequate to characterize the data within their uncertainties according to a χ2 measure, so the four-parameter fit is most appropriate. Molecular Models. Aimoli et al.19 used the transferable potentials for phase equilibria (TraPPE) models30,31 of CO2 and CH4 in their molecular simulation studies, and we employ them as well here. TraPPE CO2 is represented as three Lennard-Jones sites with partial charges on each atom center, with a rigid C−O bond of length 1.160 Å. The TraPPE potential for CH4 treats it as a single (united atom) uncharged Lennard-Jones site. (There also exists a five-site (explicit hydrogen) potential32 for CH4, but this was not used by Aimoli et al.19 because it is significantly more expensive to simulate.) All model parameters are recorded in Table 1. Standard Lorentz−Berthelot combining rules provide the cross-interaction parameters:

Tαp 2 ρκT

(7)

(4c)

Both formulas take T in Kelvin. For mixtures, the ideal gas contribution is a mole-fraction-weighted sum of these two terms. Molar Isobaric Heat Capacity, cp. The isobaric heat capacity could be given by the manipulation of the VEOS into a volumeexplicit form (V as a series in p), but we find better results by instead applying the standard thermodynamic relation:

cp = cV +

cV κTρm

where ρm is the mass density. The speed of sound is of some interest because it is accessible experimentally and in some ways is more convenient to measure than the gas density. It can be expanded as a series in density, with coefficients known as the acoustic virial coefficients. These coefficients are related27 to the pressure coefficients Bn(T) and its derivatives Bn′(T) and Bn″(T). Speed-of-sound measurements provide one of the standard approaches for the evaluation of the virial coefficients.28 Temperature Dependence. A number of the properties listed above require temperature derivatives of the virial coefficients. It is possible to develop and compute cluster integrals for these derivatives, which can be evaluated using methods similar to those for the Bij, but for the present study we evaluate them instead by differentiating a temperature fit of the coefficients. The temperature dependence in general is nontrivial, particularly in the vicinity of the critical temperature, and care must be taken when extending into the subcritical region.29 The square-well model, for which an analytical temperature dependence can be derived, suggests that the variation of the virial coefficients with temperature should be described using a polynomial in exp(a0/T), where a0 is a fitting constant. The present study is concerned with supercritical conditions, where the temperature dependence of Bij is relatively smooth. Accordingly, we employ a simple polynomial in 1/T as the basis for describing Bij(T):

+ T 2Bn′′)ρn − 1

= cVid(T ) − (2TB2′ + T 2B2′′)ρ −

cp

(6)

where H is the enthalpy. The Joule−Thomson inversion curve the locus of points in the p, T plane where μJT = 0is particularly important because it delimits the regions where an

ϵij = C

ϵiϵj , σij =

1 (σi + σj) 2

(9) DOI: 10.1021/acs.jced.6b00702 J. Chem. Eng. Data XXXX, XXX, XXX−XXX


Article

notable exception of the Joule−Thomson coefficient, for which deviations averaged about 35%. Given that this result is now established, we will not attempt to compare the VEOS to experimental data (or a surrogate thereof). The second element of the accuracy of the VEOS is in how well it describes the true thermodynamic properties for the given molecular model. This can be assessed through comparison to molecular simulation data, and we do that here. Specifically, for a property X we compute the relative difference

Table 1. Potential Parameters for Nonbonded Interactions of the TraPPE CO2/CH4 System site

ϵ/kB (K)

σ (Å)

q (e)

C O CH4

27.0 79.0 148

2.80 3.05 3.73

+0.70 −0.35 0

Computational Details. We compute all mixture virial coefficients at each of four temperatures 323.15 K, 373.15 K, 473.15 and 573.15 K, labeled T1, T2, T3, and T4, respectively. These are the same conditions studied by Aimoli et al.19 by molecular simulation, so the properties computed from our coefficient data can be compared directly to their results. Also, in consideration of our interest in the future to use the VEOS to estimate Tc of CO2, we compute pure-CO2 virial coefficients from 200 to 600 K in 10 K intervals. All coefficients are evaluated using the overlap-sampling implementation9,10 of Mayer-sampling Monte Carlo11 with a hard-sphere reference of diameter d = 4.5 Å. Wheatley’s recursive algorithm12 is adopted for calculation of the cluster sums for each configuration. For each coefficient of order n < 6, 10 independent calculations of 109 attempted Monte Carlo (MC) trials are performed, with 10 more calculations being applied for mixture virial coefficients at 323.15 K. For sixth- and seventh-order virial coefficients for pure CO2, at least 100 independent simulations of 109 MC trials are performed for better uncertainty estimates, with more calculations being applied for 323.15 and 373.15 K. For example, 318 simulations are performed for B60 at 323.15 K and 118 simulations are performed for B60 at 473.15 K. For sixth- and seventh-order mixture virial coefficients, at least 10 simulations of 109 MC trials are performed, with more sampling applied for all mixture coefficients at 323.15 and 373.15 K, and for B33, B42, B51, B43, B52, and B61 at 473.15 and 573.15 K. For example, 130 simulations are performed for B42 at 323.15 K, and 230 simulations are performed for B61 at 373.15 K. Uncertainties are obtained as 1 standard deviation of the mean, based on the variance of the independent 109-trial runs. The pure-methane coefficients B0j are just those for the Lennard-Jones model, scaled according to ⎛ ⎞ T ⎟⎟ B0j (T ) = (10−24NAvoσCH4 3)(j − 1)B*j ⎜⎜ ⎝ (ϵ/kB)CH4 ⎠

Δsim X ,n ≡ 1 −

X(VEOSn) X(sim)

(11)

where X(VEOSn) is the value of X as given by the nth-order VEOS and X(sim) is the value reported from the simulations of Aimoli et al.19 We consider VEOSn for property X to be converged with respect to simulation if |Δsim X,n | < τ for some specified tolerance τ. We require this condition not only for n but for all j ≥ n (up to the order j for which coefficient data are available); the latter requirement ensures that the agreement with simulation at order n is not fortuitous. A unique and particularly valuable feature of the VEOS in comparison to other thermodynamic models is its ability in principle to provide a self-assessment of its accuracy so that one does not need independent simulation data to determine whether it is accurate for the given model. This is accomplished by examining whether it is converged at the thermodynamic state of interest, i.e., whether the addition of more terms to the series changes the property value significantly. We also examine the accuracy of the VEOS in this respect, in a manner similar to the approach used for the comparison to simulation. We define the relative deviation from the property given by the seventhorder series (i.e., the highest for which we have coefficient data): ΔVEOS7 ≡1− X ,n

X(VEOSn) X(VEOS7)

(12)

For a given T, ρ, and y, we conclude that the series is converged | < τ for all j, n ≤ j < 7. at order n for tolerance τ if |ΔVEOS7 X,j Uncertainty. Both the comparison to simulation data and the self-assessment of convergence are clouded by stochastic uncertainties in the simulation data and the virial coefficients. Indeed, in selecting a tolerance τ we want a value that is not so loose as to make comparison meaningless but also not so tight as to lead to numerous inconclusive comparisons because it is exceeded by uncertainty. Accordingly, the tolerance τ should be selected for each property and temperature with a consideration of the size of the uncertainties in the simulation data. The uncertainties in ΔX,n (which we use to refer generally to and Δsim ΔVEOS7 X,n X,n ) are propagated from the uncertainties in the coefficients Bij(Tk), k = 1...4. Each coefficient enters in several different ways, and we must be careful to account for correlation in the intermediate quantities due to their common reliance on the same Bij values. X(VEOSn) and X(VEOS7) are correlated, for example, as are Bn, B′n, and B″n , which variously appear in the expressions for the properties X. To handle these dependences, we express the fitted function, eq 8, directly in terms of the four calculated coefficients (which is expressible this way because we are fitting the coefficients exactly)

(10)

where Bj*(T*) is the jth Lennard-Jones virial coefficient in its usual reduced units, obtained by interpolating the values reported by us recently.33 Also, NAvo is Avogadro’s number; the formula yields the coefficient in units of (cm3/mol)(j−1) when σCH4 is in Å. VEOS Accuracy Assessment. A major aim of this work is to assess the conditions for which the VEOS model of the CO2/ CH4 mixture is accurate. Ultimately, we are interested in accuracy with respect to experimental data, but this has two elements. First is the quality of the model in characterizing the true interactions of the CO2 and CH4 molecules. This assessment was made by Aimoli et al.19 when they compared their simulation results to available experimental data and extrapolation of the well-regarded GERG-2008 model34 for thermodynamic properties of the mixture. The comparison indicated that the molecular models are good: deviations of the simulation data from best estimates of experimental behavior average about 2% or less, depending on the property, with the

(

) (

T

)

Bijfit (T ) = 1 T −1 T −2 T −3 A−1 Bij (T1) Bij (T2) Bij (T3) Bij (T4)

(13) D



Article

Figure 1. Virial coefficients Bij(T) for the four temperatures examined in this work. The abscissa is the i(CO2) index scaled by n = i + j; thus the left side of each plot is for n CH4 molecules, and the right side is for n CO2 molecules. The ordinate is the coefficient, made dimensionless by ρmax = 0.02 mol/ cm3. Symbols for the Bij values indicate the coefficient order n; lines simply join the points and are provided as a guide to the eye. Filled symbols with white numerals are experimental data: B20 and B30 are from ref 21 via the REFPROP database;36 B02 and B03 are from ref 22 via REFPROP;36 and B11 at 323.15 K (not visible behind point for calculated value) is from ref 35. Except for B70 at 373.15 K and B52, B61, and B70 at 323.15 K, confidence limits (68%) are smaller than the symbol size.

where A is the matrix with elements Akl = 1/Tl−1 k . Derivatives fit and B are then obtained by differentiation of the T ′ (T) ″ (T) Bfit ij ij dependence given in the first row vector of this expression. Although Bfit ij (Tl) does not depend on Bij(Tk) for k ≠ l fit (because it is an exact fit), the derivatives Bfit ij ′(Tl) and Bij ″(Tl) do depend on Bij(Tk) for all k = 1...4. There is error in the derivatives due to the uncertainties in Bij(Tk), and these uncertainties are easily propagated via eq 13. There is also error in the derivatives due to inaccuracies in the form for Bfit ij (T), and these are more difficult to assess. To estimate them, we examined the more extensive temperature data that we have available for the pure components, which includes finely spaced data at temperatures T in the range of T1 ≤ T ≤ T4. A comparison of Bfit ij (T) to these data finds that differences on average exceed the (95%) uncertainty in Bfit ij (T) by a factor of about 1.5. This itself is not problematic because in the present study we use the fit only for the values of Tk, where it is exact. The derivatives, however, will incur errors due to the fit, even at the Tk. Comparison to derivatives evaluated by interpolation of the more finely spaced pure-species data suggests that errors in B′(T) and B″(T) may exceed their uncertainty by about 3 and 20 times, respectively. To attempt to account for this, in our analysis below we multiply our estimates of the uncertainties of the derivatives by these respective factors. The quantity ΔX,n depends on each coefficient Bij(Tk), i + j ≤ n, through its definition (eq 11 or 12), the expression for

property X given in the section above, and eq 13. The chain rule propagates the uncertainty in each Bij(Tk) to ΔX,n, and these contributions are summed (in squares) over i, j, and k to obtain the total uncertainty in ΔX,n. Uncertainty in Δsim X,n also has contributions from the uncertainties in the simulation data, and we include this in its σΔ. The MD data from Aimoli et al.19 were obtained from isobaric-ensemble simulations, so the pressure is specified and their reported densities have uncertainty. Values of X(VEOSn) entering into Δsim X,n are evaluated using the reported MD density, but its σΔ does not include contributions from the uncertainty in the density. We assess the performance of the VEOS, relative either to simulation or VEOS7, by considering the relations among the relative difference ΔX,n, its uncertainty σΔ, and the tolerance τ. There are three possible cases: 1. τ is greater than both ΔX,n and σΔ. We conclude that VEOSn for property X agrees with its reference (VEOS7 or simulation) 2. ΔX,n is the largest of the three. We conclude that VEOSn does not agree with its reference. 3. σΔ is largest, and the comparison is inconclusive. At any given state point, for the series for property X to be considered converged with respect to a reference and tolerance at order n, which is a stricter statement than asserting that it simply agrees, it must agree with the reference for all higher orders as well. Furthermore, if the comparison is inconclusive for E



Article

Figure 2. Compressibility factor, Z, as a function of density for several model CO2/CH4 mixtures. Columns correspond to different temperatures, as indicated, and rows correspond to different compositions, where the inset box shows the mole fraction of CO2. Lines are virial series VEOSn, where n is indicated by color according to the legend. (The series order may also be identified by density where the line is truncated, with higher orders extending to higher density.) For clarity, each line is truncated where it differs from VEOS7 by 10%. Points represent the simulation data of Aimoli et al.,19 showing confidence limits when larger than the symbol size.

values for the pure components and the B11 values at 323.15 K. The computed values for CH4 and the B11 value are in good agreement with experiment, and those for CO2 differ from experiment by non-neglible amounts. The computed values for B30 are consistent with those evaluated16 for other two-body potentials, which were found to require the addition of explicit three-body interactions to provide good agreement with experiment. The inaccuracies in the second and third virial coefficients for the TraPPE CO2 model are such that they will tend to compensate when summed. At first glance, the coefficients show the behavior one would like to see for a convergent series, in that their contribution generally declines with coefficient order n: roughly the magnitudes rank as |B2| ≈ |B3| > | B4| > |Bn>4|, where the n > 4 coefficients are a bit of a jumble, except for T = 573.15 K, where the progression continues. A careful examination does find some exceptions, however. The i = n (all-CO2) coefficients do not show such smooth convergence, and there is significant uncertainty in the B61 and B70 values at the lowest temperature. The n = 6 coefficients for 323.15 and 373.15 K cross from positive to negative values with increasing i, as do the n = 7 coefficients for 323.15 K. This means that for intermediate compositions these contributions will tend to cancel when forming the mixture coefficient given by eq 1b, causing this coefficient to be small. Furthermore, we note that at low T the contributions from B2 and B3 are of opposite sign, hence they will tend to cancel as well; we see this also for the CO2-rich B2 versus B3 coefficients at all temperatures. It is worth remembering at

any n, we deem the convergence assessment inconclusive overall at that state point.

■

RESULTS AND DISCUSSION Virial Coefficients. Virial coefficients for the CO2/CH4 mixture at 323.15, 373.15, 473.15, and 573.15 K are plotted in Figure 1. Values of the plotted coefficients are tabulated in machine-readable form in the Supporting Information, where additional results for the virial coefficients of pure CO2 from 200 to 600 K are also presented. All uncertainties given in the figures and tables represent 1 standard deviation of the mean (68% confidence level). Each panel in Figure 1 shows the coefficients at a different temperature, and in each plot the coefficients for all orders n = 2−7 are presented as a function of the mixture index i in Bi,n−i, rendered as i/n to allow all orders to be given on a common scale for the abscissa. The left (i = 0) side of each plot tends toward coefficients with more CH4 molecules, and the right (i = n) tends toward CO2. The plotted coefficients are made dimensionless by ρmax raised to the appropriate power, where ρmax = 0.02 mol/cm3 is selected because it is the largest density that we examine in the discussion to follow; otherwise, this choice is arbitrary. Note then that the plotted quantities, upon multiplication by the appropriate y-dependent coefficient in eq 1b, represent the contribution that each coefficient makes to the compressibility factor Z at ρmax. The computed coefficients are compared to experimental values, where available.21,22,35,36 This includes the n = 2 and 3 F



Article

Figure 3. Isothermal compressibility κT. Data are as described in Figure 2.

Figure 4. Volume expansivity, Tαp. Data are as described in Figure 2. The dashed line for y = 1 shows VEOS7 using more accurate B′(T) derived from a larger set of coefficient values.

G



Article

Figure 5. Isochoric heat capacity cV. Data are as described in Figure 2. The dashed line for y = 1 shows VEOS7 using more accurate B′(T) and B″(T) derived from a larger set of coefficient values.

and may indicate that VEOS for the property is converged at order n ∈(1,..., 6), that it is not yet converged at VEOS6, or that the test is inconclusive due to uncertainty. Likewise, the test for agreement with MD is made, and a symbol is placed at the simulation state point: numerals indicate agreement at order n (and higher), N indicates no agreement, and U indicates inconclusive agreement due to uncertainty in the VEOS and/or the MD data. The symbols are color coded using the same scheme for coloring the regions (but slightly darker to provide contrast). As mentioned above, the selection of τ in Figure 9 is made with a consideration of the uncertainties of the simulation data.19 The relative uncertainties of these data vary with the property as well as the temperature, so they are different for each plot. The relative uncertainties also vary considerably with composition, being 2 to 15 times larger for mixtures rich in CO2 compared to those with more CH4. We want, however, to select a single tolerance for a given plot, applicable for all compositions. We find that an analysis based on a tolerance that is slightly too large is more informative than one using a tolerance that is slightly too small. Accordingly, for each property and temperature we choose the tolerance τ to be twice the mean relative error reported for the corresponding simulation data but considering data only for systems having CO2 mole fraction xCO2 ≥ 0.7. The tolerance selected in this way is indicated by the value in the upper-left corner of each plot. In assessing the performance of the VEOS with respect to how quickly it converges, one should keep in mind that the absolute location of the convergence regions is dictated by the tolerance τ, which in some respects is arbitrary inasmuch as it depends on the simulation data. The tolerance value is meaningful in that it

this point that the curves in Figure 1 indicate contributions at the highest density examined, and the balance will certainly shift away from the higher-order coefficients when applied at lower density. Still, taken together, these complications make it difficult to draw conclusions about the behavior of the series from the coefficients alone. Accordingly, in the next section we turn to an examination of the properties computed from the coefficients for different state conditions. Properties and Convergence. Thermodynamic properties as given by the VEOS to various orders are presented in Figures 2 to 8, which show results for, respectively, Z, κT, αp, cV, cp, μJT, and u as a function of molar density. For each property, separate plots are given for four temperatures (323.15, 373.15, 473.15, and 573.15 K, arranged in columns) and four values of the CO2 mole fraction (1.0, 0.6, 0.3, and 0.0, arranged in rows). In each plot, we show the progression of the VEOS with increasing order of truncation, and we provide a comparison with the MD simulation data of Aimoli et al.19 Figures 2 to 8 are useful in conveying the ability of the VEOS to yield the variety of behaviors observed for the different properties and to provide a qualitative sense of their agreement with the MD data. They do not, however, provide a good quantitative description of the convergence of the VEOS or of the agreement with MD data. For this we have the plots presented in Figure 9, which maps the regions of convergence of the VEOS for each property and temperature in the density− composition plane. For each (ρ, y) point, convergence is assessed using the procedure described above, involving a comparison of the agreement measure ΔX,n, its uncertainty σΔ, and an assigned tolerance level τ. Regions are colored according to the result of the convergence test (i.e., comparison to VEOS7) H



Article

Figure 6. Isobaric heat capacity cp. Data are as described in Figure 2. The dashed line for y = 1 shows VEOS7 using more accurate B′(T) and B″(T) derived from a larger set of coefficient values.

Figure 7. Joule−Thomson coefficient μJT. Data are as described in Figure 2, except that the lines are truncated where they differ from VEOS7 by 50%.

I



Article

Figure 8. Speed of sound u. Data are as described in Figure 2.

(which is expected to occur with some small frequency). A more systematic discrepancy is seen in the self-assessment of convergence of κT at 323.15 K, where Δsim X,n shows that VEOS7 or higher is needed for agreement with simulation, in a region where the self-assessment indicates VEOS4 or VEOS5 suffices. There is no clear explanation of this discrepancy, but it may be connected to the sign change of n = 6, 7 coefficients Bij in the range where i ≈ j. This causes the magnitude of these Bn to be somewhat smaller than their constituent Bij, thereby unduly suggesting that the influence of these density orders is small. (See the discussion above in connection to Figure 1.) There is a sizable anomaly in the Z plots for 573.15 K, where ΔVEOS7 suggests that VEOS6 is converged whereas Δsim X,n X,n indicates that it is not. An especially tight tolerance is applied in this case (0.5%), and it is possible that the discrepancy reflects a systematic error in the simulation data that may be due to finitesize effects. The discrepancy does not remain when the tolerance is increased (not shown). Finally, we point out the indication by Δsim X,n of the considerable expansion of the convergence region upon addition of the seventh-order term, demonstrated in particular for the CO2-lean region of both the Z and κT maps at 373.15 K. Because this is the highest-order term for which we have coefficients, the selfassessment cannot detect this behavior. To regain perspective, a return to the direct property plots provided in Figures 2 and 3 demonstrates the size of the error in VEOS where it is not converged with respect to the tolerances used in Figure 9. We see the poor agreement with simulation in the low-temperature, CO2-rich cases and the generally good agreement everywhere else. κT especially looks good, except for the three plots forming the upper-left corner of the grid. The peak exhibited by κT for pure CO2 at 323.15 K (note the different

characterizes the precision one can expect from simulation and thus provides a gauge for the level of accuracy where VEOS is useful. Any change in τ will mainly have the effect of shifting the regions to the left (when decreasing τ) or right (when increasing τ) while adding or removing points or regions corresponding to N or U. Consistency between the measures for convergence (ΔVEOS7 ) X,n and accuracy (Δsim X,n ) can give us confidence in the self-assessment as an indicator of accuracy of the VEOS when simulation data are not available. In reference to Figure 9, this means that a numeral indicating agreement with simulation at order n should lie in the similarly colored region indicating the self-assessment of convergence at order n. Likewise, a disagreement indicated by Δsim X,7 , which is rendered as an N point in Figure 9, should lie in a indicates that convergence (gold-colored) region where ΔVEOS7 X,n is not yet reached or (a gray-colored) one where uncertainty makes the assessment inconclusive. Agreement at order n = 7 cannot be detected by the self-assessment, so we expect the 7 numerals to lie within the (gold) nonconvergence region. With all this in mind, in examining the figure we find an excellent correspondence between the self-assessment and the comparison to simulation, although there are some exceptions that we would like to understand. Focusing first on purely volumetric properties Z and κT, the and Δsim correspondence between ΔVEOS7 X,n X,n is almost perfect. Most of the disagreements are in the direction of underestimating the accuracy of the VEOS, meaning that ΔVEOS7 X,n indicates that a higher-order series compared to that shown by Δsim X,n is needed. Isolated anomalies, such as the N appearing in the VEOS3 convergence region for Z at 573.15 K, may be attributed to the simulation data differing from their respective true property values by an amount exceeding their uncertainty J



Article

Figure 9. Convergence maps for the VEOS. Each plot shows the regions (in the plane of density and CO2 mole fraction) where the VEOS is converged at a given order relative to VEOS7 and with respect to a specific tolerance (indicated by the value in the upper-left corner of each panel). Regions are colored to indicate the VEOS order of convergence, ranging from the ideal gas (VEOS1) to VEOS6. (The upper-right figure most clearly shows the progression of the order, from left to right.) Yellow regions indicate that convergence was not found up to VEOS6, and gray indicates that statistical K



Article

Figure 9. continued uncertainties in the virial coefficients made the test inconclusive. Symbols indicate the order of the VEOS where agreement is found with MD simulation data,19 where numerals indicate the order, N indicates that no agreement was found up to VEOS7, and U indicates no conclusion due to uncertainty in the data and/or virial coefficients. Rows of figures are for properties Z, κT, αP, cV, cP, μJT, and u. Columns correspond to temperatures of 323.15, 373.15, 423.15, and 573.15 K. The value in the upper-left corner of each figure indicates the tolerance τ used to generate it.

Table 2. Critical Properties of the Pure Speciesa ρc(mol/cm3)

Tc(K) CO2 CH4 a

pc(MPa)

expt

model

expt

model

expt

model

304.128(15) 190.564(12)

306.2 191(1)

0.010625(14) 0.010139(12)

0.01056 0.00997(19)

7.377(3) 4.592(2)

7.77 4.5(11)

Values are presented from experiment21,22,37,38 and from molecular simulations of the models used in this study.31,32

scale of this figure) is a vestige of its divergence at the critical point, which is in this vicinity. (See the discussion below.) Turning to the thermal properties, in reference to Figure 9 for cV at 573.15 K, there is a large region in CO2-rich compositions where simulation indicates that the VEOS is inaccurate even though the self-assessment suggests convergence at low order (VEOS3 to VEOS5). We ascribe this anomaly to an underestimation of the inaccuracy in B′n(T) and B″n (T). Inaccuracy in the temperature derivatives is expected to be worst at the end points of the fit (i.e., 323.15 and 573.15 K) because in a sense the derivative represents an extrapolation of the fitted region there. Indeed, the anomaly is not present in the maps for the lower temperatures. (MD data for 323.15 K do not extend to low enough density to determine whether the same anomaly is present there). Interestingly, this anomaly is not seen in the αp or cp convergence maps, though this difference may be due in part to the much larger tolerances used there. A look at the 573.15 K plot in Figure 5 confirms what is indicated by the convergence map, where the VEOS estimate of cV appears to converge for densities below about 0.01 mol/cm3 while disagreeing with simulation data. Figures 4 and 6 show better agreement for αp and cp as well as a greater uncertainty in the simulation data for these properties. For pure CO2, we have available a larger set of Bij(T) data, and we have used them to generate more reliable values of B′n(T) and Bn″(T) for the four temperatures of interest. We include in Figures 4 to 6 the predictions from VEOS7 for αp, cV, and cp, respectively, when using these more accurate temperature derivatives. The results are indicated with the dashed line in the figures. The effect is quite noticeable and generally moves the VEOS7 curves toward better alignment with the simulation data. This confirms that the temperature derivatives employed here are the weak link in the calculation of the thermal properties. Still, on the basis of the overall good performance of the VEOS in these applications, the four-point estimates of Bn′(T) and B″n (T) are nevertheless useful. Finally, we examine the composite properties, μJT and u. One should note the large tolerances used for the μJT convergence map given in Figure 9, reflecting the large relative uncertainties in the simulation data for this property. The Joule−Thomson coefficient differs from the other properties in that it goes through zero under the conditions examined here and consequently the relative-difference measures ΔμJT,n diverge. This behavior gives rise, for example, to the narrow line of nonconvergence cutting diagonally through the 473.15 K convergence map, which aligns well with a corresponding string of U and N points from the simulation comparison. The μJT = 0 inversion point is of interest in itself, and perhaps surprisingly we

see from Figure 7 that the VEOS performs well in estimating it. Again, performance deteriorates toward the upper-left corner of the grid, but overall the behavior of the Joule−Thomson coefficient is well described. A similar level of accuracy is seen in the description of the speed of sound, u, as shown in Figure 8. The Joule−Thomson coefficient has the interesting feature that its limiting value at ρ → 0 is not the same as the value for an ideal gas (zero) but instead depends on B2 (and B′2(T)). The limiting slope, in turn, requires B3, so the density dependence is poorly described by VEOS2, which immediately goes off in the wrong direction with increasing ρ. Once B3 is included, a good qualitative representation of the overall behavior is obtained, and the inclusion of further coefficients brings about refinement that leads to quantitative agreement with the MD data. Further interpretation of the results is aided by reference to the critical properties, which are recorded in Table 2. First, we note that the critical densities of CO2 and CH4 are almost the same, equal to about 0.010 mol/cm3; assuming that this roughly holds for the mixtures as well, the domain of the figures extends to about twice the critical density. Although all temperatures studied here are supercritical, CO2 at 323.15 K is barely above its critical temperature. (A list of the reduced temperatures of each pure component at the four temperatures examined here is provided in Table 3.) The critical point is a mathematical Table 3. Pure-Component Reduced Temperatures T/Tc T(K) CO2 CH4

323.15 1.055 1.69

373.15 1.219 1.95

473.15 1.545 2.48

573.15 1.872 3.00

singularity in the equation of state; accordingly, it affects the convergence of the VEOS for conditions in its vicinity. Thus, in particular, we observe that for this system convergence of the VEOS for all properties is severely compromised in the CO2-rich mixtures at the lower temperatures. Similar behavior was reported by Jäger in his application to argon.17 The effect of the critical point on the VEOS is particularly clear when convergence is examined in the pressure−temperature plane, presented as log p vs 1/T. Figure 10 shows this behavior for the two pure substances and puts it in the context of the vapor−liquid coexistence line. A striking feature of these figures is the apparent barrier to convergence that appears as a virtual extrapolation of the coexistence line. This suggests that the singular behavior at the critical point hampers convergence not only there but to much higher temperatures and pressures. This observation has motivated efforts to develop approximants39,40 that extend the range of validity of the VEOS by L



Article

Figure 10. Convergence maps of the compressibility factor (Z) for pure CO2 (left) and CH4 (right) in the plane of log p versus 1/T. Numerals indicates regions in which VEOS to indicated order is converged, meaning that it (and all higher available orders) is within τ = 1% of the value of Z given by the highest-order VEOS (seventh for CO2 and eighth for CH4). IG is ideal gas, and N indicates a region where no convergence is found. Vapor−liquid coexistence for the same molecular models is shown by the solid black line, ending in the filled circle indicating the critical point.

coupling it with a form that treats the singular behavior of the critical point explicitly, using the universal features exhibited there.41 Initial application to the square-well model40 shows that this is a promising approach to improving the effectiveness of the VEOS, but more work is needed to develop and test the idea, particularly for application to mixtures. Figure 10 shows a region of convergence for VEOS2, VEOS3, and VEOS4 apparently extending across the coexistence line and into the liquid phase. This behavior is in fact a manifestation of the spinodal: the VEOS up to the orders examined here do not describe the liquid correctly, and as a consequence, it is not seen as having greater stability than the vapor in this region. Instead, the equation of state is smoothly extended past the condensation point, until reaching a limit of local stability that requires an accurate characterization of only small clusters. We anticipate that if the low-order series were being compared to VEOS at much higher orders (instead of just VEOS7 or VEOS8 here) then this metastable region would shrink, perhaps to coincide with the saturation line. These considerations are discussed in more detail elsewhere.8

simulation disadvantages VEOS, given that VEOS is in fact providing much more information than what is being compared. Nevertheless, we will in this section perform a comparison of the computational effort required to generate property values for individual thermodynamic states as obtained by VEOS and simulation. A useful basis for comparison is the difficulty, which is a quantity we introduced recently42 to gauge the computational effort needed to evaluate stochastic averages. If tx is the CPU time required to evaluate an average x̅ with an uncertainty σx, then the difficulty Dx of the calculation is defined as Dx ≡ tx1/2σx

(14)

With this definition, Dx is (for sufficient sampling) independent of tx and σx and depends on them only through the group. As the name suggests, a large value of Dx indicates that more computation time is required to obtain a result of given precision, i.e., the calculation is more difficult. Typical values of Dx range over many orders of magnitude when evaluated for different properties and algorithms, so it is more convenient to work with its logarithm, which defines the difficulty index, + :

■

+ ≡ log10(Dx /x ̅ )

COMPUTATIONAL DIFFICULTY Apart from the issues relating to the accuracy of the VEOS, there is also the question of whether it offers any advantages with respect to computational efficiency. Molecular simulation provides the appropriate standard, although the fundamental differences between the VEOS and simulation complicate any attempt to compare them. The expense of the VEOS is in the calculation of the coefficients (the effort involved in computing properties via eqs 1−7 is negligible once the Bij are known), and the result provides a representation of a broad set of properties over a range of densities and mixture compositions. In contrast, all of the expense applied in molecular simulation (at least in its basic form) provides data at a single state point, and only for the properties that are decided to be measured at the outset (although if trajectories are stored they can be used to obtain other properties after the fact). A comparison of the effort to compute one property at one state point for VEOS versus

(15)

Note the division by the average x,̅ which gives a scale to the difficulty and removes any dependence on the choice of units used to describe x (and thus σx). Adopting a convention in which tx is given in seconds, values of + for typical properties and molecular models calculated by molecular simulation using present-day computers range from −3 to +3. If x is computed as a function of independent stochastic averages, which is exactly the case for the calculation of properties from the Bij, then Dx depends on how effort is allocated to the independent averages. An optimal allocation can be derived, and if this is used, then one can show42 (using variables relevant to the present application) Dxopt =

∑ ijk

M

∂x D Bij(Tk) ∂Bij (Tk)

(16) DOI: 10.1021/acs.jced.6b00702 J. Chem. Eng. Data XXXX, XXX, XXX−XXX


Article

Figure 11. Difficulty index, as defined in eq 15, describing the computational effort required to compute a property for the CO2/CH4 mixture via molecular simulation19 (points) and by the VEOS (shaded regions). N



Article

compute the difficulty, these times were scaled up by 100 to yield the (single-core) CPU time required for the 1 × 109 trials used in our production runs and then multiplied again by the number of independent runs that were used to compute the average and uncertainty for each coefficient. The time obtained this way was combined with the coefficient uncertainty according to eq 14 to yield the difficulty for the coefficient. Values so obtained for all coefficients are tabulated in the Supporting Information. It is problematic to quantify most of the issues reviewed above that contribute extraneous differences between the computational effort of VEOS and simulation, hence we do not attempt to compensate for them explicitly; rather we keep them in mind when interpreting the comparison. However, one aspectthe different number of CPU cores used for the calculationsis readily quantified: the calculation of the virial coefficients is perfectly parallelizable, so we can expect that if we had used the 16 cores available in the dual-processor cluster employed for the MD simulations (say by distributing among them the task of computing many independent averages for each coefficient), then the total time required for the calculation would be diminished by a factor of 16. The MD simulations used all 16 cores for each state point, so to render a fairer comparison we divide the CPU times for each Bij calculation by 16 or, equivalently, divide the difficulty by 4. This is done only for the comparison to follow below; the tabulated values in the Supporting Information are based on a single CPU core. All things considered, we contend that the present comparison of computational effort inherent in VEOS versus simulation is skewed in favor of simulation. With that in mind, in Figure 11 we present comparisons of the difficulty index for the calculation of each property with each approach. White regions and points show calculations requiring low computational effort for a given precision (with + = −3 being the lowest value plotted), transitioning through blue, red, and black, which indicates the greatest computational effort (+ =+ 3). In viewing the results, one should keep in mind that a unit change in the difficulty index corresponds to a 100-fold increase in the required computational effort, noting in particular that 102+ is the CPU time in hours needed to achieve a precision of 1.7%. Hence, for this level of precision the range of CPU time covered in the figure extends from a few milliseconds to 100 years. An analysis of the difficulty maps given in Figure 11 should be made in the context of the convergence maps of Figure 9. Thus, whereas molecular simulation shows its greatest advantage (blue points on a red background) under conditions close to the CO2 critical point, VEOS7 is not converged under these conditions, so the efficiency comparison is moot. Considering then only conditions where VEOS is converged, simulation seems to have an advantage for the heat capacities at low temperature for the CH4-rich phase. However, for most of the remainder of the properties and conditions, the VEOS (where converged) shows a considerable computational advantage: we observe a large range of conditions having dark-blue or red points on a light-blue or white background. The advantage diminishes but still persists under the high-density conditions for 573.15 K, where the VEOS is still converged. The Joule−Thomson coefficient is particularly interesting, where we see a lot of red-on-blue points, indicating up to a 10 000-fold reduction in computational effort for VEOS versus simulation. The prominent red-black gash running through the map at the two highest temperatures is again a manifestation of the inversion curve. The difficulty index is defined relative to the value of the property, hence it diverges where the property approaches zero. We can see the feature in

where DBij(Tk) is the difficulty for calculating coefficient Bij at temperature Tk; the sum over temperature (Tk) is needed only for properties depending on Bij′ (T) and Bij″(T), as they depend on the fit of Bij for all four temperatures. An appealing feature of this result in the present context is that we can compute VEOSbased difficulties without specifying the series order: we compute Dopt x for the highest available order, and if coefficients Bn greater than a certain n are not needed for the condition of interest, then the optimization allocates little effort to them, with the result that Dopt x is effectively unchanged whether they are included or not. The flip side of this is that the Dopt x evaluated for a particular state is meaningful only if VEOS is converged there. The evaluation of the Dx for the simulation data is comparatively straightforward. Aimoli et al.19 provide detailed data for the uncertainties of all reported property values, and they also present averages characterizing the CPU time required for each simulation; we simply combine them according to eq 14. Although timing data are not given in ref 19 for each simulated state point, a variety of averages are reported. Simulation times are most strongly affected by the mixture composition, so we use this as the key variable, and evaluate the Dx for each property and state based on the average CPU time for all simulations at the same mole fraction. Before proceeding, let us acknowledge a few more issues in addition to that mentioned at the beginning of this section. First, the data reported by Amioli et al.19 do not represent simple simulation averages but instead are the result of an analysis using the multistate Bennett’s acceptance ratio (MBAR) method. Consequently, each uncertainty benefits from input from multiple simulations at adjoining states and is necessarily smaller than what would be given by a single simulation. We do not attempt to compensate for this in our analysis. Second, the optimization leading to eq 16 would allocate different amounts of effort to the Bij(Tk) calculations for different states. In practice, the coefficients are computed without regard to which state they are being applied and hence are evaluated suboptimally. We ignore this in the analysis. Third, the MD simulations were all conducted using 500 molecules. This is a typical but arbitrary choice, and any timing measure would be different if a different system size were used. On the other hand, for a well-formulated computer code (such as LAMMPS43), the difficulty is practically independent of system size.42 The selection of system size is made to balance finite-size effects against computational cost, and to be completely thorough, a simulation study would examine two or more system sizes and extrapolate to the thermodynamic limit (N → ∞). In contrast, the VEOS directly provides results for the thermodynamic limit.44 Finally, the MD simulations were performed19 on a dual 8-core 2.3 GHz cluster using LAMMPS,43 which is a highly optimized parallel code for molecular dynamics. The virial coefficient calculations were performed in Etomica,45 which is an object-oriented code written in Java that prioritizes generality and extensibility over computational speed. Direct timing comparisons have found that MD simulations in Etomica run about 6 times slower than LAMMPS for equivalent hardware. The full set of coefficients reported here were evaluated using several different types of processors available to us. To provide consistency in the analysis of the computational effort, timings for all Bij calculations were determined from short runs (1 × 107 MC trials) performed on a 3.2 GHz Intel Core i7-4790S CPU; this is a quad-core processor, but for the purpose of the timings we used only one core for each coefficient calculation. To O



Article

both the simulation and the VEOS difficulties (both being black along this line), which is another indicator of their agreement.

computational efficiency of the VEOS as a means of computing properties. It is problematic to compare the computational effort required for property calculation via simulation versus the VEOS, but even in a comparison that is skewed to favor simulation, we observe a considerable reduction in computational effort by using the VEOS. Given a very reliable self-assessment of convergence and with the range of application expanded even further than it is now through the formulation of approximants, the VEOS presents a very appealing route to the evaluation of gas- and supercritical-fluid-phase properties of pure substances and mixtures from molecular models and, ultimately, first principles methods.

■

CONCLUDING REMARKS We have examined the performance of the virial equation of state for the calculation of seven thermodynamic properties for the CO2−CH4 mixture at four temperatures in the supercritical region. We have examined the performance both from the standpoint of accuracy in describing the properties and with respect to the computational effort required to obtain the results. All properties appear to be converged for VEOS4 or VEOS5 for densities up to the vicinity of the critical density. Also, at the highest temperature examined, corresponding to reduced temperatures of 1.9 for pure CO2 and 3.0 for pure CH4, both the self-assessment and the comparison to simulation show that all properties are converged at VEOS7 for all densities and compositions examined by the simulations, notwithstanding anomalous behavior in cV that is likely due to the poor estimation of derivatives Bij′ (T) and Bij″(T) required to compute this property from the VEOS. These derivatives were computed here in a relatively crude manner by fitting data across four temperatures. It is possible to evaluate them in a more direct and accurate manner using formulas similar to those used to compute the coefficients themselves. It would be worthwhile to give attention to methods that perform such calculations. In spite of these possible inaccuracies in the temperature derivatives, it was found that good results could be obtained for the Joule−Thomson coefficient and the speed of sound. The results suggest in particular that the VEOS7 series can be effective in mapping the Joule−Thomson inversion curve. A more focused study of these prospects is warranted, examining the performance of the high-order series for a broader range of systems. The main disadvantage of the VEOS relative to molecular simulation (apart from the former’s inability to provide transport properties) is the uncertainty in its convergence and whether the results it provides are inaccurate in this respect. In most cases, convergence can be ascertained by examining the size of the contributions made by increasing orders of the series, and in the present application, this self-assessment is largely found to be reliable. Problems may arise when there is cancellation in the Bij coefficients when they are added to yield Bn, so one should be alert to sign changes in Bi,n−i with respect to i. A remedy we have used in the past involves summing contributions in absolute value, not to obtain the property but to gauge the convergence more reliably. We have not applied such an approach to mixtures previously nor have we attempted it in this work. However, if one has particular doubt about the convergence of the VEOS for a mixture, one might consider such a test to gain more confidence in the result. The establishment of a reliable selfassessment gives a significant advantage to the VEOS over other thermodynamic models, so it is worth looking at this question carefully. Efforts to improve the convergence of the VEOS should focus on managing the critical-point singularity, attempting to remove it via the formulation of an appropriate approximant. It is understandable that this feature of the fluid-phase behavior would hamper convergence at states near the singularity, but it seems that its effect extends more broadly, suggesting that an analytical treatment of it could improve the convergence of the VEOS even under conditions quite removed from the critical point. Advances here should also yield improvement in the

■

ASSOCIATED CONTENT

* Supporting Information S

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jced.6b00702. Pure and mixture virial coefficients of the CO2/CH4 mixture using the Mayer sampling Monte Carlo method (TXT)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: kofke@buffalo.edu. Funding

This work was supported by the U.S. National Science Foundation (grant numbers CHE-1027963 and CBET1510017). Computational resources were provided by the University at Buffalo Center for Computational Research. Notes

The authors declare no competing financial interest.

■

REFERENCES

(1) Prausnitz, J. M.; Lichtenthaler, R. N.; de Azevedo, E. G. Molecular Thermodynamics of Fluid-Phase Equilibria, 3rd ed.; Prentice Hall: Upper Saddle River, NJ, 1999. (2) Donoghue, E.; Gibbs, J. H. Condensation theory for finite, closed systems. J. Chem. Phys. 1981, 74, 2975−2989. (3) Ushcats, M. Equation of State Beyond the Radius of Convergence of the Virial Expansion. Phys. Rev. Lett. 2012, 109, 040601. (4) Ushcats, M. V. Condensation of the Lennard-Jones fluid on the basis of the Gibbs single-phase approach. J. Chem. Phys. 2013, 138, 094309. (5) Ushcats, M. V. Communication: Low-temperature approximation of the virial series for the Lennard-Jones and modified Lennard-Jones models. J. Chem. Phys. 2014, 141, 101103. (6) Ushcats, M. V. Adequacy of the virial equation of state and cluster expansion. Phys. Rev. E 2013, 87, 042111. (7) Bannur, V. M. Virial expansion and condensation with a new generating function. Phys. A 2015, 419, 675−680. (8) Schultz, A. J.; Kofke, D. A. Vapor-phase metastability and condensation via the virial equation of state with extrapolated coefficients. Fluid Phase Equilib. 2016, 409, 12−18. (9) Benjamin, K. M.; Schultz, A. J.; Kofke, D. A. Gas-phase molecular clustering of TIP4P and SPC/E water models from higher-order virial coefficients. Ind. Eng. Chem. Res. 2006, 45, 5566−5573. (10) Benjamin, K. M.; Singh, J. K.; Schultz, A. J.; Kofke, D. A. Higherorder virial coefficients of water models. J. Phys. Chem. B 2007, 111, 11463−11473. (11) Singh, J. K.; Kofke, D. A. Mayer sampling: Calculation of cluster integrals using free-energy perturbation methods. Phys. Rev. Lett. 2004, 92, 220601.

P



Article

(12) Wheatley, R. J. Calculation of high-order virial coefficients with applications to hard and soft spheres. Phys. Rev. Lett. 2013, 110, 200601. (13) Kim, H. M.; Schultz, A. J.; Kofke, D. A. Second through fifth virial coefficients for model methane-ethane mixtures. Fluid Phase Equilib. 2013, 351, 69−73. (14) Yang, S.; Schultz, A. J.; Kofke, D. A.; Harvey, A. H. Interpreting Gas-Saturation Vapor-Pressure Measurements Using Virial Coefficients Derived from Molecular Models. J. Chem. Eng. Data 2014, 59, 3183− 3192. (15) Schultz, A. J.; Shaul, K. R. S.; Yang, S.; Kofke, D. A. Modeling solubility in supercritical fluids via the virial equation of state. J. Supercrit. Fluids 2010, 55, 479−484. (16) Schultz, A. J.; Kofke, D. A.; Harvey, A. H. Molecular-based virial coefficients of CO2-H2O mixtures. AIChE J. 2015, 61, 3029−3037. (17) Jäger, B. Thermodynamic Properties of Gaseous Argon from the Ab Initio Virial Equation of State. Z. Phys. Chem. 2013, 227, 303−314. (18) Jäger, B.; Hellmann, R.; Bich, E.; Vogel, E. Ab initio virial equation of state for argon using a new nonadditive three-body potential. J. Chem. Phys. 2011, 135, 084308. (19) Aimoli, C. G.; Maginn, E. J.; Abreu, C. R. Thermodynamic Properties of Supercritical Mixtures of Carbon Dioxide and Methane: A Molecular Simulation Study. J. Chem. Eng. Data 2014, 59, 3041−3054. (20) Aimoli, C. G.; Maginn, E. J.; Abreu, C. Force field comparison and thermodynamic property calculation of supercritical CO2 and CH4 using molecular dynamics simulations. Fluid Phase Equilib. 2014, 368, 80−90. (21) 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. (22) 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 100 MPa. J. Phys. Chem. Ref. Data 1991, 20, 1061−1155. (23) Castro-Marcano, F.; Olivera-Fuentes, C. G.; Colina, C. M. JouleThomson Inversion Curves and Third Virial Coefficients for Pure Fluids from Molecular-Based Models. Ind. Eng. Chem. Res. 2008, 47, 8894−8905. (24) Colina, C. M.; Olivera-Fuentes, C. Predicted Inversion Curve and Third Virial Coefficients of Carbon Dioxide at High Temperatures. Ind. Eng. Chem. Res. 2002, 41, 1064−1068. (25) Espinoza, D.; Segura, H.; Wisniak, J.; Polishuk, I. Comments on “Joule-Thomson Inversion Curves and Third Virial Coefficients for Pure Fluids from Molecular-Based Models” and “Predicted Inversion Curve and Third Virial Coefficients of Carbon Dioxide at High Temperatures. Ind. Eng. Chem. Res. 2009, 48, 6901−6903. (26) Naresh, D. J.; Singh, J. K. Virial coefficients and inversion curve of simple and associating fluids. Fluid Phase Equilib. 2009, 279, 47−55. (27) Gillis, K. A.; Moldover, M. R. Practical determination of gas densities from the speed of sound using square-well potentials. Int. J. Thermophys. 1996, 17, 1305−1324. (28) Estrada-Alexanders, A. F.; Guzmán, O.; Pérez-Vidal, B. Highprecision virial coefficients of argon and carbon dioxide from integration of speed of sound data in the pressure-temperature domain. Mol. Phys. 2012, 110, 1349−1358. (29) Schultz, A. J.; Kofke, D. A. Interpolation of virial coefficients. Mol. Phys. 2009, 107, 1431−1436. (30) 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. (31) Potoff, J. J.; Siepmann, J. I. Vapor-liquid equilibria of mixtures containing alkanes, carbon dioxide, and nitrogen. AIChE J. 2001, 47, 1676−1682. (32) Chen, B.; Siepmann, J. I. Transferable potentials for phase equilibria. 3. Explicit-hydrogen description of normal alkanes. J. Phys. Chem. B 1999, 103, 5370−5379. (33) Feng, C.; Schultz, A. J.; Chaudhary, V.; Kofke, D. A. Eighth to sixteenth virial coefficients of the Lennard-Jones model. J. Chem. Phys. 2015, 143, 044504.

(34) Kunz, O.; Wagner, W. The GERG-2008 Wide-Range Equation of State for Natural Gases and Other Mixtures: An Expansion of GERG2004. J. Chem. Eng. Data 2012, 57, 3032−3091. (35) Esper, G. J.; Bailey, D. M.; Holste, J. C.; Hall, K. R. Volumetric behavior of near-equimolar mixtures for CO2+CH4 and CO2+ N2. Fluid Phase Equilib. 1989, 49, 35−47. (36) Lemmon, E. W.; McLinden, M. O.; Huber, M. L. NIST Reference Fluid Thermodynamic and Transport Properties-REFPROP. NIST Standard Reference Database 23, Version 9.1; National Institute of Standards and Technology: Gaithersburg, MD, 2013. (37) Kleinrahm, R.; Wagner, W. Measurement and correlation of the equilibrium liquid and vapour densities and the vapour pressure along the coexistence curve of methane. J. Chem. Thermodyn. 1986, 18, 739− 760. (38) Duschek, W.; Kleinrahm, R.; Wagner, W. Measurement and correlation of the (pressure, density, temperature) relation of carbon dioxide II. Saturated-liquid and saturated-vapour densities and the vapour pressure along the entire coexistence curve. J. Chem. Thermodyn. 1990, 22, 841−864. (39) Barlow, N. S.; Schultz, A. J.; Kofke, D. A.; Weinstein, S. J. Critical isotherms from virial series using asymptotically consistent approximants. AIChE J. 2014, 60, 3336−3349. (40) Barlow, N. S.; Schultz, A. J.; Weinstein, S. J.; Kofke, D. A. Communication: Analytic continuation of the virial series through the critical point using parametric approximants. J. Chem. Phys. 2015, 143, 071103. (41) Behnejad, H.; Sengers, J.; Anisimov, M. Applied Thermodynamics of Fluids; Royal Society of Chemistry, 2010. (42) Schultz, A. J.; Kofke, D. A. Quantifying Computational Effort Required for Stochastic Averages. J. Chem. Theory Comput. 2014, 10, 5229−5234. (43) Plimpton, S. Fast parallel algorithms for short-range molecular dynamics. J. Comput. Phys. 1995, 117, 1−19. (44) Schultz, A. J.; Kofke, D. A. Sixth, seventh and eighth virial coefficients of the Lennard-Jones model. Mol. Phys. 2009, 107, 2309− 2318. (45) Schultz, A. J.; Kofke, D. A. Etomica: An object-oriented framework for molecular simulation. J. Comput. Chem. 2015, 36, 573−583.

Q


CH4 Mixtures from

Recommend Documents