Prediction of VaporâLiquid Coexistence Data for p-Cymene Using

Article pubs.acs.org/jced

Prediction of Vapor−Liquid Coexistence Data for p‑Cymene Using Equation of State Methods and Monte Carlo Simulations Madakashira Harini and Jhumpa Adhikari* Department of Chemical Engineering, Indian Institute of Technology Bombay, Mumbai 400076, India

K. Yamuna Rani Chemical Engineering Division, Indian Institute of Chemical Technology, Hyderabad 500607, India S Supporting Information *

ABSTRACT: A naturally occurring aromatic organic compound, p-cymene, finds applications as reaction intermediate, solvent in production of pharmaceuticals, fragrances, and fine chemicals. The experimental vapor liquid equilibria (VLE) data of p-cymene is available at pressures up to 537.4 kPa. In this study, the thermodynamic properties of p-cymene are determined using structure property correlations combined with equations of state (EoS) and molecular simulation techniques. Two molecular simulation techniques, Gibbs ensemble Monte Carlo (GEMC) and grand canonicaltransition matrix Monte Carlo (GC-TMMC) have been employed for prediction of VLE. The estimates of properties, including vapor pressures, heats of vaporization, coexistence densities, and critical properties have been compared with available experimental data. The thermodynamic properties predicted by molecular simulations and also that obtained from Peng− Robinson (PR) and volume translated Peng−Robinson (VTPR) EoS are generally, in broad agreement with the experimental data. Further, GEMC results for coexistence densities and saturation vapor pressures are compared with that from GC-TMMC technique. simulations10 are employed in this study. The properties estimated include the normal boiling point temperatures (Tb), the saturation temperatures, and the corresponding vapor pressures, the heats of vaporization (ΔHvap), and the coexistence densities as well as the critical state properties and the acentric factor (ω). The property prediction methods based on the structure of a molecule have found extensive applications, especially in computer-aided molecule design for specific property targets.11 In addition, these structure−property correlations can be utilized to determine the properties of molecules when experimental data are not available. The commonly used structure property correlations are the group contributionbased Joback and Reid,12 Constantinou-Gani,13 and MarreroGani (MG) methods;9 and also, the connectivity indices method.14 The group contribution method is based on the principle that “the structural aspects of individual chemical components are the same in different molecules”.15 This

1. INTRODUCTION The alkyl benzene, p-cymene, is a constituent of essential oils of thyme and oregano; and is synthetically produced by distilling turpentine oil. This compound is also obtained as a byproduct in the manufacture of sulphite paper pulp. It is used as a chemical intermediate in the manufacture of organic compounds such as p-cresol and thymol, apart from its use in the fragrance and flavor industry.1 In addition, this aromatic compound also finds application as a solvent, in carrying out reactions in a high boiling solvent,2 and in recycling polystyrene,3 where p-cymene itself is recycled after use by distillation. The knowledge of vapor liquid equilibria (VLE) data at high temperatures and pressures along with critical property information is essential, for example, to optimize the process design parameters so as to separate p-cymene efficiently. A comprehensive literature review has shown that the availability of VLE data for p-cymene is in the pressure range4,5 (0.04−537.4 kPa) including the normal boiling point; and also, the critical point.6,7 As performing experiments to determine vapor−liquid coexistence properties at high temperatures and pressures is expensive and time-consuming, theoretical methods such as equation of state methods8 combined with group contribution9 based parameters and Monte Carlo (MC) © 2014 American Chemical Society

Special Issue: Modeling and Simulation of Real Systems Received: January 29, 2014 Accepted: April 30, 2014 Published: May 8, 2014 2987

dx.doi.org/10.1021/je5001022 | J. Chem. Eng. Data 2014, 59, 2987−2994

Journal of Chemical & Engineering Data

Article

assumption makes the property prediction “inherently approximate”, although the range of its applicability is extensive. Among the group contribution-based methods, the MG method is widely used to estimate a wide range of thermodynamic and physical properties of organic compounds such as Tb, ΔHvap, critical temperature (Tc), critical pressure (Pc), critical volume (Vc), ω, viscosity, etc., with low average absolute error. We have also used the MG method along with EoS approach here to predict the properties of p-cymene. The development of various cubic equations of state (EoS) from van der Waals to volume translated Peng−Robinson (VTPR) has improved the prediction of equilibrium vapor pressure and coexistence liquid densities to a great extent. Also, the VTPR EoS predicted pure component and mixture equilibrium properties are found to be in good agreement for various molecule types such as alkanes, alcohols, and aromatics.8 These EoS require Tc, Pc, Vc, and ω as input parameters which can be predicted using either structure property correlations or molecular simulations. Alternate to the EoS approach, molecular simulations facilitate the determination of thermodynamic properties of molecules with reasonably good accuracy; provided, there exist accurate force field parameters. Siepmann and his co-workers developed a force field, TraPPE-UA (Transferable Potentials for Phase Equilibria-United Atom) to predict phase equilibrium properties, where the force field parameters for the respective UAs are transferable among different molecules. This force field has been widely used in literature to predict equilibrium properties of cyclic alkanes,16 esters,17 ketones, ethers, aldehydes,18 and aromatic compounds.19 The statistical uncertainty on critical properties of various organic compounds predicted using this force field is found to be less than 1 % for alkanes, ethers, ketones, and aldehydes;18 and within 3 % for alkyl benzenes.19 The Gibbs ensemble Monte Carlo (GEMC) technique developed by Panagiotopoulos is the most popular and widely used method to predict phase equilibria of pure compounds and mixtures.20 Conversely, as an alternative to GEMC, Errington has developed a highly efficient method to simulate VLE of pure components and mixtures; which is known as the grand canonical-transition matrix Monte Carlo (GC-TMMC) technique.21 Paluch et al. have found excellent agreement of VLE data generated using GC-TMMC with that from GEMC for pure organic molecular fluids such as ethane, n-octane, cyclohexane, 2,5-dimethyl hexane, and 1-propanol.22 The simulations performed by these authors employed the MCCCS Towhee package and TraPPE-UA force field. In general, the authors determined that the uncertainties in saturated vapor pressures and saturated vapor densities obtained by GC-TMMC were significantly lower than that from GEMC; while the uncertainties in saturated liquid densities were found to be comparable. We have adopted the same software package, including the GC-TMMC code built therein by Errington and co-workers, for our study on pcymene using both the MC simulation techniques to determine phase equilibrium properties. In view of the important applications of p-cymene, we have studied the VLE using two different approaches, EoS based on structure property prediction (MG method), and molecular simulations (GEMC and GC-TMMC with TraPPE-UA force field). The remainder of the article is organized as follows. First, we discuss the MG property estimation method and various EoS employed in this study. Details of the TraPPE-UA force field and a missing torsional parameter as predicted by density

