Improved Estimates of the Critical Point Constants ... - ACS Publications

Sep 12, 2016 - compressibility factor (Zc) for n-alkanes with chain lengths as large as C48. These are obtained for several different force field mode...
0 downloads 0 Views 1MB Size
Download PDF
Article pubs.acs.org/jced

Improved Estimates of the Critical Point Constants for Large n‑Alkanes Using Gibbs Ensemble Monte Carlo Simulations Richard A. Messerly,* Thomas A. Knotts IV, Richard L. Rowley, and W. Vincent Wilding Department of Chemical Engineering, Brigham Young University, Provo, Utah 84602, United States S Supporting Information *

ABSTRACT: In this work, we present improved estimates of the critical temperature (Tc), critical density (ρc), critical pressure (Pc), and critical compressibility factor (Zc) for n-alkanes with chain lengths as large as C48. These are obtained for several different force field models with Gibbs ensemble Monte Carlo simulations. We implement a recently proposed data analysis method designed to reduce the uncertainty in Tc, ρc, Pc, and Zc when predicted with molecular simulation. The results show a large reduction in the uncertainties compared to the simulation literature with the greatest reduction found for ρc, Pc, and Zc. Previously, even the most computationally intensive molecular simulation studies have not been able to elucidate the n-alkane Pc trend with respect to larger carbon numbers. The results of this study are significant because the uncertainty in Pc is small enough to discern between conflicting experimental data sets and prediction models for large n-alkanes. Furthermore, the results for Tc resolve a discrepancy in the simulation literature with respect to the correct Tc trend for large n-alkanes. In addition, the Zc results are reliable enough to determine the most accurate prediction trend for Zc. Finally, finite-size effects are shown to not be significant even for the relatively small system sizes required for efficient simulation of longer chain lengths.

■

INTRODUCTION The critical temperature (Tc), critical pressure (Pc), critical density (ρc), and critical compressibility factor (Zc) are important thermophysical properties in science and engineering. For example, these critical constants are required parameters for predicting many other thermophysical properties such as saturated liquid density,1 liquid heat capacity,2 second virial coefficient,3 liquid thermal conductivity,4 liquid vapor pressure,5 and liquid viscosity.6 As reliable estimates of thermophysical properties are essential for engineering design purposes, the Design Institute for Physical Properties database (DIPPR 801) requires accurate predictions of Tc, ρc, Pc, and Zc. The n-alkane family is of particular interest due to its significance in all chemical industries. Also, the n-alkane family provides insight into other families of compounds because it is traditionally assumed that the critical constants for other homologous series converge to the n-alkane family with increasing chain length.7 However, obtaining reliable critical constant values becomes difficult for larger compounds that thermally decompose at subcritical temperatures. For these reasons, many researchers have attempted to predict or develop sophisticated experimental approaches to measure the critical constants for larger n-alkanes (i.e., ranging between C20−C60). Recently, Nikitin et al. published a set of experimental Tc and Pc values for n-alkanes as large as C60.8 Unfortunately, their methodology involves indirect measurements and requires a great deal of data analysis; therefore, the reliability of the results is unclear. In fact, the results of Nikitin et al. for C36 differ significantly from their originally published values which were © 2016 American Chemical Society

obtained using a data analysis method that was only slightly different than the later approach.9 To further complicate matters, some mathematical models for predicting Tc of large nalkanes are more consistent with Nikitin’s original Tc values than with the most recent values. For example, Tsonopoulos’ model10 for Tc agrees strongly with Nikitin’s original Tc values, whereas Teja’s Tc model11 agrees more closely with Nikitin’s most recent data. The situation is even more troublesome for Pc, as some models (e.g., Tsonopoulos, Teja, and Lemmon10−12) propose a limiting, nonzero asymptote for Pc, while others (e.g., Nannoolal and Nikitin8,9,13) predict a trend with Pc converging to zero for the infinite chain length. These different models for Pc also result in very different trends for Zc. Unfortunately, experimental Zc data are only available for nalkanes as large as C18 because ρc values have not been reported for larger compounds. Furthermore, the Zc data have large uncertainties and follow a peculiar trend with increasing chain length. Therefore, experimental data alone are not sufficient to discern between different Zc models. Molecular simulation provides a viable approach to resolve these and other discrepancies for large n-alkanes. However, the molecular simulation literature is inconclusive on whether Nikitin’s original or more recent Tc values are more reliable. Furthermore, obtaining quantitatively accurate and precise values of Pc (and consequently, Zc) for larger compounds is Received: July 1, 2016 Accepted: August 31, 2016 Published: September 12, 2016 3640

DOI: 10.1021/acs.jced.6b00574 J. Chem. Eng. Data 2016, 61, 3640−3649

Journal of Chemical & Engineering Data

■

