Phase Equilibria and Critical Point Predictions of Mixtures of Molecular

Subscriber access provided by University of Idaho Library

Article

Phase Equilibria and Critical Point Predictions of Mixtures of Molecular Fluids using Grand Canonical Transition Matrix Monte Carlo Tamaghna Chakraborti, and Jhumpa Adhikari Ind. Eng. Chem. Res., Just Accepted Manuscript • Publication Date (Web): 15 May 2017 Downloaded from http://pubs.acs.org on May 22, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Industrial & Engineering Chemistry Research is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 56

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Phase Equilibria and Critical Point Predictions of Mixtures of Molecular Fluids using Grand Canonical Transition Matrix Monte Carlo Tamaghna Chakraborti, Jhumpa Adhikaria Department of Chemical Engineering, Indian Institute of Technology Bombay, Powai, Mumbai – 400476.

Revised manuscript submitted to Industrial & Engineering Chemistry Research May, 2017 a

Electronic mail: [email protected]

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract Mixtures of molecular fluids are encountered often in chemical industry and the vapor-liquid coexistence curves are required for the efficient design of separation processes. Hence, a computational scheme employing the Wang-Landau and grand canonical transition matrix Monte Carlo molecular simulation techniques is used for the first time to compute the fluid phase diagrams of molecular mixtures. Binary mixtures of n-alkanes: viz., methane-ethane, ethanepropane, propane - n-butane systems; are studied here as examples of molecular mixtures. The total molecule number probability distribution, obtained from the GC-TMMC simulations, allows us to directly calculate coexistence properties such as pressure, density and composition. On comparison with experimental data, reported in literature, it was observed that the predicted VLE properties, with the exception of pressure, are in close agreement with the corresponding experimental values. The deviation of calculated pressure from experiments could be attributed to the use of the TraPPE-UA force field where such behavior has already been reported for pure component systems. The above methodology has been coupled to the Binder cumulant intersection technique to derive the critical properties of the mixture, when one of the components is supercritical. The critical properties obtained from our simulations are a good match to those obtained from experiments. Thus, our work shows that the transition matrix methods are able to predict properties in the critical region with remarkable accuracy where other traditionally employed simulation methodologies are not very successful.

Keywords: transition matrix Monte Carlo, Wang Landau simulations, binary mixtures, n-alkanes, Binder cumulant

2 ACS Paragon Plus Environment

Page 2 of 56

Page 3 of 56

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1. Introduction Study of phase behavior of different substances is important as this enables us to understand the basic separation processes occurring in nature as well as in industry. It is important to have a molecular understanding of the phase equilibria phenomena as we can then predict the required phase behavior using parameters which describe the interactions at the level of atoms or molecules. Molecular simulation techniques have been proven to be useful to gain meaningful insights into macroscopic properties from microscopic behavior. Simulations can be performed to determine macroscopic properties using an intermolecular force field1–9 as well as estimate force field parameters from macroscopic behavior.10–17 Several different simulation techniques exist for evaluating phase equilibria, including the Gibbs Ensemble Monte Carlo (GEMC),18–21 the grand canonical transition matrix Monte Carlo (GC-TMMC)22–25 and the Gibbs-Duhem integration (GDI)26–28 techniques. The GEMC method is the most commonly used procedure primarily due to its inherent simplicity and has been used to derive phase equilibrium properties for different systems with single, binary and ternary phases using various force fields given in literature such as Transferable Potentials for Phase Equilibria (TraPPE),10–13,15–17,29 Condensedphase Optimized Molecular Potentials for Atomistic Simulation Studies (COMPASS)30,31 and Universal force field (UFF).31,32 However, the GEMC technique has several inherent limitations. The GEMC simulation method, being in the same mold as ordinary MC simulations, has relatively reduced efficiency because information from rejected moves is ignored. Being a multibox simulation algorithm, GEMC is prone to box swaps between different phases close to the critical point due to reduced barrier potential between different phases. Also, for some calculations such as the evaluation of the critical point or the surface tension values at different temperatures, there arises the requirement for construction of the entire phase space probability 3 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

distribution33–36 as also the information of scaling behavior of different intensive properties with system size. In all of the above cases, newer molecular simulation techniques that give the entire phase space probability distribution, such as the Wang-Landau (WL)37–41 algorithm and GCTMMC, are required. We note here that recently, Dinpajooh et al.42 have calculated the density probability distribution of each box near criticality to evaluate critical and near-critical phase properties via GEMC simulations. Molecular fluids have additional internal degrees of freedom and these fluids entail coupling of the previously discussed WL+GC-TMMC techniques to configurational bias Monte Carlo (CBMC)3,19,43 and expanded ensemble (EE)44–46 techniques. This increases the computational time and, hence, substantial efforts are required to improve computational efficiency by parallelization or via modifications in the algorithm to minimize the number of accessible states at the conditions of the simulation. The algorithm for the GC-TMMC method, as implemented by Shen and Errington,8 helps in reducing the number of accessible states by sampling points close to the equilibrium compositions (as opposed to the rigorous sampling of the entire phase space distribution by the algorithm implemented in other literature24). This results in significant reduction of computational time. The GC-TMMC technique has been used in the past for computing the phase envelopes for pure component monoatomic and complex molecular systems, and for mixtures of atomic systems. Errington22,23 extended the transition matrix Monte Carlo (TMMC) scheme of Smith and Bruce47 to continuous systems in the isothermal-isobaric and grand canonical ensembles and computed the phase diagram as well as the surface tension of Lennard-Jones (LJ) molecules at different temperatures. GC-TMMC has, since then, been used to compute the vapor-liquid equilibria (VLE) curves and surface tension values of square-well9 and variable range Yukawa48 fluids. Singh and Errington25 extended the above technique to compute the VLE curves of n-alkanes in 4 ACS Paragon Plus Environment

Page 4 of 56

Page 5 of 56

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


conjunction with a force field where the van der Waals interactions are modeled via the modified Buckingham exponential-6 intermolecular potential. Paluch et al.49 compared the GC-TMMC and the GEMC techniques for the case of single component molecular fluids and found convergence of intensive quantities such as the coexistence pressure to be significantly faster in case of GC-TMMC. Shen and Errington24,50 extended the technique to two-component monoatomic systems and calculated the pressure - composition (P-x-y) diagrams as well as surface tension values for a variety of symmetric, size-asymmetric mixtures with different combining rules for unlike molecules. Recently, the effect of intermolecular interaction range on binary phase diagrams of symmetric Yukawa fluids, including the effect on azeotropic to heteroazeotropic transition, has also been investigated using the same algorithm.51 GC-TMMC has also been used to study confined52,53 as well as interfacial54,55 systems. In the current work, the algorithm developed by Shen and Errington8 has been used to compute the phase diagrams of binary mixtures of n-alkanes. The original GC-TMMC method24 is not used as it includes additional (and, also, unnecessary) calculations at regions of low probability in phase space as well as the fact that different data points on the phase diagrams of these binary mixtures can be calculated in parallel using the current method. Linear alkanes are among the most widely studied organic compounds both experimentally as well as using simulations. Kay et al.56–62 studied mixtures of various hydrocarbons experimentally and determined the pressure – temperature (P – T) as well as the density– temperature (ρ – T) diagrams for various hydrocarbon mixtures at different compositions. Price and Kobayashi63 reported the vapor-liquid equilibrium (VLE) properties at low temperatures for light hydrocarbon systems. The volumetric and phase properties of methane – n-heptane system have been calculated by Reamer et al.64 Dymond and Robertson65 analyzed pure as well as 5 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

binary mixtures of n-alkanes and found their VLE properties at intermediate pressures and temperatures. Seong et al.66 computed the P-x-y diagrams of the propane and n-butane binary mixture at various temperatures and related it to values obtained from the Peng-Robinson (PR) equation of state using Wong-Sandler mixing rules for the cross terms. The equilibrium properties for asymmetric mixtures of small chain and intermediate chain alkanes have been investigated by Gardeler et al.67 and the experimental data compared with results obtained from equation of state methods. It is noted here that the equilibrium properties of the n-alkane mixtures have been extensively studied and experimental data is available in literature for the purpose of comparison with the results of the simulations performed in this study. Extensive research has been reported in literature for obtaining the equilibrium properties from equation of state methods using various mixing rules for the cross terms. Orwoll and Flory68 employed statistical mechanics to determine an equation of state for prediction of the excess thermodynamic properties of these alkane mixtures. The Predictive Soave-Redlich-Kwong (PSRK) equation of state from UNIFAC has been established by Holderbaum and Gmehling69 to compute equilibrium properties of various mixtures of n-alkanes enabling comparison with experiments. McCabe and Jackson70 used the Statistical Associating Fluid Theory – Variable Range (SAFT-VR ) to model the phase equilibrium of long chain alkanes and their mixtures with small molecules like methane and ethane, which ultimately led to the description of polymeric molecules using SAFT-VR. The critical behavior in asymmetric binary mixtures of alkanes using SAFT has been investigated by Blas and Vega71 who also reported partial miscibility phenomena in the methane – n-hexane system. Pàmies and Vega72 considered VLE and critical behavior of heavy n-alkanes and their mixtures with ethane using the soft-SAFT equation of state with transferable parameters. A universal group contribution equation of state (GCEOS) for the 6 ACS Paragon Plus Environment