functional theory (DFT) approach for p-cymene are presented in the following section, along with the simulation details of the GEMC and the GC-TMMC techniques. Further, in the subsequent section, the results obtained by employing the EoS approach are analyzed. The vapor−liquid coexistence curve, vapor pressures, and critical points for p-cymene as predicted by TraPPE-UA force field using GEMC-NVT and GC-TMMC techniques are presented. The results obtained from the MC simulations are also compared with the predictions of various EoS. The VLE data predicted by the theoretical methods and MC simulations are compared with data from literature (where available) for estimating the relative accuracy of these methods. The conclusions of this work are summarized in section 4.

2. METHODOLOGY The methodologies employed for the determination of the vapor liquid coexistence curve for p-cymene are described in this section. The structure of p-cymene with representation of various groups present in it, is shown in Scheme 1. Scheme 1. Structure of p-Cymene. 1,2,4,5 CH(aro); 3,6 C(aro); 7 CH; 8-10 CH3

2.1. Structure Property Correlations. The MG method considers the second and third order contributions in addition to first order contributions, due to which the error in the value of property estimated is significantly reduced; when compared to other methods existing in the literature, such as Joback and Reid,12 and Constantinou-Gani.13 In 2001, the first work on the MG method was published, in which a database of ∼1800 molecules were regressed for the Tb prediction.9 The database was then extended to ∼3500 molecules by Hukkerikar et al. in 2012.23 In addition, the 2012 article also reports revised and improved model parameters for prediction of pure component properties. Two sets of model parameters have been reported, one with simultaneous regression (SR) and the other with stepwise regression (SWR). Here, in this article we have employed SR parameters, as the property models analyzed using SR parameters are known to perform better than SWR parameters.23 In the p-cymene molecule, only first and second order groups are present; and the MG contributions of these groups are mentioned in Table 1. In the absence of third order groups, the property estimation model expression takes the following form, f (X ) =

∑ NC i i + w ∑ MjDj i

j

(1)

where Ci and Dj are the contributions of the first and second order groups of type i and j that occur Ni and Mj times, respectively. A binary parameter, w, indicates the presence or absence of groups of second order and is unity in our study. The function, f(X), on the left-hand side of eq 1 represents a 2988



Article

Table 1. MG-SR Contributions for Various Groups in p-Cymene23 groups

occurrence

first order

Ni

Tb1i/K

Tc1i/K

Pc1i/bar

Vc1i/(cc/kmol)

ωc1i

4 1 1 2 (Mj)

0.7332 0.0274 1.2616 0.9218 Tb2i/K

3.7648 12.144 7.9731 1.0898 Tc2i/K

0.0046 0.0183 0.0153 0.01 Pc2i/bar

45.2515 68.6439 98.8977 63.9854 Vc2i/(cc/kmol)

0.0013 0.002 0.0031 0.0017 ωc2i

1

0.1303

0.3642

−0.0043

3.7173

0.0007

CH(aro) C(aro)−CH C(aro)−CH3 CH3 second order C(aro)−CH(CH3)2

contribution

atom force field, where the hydrogen atoms that are connected to a carbon atom are treated as a single entity. The structure of p-cymene with united atom representation is shown in Scheme 1. The aromatic ring is kept rigid with a semiflexible side chain, where the bond lengths are again rigid, but with flexible bending and torsional degrees of freedom. The primary goals in the development of this force field as stated by Siepmann and co-workers are “accuracy and transferability”. For example, in nalkanes the authors have mentioned that the force field parameters for bonded interactions are obtained by fitting the experimental data or the electronic structure calculations of the small molecules (such as ethane); while for nonbonded interactions the parameters are found by fitting to the experimental vapor liquid coexistence data of the small molecules.25 The parameters of the molecule of interest here, p-cymene (as reported in Tables 2 and 3) have been obtained

simple function of the target property X, such as exponential of ratio of the target property to its adjustable parameter. The model expressions for Tb, Tc, Pc, Vc, and ω with the inclusion of adjustable parameters for MG method23 are given below (eq 2 to eq 6). Tb = 244.79 ln(∑ NT i b1i +

∑ MjTb2j)