something that has previously eluded the molecular simulation community. In previous work, we demonstrated that this is primarily due to some limitations in the traditional data analysis method used for predicting Pc.14 The traditional approach is to regress vapor pressure (Pv) data obtained from molecular simulation to the Antoine equation and calculate Pc from this regression and the predicted value of Tc.15 Due to the exponential relationship in the Antoine equation between Pv and temperature (T), this approach requires extremely accurate values of Pv and Tc. However, most force field models are not capable of predicting both Pv and Tc to this degree of accuracy. Thus, although the extrapolation of Pv to predict Pc has a strong theoretical basis, the practical use of this approach is limited due to the inadequacies of the extant force field models. For this reason, we instead implement Vetere’s method where the Rackett equation is used to predict Pc from liquid density (ρL), Tc, and ρc.16 While the traditional approach performs poorly for force field models that do not accurately predict Pv or Tc, the Vetere method is capable of overcoming these model deficiencies. Therefore, the Vetere method is ideal for the force field models utilized in this study (NERD,17 Exp-6,18 and TraPPE-UA (referred to simply as TraPPE)19) as they predict ρL more accurately than Pv. Another unique aspect of our methodology is that we use a D-optimal design to reduce the uncertainty in Tc, ρc, Pc, and Zc.14 Because Gibbs ensemble Monte Carlo (GEMC) simulation results typically have larger uncertainties than those from Grand canonical Monte Carlo (GCMC) simulations, the reduction of the critical constant uncertainties is particularly important when utilizing GEMC simulations, as in this study. The key advantage with GEMC is that finite-size scaling of the critical constant values is typically not necessary because system size dependence is generally less for GEMC than for GCMC.15 Therefore, predicting the critical constants with GEMC is an attractive approach because, as reported by Dinpajooh et al., “when resources are not available for a rigorous finite-size scaling study, GEMC simulations provide a straightforward route to determine fairly accurate critical properties using relatively small system sizes.” This is particularly important in this work because the focus is large compounds where smaller system sizes are necessary to reduce the computational cost. The purpose of the present work is fourfold. The first is to present improved estimates of Tc, ρc, Pc, and Zc using different potential models for large n-alkanes. The second is to demonstrate that, by reducing the critical constant uncertainties, particularly for Pc and Zc, we are able to elucidate the long chain length Pc and Zc trends for n-alkanes. The third goal is to resolve a conflict between the Tc values for large n-alkanes reported in the simulation literature for the NERD and Exp-6 models.17,18 The fourth and final task is to verify that finite-size effects are reasonably small for the estimated critical constant values. The outline for this document is the following. First, we provide details for the computational methods employed in this work. Next, we present results for several potential models for C6, C8, C10, C12, C16, C24, C36, and C48. This is followed by a brief consideration of possible limitations of our results. Finally, we discuss our conclusions and how they have affected recommendations for predicting the critical constant trends of large n-alkanes.

Article

COMPUTATIONAL METHODS

Simulation Protocols. Vapor−liquid orthobaric densities were obtained from GEMC simulations using the Towhee 7.0.4 package with the DX-1597-2-7 pseudorandom number generator.20 Specifically, we simulated C6, C8, C10, C12, C16, C24, C36, and C48 with the NERD17 and TraPPE19 models. Furthermore, we performed simulations with the Exp-6 model18 for C24 and C48. This was done in an attempt to quantify finite-size effects for larger compounds by comparing our results with those reported for C24 and C48 by Errington et al., who used a more precise GCMC histogram reweighting method. As it was previously demonstrated that GEMC simulations with 200 molecules for n-decane do not display significant finite-size effects,15 we limited the system size to 200 molecules for C6 to C16. For these systems, the simulation was divided into an equilibration and production period both consisting of 100 000 Monte Carlo cycles. A Monte Carlo cycle is defined as N attempted Monte Carlo moves, where N is the number of molecules used in the simulation. To assist in our analysis of finite-size effects for the longer chain lengths, we performed simulations with 200 and 1600 molecules for C24, C36, and C48. The 200 molecule simulations had an equilibration and production period of at least 200 000 Monte Carlo cycles. Due to inordinate computational costs, the equilibration and production periods for the 1600 molecule systems used fewer Monte Carlo cycles. Specifically, the number of Monte Carlo cycles for the equilibration and production periods were, respectively, 50 000 and 100 000 for C24, 30 000 and 70 000 for C36, and 20 000 and 40 000 for C48. Although we verified that the equilibration periods were sufficiently long, even for C48, the resulting uncertainties in the 1600 molecule simulations are much larger for C48 than for the shorter molecules because of the relatively short production periods. The implications of these increased uncertainties are discussed later. In this study, we utilized a novel experimental design approach known as a D-optimal design to reduce the uncertainties in Tc, ρc, Pc, and Zc.14 Only two temperatures (TR1 and TR2) are needed for the D-optimal design; however, we performed several replicate simulations at each temperature. Three replicates were performed for smaller compounds, while five replicates were required for the largest compounds to improve the statistical certainty of the results. According to the D-optimal design, it is important to perform simulations at the lowest temperature possible to minimize uncertainty in the critical constants. The lowest feasible value for TR1 was found to be near 0.7 for each chain length. The optimal value for TR2 is 0.85 when TR1 = 0.7.14 In our previous study, we also reported Tc-optimal and ρc-optimal designs that minimize uncertainties in Tc and ρc, respectively. By contrast, the D-optimal design attempts to reduce uncertainties for all of the critical constants. We did not use the Tc-optimal or ρc-optimal designs in the present work because, as demonstrated in our previous study, improvement in the uncertainty in Tc and ρc is small when compared with the D-optimal design uncertainties.14 It is important to note that we used the same Monte Carlo move probabilities as those reported in our previous publications so that the D-optimal design conditions would be valid.14,21 To be specific, the probabilities of performing a volume exchange, molecule regrowth,22 configurational bias exchange,23 trans3641

DOI: 10.1021/acs.jced.6b00574 J. Chem. Eng. Data 2016, 61, 3640−3649

Journal of Chemical & Engineering Data

Article

Table 1. Summary of Compounds, Force Field Models, Total Number of Molecules (N), State Points (T1 and T2), Number of Replicates at Each Temperature (NR), and the Number of Monte Carlo Cycles in the Equilibration (MCE) and Production (MCP) Periods compound

model

N

T1 (K)

T2 (K)

NR

MCE

MCP

n-hexane (C6H14) n-hexane (C6H14) n-octane (C8H18) n-octane (C8H18) n-decane (C10H22) n-decane (C10H22) n-dodecane (C12H26) n-dodecane (C12H26) n-hexadecane (C16H34) n-hexadecane (C16H34) n-tetracosane (C24H50) n-tetracosane (C24H50) n-tetracosane (C24H50) n-tetracosane (C24H50) n-hexatriacontane (C36H74) n-hexatriacontane (C36H74) n-hexatriacontane (C36H74) n-hexatriacontane (C36H74) n-octatetracontane (C48H98) n-octatetracontane (C48H98) n-octatetracontane (C48H98) n-octatetracontane (C48H98) n-octatetracontane (C48H98) n-octatetracontane (C48H98)

NERD TraPPE NERD TraPPE NERD TraPPE NERD TraPPE NERD TraPPE NERD TraPPE Exp-6 Exp-6 NERD NERD TraPPE TraPPE Exp-6 Exp-6 NERD NERD TraPPE TraPPE