Page 6 of 56

Page 7 of 56

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


prediction of VLE of asymmetric systems was designed to estimate the equilibrium properties of mixtures of heavy alkanes with methane and ethane.73 Both experiments and equation of state methods have been employed by Schwarz and Nieuwoudt74 to determine the properties of the propane-dotriacontane binary mixture between 378.2 K and 408.2 K. Polishuk et al.75 estimated vapor-liquid-liquid equilibria (VLLE) of methane – n-alkane and ethane – n-alkane binary mixtures using semi predictive models such as PSRK, linear combination of Vidal mixing rules (LCVM) and the global phase diagram approach (GPDA) techniques. Voutsas et al.76 evaluated the different equations of state, including the Peng-Robinson (PR) equation of state, the SanchezLacombe (SL) equation of state as well as SAFT and concluded that complexity of the equation offers no significant advantage over the simple equations of state in case of the n-alkane systems. The capacity of the PR equation of state combined with COSMO-SAC (Conductor-like screening model – segmented activity coefficients) liquid activity coefficient model using WongSandler mixing rules to predict the properties of n-alkane mixtures as well as mixtures of nalkanes with other molecular fluids like alcohols and ketones has been examined by Lee and Lin.77 Although very efficient from the point of view of computational time required, the equation of state methods do not provide insights into the molecular origins of phase transitions. These methods also perform poorly near the critical point where the presence of long-range fluctuations imply that mean-field theories can no longer be used to determine phase equilibrium properties. With the recent advances in computational power, molecular simulation techniques have taken the centre-stage in the analysis of VLE of alkane mixtures. GEMC has been employed by de Pablo and co-workers78 to compute the phase diagrams of various binary mixtures containing methane as one of the components. Martin and Siepmann79 used CBMC in the Gibbs ensemble 7 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

to compute (i) the boiling point diagram of a binary mixture of n-octane and n-dodecane at 20 kPa, (ii) the isotherm of a binary mixture of ethane and n-heptane at 366 K and (iii) the free energies of transfer of n-pentane and n-hexane from a helium vapor phase at 1atm to a n-heptane liquid phase. Mackie et al.80 reported the vapor-liquid phase diagrams for methane – n-pentane and methane – n-dodecane binary systems using different potential models and compared the resultant phase diagrams with experimental data. Nath et al.46 proposed a new force field to describe binary mixtures of long and short alkanes, and performed simulations in the constant pressure Gibbs ensemble to calculate coexistence curves and in the isothermal isobaric ensemble to calculate the Henry’s law constant. Delhommelle et al.81 compared the performances of two different types of force fields, the TraPPE-United Atom (UA) model developed by the Siepmann group10 and the anisotropic united atom (AUA) model of Toxvaerd et al.82 The vapor-liquid phase diagrams of binary and ternary mixtures containing methane, ethane and carbon dioxide using GEMC have been reported by Liu and Beck.83 Lisal et al.84 developed a new technique called reaction Gibbs ensemble Monte Carlo (RGEMC) for the computer simulation of VLE of multi-component complex fluid mixtures and employed this method to compute phase diagrams of different binary systems. Martin et al.85 used GEMC for the purpose of calculating the molecular structure (i.e., the radial distribution function) as well as the coexistence curves of a binary mixture of n-heptane and supercritical ethane. Coexistence curves have also been calculated using canonical molecular dynamics simulations and compared with existing GEMC calculations by Pàmies et al.86 In the current study, a computational scheme (using WL+GC-TMMC simulations) has been employed to determine the phase coexistence curves of mixtures of molecular fluids, viz., the methane – ethane, ethane – propane and propane – n-butane binary systems, at three different 8 ACS Paragon Plus Environment

Page 8 of 56

Page 9 of 56

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


temperatures, one of which is above the critical point of one of the components. The VLE data, thus, generated is then used in conjunction with Binder's fourth order cumulant intersection method to predict the critical pressure, critical density and critical composition. This paper is organized as follows. In Section 2, the TraPPE-UA force field, which is employed to describe the interactions present in the alkane mixtures, is discussed. In Section 3, the simulation methodology adopted in this paper is reviewed and the details of the simulations used to compute the results, are reported. This is followed by analysis of the results in Section 4. The conclusions of this study are summarized in Section 5. 2. Force Field Model The first and most important step of any molecular simulation is the choice of the force field or the interatomic potential as the correctness of the properties predicted by any molecular simulation technique depends on the accuracy of the force field selected to represent the molecular interactions in the mixture. The TraPPE-UA force field, developed by Martin and Siepmann,10 is employed in the current study because this force field is known to give excellent predictions for VLE properties of pure component alkanes31 and the parameters are known to be highly transferable specifically for prediction of phase equilibria. One of the main features of the TraPPE-UA force field is that the methyl (CH3-) and the methylene (-CH2-) groups are modeled as united atoms and the non-bonded interactions involving the united atom are represented by only a single set of values of the interaction parameters: collision diameter (σ) and well-depth (ε) of the Lennard-Jones (LJ) interaction which is given by,

 σ ij 12  σ ij 6  uij (r ) = 4ε ij   −     r   r  


(1)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 56

The non-bonded interaction parameters for the TraPPE force field used in this study are listed in Table 1. The Lorentz-Berthelot (LB) mixing rules are employed to describe the interaction parameters between unlike united atoms. Electrostatic interactions are absent in the systems studied due to the absence of partial charges on the united atoms in this work. The bonded interactions in TraPPE-UA are semi-flexible as bond lengths are rigid while bond angles and torsions are flexible. The bond lengths considered in this work are listed in Table 2(a). The bond bending energies (Eq. 2) are assumed to vary harmonically with bond angle, being dependent on two interaction parameters, the equilibrium bond angle (θ0) and the force constant (kθ). The bond bending energy is as shown below: -

1 2 uθ (θ ) = kθ (θ − θ0 ) 2

(2)

The values of these two interaction parameters employed in this work are reported in Table 2(b). The dihedral angle interaction energies in TraPPE are assumed to follow a cosine series function of the form, uφ (φ ) = c0 + c1 (1 + cos φ ) + c2 (1 − cos 2φ ) + c3 (1 + cos 3φ )

(3)

The parameters for the torsional energy potential used in this work are given in Table 2(c). As intramolecular non-bonded interactions appear only for those atoms which are separated by more than three bonds and only molecules up to n-butane in the n-alkane homologous series have been considered, intramolecular non-bonded interactions have not been considered. We note here that the complexity of the bonded interactions increases from none in methane (modeled as one UA in the TraPPE-UA force field) to the presence of a rigid bond in ethane, bond bending as well as rigid bonds in propane and additional torsional interactions in butane. 10 ACS Paragon Plus Environment

Page 11 of 56

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


3. Simulation Methodology and Details 3.1 WL + GC-TMMC Simulation Methodology The calculation of phase coexistence properties requires the equilibration of temperature, pressure and chemical potential between the different coexisting phases. In GEMC, this is incorporated into the simulations by virtue of the different simulation moves like particle displacement, volume change and particle exchange between the different boxes. In contrast, GC-TMMC simulations are performed in a single box at constant temperature and chemical potential, and only the pressure is equilibrated post simulations using histogram reweighting. A brief background of the GC-TMMC methodology is presented below. The GC-TMMC method implemented here is the flat histogram method for evaluating multicomponent phase equilibria as described by Shen and Errington.8 The method adopted in this paper is unlike the more rigorous technique given by the same authors in a previous publication24 and implemented by us earlier,51 where the probability distribution of each state point (N1, N2) in the two component particle number space was calculated and based on the location of the maxima and the corresponding minima, regions belonging to different phases were ascertained. The reason for the selection of this new simulation technique here is that the previous method was computationally intensive and could be applied in practice only to binary mixtures of atomic fluids. In the present implementation, the ratio of the fugacity of the different species constitutes a fixed variable and the total number of molecules of both species together is taken as a macrostate. This reduces the previously encountered two-dimensional problem to a onedimensional one. The details of the methodology are given in literature87 and the important steps are described briefly below.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 56

A collection matrix C is maintained and updated where transition from macrostate S to macrostate T (in our case, the macrostate being the total particle number) is due to transfer from microstate s to microstate t, C S ,T = C S , T + b s , t C S , S = C S , S + 1 − bs , t

(4)

π bs,t = min 1 , t  πs  

(5)

where

πs being the probability of microstate s which is defined in the grand canonical ensemble as,

e βµ1N1 e βµ2 N2 q1 (T ) N1 q2 (T ) N 2 V N1 + N 2 e πs = N1 ! N 2 !Ξ ( µ1 , µ2 ,V , T )

(

− βU s N1+ N 2 ;V

)

where,

N i = Number of molecules of species i with chemical potential µi , V = Volume of the system, T = Temperature, = 1/, U = Internal energy Ξ = Grand canonical partition function, and, qi(T) = Partition function of single molecule of species i. 12 ACS Paragon Plus Environment

(6)