i

(2)

j

Tc = 181.67 ln(∑ NT i c1i +

∑ MjTc2j)

i

(3)

j

(Pc − 0.0519)−0.5 − 0.1155 = (∑ NP i c1i + i

∑ MjPc2j) j

(4)

(Vc − 14.62) = (∑ NV i c1i + i

∑ MjVc2j) (5)

j

⎡ ⎛ ω ⎞⎤0.0447 ⎟⎥ − 1.0039 = (∑ Niω1i + ⎢exp⎜ ⎣ ⎝ 0.9132 ⎠⎦ i

Table 2. TraPPE-UA Parameters for Nonbonded LJ Interactions for p-Cymene

∑ M j ω 2j ) j

(6)

2.2. Equations of State. The ease of usage and extensive applicability of different EoS, made these methods more prevalent to predict VLE of various organic compounds. Peng− Robinson (PR) EoS has been widely used to predict VLE of pure components and mixtures. Our earlier work24 on the prediction of VLE of pharmaceutical intermediate, phenylacetylcarbinol, lists the advantages and drawbacks of different EoS along with the parameters employed. Hence, we do not repeat the detailed discussion on the EoS approach to determine VLE here. The input parameters required to determine VLE using EoS are Tc, Pc, Vc, and ω. In the classical EoS approach, these parameters are obtained from experimental data in literature.6,7 However, we have predicted the input parameters using the MG method which are necessarily less accurate than the experimental values. Thus, when MG method predicted input parameters are employed, the errors in the property estimates using these MG parameters will propagate to predictions by EoS in addition to the error due to inherent approximations in the EoS itself. In this study, the vapor pressures and coexistence liquid and vapor densities at various temperatures have been computed using Soave−Redlich−Kwong (SRK), PR, and VTPR EoS. Further, the predicted properties from EoS are compared with the MC simulation method predictions described in section 3. 2.3. Simulation Details. 2.3.1. Force Field. The TraPPEUA force field, developed by the Siepmann group at the University of Minnesota has been employed in this study to predict VLE of p-cymene. As the name states, it is a united

group

(ε/kB)/K

(Å)

ref

CH3 CH C(aro) CH(aro)

98.0 10.0 21.0 50.5

3.75 4.68 3.88 3.695

18 18 19 19

Table 3. TraPPE-UA Force Field Parameters for Bonded Interactions of p-Cymene vibration

bond length r0/Å

ref

CH(aro)−CH(aro) CH(aro)−C(aro) C(aro)−CH C(aro)−CH3 CH−CH3 bending

1.40

19

1.54

19

CH(aro)−CH(aro)−C(aro) CH(aro)−C(aro)−CH(aro) C(aro)−CH−CH3 CH(aro)−CH−CH3 torsion CH(aro)−C(aro)−CH(aro)− CH(aro) CH(aro)−CH(aro)−C(aro)− CH(aro) C(aro)−CH(aro)−CH(aro)− C(aro) CH(aro)−CH(aro)−C(aro)−CH CH(aro)−CH(aro)−C(aro)− CH3 CH(aro)−C(aro)−CH−CH3

2989

kθ/2kB/(K/rad−2)

θ0/deg

ref

rigid

120.0

19

31250.0

112.0

19

constants (K)

equation type

ref

rigid

9

19

C0/kB: 583.7 C1/kB: 775.0

10

this work



Article

basis set using the Gaussian 03 software.29 Scans are performed at 10° interval from −180° to 180° for CH(aro)−C(aro)− CH−CH3 dihedrals. The Fourier coefficients are determined by fitting eq 10, a TraPPE simple cosine series, to the DFT derived potential energy surfaces with the curve fitting toolbox (cftool) in Matlab. Figure 1 shows the torsional energy from potential

from molecules other than p-cymene by assuming transferability of parameters for the same UAs present in different molecules. In addition, this force field has been proven to predict the vapor liquid coexistence densities and other thermodynamic properties of various organic compounds, including aromatics, accurately.19 Only a brief description of this force field is provided in this article, as this force field has been widely discussed in detail elsewhere in the literature.16−19 In TraPPE-UA force field, the nonbonded van der Waals interactions are governed by Lennard-Jones (LJ) given as follows: ⎡⎛ ⎞12 ⎛ ⎞6 ⎤ σij σij u(rij) = 4εij⎢⎢⎜⎜ ⎟⎟ − ⎜⎜ ⎟⎟ ⎥⎥ r ⎝ rij ⎠ ⎦ ⎣⎝ ij ⎠

(7)

where rij, εij, and σij are the separation distance between interaction sites i and j, the LJ well depth, and the LJ diameter, respectively. The LJ parameters for p-cymene are summarized in Table 2. Lorentz−Berthelot combining rules are employed to determine parameters for LJ unlike interactions between sites. The LJ interactions are truncated at a cutoff distance of 14 Å, and analytical tail corrections are used. There are no partial charges present on any of the UAs of p-cymene. The bonded interactions in the aromatic ring are held rigid. However, for the side chain bonded interactions, the bond lengths are kept rigid and bond bending is modeled by a harmonic potential specified as

Figure 1. Torsional energy as a function of dihedral angle of p-cymene for CH(aro)−C(aro)−CH−CH3. Potential energy scans at B3PW91/ 6-31++G(d,p) are shown as open circles. The solid line is the fit to the cosine series represented by eq 10 to the DFT data.