200 200 200 200 200 200 200 200 200 200 200 200 200 1600 200 1600 200 1600 200 1600 200 1600 200 1600

360 360 390 390 435 435 465 465 510 510 580 580 560 560 650 650 650 650 700 700 730 730 730 730

440 440 490 490 535 535 575 575 625 625 700 700 680 680 765 765 765 765 800 800 830 830 830 830

3 3 3 3 3 3 3 3 3 3 3 3 5 5 5 5 5 5 5 5 5 5 5 5

100 000 100 000 100 000 100 000 100 000 100 000 100 000 100 000 100 000 100 000 100 000 100 000 200 000 50 000 200 000 30 000 200 000 30 000 200 000 20 000 200 000 20 000 200 000 20 000

100 000 100 000 100 000 100 000 100 000 100 000 100 000 100 000 100 000 100 000 100 000 100 000 200 000 100 000 400 000 70 000 400 000 70 000 300 000 40 000 300 000 40 000 300 000 40 000

regressing the GEMC orthobaric densities to the law of rectilinear diameters25 and density scaling law.26 The law of rectilinear diameters can be expressed as

lation, and rotation move were 1, 13, 20, 33, and 33%, respectively. An initial estimate for Tc is required to calculate T1 = 0.7 × Tc and T2 = 0.85 × Tc. The initial estimate for Tc was obtained from previous studies for the NERD, Exp-6, and TraPPE models;17−19 however, because reliable simulations of C24, C36, and C48 have not been performed for the TraPPE model, we estimated the TraPPE Tc to be the same as that found in the literature for the NERD model. These initial approximations for Tc do not affect the final values reported for Tc, ρc, Pc, and Zc. The TraPPE and Exp-6 nonbonded interactions were truncated at a cutoff distance of 14 Å, while the NERD potential used a value of 13.8 Å. The cutoff distance for the TraPPE and NERD models were chosen to be consistent with the values used in the development of the respective force field. However, because Errington et al. did not report a long-range cutoff distance for the Exp-6 model, we chose the larger value of 14 Å. For each model, we implemented standard long-range tail corrections.24 (For an important discussion on the use of this approach with polymeric molecules, see Note 18 in ref 17.) In addition, we employed a hard inner cutoff value of 1 Å to improve computational efficiency. We provide a summary of the simulation specifications in Table 1. In this table, we include each compound that was simulated, the force field model used, the total number of molecules between the two simulation boxes (approximately 20% in the vapor phase), the two temperatures (D-optimal conditions), the number of replicates at each temperature, and the number of Monte Carlo cycles in the equilibration and production periods. Data Analysis Methodology. Having discussed the simulation protocol, we now focus our attention on the data analysis methodology. Estimates for Tc and ρc were obtained by

ρr ≡

ρL + ρv 2

= ρc + A(Tc − T )

(1)

where ρr is the rectilinear density, ρv is the vapor density, and A is a fitting parameter. The density scaling law is ρs ≡ ρL − ρv = B(Tc − T )β

(2)

where ρs is the scaling density, B is a fitting parameter, and β is a constant, typically 0.326.15 The estimate for Pc is calculated according to Vetere’s method by rearranging the Rackett equation14,27 as (1 − Ti / Tc)−γ

R gTcρc ⎛ ρc ⎞ ⎜⎜ ⎟⎟ Pc = M w ⎝ ρL ⎠

(3)

where Rg is the universal gas constant, Mw is the molecular weight, Ti is the temperature corresponding to ρL, and γ is an empirical parameter estimated to be around 2 . The critical 7 compressibility factor is defined as

Zc ≡

PcM w R gTcρc

(4)

Finally, the uncertainties in Tc, ρc, Pc, and Zc were obtained using the nonlinear approach explained in our previous studies.14,21 It should be reiterated that a primary goal of this paper is to validate the use of eq 3 with a D-optimal design to improve the predictability of Pc and Zc from GEMC simulations. 3642

DOI: 10.1021/acs.jced.6b00574 J. Chem. Eng. Data 2016, 61, 3640−3649

Journal of Chemical & Engineering Data

Article

Table 2. Predicted Tc, ρc, Pc, and Zc Values for the Different Compounds, Force Field Models, and System Sizes Simulated in This Studya

a

compound

model

N

n-hexane (C6H14) n-hexane (C6H14) n-octane (C8H18) n-octane (C8H18) n-decane (C10H22) n-decane (C10H22) n-dodecane (C12H26) n-dodecane (C12H26) n-hexadecane (C16H34) n-hexadecane (C16H34) n-tetracosane (C24H50) n-tetracosane (C24H50) n-tetracosane (C24H50) n-tetracosane (C24H50) n-hexatriacontane (C36H74) n-hexatriacontane (C36H74) n-hexatriacontane (C36H74) n-hexatriacontane (C36H74) n-octatetracontane (C48H98) n-octatetracontane (C48H98) n-octatetracontane (C48H98) n-octatetracontane (C48H98) n-octatetracontane (C48H98) n-octatetracontane (C48H98)

NERD TraPPE NERD TraPPE NERD TraPPE NERD TraPPE NERD TraPPE NERD TraPPE Exp-6 Exp-6 NERD NERD TraPPE TraPPE Exp-6 Exp-6 NERD NERD TraPPE TraPPE

200 200 200 200 200 200 200 200 200 200 200 200 200 1600 200 1600 200 1600 200 1600 200 1600 200 1600

Tc (K) 519 505 579 571 626 622 664 662 725 726 813 817 813 815 887 890 901 902 946 948 945 947 959 961

± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ±

1 2 1 1 2 1 2 1 3 2 2 2 3 2 2 2 1 2 5 7 3 5 3 4

ρc (kg/L) 0.229 0.240 0.230 0.238 0.230 0.236 0.227 0.233 0.221 0.224 0.213 0.213 0.207 0.206 0.199 0.196 0.200 0.197 0.189 0.180 0.185 0.186 0.186 0.187

± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ±

0.002 0.005 0.001 0.001 0.003 0.002 0.003 0.002 0.004 0.003 0.003 0.002 0.003 0.002 0.002 0.003 0.002 0.002 0.004 0.006 0.005 0.006 0.004 0.005

Pc (MPa) 3.10 3.2 2.52 2.60 2.15 2.18 1.81 1.87 1.38 1.41 0.96 0.95 0.89 0.89 0.60 0.59 0.61 0.60 0.44 0.42 0.43 0.42 0.43 0.45

± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ±

0.05 0.1 0.02 0.03 0.06 0.04 0.04 0.03 0.06 0.04 0.03 0.02 0.02 0.02 0.01 0.02 0.01 0.02 0.02 0.03 0.02 0.05 0.02 0.02

Zc 0.270 0.271 0.260 0.262 0.255 0.255 0.246 0.248 0.234 0.235 0.225 0.222 0.217 0.218 0.206 0.205 0.208 0.205 0.202 0.198 0.197 0.198 0.196 0.200

± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ±

0.003 0.006 0.001 0.002 0.004 0.003 0.004 0.002 0.006 0.004 0.003 0.002 0.003 0.002 0.003 0.004 0.002 0.003 0.005 0.008 0.004 0.009 0.003 0.004

N is the total number of molecules; uncertainties are reported at the 95% confidence level.

■

RESULTS Table 2 contains the results of the simulations for each condition outlined in Table 1. The predicted critical constant values with their corresponding 95% confidence intervals are listed. A detailed discussion of these results is provided in the following sections. For reference, see Supporting Information for the vapor−liquid coexistence data of each system under consideration. In the discussion that follows, we compare our simulation results with (1) experimental data, (2) prediction models, and (3) existing simulation literature. We first focus on the smaller n-alkanes for which existing simulation results can be found in the literature for the force field models in question. We then direct our attention to the larger compounds which have fewer reliable simulation data. The majority of the discussion will focus on these larger systems because a primary goal of this study is to obtain more accurate estimates of the critical constants for large n-alkanes. Validation of D-Optimal Design. The primary purpose of this section is to demonstrate that the predicted Tc, ρc, and Pc values obtained with the D-optimal design are reliable (i.e., accurate with small uncertainties). For this purpose, we simulated n-alkanes with a carbon number (CN) range between 6 and 24 with the NERD and TraPPE models and compared these results with those found in the literature. Specifically, our results for Tc and ρc can be directly compared to those reported by Nath et al. and Martin et al. because they were obtained through the same data analysis method, namely by using eqs 1 and 2. However, we analyzed the literature data using eq 3 to obtain values of Pc to use for comparison. These estimates for Pc are referred to as Nath-Vetere and Siepmann-Vetere, respectively (see our previous work for a more complete explanation14). By using the same data analysis method for Pc

on both our data and those found in the literature, we were able to solely investigate the performance of the D-optimal design. Figures 1 and 2 provide a comparison between our Tc and ρc results and those found in the literature. In each of these

Figure 1. Trend for Tc with respect to CN for lower CN. Three experimental data sets are included. The values reported by Nath et al. and Martin et al. are provided for comparison with our simulation results using the NERD and TraPPE models with a D-optimal design. Error bars correspond to 95% confidence intervals. 3643

DOI: 10.1021/acs.jced.6b00574 J. Chem. Eng. Data 2016, 61, 3640−3649

Journal of Chemical & Engineering Data

Article

Figure 2. Trend for ρc with respect to CN for lower CN. Included are simulation results for the NERD and TraPPE models from the literature and this work using the D-optimal design. Error bars correspond to 95% confidence intervals.

Figure 3. Trend for Pc with respect to CN for lower CN. Three conflicting experimental data sets are included. The Nath-Vetere and Siepmann-Vetere results were obtained by analyzing the values from Nath et al. and Martin et al. with eq 3. Included are simulation results from this work using the NERD and TraPPE models with a D-optimal design. Error bars correspond to 95% confidence intervals.

graphs, the critical property is plotted as a function of carbon number. The uncertainties for the literature values of the TraPPE model are those reported by Martin et al. However, the uncertainties for Nath et al. are those obtained from a rigorous nonlinear analysis found in a previous study.21 The key conclusion of Figure 1 is that our Tc values are consistent with those reported in the literature. As we discussed in our previous study, the uncertainty in Tc is not significantly reduced with the D-optimal design. However, in Figure 2 we see that, in general, the D-optimal design greatly reduces the uncertainty for ρc, as predicted in our previous work.14 Furthermore, our ρc results are consistent with the literature values for most of the smaller n-alkanes. In Figure 3, we compare our Pc results for the NERD and TraPPE models with the Nath-Vetere and Siepmann-Vetere values. Again, this comparison provides insight into the performance of the D-optimal design because the literature values were obtained from the same potential models and used the same data analysis method as our results, i.e., eq 3. There are three significant conclusions that can be drawn from Figure 3. First, our estimates for Pc agree well with the experimental data. Second, our results for the NERD model agree with the Nath-Vetere values, and our results for the TraPPE model agree well with the Siepmann-Vetere values. Finally, by using a Doptimal design, the corresponding uncertainties are small enough to provide quantitative insight (which will be important for the discussion of large n-alkanes). Notice that our uncertainties are typically much smaller than those for NathVetere and Siepmann-Vetere. Results for C24, C36, and C48. Having demonstrated that the D-optimal design predicts properties well for smaller nalkanes, we focus the remainder of the discussion on larger nalkanes. First, we analyzed Tc and ρc as this provides a useful comparison with existing simulation results because our data analysis method for these two properties is the same as that found in the literature (eqs 1 and 2). To provide further insight,