Page 13 of 56

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


After completion of the simulation, the transition probabilities between two macrostates are obtained from the C-matrix as, C S ,T

PS ,T =

∑C

S ,T '

(7)

T'

The transition probabilities are subsequently used to obtain the actual macrostate probabilities using detailed balance. ΠS PS ,T = ΠT PT ,S

(8)

The total particle number probability distribution (which will be a bimodal distribution for a twophase system), will be used to calculate the pressure of each phase using the well-known statistical thermodynamic equation given below: -



∑ Π ( N ; µ , ∆µ  α

β p (α )V = ln 

{N }∈

1

12

 ,V , T )  − ln Π ( 0; µ1 , ∆µ12 ,V , T ) 

(9)

The separation between the phases is achieved by considering the two points of inflection as the limits of summation over the vapor and liquid phases respectively. This is done because the data within the points of inflection are phase points consisting of both the vapor and liquid phases. Ignoring this detail leads to little error when calculating VLE data away from the critical point but may contribute significantly near the critical point. The difference of the pressures calculated for the different phases will be non-zero for all values of the chemical potential, µ1, except at coexistence, µ1coex. The probability distributions at this value of the coexistence chemical potential are calculated using histogram reweighting. 13 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ln Π ( N ; µ1coex , ∆µ12 , V , T ) = ln Π ( N ; µ1 , ∆µ12 , V , T ) + β N ( µ1coex − µ1 )

Page 14 of 56

(10)

The probability distributions can be used to calculate the equilibrium properties and obtain the phase diagram. It is noted here that standard Monte Carlo averaging has been performed to obtain the individual molecule numbers of each species (N1 or N2) at each value of the total molecule number (N = N1 + N2) and separate book-keeping steps have been added in the code for this. The species molecule number average in a particular phase is obtained using the following relation,

N

(α ) i

∑ N Π(N) = ∑α Π ( N ) i

N ∈α

N∈

(11)

where, Ni = Standard MC average number of molecules of species i as obtained in the semi-grand

ensemble. The intrinsic properties of the mixture such as mole fraction of species i in the vapour phase (yi) or liquid phase (xi), and density ( ρ mix ) have been calculated using the following equations.

xi or yi =

Ni(α )

∑ Nα

( ) i

i

ρ mix =

∑

(12)

N i(α )

i

V


(13)

Page 15 of 56

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The success of a TMMC simulation depends on proper sampling of the entire domain of the phase space. Multicanonical sampling allows the acceptance ratio to be biased using weighting functions such that all regions of phase space are sampled properly. It has been observed88 that having an approximate idea of the weighting function before a simulation run reduces time required for convergence of probability distribution in the GC-TMMC run. WL simulations provide an alternative technique for the calculation of the probability distributions at different macrostates. In WL simulations,89 the acceptance probabilities, of moves that change from one macrostate to the other, are such that all macrostates are visited with equal probability. However, these acceptance probabilities are unknown and have to be determined by updating a weighting function while the simulation is in progress. The weighting function corresponding to a particular macrostate is updated by a constant factor every time the system exists in that particular macrostate, S.90

ln Π(S) = ln Π(S) + ln f

(14)

where f is an arbitrary predefined constant, fixed during the simulation. We also update a histogram, H(S), containing information about the number of steps for which the system has been in a particular macrostate, when the system occurs in the particular macrostate.

H(S) = H(S) +1

(15)

After the above histogram has achieved a degree of flatness, which is defined to be the difference between the maximum and the minimum values of the histogram being equal to 40% of the

(

minimum value, the value of f is reduced f =

)

f , f > 1 and the simulation restarted with the

reduced value of f. This is continued till f reaches the prescribed minimum value depending on



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 56

the time taken for simulations. The weighting functions required for input to GC-TMMC simulations are given by Π(S). It is noted here that the above formulation of the WL+GC-TMMC technique for calculation of the vapor-liquid coexistence properties has not been used previously, to the best of our knowledge, for the computation of phase equilibria properties of molecular mixtures. 3.2 Critical Point The calculation of the entire phase space probability distribution in the WL+GC-TMMC method allows for the evaluation of critical properties. Different techniques are available in literature for the calculation of the critical properties from the phase space probability distributions. The most comprehensive method to determine the critical point has been given by Bruce and Wilding36,91 who calculated the entire probability distribution in the (N1, N2, E) space and obtained the critical parameters from mixed field finite size scaling analysis. This technique has been further extended to two-component systems by Potoff and Panagiotopoulos.34 The more widely used and easily understandable method is by using the Wegner scaling law and the law of rectilinear diameters to determine the critical temperature and density respectively. The critical pressure is subsequently obtained by invoking a Clausius-Clapeyron-type or an equivalent equation such as the Antoine equation for single component systems. This technique has been used widely by many authors to calculate the critical properties10–17 for the VLE of pure component systems. For the determination of the critical point in the pressure – composition (P-x-y) diagrams for binary systems, suitable scaling laws to be employed for the computations of critical pressure and critical composition are yet to be developed although several attempts have been made to calculate the critical composition in the constant temperature plane.92 The technique we have


Page 17 of 56

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


employed for the computation of critical properties is by use of Binder’s fourth order cumulant.93–95 This method requires as input the probability distribution of the total molecule number in molecule number phase space which is given by the WL+GC-TMMC technique. The cumulant is given by,

UL = 1−

M4 3 M

L 2 2

(16)

L

where M is the order parameter for the determination of the critical point of VLE and is defined as, M =ρ− ρ

(17)

At the critical point, the probability distribution is, after suitable scaling, the same for all systems belonging to the Ising universality class96 and therefore, the Binder cumulant, which reflects the properties of the probability distribution, tends to a universal constant value of approximately 0.4655 for all systems belonging to the Ising universality class with similar boundary conditions. For temperatures above the critical temperature, the probability distribution tends to a Gaussian distribution for infinite system size and the Binder cumulant becomes zero. For temperatures below the critical temperature, the probability distribution becomes bimodal and at infinite system size, it tends to become confined to two very narrow regions centered on two values corresponding to the coexisting vapor and liquid phases. The value of the Binder cumulant then tends to a value of 2 3 . For finite systems, the value of the Binder cumulant varies continuously as we pass through the point of criticality and the point of intersection of the different cumulant curves at different system sizes of the simulation box give an estimate of the critical point.94 In case of binary systems, for phase diagrams being drawn at constant temperature, the Binder 17 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 56

cumulant technique is used to obtain an estimate of the critical pressure (PC), critical composition (xC) and critical density (ρC) of the different mixtures at the given value of temperature. 3.3 Simulation Details The alkane mixtures studied are methane-ethane, ethane-propane and propane-butane. The GCTMMC simulations have been used to verify the accuracy of the TraPPE-UA force field by setting the simulation conditions according to the data available in literature. For propane-butane mixtures, the simulations have been run at 320 K, 353 K and 373.15 K. The temperatures 320 K and 353 K are below the critical point of propane (368 K according to pure component TraPPEUA predictions10) while 373.15 K is above 368 K. The reason behind the choice of these temperatures is the availability of experimental data at these temperatures which allows us to verify the applicability of TraPPE-UA for fluid phase equilibria of n-alkane mixtures. One of the temperatures has been selected above the critical point of the more volatile component to investigate the proficiency of GC-TMMC in capturing closed-loop phase diagrams and in the estimation of the critical point. For the ethane-propane system (critical point of ethane modeled using the TraPPE-UA force field is 304 K10), the temperatures at which the coexistence curve have been computed, are 255.37 K, 283.15 K and 329.82 K. For the methane-ethane system (critical point of methane described by the TraPPE-UA force field is 191.4 K10), the temperatures are 172.04 K, 186.09 K and 199.93 K. The box-lengths for the different simulation runs are given in Table 3. The data has been generated at 22 different values of the fugacity ratio which are so selected such that the entire phase diagram is computed. All simulations are done at arbitrary chemical potentials of the smaller molecule and the probability distributions obtained are subsequently modified using the technique of histogram reweighting to coexistence values. 18 ACS Paragon Plus Environment

Page 19 of 56

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The standard cut-off radius of 1.4 nm as prescribed by Martin and Siepmann,10 is used for the purpose of potential energy truncation. Standard tail corrections are also applied. The minimum value of the total molecule number is 0 and the maximum value is a few molecules greater than the molecule number corresponding to the larger pure component liquid density. Initial estimates for the weighting functions are evaluated using WL simulations. All simulations are terminated once all macrostates have been visited at least 1,000,000 times. For the estimation of the critical point parameters, the point of intersection of the Binder cumulant at different lengths should be taken as the critical point. However, due to the presence of finite size effects, the points of intersection corresponding to different lengths are close but do not coincide. This is consistent with the observation of system size dependence reported by Pérez-Pellitero et al.93 Hence, a simple averaging is performed to obtain the point of criticality. The WL simulations, which have been run to get an estimate of the initial weighting functions for the GC-TMMC runs, are started with an initial value of ln f = 1 and the runs were terminated when the values of ln f reached values of the order of 10-5. The value of f was reduced to

f

when the difference of the maximum and the minimum values of H(S) reached 40% of the minimum value of H(S).