energy scans as a function of dihedral angle. Figure 1 depicts that the fit of eq 10 is in good agreement with the DFT calculations. The discrepancies in the fit at points other than maxima and minima have been observed in the literature for DFT calculated torsional energies of hemiacetals and glyoxal derivatives.30 2.3.2. Molecular Simulation. All the simulations in this study are carried out using MCCCS Towhee package,31 version 7.4. Constant (total) volume Gibbs ensemble Monte Carlo (GEMC-NVT) and grand canonical transition matrix Monte Carlo (GC-TMMC) simulations have been performed over a wide temperature range. Both GEMC and GC-TMMC offer specific advantages, and therefore, it is hard to comment about the superiority of either method. Nevertheless, GC-TMMC has the advantage of performing the simulation in a single box, as compared to GEMC where two box simulations at each coexistence point are necessary. Additional details on the comparison of these two methods can be found in Paluch et al.22 a. Gibbs Ensemble Monte Carlo Simulation. The pure component VLE simulations have been carried out in the temperature range of 400 K to 600 K using GEMC-NVT, where a total of 300 molecules have been considered. There are two simulation boxes with simple cubic initial configurations: one with 260 molecules, representing the liquid phase and the other with 40 molecules, representing a vapor phase. The two boxes are in thermodynamic contact but do not share an interface. Periodic boundary conditions are implemented in all three directions of the simulation boxes. The criteria for VLE in GEMC-NVT simulations has been achieved by performing volume exchange moves between two boxes for pressure equilibration, configuration biased Monte Carlo (CBMC) interbox particle transfer moves for attaining chemical equilibrium. Intra-box CBMC particle transfer moves, aggregation volume biased moves, center of mass displacement, and rotational moves are performed to attain thermal equilibrium, and to sample the remaining degrees of freedom. The frequency ratio in which the above-mentioned moves are performed is 10:10:5:15:10:25:25 for volume exchange, inter-box swap CBMC, intra-box swap CBMC, aggregation volume biased,

kθ (θ − θeq)2 (8) 2 where kθ, θ, and θeq are the bending force constant, the bond angle, and the equilibrium bond angle, respectively. The bond length and bond angle parameters for p-cymene are reported in Table 3. The bond rotations are described by a torsional potential for united atoms separated by three bonds. The following torsional form (eq 9) has been employed for all rigid torsions of p-cymene; where x varies from 1 to “ntorloop”, the number of torsion loops, which determines the number of stable, equilibrium values for the torsion angle. The value of ntorloop is 2. ubend =

⎧∞ if |ϕ − C(x)| > C0 utorsion(ϕ) = ⎨ ⎩ 0 otherwise

(9)

where ϕ is the torsion angle, C(x) is the equilibrium value of torsion angle, and tolerance (C0) is taken to be 10−5 radians. The equilibrium values of torsion angles, C(0) and C(1) are 3.1459 and zero radians, respectively. For the torsion involving the UAs, CH(aro)−C(aro)−CH− CH3, the expression given below is employed where Co and C1 are the constants and ϕ is the torsional angle. utorsion(ϕ) = Co(1 − cos(2[ϕ − π ])) + C1 (10) The torsion parameters (Co and C1 in eq 10) for the group mentioned above are found missing in TraPPE-UA force field literature, and its parameters are essential to model the conformational behavior of p-cymene accurately. The parameters for missing torsional potential have been determined using the DFT approach. This methodology has been successfully implemented earlier in literature to determine missing torsions in various other organic compounds including aromatics.16,26,27 The potential energy surface scan has been carried out at the B3PW91 level of theory28 with 6-31++G(d,p) 2990



Article

The plots of estimated vapor pressure (ln Psat) from various EoS as a function of temperature (1/T) are shown in Figure 2;

CBMC regrowth, centre of mass translational and rotational moves, respectively. The system is allowed to equilibrate over 200 000 MC cycles with averages computed at the end of 50 000 cycles of the production run. The statistical uncertainties (standard deviations) are computed during the production runs by taking a block of 1000 MC cycles. b. Transition Matrix Monte Carlo Simulation. GC-TMMC simulation requires the knowledge of chemical potential (μ), system volume (V), temperature (T), and the maximum number of molecules (Nmax) allowed in the system as input values. A simulation box of volume 64 000 Å3 and Nmax of 1000 has been used to determine VLE of p-cymene in a temperature range of 350 K to 575 K. CBMC of the type couple to prenonbond formulation has been employed along with aggregation volume biased MC moves to enhance the sampling further. The fraction of moves attempted to attain equilibrium are 40 % particle insertions/deletions, 15 % aggregation volume biased, 10 % configuration biased regrowths, 15 % center of mass translation, and 20 % rotation. The statistical uncertainties (standard deviations) are recorded by performing three independent sets of simulations. The vapor and liquid compositions and pressure at coexistence has been determined from the probability distributions obtained at the end of the simulation run.22 The VLE data predicted using GC-TMMC have been compared to that from GEMC-NVT simulations.

Figure 2. Comparison between vapor pressures obtained from EoS approach, GEMC-NVT simulations, GC-TMMC simulations from this work, and literature reported experimental4,5 data. The error bars for simulation data are smaller than the size of the markers.