we also implemented the rigorous uncertainty analysis developed in our previous work to obtain a joint-confidence region for Tc and ρc.21 Our Tc results resolve a conflict in the simulation literature as to which data and prediction models are most reliable for large n-alkanes. We will not discuss in detail our results for ρc as they are consistent with those previously reported in the simulation literature.17,18 In short, Teja’s ρc model is more accurate than that of Tsonopoulos for large nalkanes.10,11 After the description of the results for Tc and ρc, we shift our attention to Pc to demonstrate the strength of using the Rackett equation and a D-optimal design for predicting Pc both accurately and precisely. The results suggest a definitive answer as to which experimental data and prediction methods are most reliable for Pc of large n-alkanes. Finally, we present our Zc values with respect to CN to validate the internal consistency of our Tc, ρc, and Pc results. Our Zc results also help discern between the different prediction models and their trends with increasing chain length. Tc and ρc Results for Larger n-Alkanes. Figure 4 contains an uncertainty analysis from several different methods for estimating Tc and ρc. Panel a is for C24, panel b is for C36, and panel c is for C48. The ordinate is the critical temperature, while the abscissa is the critical density. The results are presented as uncertainty regions at the 95% confidence level. The confidence region for the experimental data is presented as two lines rather than boxes or ovals (denoting joint-confidence regions) because experimental data exist for only Tc for these compounds. The lines correspond to the 95% confidence intervals in Tc reported by Nikitin et al., where we included both the values reported in 1997 and 2014.8,9 However, notice that no experimental data exist for C48. We also provided uncertainty regions for two different prediction models, that of Teja and Tsonopoulos.10,11 These joint-confidence regions are 3644

DOI: 10.1021/acs.jced.6b00574 J. Chem. Eng. Data 2016, 61, 3640−3649

Journal of Chemical & Engineering Data

Article

Figure 4. Comparison between uncertainty regions of experimental data from Nikitin, two different prediction models from Tsonopoulos and Teja, literature simulation results from Errington (with Tc and ρc shifted by 1% to correct for finite-size effects) and Nath, and simulation results from this work using the Exp-6, NERD, and TraPPE models. Panels a, b, and c correspond to C24, C36, and C48, respectively. Each region represents the 95% joint-confidence region for Tc and ρc.

represented as boxes because the Tc and ρc prediction models were regressed independently and, therefore, the correlation should be zero. Although some correlation is expected, the molecular simulation results from Errington et al. are also reported as boxes because the author only reported the standard error for Tc and ρc. The oval shaped joint-confidence regions were obtained using the algorithm outlined in our previous study.21 This algorithm rigorously accounts for the correlation between regressing Tc and ρc from eqs 1 and 2 and was performed not only on the data from this work but also on the data of Nath et al.17 The results presented in Figure 4 support several key conclusions. First, the confidence regions from the D-optimal design (denoted as “This Work” in Figure 4) are very small, considering only 10 GEMC simulations were performed. Notice that the relatively simple D-optimal design yielded a smaller confidence region than that of the computationally intensive GCMC histogram reweighting approach employed by Errington et al. (Figures 4a and c). This is significant when trying to assess which prediction model and experimental data are most reliable. Elucidating Correct Tc Trend. A significant observation from Figure 4c is that the C48 results from this work for the NERD and Exp-6 models strongly agree (dashed/solid dark blue and dashed/solid dark green ovals). By contrast, the confidence regions for Nath et al. using the NERD model (dashed-dotted light blue oval) and Errington et al. using the Exp-6 model (dashed-dotted purple rectangle) differ considerably, particularly for Tc. The importance of this observation is best demonstrated by comparing the different methods of estimating Tc with respect to carbon number for larger nalkanes. In Figure 5, we replotted the experimental data, prediction models, and simulation results as found in Figure 4 but focused on the dependence of Tc on CN for longer chain lengths. As depicted in Figure 5, the Tc value for C48 reported by Nath et al. for the NERD model was much lower than that reported by Errington et al. for the Exp-6 model. This is important because

Figure 5. Trend for Tc with respect to CN for higher CN. Included are two different sets of experimental data from Nikitin, two different prediction models from Tsonopoulos and Teja, literature simulation results from Errington et al., Nath et al., and Siepmann’s group (Zhuravlev et al.), and simulation results from this work using the Exp6, NERD, and TraPPE models. Results from only the 200 molecule simulations were included for clarity. The 1600 molecule simulation results overlap considerably with those shown for 200 molecules. Error bars correspond to 95% confidence intervals. For clarity, only a single representative error bar is presented for Nikitin’s experimental data.

the results from Nath et al. agree more with Tsonopoulos’ Tc model, whereas the results from Errington et al. support Teja’s Tc model. Our results resolve the conflict between the NERD and Exp-6 models in that both conclude that Teja’s Tc model is more accurate than Tsonopoulos’ Tc model for larger n-alkanes. In addition, although the TraPPE model is known to over predict Tc for larger n-alkanes, the TraPPE model also agrees more closely with Teja’s model than with Tsonopoulos’. 3645

DOI: 10.1021/acs.jced.6b00574 J. Chem. Eng. Data 2016, 61, 3640−3649

Journal of Chemical & Engineering Data

Article

GEMC simulations is significantly smaller than those from grand-canonical ensemble simulations.”15 Pc Results for Larger n-Alkanes. Having compared our results with the literature for Tc and ρc, we now turn our attention to Pc. Figure 6 depicts the dependence of Pc upon