4. Results and Discussion We have employed the WL + GC-TMMC simulations to predict the coexistence of n-alkane mixtures (studied here as examples of molecular fluids). The coexistence properties of the binary molecular systems are directly calculated from the complete probability distribution of the molecular number obtained from the GC-TMMC simulations.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 56

4.1 Methane-Ethane System Figure 1(a) shows the P-x-y diagram (with x and y representing the mole fractions of the more volatile component, methane) at temperatures of 172.04 K and 186.09 K. Both temperatures are below 191.4 K, the critical point of methane when methane is modeled as TraPPE-UA pseudoatom.10 The predictions obtained from the WL+GC-TMMC simulations using the TraPPE-UA force field with LB combining rules are in good agreement with experiment for the saturated liquid phase; though significant deviations are observed in the coexisting vapor phase composition. The agreement between simulation and experiment in the vapor phase is seen to improve in the pure methane limit where TraPPE is reasonably successful in predicting the phase behavior of the pure component.10 There are, however, considerable deviations in the calculated coexistence pressures in the pure ethane limit which are not visible in Figure 1(a) due to the large pressure scale used to represent the data. The predictions of TraPPE are higher by about 50% than the experimental values of Wichterle and Kobayashi97 in the pure ethane limit, and this results in the phase equilibria curve shifting upwards. The variations of pressure and composition with the density of the molecular mixture at VLE are illustrated in Figures 1(b) and 1(c). The density of the mixture is seen to vary non-linearly with pressure and composition; and passes through a maximum for the liquid phase as one goes from the pure ethane limit to the pure methane limit. There are no experimental data for comparison here. Figures S1 and S2 (in the supporting information) show the y-x plot of the methane-ethane system at these two temperatures where slight differences with experimental data are observed. These slight differences have been attributed by some authors98 to the unsuitability of the LB combining rules in describing the force field parameters for the unlike pseudo-atom interactions. The agreement of the simulation predicted separation factors with experimental data has been reported for 20 ACS Paragon Plus Environment

Page 21 of 56

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


mixtures of organic molecules modeled using the TraPPE-UA force field by Wick et al.12 for the n-pentane – benzene system at 318.15 K and by Stubbs et al.13 for a mixture of diethyl ether and ethanol at 283.15 K. It is noted here that the parameterization of the TraPPE-UA force field and the combining rules has not utilized any experimental binary VLE data. The P-x-y diagram in Figure 2(a) is for the methane-ethane system at a temperature of 199.93 K which is above the critical temperature of methane. A closed loop behavior is expected at this temperature with the phase diagram ending in a critical point. A good match of the liquid phase data with experimental values is again observed, though discrepancies exist in the vapor phase. The saturation vapor pressure predicted for ethane using the TraPPE-UA force field at this temperature is again found to be higher than actual values as previously noted (Table 4). Figures 2(b) and 2(c) display the pressure-density and the composition - density plots of the methaneethane system at 199.93 K. The densities obtained by our simulation methodology show satisfactory agreement with the equation of state predictions of Muller et al.99 The y-x data obtained from simulations and displayed in Figure S3, also show reasonable match with experiments. This is important since these diagrams are useful in the design of separation equipment in the chemical industry. The calculation of the critical pressure requires the computation of the Binder cumulant which is plotted as a function of mean composition and pressure along the coexistence curve in Figures 3(a) and 3(b) respectively. It is seen that the cumulant curves intersect each other around a value of 0.48 which is close to the value of 0.4655 for systems belonging to the Ising universality class. The value of the critical pressure obtained is 50.71±0.07 bars, which is in agreement with the experimentally observed value of 51.4 bar.97 This is expected since TraPPE-UA gives excellent prediction of the VLE behavior in the pure methane limit. Predictions based on Backone family of equations of state by Müller et al.99 21 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 56

estimate the critical pressure to be 53.28 bars which is higher than both our value and that from experiment. The values of the critical molar composition and the critical density as predicted by our method are 0.9308±0.0011 and 10301±3mol/m3 respectively, which match the molar composition of 0.9553 obtained experimentally and the density of 10538 mol/m3 obtained using equation of state techniques. Also, the densities near the critical region, as plotted in Figures 2(b) and 2(c), show negligible finite size effects for the values of the molar compositions explored in our work. This is further demonstrated in Figure S10(a) where the density difference between the liquid and the vapor phases is plotted versus the difference between the critical and mean mole fractions of the mixture. This might be due to the fact that the methane molecule is not a very large molecule and the simulation box-lengths chosen are sufficiently large to display correct VLE behavior for the chosen values of the molar compositions. This is confirmed by the fact that all the values of the Binder cumulant below the critical composition are close to the value of 0.67 which is the expected value for a system displaying two phase behavior. 4.2 Ethane-Propane System Figure 4(a) shows the pressure-composition diagram for the ethane-propane mixture at 255.37 K and 283.15 K (which are below the critical point of ethane estimated using TraPPE-UA at 304 K10) with the composition reported in mole fractions of the more volatile component ethane. The limitations of the TraPPE-UA force field in correctly estimating the vapor pressures of the pure components (Table 4) and the shortcomings of the LB combining rules (the TraPPE-UA parameters and LB combining rules have not been fitted to any binary experimental VLE data as mentioned previously) result in a P-x-y diagram where the pressure is significantly overestimated with respect to the experimental data reported by Matschke and Thodos.100 Our simulations are only able to qualitatively capture the features of the phase behavior and significant deviations 22 ACS Paragon Plus Environment

Page 23 of 56

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


occur both in the liquid and the vapor phases. The TraPPE-UA estimates of the pure component vapor pressures in both the pure ethane and pure propane limit are found to be significantly higher than the values reported in NIST101 by around 15-20% for ethane and 30-35% for propane; and these deviations decrease with an increase in temperature; an observation also noted by Martin and Siepmann.10 However, the y-x plots, as illustrated in Figures S4 and S5, show reasonable match with experimental data and can be utilized for designing separation equipment. Figures 4(b) and 4(c) show the P-ρ and the x-y-ρ diagrams for the same mixture at the same temperature. We note that, unlike the methane-ethane system, the liquid densities in Figure 4(c) at 255.37 K do not show a non-linear behavior and the values increase monotonically from the pure propane to the pure ethane limit possibly because the difference in molecular sizes for the ethane-propane system is less than that of the methane-ethane system. The P-x-y diagram in Figure 5(a) for the ethane-propane system at a temperature of 329.82 K (which is above the critical temperature of ethane) displays a closed loop phase behavior as expected. The inability of the TraPPE-UA force field to predict the phase behavior with quantitative accuracy is again noted although the qualitative features are captured. The deviation of the vapor pressures in the pure propane limit is 18.4% and this leads to the phase envelope moving upwards to higher pressure. Finite size effects are observed in the vapor phase in Figure 5(a) although finite size effects in the liquid phase are insignificant. Similar finite size effects are also observed in the y-x diagram of the same system in Figure S6 of the supporting information. Distinct finite size effects are, however, observed in pressure-density and the densitycomposition plots of Figures 5(b) and 5(c) near the critical point where significant differences in coexisting vapor and liquid densities are noted with variation in simulation box length. The decrease in simulation box length appears to shift the liquid compositions and densities to higher 23 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 56

values to the right of the critical composition and density respectively, and the vapor compositions and densities to lower values to the left. Thus, the critical pressure appears to shift to a higher value with decrease in system size. The result is further emphasized in Figure S10(b). This effect was also observed by Mon and Binder for two-dimensional Ising lattice-gas systems.102 Finite-size effects are observed near criticality due to the divergence of the correlation length at criticality. The finite size of the simulation box places a limit on the maximum correlation length within a simulation run and causes the free energy to not display the singular behavior that should be observed near criticality in the thermodynamic limit. The result is that the coexistence compositions and densities tend to display a mean-field effect with the singularities being tapered off by the finite-sized simulation box. The value of the critical pressure, thus, determined from simulations in the box of finite length increases with decreasing system size. These simulations can however be used to predict the correct values of the critical point using the Binder cumulant intersection method as shown in Figures 6(a) and 6(b). The intersection of the cumulant curves at different simulation box-lengths occur at a cumulant value of approximately 0.4 which is less than the value of 0.4655 for the Ising criticality class. This suggests that significant finite size corrections persist in the calculated values of the Binder cumulant. The critical pressure calculated using the above methodology is 57.19±0.78 bars while the experimental critical pressure is given by 52.39 bars.100 The overestimation of the critical pressure is consistent with the observed trend of overprediction of pressure while using the TraPPE-UA force field. The critical composition is estimated using our simulations to be 0.7098±0.0103 which is close to the experimental value of 0.7029.100 This indicates that although TraPPE-UA yields incorrect values of the critical pressure, it is far more useful for the purpose of calculating the critical composition. The critical density is determined to be


Page 25 of 56

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