where, the input parameters for EoS are predicted using the MG method. The vapor pressure predictions from VTPR and PR EoS coincide and are also in good agreement with experimental4,5 data. The SRK EoS consistently overestimates the vapor pressures as compared to the values from experiment and also the predictions by PR and VTPR EoS. The deviations in the prediction of vapor pressure using SRK EoS from that of PR and VTPR EoS decrease with increase in temperature. This same trend is also observed for SRK EoS predicted deviations when compared with the experimental data which is available until 523.15 K. Furthermore, the ΔHvap data of p-cymene determined using the Claperyon equation for various EoS at different temperatures is shown in Figure 3. The literature

3. RESULTS AND DISCUSSION The predicted VLE data for p-cymene using the EoS approach and the molecular simulation techniques are discussed in this section. The saturated vapor pressures, enthalpies of vaporization, and coexistence vapor and liquid densities obtained from these methods are compared. In addition, the equilibrium vapor pressures and enthalpies of vaporization obtained from the theoretical and simulation methods are also compared with the existing experimental data.4−7 3.1. EoS Approach with Structure Property Correlations. To begin with, the Tb, Tc, and Pc have been predicted using the MG method to facilitate the comparison between property values estimated using structure property correlations and the experimental data available in the literature.4−7 As shown in Table 4, the property values predicted from the MG Table 4. Property Data As Predicted by GEMC-NVT and GC-TMMC Simulations and MG Method for p-Cymene property

GEMC-NVT

GC-TMMC

MG

experiment4−7

Tb (K) Tc (K) Pc (bar) Vc (cc/mol) ω

447(3) 652(4) 27(2) 469(4) 0.2519

449(2) 653(2) 29(2) 508(3) 0.2739

446(6) 659(8) 30(1) 495(8) 0.37(5)

450(1) 652 28 505 0.3703

Figure 3. Comparison of enthalpy of vaporization for p-cymene obtained from various EoS, GEMC-NVT simulation, experimental32 data, and the Riedel equation estimation at Tb. The error bars for simulation data are smaller than the size of the markers.

reported experimental data32 of ΔHvap in the temperature range of 273.15 K to 450 K is also compared with EoS predictions in Figure 3. The VTPR predicted ΔHvap is in good agreement with the experimental data32 at temperatures below 450 K. The PR EoS predictions slightly underestimate VTPR predictions (and also the experimental results) below 450 K. However, the estimates from PR and VTPR EoS are in agreement with increase in temperature; while the SRK EoS results consistently underestimate both PR and VTPR predictions over the entire temperature range. In addition, the vapor and liquid coexistence densities have been determined using various EoS by utilizing the critical property data and ω value (as and when required), estimated by means of the MG method. The vapor phase

method are in accordance with the literature reported experimental data.4−7 The Tc, Pc, Vc, and Tb are found to be ∼7 K higher, ∼2 bar higher, 10 cc/mol, and ∼4 K lower, respectively, than the corresponding experimental values. Further, the MG method has also been utilized to determine the ω value, which may be required as input parameter to each EoS along with Tc, Pc, and Vc, and are reported in Table 4. It is to be noted that the uncertainties in Tc, Pc, Vc, and ω values determined by the MG method as reported in Table 4 will propagate into the calculations of all the reported properties using this EoS approach. 2991



Article

6.3 % at 600 K), though comparable to that from PR EoS with a maximum deviation of 11 % at 350 K (and a minimum deviation of 6.1 % at 600 K); followed by SRK EoS with a maximum deviation of 8.9 % also at 350 K (and a minimum of 5 % at 600 K). On the other hand, the molecular simulation approach discussed below does not require Tc, Pc, Vc, and ω as input parameters. Hence, the accuracy of the results from GEMCNVT and GC-TMMC simulations are not influenced by the uncertainties in the estimates of Tc, Pc, Vc, and ω from structure−property correlations (or in the experimental measurements). 3.2. Molecular Simulation. Isobaric isothermal single-box MC simulation has been carried out at a temperature of 300 K and a pressure of 1 atm, to compute the density of p-cymene. The density obtained from simulation (0.854 g/cm3) is found to be in good agreement with the experimental value (0.857 g/ cm3) reported in literature.33 Therefore, we believe the TraPPE-UA force field is sufficiently accurate and is used further to determine the VLE of p-cymene. Both GEMC-NVT and GC-TMMC simulations have been performed to estimate the coexistence properties, which include saturation pressures, coexistence liquid, and vapor densities and Tb. Both GEMCNVT and GC-TMMC simulation generated results have been employed to calculate the critical state properties (Tc, Pc, Vc) and ω, the values of these properties are reported in Table 4. Only GEMC-NVT simulations have been employed to measure average ΔHvap values, and a comparison of these values with various EoS predictions and the available experimental data32 is shown in Figure 3. The comparison of predicted vapor pressures from various EoS, GEMC-NVT, and GC-TMMC techniques along with the available experimental data4,5 in the literature are shown in Figure 2. As is well established, VTPR EoS provides better vapor pressure estimates than SRK EoS8, and this is found to be true in this case also on comparison with experimental results. The vapor pressure predictions by GEMC-NVT are in agreement with the predictions of SRK EoS at temperatures above 400 K and overestimate the literature reported experimental data. The GC-TMMC predicted vapor pressures concur with the predictions of VTPR and PR EoS at temperatures above 375 K, and the predictions are also in good agreement with the available experimental4,5 data. The ΔHvap data as a function of temperature have been estimated using GEMC-NVT simulations. As illustrated in Figure 3, the ΔHvap values from all the methods employed here decrease with increase in temperature. The ΔHvap value predicted by the Riedel equation at Tb using MG parameters is in excellent agreement with experimental data while showing a deviation of 4.8 % from GEMC-NVT predictions. The VTPR and PR EoS predicted ΔHvap values are consistently overestimated, whereas SRK EoS predictions always underestimate GEMC predictions; though the deviations between various EoS and GEMC predictions decrease with increase in temperature. The available experimental data32 (below 450 K) almost coincides with VTPR predictions (except at Tb), while the GEMC predictions are significantly lower than the experimental32 results. The correctness of property values obtained by molecular simulation techniques is limited by the accuracy of the TraPPEUA force field. Nevertheless, as it is difficult to perform experiments at high temperatures and high pressures, these simulation techniques allow us to estimate the VLE of pcymene without requiring Tc, Pc, Vc, and ω as inputs, unlike the