Our simulation data presented in Figure 5 also help discern between the conflicting experimental Tc data. Notice that our simulation results for C24, C36, and C48 suggest that the more recent experimental values from Nikitin et al. are more reliable than their original values. Interestingly, the Tc values from Nath et al. are in better agreement with the original values reported by Nikitin et al. As explained in the next section, we attribute the discrepancy between our results and those of Nath et al. to finite-size effects. Quantifying Finite-Size Effects. Finite-size effects can become significant when simulating large compounds because the number of molecules typically simulated are relatively small. Unfortunately, quantifying finite-size effects with GEMC is not a straightforward task.15 We implemented two different methods to quantify the finite-size effects found in the results presented in Figure 4. First, we employed a common approach found in the literature of simply comparing the results from simulations using 200 and 1600 molecules.15 Second, we compared our results for the Exp-6 model to those of Errington et al. Both approaches are now explained with the former being presented first. Recently, it was shown that, for GEMC systems as large as ndecane, finite-size effects for Tc and ρc were smaller than the corresponding statistical uncertainty.15 Following this methodology, we arrive at the same conclusion for the large n-alkanes simulated in this study. Notice in Figure 4c that even for C48 the 95% confidence regions overlap between the 200 and 1600 molecule systems sizes for a given force field model. On the other hand, our results for the NERD model differ considerably from the results of Nath et al. for C36 and C48 (see Figures 4b and c). In the prior work, Nath et al. state that due to finite-size effects, their “results should therefore be regarded with caution” as they were obtained with a system size of only 120 molecules. The results in Figures 4b and c seem to support this statement. In short, the fact that the results for 120 and 200 molecules are significantly different while those for 200 and 1600 molecules are comparable suggests that the 200 molecule simulations have converged to the infinite system value. This conclusion is consistent with what was previously observed for GEMC simulations of n-decane.15 While the previous comparison of 200 and 1600 molecule systems is promising, another approach to quantify finite-size effects is to compare our results for the Exp-6 model with those originally reported by Errington et al. This comparison was performed for C24 and C48, which are the two largest molecules previously simulated with the Exp-6 model. Errington et al. demonstrated that finite-size effects resulted in a systematic under prediction of Tc and ρc of approximately 1% for ethane and n-octane with the system sizes used in developing the Exp6 model.18 Because Errington et al. implemented the same system size for the larger n-alkanes, it is assumed that C24 and C48 also suffer from this systematic deviation. To account for finite-size effects and to allow for a useful comparison with our results, we increased the values that Errington et al. reported for Tc and ρc by 1% (Figure 4). As Figures 4a and c depict, the results from this study for C24 and C48 overlap considerably at the 95% confidence level with the shifted uncertainty regions for the Errington et al. values. It is therefore reasonable to assume that finite-size effects do not affect our results to a large degree. In brief, the results for long chain length molecules agree with the conclusion reported by Dinpajooh et al. that “the finite-size dependence of the critical properties obtained from

Figure 6. Trend for Pc with respect to CN for higher CN. Two conflicting experimental data sets are included as well as three prediction models regressed to experimental data. Simulation results for the TraPPE and Exp-6 models are as reported by Zhuravlev et al. and Errington et al., respectively. The Nath-Vetere and SiepmannVetere results are also shown. Included are simulation results from this work using the Exp-6, NERD, and TraPPE models with a D-optimal experimental design. Results from only the 200 molecule simulations were included for clarity. The 1600 molecule simulation results overlap considerably with those shown for 200 molecules. Error bars correspond to 95% confidence intervals. For clarity, only a single representative error bar is presented for Nikitin’s experimental data.

carbon number for large n-alkanes. Included in Figure 6 are both sets of experimental values reported by Nikitin et al.,8,9 three different prediction models,8,12,13 molecular simulation results from three different sources,18,19,28 and the results we obtained for the TraPPE, NERD, and Exp-6 models. Remember that our results utilize Vetere’s method to predict Pc from ρL, while those from the literature used Pv to predict Pc. However, we also implemented the Vetere method with the orthobaric density results reported by Nath et al. for the NERD model and those reported by Zhuravlev et al. for the TraPPE model. These Pc values are referred to as Nath-Vetere and Siepmann-Vetere in Figure 6. The uncertainties for NathVetere, Siepmann-Vetere, and the results from this work were obtained with the rigorous nonlinear analysis,14 while the uncertainties for the simulation literature values are those reported by Errington et al.18 and Zhuravlev et al.28 Also, we 3646

DOI: 10.1021/acs.jced.6b00574 J. Chem. Eng. Data 2016, 61, 3640−3649

Journal of Chemical & Engineering Data

Article

included a single error bar for each set of experimental data that is representative of the 95% confidence intervals that were reported by Nikitin et al.8,9 There are at least two important conclusions that can be drawn from the results presented in Figure 6. First, our simulation results for all three potential models (Exp-6, NERD, and TraPPE) agree with Nikitin’s most recently reported values for Pc. Again, although Nath-Vetere agrees most with Nikitin’s older set of data and with Nannoolal’s model, we attribute this to finite-size effects. The deviation in Errington’s results is likely due to the traditional approach of obtaining Pc from Pv. Even though the Exp-6 model predicts Pv and Tc fairly accurately, the over prediction of Pc is probably “the result of compound errors in the predictions of critical temperatures and vapor pressure.”29 By contrast, we were able to obtain consistent predictions of Pc for the simple Lennard-Jones 12-6 models (NERD and TraPPE) and the more advanced Exp-6 model by utilizing the Vetere method. Second, by using a D-optimal design with the Vetere method, the uncertainties in our results are small enough to elucidate the correct long chain length trend. Notice that the uncertainties in Pc from simulation literature data are typically very large even when the Vetere method is employed. We believe that these are the first quantitatively accurate and precise estimates of Pc from molecular simulation for large nalkanes. These results are significant because they demonstrate that Nikitin’s most recent Pc data should be accepted over his previous data. Our results also suggest that the models developed by Lemmon et al. and Nannoolal et al. are not as reliable as the model reported by Nikitin et al. for predicting Pc of long chain length n-alkanes. Zc Results. The critical compressibility factor is a key constant when assessing the reliability of Tc, ρc, and Pc values collectively. Figure 7 provides the Zc trend with respect to CN for the entire range of compounds considered in this work. Included in Figure 7 are the experimental values reported by Teja et al.11 (because Nikitin et al. did not provide any ρc values, Zc cannot be calculated for their data), four different prediction models,8,10,11,13,30 molecular simulation results from three different sources,18,19,28 the Siepmann-Vetere and NathVetere results, and the results we obtained for the TraPPE, NERD, and Exp-6 models. The model labeled “Nikitin” is a composite of two different studies. Specifically, the Tc and Pc models are from the 2014 publication,8 while the ρc model is from a 1998 publication,30 which is the most recent ρc model reported by this group. We included a single representative experimental uncertainty for C11 that was obtained from the DIPPR 801 database.31 The uncertainties corresponding to the literature simulation values for the Exp-6 and TraPPE models were those reported by Errington et al. and Zhuravlev et al., respectively.18,28 Notice that for the Siepmann values, we provided only a single error bar for clarity because the uncertainties reported by Martin et al. for C8 and C12 are of the same size as those reported by Zhuravlev et al. It should be mentioned that the propagation of error derivation used to calculate the literature uncertainties assumes that the uncertainties in Tc, ρc, and Pc are not correlated, which is clearly not a valid assumption. However, this is the standard approach found in the literature. The uncertainties for Siepmann-Vetere, Nath-Vetere, and our results were obtained using a statistically rigorous algorithm that accounts for this correlation.14