6477.72±19.15 mol/m3 from our simulations; there is, however, no corresponding experimental data reported in literature. 4.3 Propane-Butane System Figure 7(a) plots the phase behavior of the propane-butane binary mixture in the P-x-y plane for 320 K and 353 K, which are again below the critical temperature of propane (368 K) as reported in TraPPE-UA literature.10 A comparison with the experimental data of Kay56 has been provided at 353 K and this shows that the simulation data differs significantly from experiment. The reason for this deviation, as seen from the data in Table 4, is the deviation of pure component vapor pressures at the said temperatures. Unlike the other two mixtures, Kay’s experimental results56 for pressure vs. density and composition vs. density plots are also available. Figures 7(b) and 7(c) give pressure-density and composition-density diagrams of the propane-butane system at the two temperatures. It is observed from the above graphs that, over-estimation of vapor pressures notwithstanding, the TraPPE-UA force field performs remarkably well with regard to calculation of coexistence compositions and densities. The liquid densities calculated using TraPPE-UA force field by our simulation methodology show a reasonable agreement with experimental liquid densities of Kay56 while the vapor densities deviate slightly from the experimental values in Figure 7(c). Comparison of simulation data with experiment on the y-x diagram of the propane-butane system at 353 K in Figure S8 of supporting information also demonstrates that the molar compositions are well predicted using TraPPE-UA. This reestablishes the fact that TraPPE-UA is one of the best potential models for studying phase equilibria using molecular simulations.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 56

The P-x-y diagram in Figure 8(a) is for the propane-butane system at 373.15 K which is above the critical point of propane and a closed loop phase behavior is hence observed. The deviation of pure component vapor pressure of butane (Table 4) is around 21.6% and the phase envelope shifts upwards to higher pressures as seen in the previous cases. Since data on liquid density is available, it is possible to make an effective comparison between experiments and simulations in this case study. Figures 8(b) and 8(c) give the pressure-density and the density-composition diagrams respectively. It is again observed that TraPPE-UA, despite its simplicity and inability to quantitatively predict the pressures at equilibrium, accurately predicts the equilibrium densities and compositions. Further, in Figure 8(b), we note that there are slight deviations of the TraPPE-UA predictions of the vapor densities away from criticality while liquid densities deviate from experimental values near criticality. Finite size effects are distinctly observed in Figures 8(b) and 8(c) and the composition-density experimental data are best reproduced by the largest simulation box length of 6 nm. The cumulant curves in Figures 9(a) and 9(b) at different lengths intersect at a cumulant value of 0.4 which is similar to the intersection obtained for the ethanepropane system at 329.82 K. The critical pressure and composition are 48.65±0.26 bars and 0.9663±0.0026 respectively, while experiments56 report the two values as 42.74 bars and 0.9572. This lends further credence to the use of this technique for effective computation of critical composition. The difference between the estimated critical pressure and the experimentally determined one is again an artifact of the TraPPE-UA force field and its inability to compute pressure accurately. Figures 8(b) and 8(c) give the P-ρ and the x-y-ρ diagrams for the above system and the predicted composition and densities show remarkable agreement with experimentally determined data in Figure 8(c). Figure 8(c) also shows considerable finite size effects and establishes the fact that increasing the system size results in better semblance of


Page 27 of 56

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


simulation results to experiments. A comparison of the critical density obtained via simulations with that from experiments is also seen in Figures 8(b) and 8(c). The predicted value of the critical density at the given temperature using TraPPE-UA force field is 5142±3mol/m3, which is near the experimentally56 obtained critical density of 4927.67 mol/m3. This shows that our simulations are also capable of predicting the critical densities of n-alkane mixtures.

5. Conclusions Real systems encountered in industry generally involve mixtures of molecular fluids and the vapor-liquid coexistence data of these mixtures are required for the efficient design of separation processes and equipment. Towards achieving this goal, the WL+GC-TMMC technique, which has been successfully used to predict the coexistence properties of pure molecular fluids in literature,25,49 has been extended to molecular mixtures (viz., methane-ethane, ethane- propane, propane – n-butane systems) for the first time, to the best of our knowledge, in this work. This computational scheme has been used to determine the entire molecule number probability distribution at coexistence, where the system exhibits a bimodal distribution in total molecule number. Using this bimodal distribution, the different properties of the phase envelope of a binary mixture such as density, composition and pressure are directly computed. It is observed that the methodology employed is able to reproduce the experimental composition and density data with sufficient accuracy for the n-alkane mixtures, which have been studied as examples of mixtures of molecular fluids. Discrepancies are observed in the calculation of pressure which may be attributed to the force field used in the simulations, namely TraPPE-UA. The method is also able to capture the finite size effects near the critical point, when one of the components is supercritical. Since the WL+GC-TMMC methodology entails the calculation of the entire molecule number probability distribution, it is coupled to the Binder cumulant intersection 27 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 56

methodology to determine the critical properties of the system. The deviation from the expected value of the Binder cumulant (0.4655) is due to the presence of finite size effects and this has been observed in literature also.93,103 We find that, within the limitations of the TraPPE-UA force field, the methodology was able to predict the critical properties of the binary molecular mixtures with sufficient accuracy. However, a trade-off with computational efficiency should also be considered as seen in Table S1 of supporting information where the different simulation times taken by the WL and GC-TMMC technique separately on an Intel(R) Xeon (R) CPU E5440 @ 2.83 GHz are reported. It is noted that the WL+GC-TMMC technique needs considerable computational resources for its execution. Thus, if computational efficiency is a primary criteria for determination of the choice of simulation methodology, the technique should only be applied to evaluate the phase coexistence properties near the critical region where multiple box simulation methodologies like the Gibbs ensemble Monte Carlo of Panagiotopoulos18 and the activity fraction expanded ensemble (AFEE) technique of the Errington group55 are not very successful.