coexistence densities as predicted by PR and VTPR EoS coincide, and SRK EoS slightly overestimates the values as compared to predictions from both PR and VTPR EoS (Figure 4). On the contrary, SRK EoS predicted coexistence liquid

Figure 4. Comparison of predicted densities of p-cymene from various EoS (using inputs from MG method), GEMC-NVT simulation, and GC-TMMC simulation. The error bars for simulation data are smaller than the size of the markers.

densities are significantly lower than that from PR and VTPR EoS predictions (Figure 4). The agreement between PR and VTPR EoS predicted coexistence liquid densities improve with increase in temperature. We note here there is no experimental data available in literature for comparison of coexistence densities. Further, the experimental values6,7 of Tc, Pc, Vc, and ω are set as input to the various EoS in the classical approach and the deviations of these results from the property predictions of the EoS with MG parameters as well as experimental data are discussed here. The discrepancies in the natural logarithm of vapor pressure values predicted from the classical EoS and EoS with MG parameters approach decrease with increase in temperature and are within a maximum deviation of 9.5 %, 8.7 % and 4.4 % for VTPR, PR, and SRK EoS, respectively. The saturated pressure (ln Psat) vs temperature (1/T) calculated from the classical EoS show the same trends as the EoS with MG parameters. The ln Psat data calculated using the classical approach show greater discrepancy from experimental data4,5 for VTPR and PR EoS unlike the EoS with MG parameters approach, while in the classical SRK EoS calculations the discrepancy from experiment is further magnified as compared to the deviation of SRK EoS with MG parameters from experimental data.4,5 The ΔHvap values as determined using the classical approach also differ from that calculated using EoS with MG parameters (reported in Figure 3). The variations in the ΔHvap values obtained from two methods show maximum deviations of 7.1 %, 22.2 %, and 22.8 % for SRK, PR, and VTPR EoS, respectively. As discussed previously, the VTPR EoS with MG parameters predicts ΔHvap values which are in good agreement with experimental data.32 However, the classical VTPR EoS approach shows significant deviations from the experimental data.32 Besides, we also find that the differences in the classical approach results from that of the EoS with MG parameters in the case of liquid coexistence densities increase with increase in temperature for all EoS considered here, and a maximum deviation of 9 % is noted. On the contrary, in the case of vapor coexistence densities, deviations in the values calculated using the two EoS approaches decrease with increase in temperature for all EoS considered. The VTPR EoS has the highest deviation of 11.4 % at 350 K (as against minimum deviation of 2992



Article

data6,7 by 3.6 %, 3.6 %, and 7.1 %, respectively. When compared with MG predictions and GEMC-NVT simulations, the estimated Vc from GC-TMMC is in good agreement with experimental6,7 data showing a deviation of 0.6 % only. The estimated Vc values from GEMC-NVT and MG prediction are lower than the experimental6,7 data with deviations of 7.1 % and 2.0 %, respectively. It is to be noted that the ω value as predicted by the MG method is in good agreement with the experimental data,5 but is ∼32 % and ∼26 % higher than GEMC-NVT and GC-TMMC predicted values, respectively. It may be noted here that the deviations are measured as the percentage of the absolute value of the difference in the GEMC-NVT or GC-TMMC or MG prediction, and the corresponding experimental data relative to the experimentally determined values and only the maximum deviations are reported here.

EoS approach. The VLE data generated using molecular simulation techniques can also be further used in process design, for example, to separate p-cymene efficiently by distillation. The GEMC-NVT and GC-TMMC determined Tb is around 447 K and 449 K, respectively, where both values are close to the experimental6,7 value (450 K) and MG method predictions (446 K) as reported in Table 4. The TraPPE-UA predicted Tb underestimates experimental value by a maximum of 0.67 %. A maximum deviation of 3 % has been observed in literature for aromatic hydrocarbons19 when the TraPPE-UA force field has been employed; therefore, the Tb for p-cymene obtained in this study is a close match with that from experiment. Figure 4 illustrates the coexistence liquid and vapor densities for p-cymene using GEMC-NVT, GC-TMMC, and various EoS. The vapor phase coexistence densities as predicted by PR and VTPR EoS are in good agreement with the GEMC-NVT and GC-TMMC predictions at low temperatures (below 450 K). However, with increase in temperature (above 450 K), the GEMC-NVT predicted vapor densities are less than that predicted by EoS and GC-TMMC methods. It is observed that GC-TMMC predictions coincide with SRK EoS predictions. For the liquid phase, the GEMC-NVT and GC-TMMC predicted liquid densities agree in the temperature range considered. The SRK and VTPR EoS predictions are consistently lower than the GEMC-NVT and GC-TMMC predicted liquid densities, while the PR EoS predictions are always higher until 475 K. In the high temperature range (above 475 K) the predictions from the two simulation techniques and PR EoS are in agreement. Again, as stated previously, the VTPR EoS is believed to improve the prediction of liquid densities as compared to the PR EoS.8 However, it should be noted that the VTPR predictions employ the generalized Twu parameters,34 as parameters specific to pcymene are not available in literature. We note here there is no experimental data available in literature for comparison of coexistence densities. The GEMC-NVT and GC-TMMC generated data could not be obtained beyond temperatures approaching 0.92Tc (∼600 K) and 0.88Tc (∼575 K) as molecular simulations cannot be performed at temperatures very close to Tc.10 It is observed that at temperatures above 600 K in GEMC-NVT, the fluctuations in property values increase, and densities in both the boxes become identical and the boxes switch identities. For GC-TMMC simulations above 575 K, we are unable to distinguish two distinct maxima in the particle number probability distribution plot. For GEMC-NVT and GC-TMMC simulations, the Tc and critical density (ρc) have been estimated by fitting the VLE data using the saturated density scaling law and the law of rectilinear diameters10 with a scaling exponent, β = 0.325. The Tc obtained from both methods are in close agreement with experimental6,7 data and are reported in Table 4. The ρc values are found to be 0.2857 g/cm3 and 0.2641 g/cm3 from GEMC-NVT and GCTMMC simulations, respectively. The Pc has been estimated by extrapolating to Tc, the fit of the simulation predicted vapor pressure data to Antoine equation. The acentric factor, ω, is calculated from vapor pressure data corresponding to the reduced temperature, Tr = (T/TC) of 0.7. The estimated values of Pc, Vc, and ω using the MG method, GEMC-NVT, and GC-TMMC techniques along with the experimental data obtained from literature6,7 are also reported in Table 4. The Pc values calculated from GEMC-NVT, GCTMMC data, and MG predictions differ from experimental

4. CONCLUSIONS The thermodynamic properties of p-cymene as determined by MG method combined with various EoS, GEMC-NVT, and GC-TMMC techniques, in the absence of experimental data at high pressures (above 537.4 kPa) are reported here. The VLE is determined by using EoS (SRK, PR, and VTPR) by making use of critical property data and acentric factor value as predicted by the MG method. Complementarily, GEMC-NVT and GCTMMC simulations are performed to determine VLE using the TraPPE-UA force field, in which an unknown torsion has been determined from the DFT approach. The predicted vapor pressures from GC-TMMC simulations and PR and VTPR EoS are in good agreement with the reported experimental data (below 537.4 kPa). The GEMC-NVT predicted vapor pressures match well with the predictions of VTPR and PR EoS with increase in temperature. The VTPR EoS predicted ΔH vap data is in good agreement with the reported experimental32 results, while the GEMC-NVT predicted ΔHvap data lie between the predictions of VTPR and SRK EoS. The estimated coexistence liquid densities from GEMCNVT and GC-TMMC coincide and lie between the predictions of PR and VTPR EoS at low temperatures (until 475 K) and are in better agreement with PR EoS predictions as temperature increases. However, the predicted coexistence vapor densities from GEMC-NVT and GC-TMMC are not in agreement. The GC-TMMC predicted vapor densities match well with the predictions of SRK EoS. The GEMC-NVT results underestimate the vapor densities when compared to GC-TMMC, VTPR, PR, and SRK EoS at high temperatures, and the vapor density values are close to the PR EoS predictions at temperatures below 475 K. The estimated critical properties from GC-TMMC and GEMC-NVT simulations are in good agreement with experimental6,7 data with the exception of critical volume where only GC-TMMC predictions are close.

■

ASSOCIATED CONTENT

S Supporting Information *

Tables listing the deviations in property (pressure, enthalpy of vaporization, coexistence densities) values estimated from the classical EoS approach with the MG method predicted properties as input to EoS and the uncertainty in the MG predicted properties as input to EoS. This material is available free of charge via the Internet at http://pubs.acs.org/. 2993



■

Article

(19) 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. (20) Panagiotopoulos, A. Z. Direct Determination of Phase Coexistence Properties of Fluids by Monte-Carlo Simulation in a New Ensemble. Mol. Phys. 1987, 61, 813−826. (21) Errington, J. R. Direct Calculation of Liquid−Vapor Phase Equilibria from Transition Matrix Monte Carlo Simulation. J. Chem. Phys. 2003, 118, 9915−9925. (22) 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. (23) Hukkerikar, A. S.; Sarup, B.; Ten Kate, A.; Abildskov, J.; Sin, G.; Gani, R. Group-Contribution+ (GC+) Based Estimation of Properties of Pure Components: Improved Property Estimation and Uncertainty Analysis. Fluid Phase Equilib. 2012, 321, 25−43. (24) Harini, M.; Adhikari, J.; Rani, K. Y. Prediction of Vapour−liquid Coexistence Data of Phenylacetylcarbinol. Fluid Phase Equilib. 2014, 364, 6−14. (25) Eggimann, B. L.; Sunnarborg, A. J.; Stern, H. D.; Bliss, A. P.; Siepmann, J. I. An Online Parameter and Property Database for the TraPPE Force Field. Mol. Simul. 2014, 40, 101−105. (26) Hao, M.; Haq, O.; Muegge, I. Torsion Angle Preference and Energetics of Small-Molecule Ligands Bound to Proteins. J. Chem. Inf. Model. 2007, 47, 2242−2252. (27) Maerzke, K. A.; Schultz, N. E.; Ross, R. B.; Siepmann, J. I. TraPPE-UA Force Field for Acrylates and Monte Carlo Simulations for Their Mixtures. J. Phys. Chem. B 2009, 113, 6415−6425. (28) Perdew, J. P.; Chevary, J. A.; Vosko, S. H.; Jackson, K. A.; Pederson, M. R.; Singh, D. J.; Fiolhais, C. Atoms, Molecules, Solids, And Surfaces: Applications of the Generalized Gradient Approximation for Exchange and Correlation. Phys. Rev. B 1992, 46, 6671− 6687. (29) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, Jr., J. A.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03, revision C.01; Gaussian Inc.: Wallingford, CT, 2003. (30) Kahn, K.; Bruice, T. C. Parameterization of OPLSAA Force Field for the Conformational Analysis of Macrocyclic Polyketides. J. Comput. Chem. 2002, 23, 977−996. (31) MCCCS Towhee. Available at: http://towhee.sourceforge.net. (Accessed March 2014). (32) Timmermans, J. Physico-Chemical Constants of Pure Organic Substances; Elsevier, New York, 1965; Vol II. (33) ChemSpider. http://www.chemspider.com/Chemical-Structure. 7183.html. (Accessed Jan 2014). (34) Twu, C. H.; Bluck, D.; Cunningham, J. R.; Coon, J. E. A Cubic Equation of State with a New Alpha Function and a New Mixing Rule. Fluid Phase Equilib. 1991, 69, 33−50.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Fax: +91 (22) 2572 6895. Tel.: +91 (22) 2576 7245. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The author M. Harini is grateful to Council of Scientific and Industrial Research (CSIR), New Delhi for financial support. We, also, thank Computer Centre at Indian Institute of Technology, Bombay and Center for Development of Advance Computing (CDAC), Pune head quarters for providing the computational resources.

■

REFERENCES

(1) Kirk-Othmers Encyclopdeia of Chemical Technology, 4th ed.; John Wiley and Sons: NewYork, 1996. (2) Mahato, S. B.; Banerjee, S. K.; Chakravarti, R. N. Reactions in High Boiling Solvents. II. Effect of Raney Nickel on Triterpenoids. Tetrahedron 1971, 27, 177−186. (3) García, M. T.; Gracia, I.; Duque, G.; Lucas, A.; De Rodríguez, J. F. Study of the Solubility and Stability of Polystyrene Wastes in a Dissolution Recycling Process. Waste Manag. 2009, 29, 1814−1818. (4) Mcdonald, R. A.; Shrader, S. A.; Stull, D. R. Vapor Pressures and Freezing Points of 30 Organics. J. Chem. Eng. Data 1959, 4, 311−313. (5) Kobe, K. A.; Okabe, T. S.; Ramstad, M. T.; Huemmer, P. M. pCymene Studies. VI. Vapor Pressure of p-Cymene, Some of Its Derivatives and Related Compounds. J. Am. Chem. Soc. 1941, 63, 264− 265. (6) Tsonopoulos, C.; Ambrose, D. Vapor−liquid Critical Properties of Elements and Compounds. 3. Aromatic Hydrocarbons. J. Chem. Eng. Data 1995, 40, 547−558. (7) Riddick, J. A., Bunger, W. B., Sakano, T. K. Organic Solvents: Physical Properties and Methods of Purification, 4th ed.; Wiley Interscience: New York, 1986. (8) Schmid, B.; Gmehling, J. From van Der Waals to VTPR: The Systematic Improvement of the Van Der Waals Equation of State. J. Supercrit. Fluids 2010, 55, 438−447. (9) Marrero, J.; Gani, R. Group-Contribution Based Estimation of Pure Component Properties. Fluid Phase Equilib. 2001, 183−184, 183−208. (10) Frenkel, D.; Smit, B. Understanding Molecular Simulation: From Algorithm to Applications; Academic Press: San Diego, 1996. (11) Harper, P. M.; Gani, R. A Multistep and Multilevel Approach for Computer Aided Molecular Design. Comput. Chem. Eng. 2000, 24, 677−683. (12) Joback, K. G.; Reid, R. C. Estimation of Pure-Component Properties from Group-Contributions. Chem. Eng. Commun. 1987, 57, 233−243. (13) Constantinou, L.; Gani, R. New Group Contribution Method for Estimating Properties of Pure Compounds. AIChE J. 1994, 40, 1697−1710. (14) Randic, M. The Connectivity Index 25 Years After. J. Mol. Graph. Model. 2001, 20, 19−35. (15) Sandler, S. I. Chemical and Engineering Thermodynamics; John Wiley and Sons: New York, 1999. (16) 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. (17) Ferrando, N.; Lachet, V.; Boutin, A. Transferable Force Field for Carboxylate Esters: Application to Fatty Acid Methylic Ester Phase Equilibria Prediction. J. Phys. Chem. B 2012, 116, 3239−3248. (18) 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. 2994


Prediction of VaporâLiquid Coexistence Data for p-Cymene Using

Recommend Documents

Prediction of VaporâLiquid Coexistence Data for p-Cymene Using