Figure 7. Trend for Zc with respect to CN. Included are experimental data, four prediction models, and simulation results from the literature and from this work using the Exp-6, NERD, and TraPPE models with a D-optimal experimental design. Results from only the 200 molecule simulations were included for clarity. The 1600 molecule simulation results overlap considerably with those shown for 200 molecules. Error bars correspond to 95% confidence intervals. For clarity, only a single representative error bar is presented for the Siepmann values and experimental data.

There are several key observations for the results presented in Figure 7. First, our results appear to be the most internally consistent set. Notice that the values from Errington et al., Siepmann, and Nath-Vetere have large oscillations between different CN. The Siepmann-Vetere results agree surprisingly well with the Nikitin model, however, the uncertainties are too large to provide much confidence. Because Zc depends on all three of the critical constants, the uncertainties in Zc are typically quite large. By contrast, the uncertainties for our results are small enough to discern between the different prediction models. It is also worth noting that our simulation results appear to be approaching an asymptotic value around 0.2. Recently, the statistical associating fluid theory (SAFT) was 1 used to predict that Zc → 5 for the infinite chain length nalkane.32 Notice that the Nath-Vetere results are much lower than this value, and those of Errington et al. appear to be approaching a much higher asymptote. Furthermore, the models from Teja et al. and Tsonopolous et al. diverge because Zc → ∞, while the Nannoolal model predicts that Zc → 0.10,11,13 It is worth mentioning that the asymptotic Zc value for the Nikitin composite model is 0.15. Although the Nikitin model disagrees with the limiting value obtained from the SAFT equation of state, it does agree well with our simulation results. Therefore, we believe that the Nikitin composite model 3647

DOI: 10.1021/acs.jced.6b00574 J. Chem. Eng. Data 2016, 61, 3640−3649

Journal of Chemical & Engineering Data

Article

temperature values (TR ≈ 0.7 and 0.85). The other temperatures shown in Figure 8 are simply to demonstrate the validity of eqs 1 and 2 over a large temperature range. The same simulation specifications and number of replicates were used at each temperature.

is the most reliable approach for predicting Zc for larger nalkanes, at least as large as C48.

■

LIMITATIONS As mentioned in our previous study,14 some limitations exist for using the Vetere approach with the D-optimal design. One such limitation is whether or not the law of rectilinear diameters and the density scaling law are valid at temperatures far removed from Tc. For this reason, in Figure 8, we present the linearized forms of eqs 1 and 2 with respect to TR for C24, C36, and C48. For clarity, the results for only the 200 molecule systems are shown in Figure 8. The main point of this figure is that the low TR values are consistent with a straight-line regression. It is worth clarifying that the results presented previously in Figures 4−7 were obtained from only two

■