7. Acknowledgements Computational resources have been provided by the Computer Centre at Indian Institute of Technology, Bombay and funding for this project has been provided by the Department of Science and Technology, India (Reference #: SB/S3/CE/054/2014).The authors would also like to acknowledge the computational resources provided by the Center for Development of Advanced Computing (CDAC), Pune, India.


Page 29 of 56

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


8. Supporting Information The supporting information includes the following tables and figures: Table S1 with simulation times required for WL and GC-TMMC simulations, listed for each box length considered in the methane - ethane, ethane - propane and propane - butane systems; Table S2 with information on the total number of molecules of each system at the critical point; Figures S1 to S9 showing the y vs x data obtained from our simulations (and comparison with experimental results or equation of state predictions where available in literature) for each temperature considered in the methane ethane, ethane - propane and propane - butane systems and Figure S10 demonstrating the finite size effects in the plots of the density difference between the liquid and the vapor phases vs the difference between the critical and mean mole fractions of the mixture for the methane-ethane system at 199.93 K, the ethane-propane system at 329.82 K and the propane-butane system at 373.15 K for the different simulation box-lengths used in this study.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

References (1)

Ibergay, C.; Ghoufi, A.; Goujon, F.; Ungerer, P.; Boutin, A.; Rousseau, B.; Malfreyt, P. Molecular simulations of the n -alkane liquid-vapor interface: Interfacial properties and their long range corrections. Phys. Rev. E 2007, 75, 51602.

(2)

Eike, D. M. ; Maginn, E. J. Atomistic simulation of solid-liquid coexistence for molecular systems: Application to triazole and benzene. J. Chem. Phys. 2006, 124, 164503 .

(3)

Smit, B.; Karaborni, S.; Siepmann, J. I. Computer simulations of vapor–liquid phase equilibria of n-alkanes. J. Chem. Phys. 1995, 102, 2126 .

(4)

Harini, M.; Adhikari, J.; Rani, K. Y. Prediction of Vapor − Liquid Coexistence Data for p

‑ Cymene Using Equation of State Methods and Monte Carlo Simulations. J. Chem. Eng. Data 2014, 59, 2987–2994 . (5)

Fan, C. F.; Olafson, B. D.; Blanco, M.; Hsu, S. L. Application of molecular simulation to derive phase diagrams of binary mixtures. Macromolecules 1992, 25, 3667–3676 .

(6)

Delhommelle, J.; Boutin, A.; Fuchs, A. H. Molecular Simulation of Vapour-Liquid Coexistence Curves for Hydrogen Sulfide-Alkane and Carbon Dioxide-Alkane Mixtures. Mol. Simul. 1999, 22, 351–368.

(7)

Pàmies, J. C.; McCabe, C.; Cummings, P. T.; Vega, L. F. Coexistence Densities of Methane and Propane by Canonical Molecular Dynamics and Gibbs Ensemble Monte Carlo Simulations. Mol. Simul. 2003, 29, 463–470 .

(8)

Errington, J. R.; Shen, V. K. Direct evaluation of multicomponent phase equilibria using flat-histogram methods. J. Chem. Phys. 2005, 123, 164103 .

(9)

Singh, J. K.; Kofke, D. A.; Errington, J. R. Surface tension and vapor–liquid phase coexistence of the square-well fluid. J. Chem. Phys. 2003, 119, 3405.

(10)

Martin, M. G.; Siepmann, J. I. Transferable Potentials for Phase Equilibria. 1. UnitedAtom Description of n -Alkanes. J. Phys. Chem. B 1998, 102, 2569–2577.

(11)

Chen, B.; Siepmann, J. I. Transferable Potentials for Phase Equilibria. 3. Explicit30 ACS Paragon Plus Environment

Page 30 of 56

Page 31 of 56

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Hydrogen Description of Normal Alkanes. J. Phys. Chem. B 1999, 103, 5370–5379 . (12)

Wick, C. D.; Martin, M. G.; Siepmann, J. I. Transferable Potentials for Phase Equilibria. 4. United-Atom Description of Linear and Branched Alkenes and Alkylbenzenes. J. Phys. Chem. B 2000, 104, 8008–8016 .

(13)

Stubbs, J. M.; Potoff, J. J.; Siepmann, J. I. Transferable Potentials for Phase Equilibria. 6. United-Atom Description for Ethers, Glycols, Ketones, and Aldehydes. J. Phys. Chem. B 2004, 108, 17596–17605.

(14)

Keasler, S. J.; Charan, S. M.; Wick, C. D.; Economou, I. G.; Siepmann, J. I. Transferable potentials for phase equilibria-united atom description of five- and six-membered cyclic alkanes and ethers. J. Phys. Chem. B 2012, 116, 11234–11246.

(15)

Rai, N.; Siepmann, J. I. Transferable potentials for phase equilibria. 10. Explicit-hydrogen description of substituted benzenes and polycyclic aromatic compounds. J. Phys. Chem. B 2013, 117, 273–288.

(16)

Lubna, N.; Kamath, G.; Potoff, J. J.; Rai, N.; Siepmann, J. I. Transferable potentials for phase equilibria. 8. United-atom description for thiols, sulfides, disulfides, and thiophene. J. Phys. Chem. B 2005, 109, 24100–24107.

(17)

Rai, N.; Siepmann, J. I. Transferable potentials for phase equilibria. 9. Explicit hydrogen description of benzene and five-membered and six-membered heterocyclic aromatic compounds. J. Phys. Chem. B 2007, 111, 10790–10799.

(18)

Panagiotopoulos, A. Z. Monte Carlo methods for phase equilibria of fluids. J. Phys.: Condens. Matter 2000, 12, R25–R52.

(19)

Cui, S. T.; Cummings, P. T.; Cochran, H. D. Configurational bias Gibbs ensemble Monte Carlo simulation of vapor-liquid equilibria of linear and short-branched alkanes. Fluid Phase Equilib. 1997, 141, 45–61.

(20)

Vega, L.; de Miguel, E.; Rull, L. F.; Jackson, G.; McLure, I. A. Phase equilibria and critical behavior of square-well fluids of variable width by Gibbs ensemble Monte Carlo simulation. J. Chem. Phys. 1992, 96, 2296. 31 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(21)

Cortés Morales, A. D.; Economou, I. G.; Peters, C. J.; Ilja Siepmann, J. Influence of simulation protocols on the efficiency of Gibbs ensemble Monte Carlo simulations. Mol. Simul. 2013, 39, 1135–1142.

(22)

Errington, J. R. Direct calculation of liquid–vapor phase equilibria from transition matrix Monte Carlo simulation. J. Chem. Phys. 2003, 118, 9915.

(23)

Errington, J. R. Evaluating surface tension using grand-canonical transition-matrix Monte Carlo simulation and finite-size scaling. Phys. Rev. E 2003, 67, 12102.

(24)

Shen, V. K.; Errington, J. R. Determination of fluid-phase behavior using transitionmatrix Monte Carlo: Binary Lennard-Jones mixtures. J. Chem. Phys. 2005, 122, 64508.

(25)

Singh, J. K.; Errington, J. R. Calculation of phase coexistence properties and surface tensions of n-alkanes with grand-canonical transition-matrix Monte Carlo simulation and finite-size scaling. J. Phys. Chem. B 2006, 110, 1369–1376.

(26)

Kofke, D. A. Gibbs-Duhem integration: a new method for direct evaluation of phase coexistence by molecular simulation. Mol. Phys. 1993, 78, 1331–1336.

(27)

Lísal, M.; Vacek, V. Direct Evaluation of Solid–Liquid Equilibria by Molecular Dynamics Using Gibbs-Duhem Integration. Mol. Sim. 1997,19, 43–61.

(28)

Lísal, M.; Vacek, V. Direct Evaluation of Vapour-Liquid Equilibria by Molecular Dynamics using Gibbs-Duhem Integration. Mol. Sim. 1996, 17, 27–39.

(29)

Chen, B.; Potoff, J. J.; Siepmann, J. I. Monte Carlo calculations for alcohols and their mixtures with alkanes. Transferable potentials for phase equilibria. 5. United-atom description of primary, secondary, and tertiary alcohols. J. Phys. Chem. B 2002, 105, 3093–3104.

(30)

Sun, H. COMPASS: An ab initio Force-Field Optimized for Condensed-Phase Applications - Overview with Details on Alkane and Benzene Compounds. J. Phys. Chem. B 1998, 5647, 7338–7364.

(31)

Martin, M. G. Comparison of the AMBER , CHARMM , COMPASS , GROMOS , OPLS 32 ACS Paragon Plus Environment

Page 32 of 56

Page 33 of 56

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


, TraPPE and UFF force fields for prediction of vapor – liquid coexistence curves and liquid densities. Fluid Phase Equilib. 2006, 248, 50–55. (32)

Rappe, A. K.; Casewit, C. J.; Colwell, K. S.; Goddard III, W. A.; Skiff, W. M. UFF, a Full Periodic Table Force Field for Molecular Mechanics and Molecular Dynamics Simulations. J. Am. Chem. Soc. 1992, 114, 10024–10035.

(33)

Caillol, J. M. Critical-point of the Lennard-Jones fluid : A finite-size scaling study. J. Chem. Phys. 2014, 4885, 4885–4893.

(34)

Potoff, J. J.; Panagiotopoulos, A. Z. Critical point and phase behavior of the pure fluid and a Lennard-Jones mixture. J. Chem. Phys. 1998, 109, 10914.

(35)

Wilding, N. B. Critical-point and coexistence-curve properties of the Lennard-Jones fluid: A finite-size scaling study. Phys. Rev. E 1995, 52, 602–611.

(36)

Bruce, A. D.; Wilding, N. B. Scaling Fields and Universality of the Liquid-Gas Critical Point. Phys. Rev. Lett. 1992, 68, 193–196.

(37)

Shell, M.; Debenedetti, P.; Panagiotopoulos, A. Z. Generalization of the Wang-Landau method for off-lattice simulations. Phys. Rev. E 2002, 66, 56703.

(38)

Ngale, K. N.; Desgranges, C.; Delhommelle, J. Wang – Landau configurational bias Monte Carlo simulations : vapour – liquid equilibria of alkenes. Mol. Simul. 2012, 38, 653–658.

(39)

Poulain, P.; Calvo, F.; Antoine, R.; Broyer, M.; Dugourd, P. Performances of WangLandau algorithms for continuous systems. Phys. Rev. E 2006, 73, 56704.

(40)

Maerzke, K. A.; Gai, L.; Cummings, P. T.; McCabe, C. Simulating Phase Equilibria using Wang-Landau-Transition Matrix Monte Carlo. J. Phys.: Conf. Ser. 2014, 487, 012002.

(41)

Wang, F.; Landau, D. P. Efficient, Multiple-Range Random Walk Algorithm to Calculate the Density of States. Phys. Rev. Lett. 2001, 86, 2050–2053.

(42)

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, 33 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

143. (43)

Martin, M. G.; Siepmann, J. I. Novel Configurational-Bias Monte Carlo Method for Branched Molecules. Transferable Potentials for Phase Equilibria. 2. United-Atom Description of Branched Alkanes. J. Phys. Chem. B 1999, 103, 4508–4517.

(44)

Chang, J.; Sandler, S. I. Determination of liquid–solid transition using histogram reweighting method and expanded ensemble simulations. J. Chem. Phys. , 2003, 118, 8390.

(45)

Escobedo, F. A.; de Pablo, J. J. Monte Carlo simulation of the chemical potential of polymers in an expanded ensemble. J. Chem. Phys. 1995, 103, 2703.

(46)

Nath, S. K.; Escobedo, F. A.; de Pablo, J. J.; Patramai, I. Simulation of Vapor - Liquid Equilibria for Alkane Mixtures. Ind. Eng. Chem. Res. 1998, 37, 3195–3202.

(47)

Smith, G. R.; Bruce, A. D. A study of the multi-canonical Monte Carlo method. J. Phys. A:Math. Gen. 1995, 28, 6623.

(48)

Singh, J. K. Surface tension and vapour–liquid phase coexistence of variable-range hardcore attractive Yukawa fluids. Mol. Simul. 2009, 35, 880–887.

(49)

Paluch, A. S.; Shen, V. K.; Errington, J. R. Comparing the Use of Gibbs Ensemble and Grand-Canonical Transition-Matrix Monte Carlo Methods to Determine Phase Equilibria. Ind. Eng. Chem. Res. 2008, 47, 4533–4541.

(50)

Shen, V. K.; Errington, J. R. Determination of surface tension in binary mixtures using transition-matrix Monte Carlo. J. Chem. Phys. 2006, 124, 24721.

(51)

Chakraborti, T.; Adhikari, J. Prediction of fluid-phase behavior of symmetrical binary Yukawa fluids using transition matrix Monte Carlo. Fluid Phase Equilib. 2016, 415, 64– 74.

(52)

Siderius, D. W.; Shen, V. K. Use of the Grand Canonical Transition-Matrix Monte Carlo Method to Model Gas Adsorption in Porous Materials. J. Phys. Chem. C 2013, 117, 5861– 5872. 34 ACS Paragon Plus Environment

Page 34 of 56

Page 35 of 56

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(53)

Singh, J. K.; Kwak, S. K. Surface tension and vapor-liquid phase coexistence of confined square-well fluid. J. Chem. Phys. 2007, 126, 024702.

(54)

Grzelak, E. M.; Errington, J. R. Computation of interfacial properties via grand canonical transition matrix Monte Carlo simulation. J. Chem. Phys. 2008, 128, 14710 .

(55)

Kumar, V.; Errington, J. R. Monte Carlo simulation strategies to compute interfacial and bulk properties of binary fluid mixtures. J. Chem. Phys. 2013, 138, 174112 .

(56)

Kay, W. B. Vapor-Liquid Equilibrium Relations of Binary Systems. n-Butane and nPentane. J. Chem. Eng. Data 1970, 15, 46–52.

(57)

Kay, W. B. The Critical Locus Curve and the Phase Behavior of Mixtures. Acc. Chem. Res. 1968, 1, 344–351.

(58)

Hissong, D. W.; Kay, W. B. The Calculation of the Critical Locus Curve of a Binary Hydrocarbon System. AIChE J. 1970, 16, 580–587.

(59)

Pak, S. C.; Kay, W. B. The Critical Properties of Binary Hydrocarbon Systems. Ind. Eng. Chem. Fundam. 1972, 11, 255–267.

(60)

Kay, W. B. Liquid-Vapor Phase Equilibrium Relations in the Ethane- n-Heptane System. Ind. Eng. Chem. 1938, 30, 459–464.

(61)

Kreglewski, A.; Kay, W. B. The Critical Constants of Conformal Mixtures. J. Phys. Chem. 1969, 73, 3359–3366.

(62)

Kay, W. B. Vapor-liquid Equilibrium Relationships of Binary Systems. J. Chem. Eng. Data 1971, 16, 137–140 .

(63)

Price, A. R.; Kobayashi, R. Low Temperature Vapor-Liquid Equilibrium in Light Hydrocarbon Mixtures: Methane-Ethane-Propane System. J. Chem. Eng. Data 1959, 4, 40–52.

(64)

Reamer, H. H.; Sage, B. H.; Lacey, W. N. Volumetric and Phase Behavior of the Methane-n-Heptane System. Ind. Eng. Chem. 1956, 1, 29–42.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(65)

Dymond, J. H.; Robertson, J. (p , ρ,T ) of some pure n-alkanes and binary mixtures of nalkanes in the range 298 to 373 K and 0.1 to 500 MPa. J. Chem. Thermodyn. 1982, 14, 51–59.

(66)

Seong, G.; Yoo, K.; Lim, J. S. Vapor-Liquid Equilibria for Propane (R290) + n-Butane (R600) at Various Temperatures. J. Chem. Eng. Data 2008, 53, 2783–2786.

(67)

Gardeler, H.; Fischer, K.; Gmehling, J. Experimental Determination of Vapor - Liquid Equilibrium Data for Asymmetric Systems. Ind. Eng. Chem. Res. 2002, 41, 1051–1056.

(68)

Orwoll, R. A.; Flory, P. J. Thermodynamic Properties of Binary Mixtures of n-Alkanes. J. Am. Chem. Soc. 1967, 89, 6822–6829.

(69)

Holderbaum, T.; Gmehling, J. PSRK: A Group Contribution Equation of State Based on UNIFAC. Fluid Phase Equilib. 1991, 70, 251–265.

(70)

McCabe, C.; Jackson, G. SAFT-VR modelling of the phase equilibrium of long-chain nalkanes. Phys. Chem. Chem. Phys. 1999, 1, 2057–2064.

(71)

Blas, F. J.; Vega, L. F. Critical behavior and partial miscibility phenomena in binary mixtures of hydrocarbons by the statistical associating fluid theory. J. Chem. Phys. 1998,

109, 7405. (72)

Pamies, J. C.; Vega, L. F. Vapor - Liquid Equilibria and Critical Behavior of Heavy n Alkanes Using Transferable Parameters from the Soft-SAFT Equation of State. Ind. Eng. Chem. Res. 2001, 40, 2532–2543.

(73)

Ahlers, J.; Gmehling, J. Development of a Universal Group Contribution Equation of State . 2 . Prediction of Vapor - Liquid Equilibria for Asymmetric Systems. Ind. Eng. Chem. Res. 2002, 41, 3489–3498.

(74)

Schwarz, C. E.; Nieuwoudt, I. Phase equilibrium of propane and alkanes Part I. Experimental procedures , dotriacontane equilibrium and EOS modelling. J. Supercrit. Fluids 2003, 27, 133–144.

(75)

Polishuk, I.; Wisniak, J.; Segura, H. Estimation of Liquid-Liquid-Vapor Equilibria in 36 ACS Paragon Plus Environment

Page 36 of 56

Page 37 of 56

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Binary Mixtures of n-Alkanes. Ind. Eng. Chem. Res. 2004, 43, 5957–5964. (76)

Voutsas, E. C.; Pappa, G. D.; Magoulas, K.; Tassios, D. P. Vapor liquid equilibrium modeling of alkane systems with Equations of State: ‘simplicity versus complexity’. Fluid Phase Equilib. 2006, 240, 127–139.

(77)

Lee, M.; Lin, S. Prediction of mixture vapor–liquid equilibrium from the combined use of Peng–Robinson equation of state and COSMO-SAC activity coefficient model through the Wong–Sandler mixing rule. Fluid Phase Equilib. 2007, 254, 28–34.

(78)

de Pablo, J. J.; Bonnin, M.; Prausnitz, J. M. Vapor-liquid equilibria for polyatomic fluids from site-site computer simulations: pure hydrocarbons and binary mixtures containing methane. Fluid Phase Equilib. 1992, 73, 187–210.

(79)

Martin, M. G.; Siepmann, J. I. Predicting Multicomponent Phase Equilibria and Free Energies of Transfer for Alkanes by Molecular Simulation. J. Am. Chem. Soc. 1997, 119, 8921–8924.

(80)

Mackie, A. D., Tavitian, B., Boutin, A.; Fuchs, A. H. Vapour-Liquid Phase Equilibria Predictions of Methane–Alkane Mixtures by Monte Carlo Simulation. Mol. Simul. 1997,

19, 1–15. (81)

Delhommelle, J.; Boutin, A.; Tavitian, B.; Mackie, A. D.; Fuchs, A. H. Vapour-liquid coexistence curves of the united-atom and anisotropic united-atom force fields for alkane mixtures. Mol. Phys. 1999, 96, 1517–1524.

(82)

Nieto-Draghi, C.; Ungerer, P.; Rousseau, B. Optimization of the anisotropic united atoms intermolecular potential for n-alkanes: Improvement of transport properties. J. Chem. Phys. 2006, 125, 044517.

(83)

Liu, A.; Beck, T. L. Vapor-Liquid Equilibria of Binary and Ternary Mixtures Containing Methane, Ethane, and Carbon Dioxide from Gibbs Ensemble Simulations. J. Phys. Chem. B 102, 1998, 7627–7631.

(84)

Lisal, M.; Smith, W. R.; Nezbeda, I. Accurate Computer Simulation of Phase Equilibrium for Complex Fluid Mixtures. Application to Binaries Involving Isobutene, Methanol, 37 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Methyl tert -Butyl Ether, and n -Butane. J. Phys. Chem. B 1999, 103, 10496–10505 . (85)

Martin, M. G.; Chen, B.; Siepmann, J. I. Molecular Structure and Phase Diagram of the Binary Mixture of n -Heptane and Supercritical Ethane: A Gibbs Ensemble Monte Carlo Study. J. Phys. Chem. B 2000, 104, 2415–2423.

(86)

Pàmies, J. C.; McCabe, C.; Cummings, P. T.; Vega, L. F. Coexistence Densities of Methane and Propane by Canonical Molecular Dynamics and Gibbs Ensemble Monte Carlo Simulations. Mol. Simul. 2003, 29, 463–470.

(87)

Errington, J. R.; Shen, V. K. Direct evaluation of multicomponent phase equilibria using flat-histogram methods. J. Chem. Phys. 2005, 123, 164103.

(88)

Rane, K. S.; Murali, S.; Errington, J. R. Monte Carlo Simulation Methods for Computing Liquid − Vapor Saturation Properties of Model Systems. J. Chem. Theory Comput. 2013,

9, 2552–2566. (89)

Ganzenmüller, G.; Camp, P. J. Applications of Wang-Landau sampling to determine phase equilibria in complex fluids. J. Chem. Phys. 2007, 127, 154504.

(90)

Wang, F.; Landau, D. P. Determining the density of states for classical statistical models: A random walk algorithm to produce a flat histogram. Phys. Rev. E 2001, 64, 56101.

(91)

Wilding, N. B.; Bruce, A. D. Density fluctuations and field mixing in the critical fluid. J. Phys.: Condens. Matter 1992, 4, 3087–3108.

(92)

Ewing, M. B.; Johnson, K. A.; McGlashan, M. L. The (liquid + liquid) critical state of (cyxlohexane + methanol) IV. (T, x)p coexistence curve and the slope of the critical line. J. Chem. Thermodyn. 1988, 20, 49–62.

(93)

Pérez-Pellitero, J.; Ungerer, P.; Orkoulas, G.; Mackie, A. D. Critical point estimation of the Lennard-Jones pure fluid and binary mixtures. J. Chem. Phys. 2006, 125, 54515.

(94)

Singh, S. K.; Singh, J. K. A comparative study of critical temperature estimation of atomic fluid and chain molecules using fourth-order Binder cumulant and simplified scaling laws. Mol. Simul. 2013, 39, 154–159. 38 ACS Paragon Plus Environment

Page 38 of 56

Page 39 of 56

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(95)

Binder, K. Finite Size Scaling Analysis of Ising Model Block Distribution Functions. Z. Phys. B - Condens. Matter 1981, 43, 119–140 .

(96)

Bruce, A. D. Probability density functions for collective coordinates in Ising-like systems. J. Phys. C: Solid State Phys. 1981, 14, 3667–3688.

(97)

Wichterle, V.; Kobayashi, R. Vapor-Liquid Equilibrium of Methane-Ethane System at Low Temperatures and High Pressures. J. Chem. Eng. Data 1972, 17, 9–12.

(98)

Vrabec, J.; Fisher, J. Vapor-Liquid Equilibria of Binary Mixtures Containing Methane , Ethane , and Carbon Dioxide from Molecular Simulation. Int. J. Thermophys. 1996, 17, 889–908.

(99)

Muller, A.; Winkelmann, J.; Fisher, J. Backone Family of Equations of State : 1. Nonpolar and Polar Pure Fluids. AIChE J. 1996, 42, 1116–1126.

(100) Matschke, D. E.; Thodos, G. Vapor-Liquid Equilibria for the Ethane-Propane System. J. Chem. Eng. Data 1962, 7, 232–234. (101) Linstrom, P. J.; Mallard, W. G. in NIST Chemistry Webbook,NIST Standard Reference Database no. 69 Gaithersburg MD 20899 (102) Mon, K. K.; Binder, K. Finite size effects for the simulation of phase coexistence in the Gibbs ensemble near the critical point. J. Chem. Phys. 1992, 96, 6989–6995. (103) Luijten, E.; Fisher, M.; Panagiotopoulos, A. Universality Class of Criticality in the Restricted Primitive Model Electrolyte. Phys. Rev. Lett. 2002, 88, 185701.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Table 1: Interaction potential parameters for the non-bonded interactions in TraPPE-UA model for this work. Lorentz-Berthelot combining rules are used for the cross terms. Data taken from Martin and Siepmann.10

United Atom

ε/kB (K)

σ (nm)

CH4

148

0.373

-CH3

98

0.375

-CH2-

46

0.395


Page 40 of 56

Page 41 of 56

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 2:Interaction potential parameters for the bonded interactions of the TraPPE-UA model used in this work. Data taken from Martin and Siepmann.10

(a) Bond Length Parameters United Atoms

Bond length (nm)

CH3-CH3

0.154

-CH2-CH2-

0.154

CH2-CH3

0.154

(b) Bond Angle Parameters United Atoms

kθ/kB(K rad-2)

θ0 (◦)

CH3-CH2-CH3

62500

114

-CH2-CH2-CH3

62500

114

(c) Torsional Angle Parameters United Atoms

c0/kB (K)

c1/kB (K)

c2/kB (K)

c3/kB (K)

CH3-CH2-CH2-CH3

0.0

355.03

-68.19

791.32



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 42 of 56

Table 3:Summary of the different simulation box-lengths used in the current study System

Temperature (K)

Box Length (nm)

320

3.4

353

3.4

Propane - Butane

3.4 373.15

4.5 6.0

255.37

3.0

283.15

3.0

Ethane - Propane

3.0 329.82

4.2 5.5

172.04

2.85

186.09

2.85

Methane - Ethane

2.85 199.93

3.4 4.1


Page 43 of 56

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 4: Comparison of pure component vapor pressures, obtained via simulations using the TraPPE force field and using experiments, of the different species studied in the current work. Actual* data reported are from the NIST website.101 The numbers reported in parenthesis are the statistical uncertainties (standard deviations obtained by performing three independent sets of simulations) in the last digit(s).

Species

CH4

CH3-CH3

CH3-CH2-CH3

CH3-CH2-CH2-CH3

Temperature

Simulation pressure

Actual pressure

(K)

(bar)

(bar)

172.04

25.18(1)

25.04

0.55

186.09

39.54(9)

39.97

-1.09

172.04

0.77(1)

0.49

56.93

186.09

1.60(1)

1.10

45.67

199.93

3.01(1)

2.17

38.76

255.37

18.09(1)

15.12

19.61

283.15

34.49(1)

30.18

14.28

255.37

3.65(1)

2.65

37.87

283.15

8.30(13)

6.37

30.44

320.00

19.27(4)

15.92

21.00

329.82

23.40(5)

19.08

18.41

353.00

35.91(2)

31.24

14.95

320.00

5.97(6)

4.56

30.86

353.00

12.61(10)

10.08

25.08

373.15

18.55(3)

15.26

21.58

*Actual value refers to data estimated from equations fitted by NIST to experimental data.

.


% Deviation


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 44 of 56

Figure Captions Figure 1: Phase diagrams of the methane(1)-ethane(2) mixture at two temperatures below the critical point of the more volatile component methane (191.4 K) in (a) the P-x-y plane, (b) the P-ρ plane and (c) the x-y-ρ plane. The experimental results shown are that of Wichterle and Kobayashi.97 The lines are guides to the eye.

Figure 2: Phase diagrams of the methane(1)-ethane(2) mixture at 199.93 K above the critical point of methane plotted for three different simulation box lengths. The graphs are plotted in (a) the P-x-y plane, (b) the P-ρ plane and (c) the x-y-ρ plane. The experimental results are that of Wichterle and Kobayashi.97 The equation of state results are taken from Müller et al.99 The lines are guides to the eye.

Figure 3: Binder fourth order cumulant values of the methane(1)-ethane(2) mixture plotted as a function of (a) mean composition in mole fractions and (b) pressure for three different simulation box lengths at a temperature of 199.93 K which is above the critical point of methane. The lines are guides to the eye.

Figure 4: Phase diagrams of the ethane(1)-propane(2) mixture at two temperatures below the critical point of ethane (304 K) in (a) the P-x-y plane, (b) the P-ρ plane and (c) thex-y-ρ plane. The experimental results are obtained from Matschke and Thodos.100 The lines are guides to the eye. 44 ACS Paragon Plus Environment

Page 45 of 56

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 5: Phase diagrams for the ethane(1)-propane(2) mixture at a temperature of 329.82 K, which is above the critical temperature of ethane, for three different simulation box-lengths. The figures are plotted in (a) the P-x-y plane, (b) the P-ρ plane and (c) thex-y-ρ plane. The experimental results shown are taken from Matschke and Thodos.100 The lines are guides to the eye.

Figure 6: Binder fourth order cumulant values of the ethane(1)-propane(2) mixture plotted as a function of (a) mean composition in mole fractions and (b) pressure for three different simulation box lengths at a temperature of 329.82 K which is above the critical temperature of ethane. The lines are guides to the eye.

Figure 7: Phase diagrams of the propane(1)-butane(2) mix at two temperatures below the critical point of propane (368 K) in (a) the P-x-y plane, (b) the P-ρ plane and (c) the x-y-ρ plane. Experimental data have been taken from Kay.56 The lines are guides to the eye.

Figure 8: Phase diagrams for the propane(1)-butane(2) mixture at 373.15 K, which is above the critical temperature of propane, for three different simulation box lengths. The plots are in (a) the P-x-y plane, (b) the P-ρ plane and (c) the x-y-ρ plane. The experimental data have been taken from Kay.56 The lines drawn are guides to the eye.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 9: Binder fourth order cumulant values of the propane(1)-butane(2) mixture plotted as a function of (a) mean composition in mole fractions and (b) pressure for three different simulation box lengths at a temperature of 373.15 K which is above the critical temperature of propane. The lines drawn are guides to the eye.


Page 46 of 56

Page 47 of 56

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 1:



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 2:


Page 48 of 56

Page 49 of 56

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 3:



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 4:


Page 50 of 56

Page 51 of 56

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 5:



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 6:


Page 52 of 56

Page 53 of 56

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 7:



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 8:


Page 54 of 56

Page 55 of 56

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 9:



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

TOC Figure 164x170mm (96 x 96 DPI)

ACS Paragon Plus Environment

Page 56 of 56

Phase Equilibria and Critical Point Predictions of Mixtures of Molecular

Recommend Documents