CONCLUSIONS In this work we have demonstrated that the uncertainties of critical constants obtained from GEMC simulation can be reduced by utilizing a D-optimal design. Specifically, the uncertainty in ρc is greatly reduced by performing simulations at reduced temperature values near 0.7 and 0.85. Furthermore, by reducing the uncertainty in ρc, the Rackett equation yields reliable estimates for Pc. A reduction in the uncertainties for ρc and Pc leads to significantly smaller uncertainties in Zc than typically found in the literature. Finally, we verified that finitesize effects are of approximately the same magnitude as the statistical uncertainty. There are three main implications from our simulation results. First, our results rectified the discrepancy in the simulation literature between the NERD and Exp-6 models for the Tc trend of large n-alkanes. Both models support Teja’s Tc model over Tsonopoulos’ model. Second, we obtained accurate and precise results of Pc using the TraPPE, NERD, and Exp-6 models for large n-alkanes. These results helped us evaluate which experimental data were most reliable for the critical constants. Specifically, it was found that the most recent Tc and Pc values are indeed more reliable than the original values reported by Nikitin et al. Third, our results support the limiting Zc value of 0.2 predicted by the SAFT model. This suggests that an acceptable set of Tc, ρc, and Pc prediction models should result in a noninfinite and nonzero Zc value for the infinite chain length n-alkane. In total, this study demonstrates how reducing the uncertainty in molecular simulation results can assist in analyzing the validity of experimental data and prediction methods.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jced.6b00574. Intramolecular and intermolecular potential parameters for the TraPPE, NERD, and Exp-6 force fields and vapor−liquid coexistence simulation data for each compound, force field, system size, and temperature investigated in this study (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Funding

The authors are grateful to the Design Institute for Physical Properties (DIPPR 801) for funding. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors acknowledge the Brigham Young University Fulton Supercomputing Lab for providing computational resources.

Figure 8. Validation of the regression equations over a large range in TR for C24, C36, and C48 using 200 molecules. In panels a and b we see follow a linear trend over this entire TR that, respectively, ρr and ρ1/β s range. 3648

DOI: 10.1021/acs.jced.6b00574 J. Chem. Eng. Data 2016, 61, 3640−3649

Journal of Chemical & Engineering Data

■

Article

segments of chain molecules with strong intramolecular interactions. Macromolecules 2000, 33, 7207−7218. (24) Allen, M. P.; Tildesley, D. J. Computer simulation of liquids; Clarendon Press: Oxford, U.K., 1987; pp 64−65. (25) Biswas, A. C. The law of rectilinear diameter for the liquid-gas phase transition. Pramana 1973, 1, 109−111. (26) Rowlinson, J. S. Liquids and liquid mixtures: Butterworths monographs in chemistry, 3rd ed.; Butterworth Scientific: London, 1982; p 328. (27) Rackett, H. G. Equation of State for Saturated Liquids. J. Chem. Eng. Data 1970, 15, 514−517. (28) Zhuravlev, N. D.; Martin, M. G.; Siepmann, J. I. Vapor-liquid phase equilibria of triacontane isomers: Deviations from the principle of corresponding states. Fluid Phase Equilib. 2002, 202, 307−324. (29) Potoff, J. J.; Bernard-Brunel, D. A. Mie Potentials for Phase Equilibria Calculations: Applications to Alkanes and Perfluoroalkanes. J. Phys. Chem. B 2009, 113, 14725−14731. (30) Nikitin, E. D. The critical properties of thermally unstable substances: Measurement methods, some results and correlations. High Temp. 1998, 36, 305−318. (31) Rowley, R. L.; Wilding, W. V.; Oscarson, J. L.; Knotts, T. A.; Giles, N. F. DIPPR Data Compilation of Pure Chemical Properties; Design Institute for Physical Properties, AIChE: New York, 2013. (32) Pamies, J. C.; Vega, L. F. Critical properties of homopolymer fluids studied by a Lennard-Jones statistical associating fluid theory. Mol. Phys. 2002, 100, 2519−2529.

REFERENCES

(1) Bhirud, V. L. Saturated Liquid Densities of Normal Fluids. AIChE J. 1978, 24, 1127−1131. (2) Chueh, C. F.; Swanson, A. C. Estimating Liquid Heat-Capacity. Chem. Eng. Prog. 1973, 69, 83−85. (3) Tsonopoulos, C. Empirical Correlation of Second VirialCoefficients. AIChE J. 1974, 20, 263−272. (4) Sastri, S. R. S.; Rao, K. K. A new temperature-thermal conductivity relationship for predicting saturated liquid thermal conductivity. Chem. Eng. J. 1999, 74, 161−169. (5) Lee, B. I.; Kesler, M. G. Generalized Thermodynamic Correlation Based on 3-Parameter Corresponding States. AIChE J. 1975, 21, 510− 527. (6) Przedziecki, J. W.; Sridhar, T. Prediction of Liquid Viscosities. AIChE J. 1985, 31, 333−335. (7) Chickos, J.; Wang, T.; Sharma, E. Hypothetical thermodynamic properties: Vapor pressures and vaporization enthalpies of the even nalkanes from C-40 to C-76 at T = 298.15 K by correlation-gas chromatography. Are the vaporization enthalpies a linear function of carbon number? J. Chem. Eng. Data 2008, 53, 481−491. (8) Nikitin, E. D.; Popov, A. P. Critical temperatures and pressures of C-40, C-44, and C-60 normal alkanes measured by the pulse-heating technique. Fluid Phase Equilib. 2014, 379, 191−195. (9) Nikitin, E. D.; Pavlov, P. A.; Popov, A. P. Vapour-liquid critical temperatures and pressures of normal alkanes with from 19 to 36 carbon atoms, naphthalene and m-terphenyl determined by the pulseheating technique. Fluid Phase Equilib. 1997, 141, 155−164. (10) Tsonopoulos, C.; Tan, Z. M. The Critical Constants of Normal Alkanes from Methane to Polyethylene 0.2. Application of the Flory Theory. Fluid Phase Equilib. 1993, 83, 127−138. (11) Teja, A. S.; Lee, R. J.; Rosenthal, D.; Anselme, M. Correlation of the Critical Properties of Alkanes and Alkanols. Fluid Phase Equilib. 1990, 56, 153−169. (12) Lemmon, E. W.; Goodwin, A. R. H. Critical properties and vapor pressure equation for alkanes CnH2n+2: Normal alkanes with n ≤ 36 and isomers for n = 4 through n = 9. J. Phys. Chem. Ref. Data 2000, 29, 1−39. (13) Nannoolal, Y.; Rarey, J.; Ramjugernath, D. Estimation of pure component properties Part 2. Estimation of critical property data by group contribution. Fluid Phase Equilib. 2007, 252, 1−27. (14) Messerly, R. A.; Knotts, T. A.; Rowley, R. L.; Wilding, W. V. An improved approach for predicting the critical constants of large molecules with Gibbs Ensemble Monte Carlo simulation. Fluid Phase Equilib. 2016, 425, 432−442. (15) Dinpajooh, M.; Bai, P.; Allan, D. A.; Siepmann, J. I. Accurate and precise determination of critical properties from Gibbs ensemble Monte Carlo simulations. J. Chem. Phys. 2015, 143, 114113. (16) Vetere, A. Estimation of Critical Pressures by the Rackett Equation. Chem. Eng. Sci. 1989, 44, 791−795. (17) Nath, S. K.; Escobedo, F. A.; de Pablo, J. J. On the simulation of vapor-liquid equilibria for alkanes. J. Chem. Phys. 1998, 108, 9905− 9911. (18) Errington, J. R.; Panagiotopoulos, A. Z. A new intermolecular potential model for the n-alkane homologous series. J. Phys. Chem. B 1999, 103, 6314−6322. (19) 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. (20) Martin, M. G. MCCCS Towhee. http://towhee.sourceforge.net (accessed June 12, 2013). (21) Messerly, R. A.; Rowley, R. L.; Knotts, T. A.; Wilding, W. V. An improved statistical analysis for predicting the critical temperature and critical density with Gibbs ensemble Monte Carlo simulation. J. Chem. Phys. 2015, 143, 104101. (22) Siepmann, J. I.; Frenkel, D. Configurational Bias Monte-Carlo a New Sampling Scheme for Flexible Chains. Mol. Phys. 1992, 75, 59− 70. (23) Wick, C. D.; Siepmann, J. I. Self-adapting fixed-end-point configurational-bias Monte Carlo method for the regrowth of interior 3649

DOI: 10.1021/acs.jced.6b00574 J. Chem. Eng. Data 2016, 61, 3640−3649