Vapor–liquid equilibria for the binary systems ethylene + water, ethylene + ethanol, and ethanol + water, and the ternary system ethylene + water + ethanol from Gibbs-ensemble molecular simulation

Accepted Manuscript Title: Vapor-liquid equilibria for the binary systems ethylene + water, ethylene + ethanol, and ethanol + water, and the ternary system ethylene + water + ethanol from Gibbs-ensemble molecular simulation Author: Y. Mauricio Mu˜noz-Mu˜noz Mario Llano-Restrepo PII: DOI: Reference:

S0378-3812(15)00109-0 http://dx.doi.org/doi:10.1016/j.fluid.2015.03.007 FLUID 10477

To appear in:

Fluid Phase Equilibria

Received date: Revised date: Accepted date:

15-1-2015 26-2-2015 2-3-2015

Please cite this article as: Y.Mauricio Mu˜noz-Mu˜noz, Mario Llano-Restrepo, Vapor-liquid equilibria for the binary systems ethylene + water, ethylene + ethanol, and ethanol + water, and the ternary system ethylene + water + ethanol from Gibbs-ensemble molecular simulation, Fluid Phase Equilibria http://dx.doi.org/10.1016/j.fluid.2015.03.007 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Vapor-liquid equilibria for the binary systems ethylene + water, ethylene + ethanol, and ethanol + water, and the ternary system ethylene + water + ethanol from Gibbs-ensemble molecular simulation

Y. Mauricio Muñoz-Muñoz, Mario Llano-Restrepo*[email protected]

School of Chemical Engineering, Universidad del Valle, Ciudad Universitaria Melendez, Building 336, Apartado 25360, Cali, Colombia

*

Corresponding author. Tel: + 57-2-3312935; fax: + 57-2-3392335

2

Highlights  We corroborated the validity of a set of Lennard-Jones plus point-charge intermolecular potential models for ethylene, water, and ethanol.  These models are capable of predicting the experimental VLE phase diagrams of the binary systems ethylene + water, ethylene + ethanol, and ethanol + water.  Molecular simulation predictions for the VLE of the ternary system ethylene + water + ethanol are in very good agreement with predictions made by means of the PRSV2-WSUNIQUAC thermodynamic model.

3

Abstract The conceptual design of a reactive separation process for the hydration of ethylene to ethanol requires reliable vapor-liquid equilibrium (VLE) data for the ternary system ethylene + water + ethanol. Due to the paucity of experimental data points in the VLE phase diagrams that have been reported for that system, molecular simulation looks appealing in order to predict such data. In this work, the Gibbs-ensemble Monte Carlo (GEMC) method was used to calculate the VLE of the pure components (ethylene, water, and ethanol), the binary subsystems (ethylene + water, ethylene + ethanol, and ethanol + water), and the ternary system (ethylene + water + ethanol). A set of previously validated Lennard-Jones plus point-charge potential models were chosen for the pure components, and the validity of these models was corroborated from the good agreement of the GEMC simulation results for the vapor pressure and the VLE phase diagrams of those components with respect to calculations carried out by means of the most accurate (reference) multiparameter equations of state currently available for ethylene, water, and ethanol. These potential models were found to be capable of predicting the available VLE phase diagrams of the binary subsystems: ethylene + water at 200 and 250°C, ethylene + ethanol at 150, 170, 190, 200, and 220°C, and ethanol + water at 200, 250, 275, and 300°C. Molecular simulation predictions for the VLE phase diagrams of the ternary system at 200°C and pressures of 30, 40, 50, 60, 80, and 100 atm, were found to be in very good agreement with predictions previously made by use of a thermodynamic model that combines the Peng-Robinson-Stryjek-Vera equation of state, the Wong-Sandler mixing rules, and the UNIQUAC activity coefficient model. The agreement between the predictions of these two independent approaches gives confidence for the subsequent use of molecular simulation to predict the combined phase and chemical equilibrium of the ternary system and check the validity of predictions previously made by means of the thermodynamic model.

Keywords: Vapor-liquid equilibrium, phase diagrams, ternary systems, ethylene hydration, petrochemical ethanol.

4

Nomenclature AAD

average absolute deviation

ARD

average relative deviation

kB

Boltzmann constant

k ij

interaction parameter in the Wong-Sandler mixing rule

N

number of molecules

P

pressure

P sat

vapor pressure

q

positive electric charge

q

negative electric charge

T

absolute temperature

U

configurational energy

Greek symbols



correction factor for the Lorentz combining rule

H v

heat of vaporization

u ij

characteristic energy parameter of the UNIQUAC model

U T

change of total configurational energy



energy parameter of the LJ potential model



molar density



size parameter of the LJ potential model

Subscripts Coul

Coulombic

LJ

Lennard-Jones

5

T

total

Superscripts real

real-space contribution

recip

reciprocal-space contribution

self

self-interaction contribution

6

1. Introduction The economics of the production of petrochemical ethanol by the direct hydration of ethylene might improve by the application of the concept of process intensification, i.e., by making the reaction and the separation of the product (ethanol) and the reactants (ethylene and water) to occur simultaneously in the same piece of equipment. In the intensified process, ethanol would be produced from ethylene by means of a vapor-liquid mixed-phase hydration, by feeding gaseous ethylene and liquid water into a reactive separation column [1]. In a previous study [1], we used a thermodynamic model that combines the Peng-RobinsonStryjek-Vera (PRSV2) equation of state [2, 3], the Wong-Sandler (WS) mixing rules [4], and the UNIQUAC activity coefficient model [5] in order to correlate the available vapor-liquid equilibrium (VLE) experimental data for the binary subsystems (ethylene + water, ethylene + ethanol, and ethanol + water) at 200°C. From the optimum values obtained for the adjustable parameters of the PRSV2-WS-UNIQUAC thermodynamic model, both the VLE and the combined phase and chemical equilibrium (CPE) of the ternary system (ethylene + water + ethanol) were predicted at 200°C and various pressures. The thermodynamic model predicts that for many values of the ethylene to water feed mole ratio, the vapor-liquid mixed-phase hydration of ethylene achieves equilibrium conversions much higher than those calculated for a vaporphase reaction that would hypothetically occur at the same conditions of temperature, pressure, and ethylene to water feed mole ratio. Another prediction of the thermodynamic model is that the reactive phase diagram of the ternary system exhibits a critical point at 200°C and 155 atm. Due to the paucity of experimental data points in the VLE phase diagrams reported by Tsiklis et al. [6] for the ternary system (in the absence of chemical reaction), it is difficult to make a direct assessment of the validity of the predictions made by the thermodynamic model for both the VLE and the CPE of the ternary system. However, an indirect assessment of those predictions is

7 possible by use of another method of predictive nature and entirely independent from the thermodynamic model, such as molecular simulation. The Gibbs-ensemble Monte Carlo (GEMC) method of molecular simulation [7-14] has been used widely and successfully by many authors for the computation of VLE phase diagrams of binary and ternary systems of industrial interest. In this work, the GEMC method was carefully implemented in order to compute the VLE of the pure components (ethylene, water, and ethanol), the binary subsystems (ethylene + water, ethylene + ethanol, and ethanol + water), and the ternary system (ethylene + water + ethanol). A set of previously validated Lennard-Jones plus pointcharge potential models [15-17] were chosen for ethylene, water, and ethanol. The outline of the paper is as follows. The intermolecular potential models for the pure components are presented in Section 2. The simulation methods are described in Section 3. In Section 4, simulation results for the pure components are reported, discussed, and compared with calculations carried out by means of the most accurate (reference) equations of state currently available for ethylene, water, and ethanol [18-20]. In Section 5, simulation results for the VLE phase diagrams of the binary subsystems are reported, discussed, and compared with both the available experimental data [2123] and the correlation results obtained by use of the PRSV2-WS-UNIQUAC thermodynamic model described in our previous study [1]. In Section 6, simulation results for the VLE phase diagrams of the ternary system are reported, discussed, and compared with the few experimental data available [6] and, more importantly, with the predictions previously made [1] by means of the thermodynamic model.

2. Intermolecular potential models As explained in more detail below, the Lennard-Jones (LJ) plus point-charge intermolecular potential models recently devised by Weitz and Potoff [15] for ethylene, Huang et al. [16] for water, and Schnabel et al. [17] for ethanol were chosen for the present work due to their greater accuracy relative to other models available for those components. These potential models are planar and rigid (i.e., they have no internal degrees of freedom).

8 The potential model for ethylene devised by Weitz and Potoff [15] comprises two LJ sites for the methylene (CH2) united-atom groups, each site with a positive point charge q  0.85 e , where e is the electron charge, and separated by a distance of 1.33 Å. A negative point charge

q  1.70 e is located in the middle point (center of mass, COM) between the two LJ sites. The point charges and the distance between the CH2 groups reproduce the gas-phase quadrupole moment of ethylene. The local coordinates and charges of the sites, and the size and energy parameters of the LJ potential are given in Table 1. With this model, Weitz and Potoff used the histogram-reweighting technique [14, 24] in combination with grand canonical ensemble Monte Carlo simulations to compute the VLE of pure ethylene and the binary mixtures ethylene + carbon dioxide, ethylene + xenon, and ethylene + n-butane. Saturated liquid densities of pure ethylene were predicted in very close agreement with experiment whereas its vapor pressure was reproduced with a deviation of 2%. By comparing their potential model with other two unitedatom models (NERD and TraPPE-UA) available for ethylene, Weitz and Potoff [15] found that their proposed model yields more accurate predictions for the saturated liquid densities, vapor pressures, critical density, and normal boiling point of ethylene. However, no predictions for the heat of vaporization of ethylene were reported by Weitz and Potoff. The potential model for water devised by Huang et al. [16] is an optimized version of the TIP4P model for water proposed 30 years ago by Jorgensen et al [25]. The optimized TIP4P model by Huang et al. comprises three point charges and one LJ site. One negative point charge (NPC)

q  0.8391 e is located at a distance of 0.20482 Å from the oxygen (O) atom, on the line bisecting the HOH angle, and two positive point charges, each with a value q  0.41955 e , are located on each hydrogen (H) atom. Each hydrogen atom is located at a distance of 1.1549 Å from the oxygen atom. The only LJ site of the model is located on the oxygen atom. The local coordinates and charges of the sites and the size and energy parameters of the LJ potential are given in Table 1. With this model, Huang et al. used the grand equilibrium method [26] in combination with NPT-ensemble simulations and the gradual insertion method [27, 28] to compute the VLE of pure water and the binary systems water + ethylene oxide and water + ethylene glycol. The vapor pressure, the saturated liquid densities, and the heat of vaporization of pure water were reproduced with deviations of 7.2%, 1.1%, and 2.8%, respectively. By

9 comparing their potential model with other four united-atom models (TIP4P, TIP4P-Ew, TIP4P / 2005, and SPC / E) available for water, Huang et al. [16] found that their proposed model yields more accurate predictions for the saturated liquid densities, vapor pressures, and heat of vaporization of water. The potential model for ethanol devised by Schnabel et al. [17] comprises three LJ sites and three point charges. The LJ sites are located on the methyl (CH3) and methylene (CH2) united-atom groups, and on the oxygen atom of the hydroxyl (OH) group. Positive point charges q  ,1  0.2556 e and q  , 2  0.44151 e are located on the CH2 group and the hydrogen atom of the

OH group, respectively. A negative charge q  0.69711 e is located on the oxygen atom of the OH group. The CH2 group is located at a distance of 1.9842 Å from the CH3 group. The oxygen atom of the OH group is located at a distance of 1.71581 Å from the CH2 group. The hydrogen and oxygen atoms of the OH group are separated by a distance of 0.95053 Å. The (CH3)(CH2)O and (CH2)OH angles are 90.950° and 106.368°, respectively. The local coordinates and charges of the sites and the size and energy parameters of the LJ potential are given in Table 1. With this model, Schnabel et al. used the NPT + test particle method of Möller and Fischer [29, 30] in combination with the gradual insertion method [27, 28] to compute the VLE of pure ethanol, and the grand equilibrium method [26] to compute the VLE of the binary system ethanol + carbon dioxide. The vapor pressure, the saturated liquid densities, and the heat of vaporization of pure ethanol were reproduced with deviations of 3.7%, 0.3%, and 0.9%, respectively. By comparing their potential model with other two united-atom models (OPLS-UA and TraPPE-UA) available for ethanol, Schnabel et al. [17] found that their proposed model yields more accurate predictions for the saturated liquid densities, vapor pressures, and heat of vaporization of ethanol.

3. Simulation methods In contrast to the methods used in the works by Weitz and Potoff [15], Huang et al. [16], and Schnabel et al. [17], the Gibbs ensemble Monte Carlo (GEMC) method [7-14] was used in the present work, in order to compute the VLE of the pure components ethylene, water, and ethanol, the binary subsystems (ethylene + water, ethylene + ethanol, and ethanol + water) and the ternary mixture (ethylene + water + ethanol). The canonical version (GEMC-NVT) of the GEMC method

10 was applied to the pure components and its isothermal-isobaric version (GEMC-NPT) was applied to the binary and ternary systems. Our two previous implementations of the GEMC method [31, 32] were taken as the starting point to develop the computational strategies explained in Sections 3.1 and 3.2. There are three types of random moves for the GEMC method when dealing with rigid molecules: translational displacement and rotation of molecules inside each of the two simulation boxes, volume changes for the boxes, and transfer of molecules between the two boxes (i.e, simultaneous removal and insertion moves). These moves are accepted or rejected in accordance with a particular probability recipe that involves the calculation of the total intermolecular potential energy change U T  U Tn  U To between the new (trial) ( U Tn ) and the old ( U To ) configurations. The probability formulas for the acceptance or rejection of the three kinds of moves in the GEMC method have been discussed in various papers [7-14] and also in the textbook by Frenkel and Smit [33]. The total intermolecular energy was computed as the sum of non-electrostatic and electrostatic contributions. The non-electrostatic contribution ( U LJ ) was calculated by means of the Lennard-Jones pair potential model and the electrostatic contribution ( U Coul ) was calculated by means of the Ewald summation method [33-37] for the Coulombic potential, as follows: trunc corr real recip self U T  U LJ  U Coul  U LJ  U LJ  U Coul  U Coul  U Coul

(1)

trunc corr is the truncated LJ potential, U LJ is the corresponding LJ long-range correction where U LJ

real recip [33, 34, 38], and U Coul , U Coul , and U LJself are the real space, reciprocal space, and self-interaction

terms of the Ewald sum, respectively [33-37]. Periodic boundary conditions and the minimum image convention, as explained in detail by Allen and Tildesley [34], were applied to the calculation of the total intermolecular potential energy. For the rotational moves required for water and ethanol, the orientational displacements followed the scheme [34] that chooses random values for the Euler angles in the rotation matrix, and employs the internal coordinates of the sites of the molecule (see Table 1) to calculate their simulation-box coordinates. For all

11 simulations in this work, a spherical cutoff distance rc  8.5 Å was used to truncate the LJ potential and calculate its long-range correction from the analytical expression given in Refs. [33, 34] and extended to mixtures by de Pablo and Prausnitz [38]. The value of 8.5 Å for the cutoff distance is a compromise between the values obtained by applying the typical criterion rc  2.5 to the LJ contributions to the potentials models of ethylene, water, and ethanol (see Table 1), and was always less than half the simulation box length so that consistency with the minimum image convention was preserved. Since the Ewald summation method was used in the present work to compute the Coulombic interactions between point charges, then there was no need to truncate those interactions and, accordingly, no need either for such a long cutoff distance (e.g., 17.5 Å) as used in other works [16, 17] for the calculation of the Coulombic interactions by truncation and correction by means of the reaction field method. The time-saving strategy devised by Fartaria et al. [39] was carefully implemented in the present work to obtain a significant decrease in the computer time required for the calculation of the Ewald sum. In that strategy, two repository matrices are used to store pairwise potential energies like the real space term of the Ewald sum and the truncated LJ potential, which can be retrieved when needed, thereby saving time for the calculation of the intermolecular energy associated to trial configurations. Whenever a translational displacement or rotation of a molecule is carried out, the only change experienced by the repository matrices is for the row and column associated to the molecule being moved. Therefore, the change of total configurational energy is given by the summation (over all of the other molecules) of the difference between the trial and the old pair energies of each of those molecules and the molecule being moved. Whenever a trial configuration is accepted, the corresponding pairwise potential energies are recorded in the repository matrices, to be retrieved whenever needed along the simulation run. Fartaria et al. [39] showed that time-savings of up to 52% are possible when this strategy is implemented for molecular models comprising point charges and the Ewald sum is used to compute the Coulombic interactions. Besides the two repository matrices already mentioned, an additional repository matrix was used in the present work, in order to store and retrieve the complexvariable summation [33-37] that runs over molecular sites for each reciprocal space vector, allowing an efficient calculation of changes in the reciprocal space term of the Ewald sum.

12

3.1 Simulation method for pure components For the simulation of VLE for the pure components, the following five-stage strategy was implemented. In the first stage, estimates of the saturated liquid and vapor densities of the pure components at a specified temperature T were obtained by means of reference multiparameter equations of state [18-20]. These equations of state were solved by means of the calculation approach proposed by Mao et al. [40] and Shi and Mao [41] for water and steam, and extending its application to ethylene and ethanol. In the second stage, for each phase of the pure component, an NVT-ensemble simulation (with N  400 molecules), with the corresponding density fixed at the value estimated from the equation of state, was carried out for a total number of 2  10 7 moves (molecular displacements and rotations), 60% of which were typically used to

equilibrate the configurational energy. In the third stage, starting from the final configuration obtained after the NVT runs, an NPT-ensemble simulation, with the pressure fixed at the value of the vapor pressure estimated from the equation of state, was carried out for each phase for a total number of 3 10 7 moves (using a ratio of one box-volume change to N molecular displacements and rotations), 60% of which were typically used to equilibrate the molecular density and the configurational energy. In the fourth stage, starting from the final configurations obtained after the NPT runs, a GEMC-NVT simulation was carried out for the set of two boxes (each of them with an initial number of N  400 molecules), for a total number of at least 1 10 5 moves (using a ratio of one volume change to 2N molecular displacements and rotations and 2N molecular transfers between the boxes), 60% of which were typically used to equilibrate the properties being averaged (molecular densities and configurational energies of the two phases). The attainment of phase equilibrium was ascertained from the statistical equality of the chemical potentials for the two phases, calculated by means of the particle insertion method, and the length of the simulation run was extended until that equality was achieved. Since the equality of the chemical potentials is achieved precisely by means of the transfer moves, the achievement of that equality in the simulation runs indicates that the number of successful transfer moves was sufficient. Statistical uncertainties associated to the GEMC-NVT ensemble-averages were calculated by means of the method by Flyvbjerg and Petersen [42]. In the fifth stage, an

13 additional NVT-ensemble simulation was carried out for the final vapor-phase box obtained after the GEMC-NVT run, in order to compute the pure-component vapor pressure by means of the method of Harismiadis et al. [43]. 3.2 Simulation method for mixtures For the simulation of VLE for the binary and ternary mixtures, the following four-stage strategy was implemented. In the first stage, by specifying both the temperature and the pressure of the system, estimates of the molar compositions and the densities of the coexisting phases were obtained from a T-P flash calculation with the help of the PRSV2-WS-UNIQUAC thermodynamic model described in our previous study [1]. In the second stage, for each phase of the mixture, an NVT-ensemble simulation (with N  400 molecules), with the density and mole fractions fixed at the values estimated from the thermodynamic model, was carried out for a total number of 2  10 7 moves (molecular displacements and rotations), 60% of which were used to equilibrate the configurational energy. In the third stage, starting from the final configuration obtained after the NVT runs, an NPT-ensemble simulation was carried out for each phase for a total number of 3 10 7 moves (using a ratio of one box-volume change to N molecular displacements and rotations), 60% of which were used to equilibrate the molecular density and the configurational energy. In the fourth stage, starting from the final configurations obtained after the NPT runs, a GEMC-NPT simulation was carried out for the set of two boxes (each of them with an initial number of N  400 molecules), for a total number of moves in the range from 1 10 7 to 1.5  10 7 (using a ratio of one volume change for each box to 2N molecular displacements and rotations and 2N molecular transfers between the boxes). Properties of the coexisting phases were sampled every 5  10 5 moves, and running averages were recalculated until the statistical equality for the chemical potentials of each component in the two phases was attained. Statistical uncertainties associated to the GEMC-NPT ensemble-averages were also calculated by means of the method by Flyvbjerg and Petersen [42]. For all VLE simulations for the binary and ternary mixtures, the Lorentz-Berthelot combining rules were used to calculate the size and energy parameters of the LJ potential for the unlike interactions.

4. Simulation results for the pure components

14

By following the five-stage strategy explained in Section 3.1, pure-component VLE molecular simulations were carried out by using the potential models of ethylene, water, and ethanol described in Section 2. The VLE phase diagrams of the pure components were computed over the temperature intervals [160 K, 265 K] for ethylene, [300 K, 600 K] for water, and [270 K, 493 K] for ethanol. These phase diagrams are shown in Figs. 1-3 as graphs of absolute temperature vs. molar density of the two phases, where the empty squares correspond to the molecular simulation results and the solid lines correspond to the results obtained by means of the most accurate (reference) multiparameter equations of state currently available: the equation of Smukala et al. [18] for ethylene, the equation of Wagner and Pruß [19] for water, and the equation of Dillon and Penoncello [20] for ethanol. A comprehensive account of multiparameter equations of state as an accurate source of thermodynamic property data has been given in the book by Span [44]. These equations are regarded as fundamental equations of state because they are based on the Helmholtz free energy. The reference character of these multiparameter equations of state is supported by the very low deviations that are associated to the thermodynamic properties calculated from those equations. Indeed, the reference equation of state by Smukala et al. [18] for ethylene reproduces the experimental data for the saturated liquid and vapor densities with deviations within 0.004% and 0.02%, respectively, and the data for the vapor pressure with a deviation less than 0.01%. The reference equation of state by Wagner and Pruß [19] for water, which is known as the IAPWS-95 formulation, reproduces the experimental data for the saturated liquid and vapor densities with deviations within 0.0025% and 0.1%, respectively, and the data for the vapor pressure with a deviation within 0.025%. The reference equation of state by Dillon and Penoncello [20] for ethanol reproduces the experimental data for the saturated liquid and vapor densities with deviations less than 0.2% and 0.5%, respectively, and the data for the vapor pressure with a deviation less than 0.4%. These reference equations of state were compared by their proponents [18-20] with previous equations and were found to be much more accurate. The average relative deviation (ARD) of the simulation results for the saturated liquid densities with respect to the values obtained from the corresponding reference equation of state is 2.5% for

15 ethylene, 1.1% for water, and 1.8% for ethanol. An assessment of the deviation for ethylene cannot be made because no corresponding deviation value was reported by Weitz and Potoff [15], using the same potential model for ethylene. The deviation for water agrees perfectly with the value (1.1%) reported by Huang et al [16], using the same potential model for water but a different simulation method. In contrast, the deviation for ethanol (1.8%) turns out to be larger than the value (0.3%) reported by Schnabel et al. [17], using the same potential model for ethanol but a different simulation method. As mentioned in Section 3.1, the vapor pressure of the pure components was computed by applying the method of Harismiadis et al. [43] to the final vapor-phase simulation box obtained after each GEMC-NVT simulation run. The resulting Clausius-Clapeyron diagrams (i.e. graphs of

ln P sat vs. 1 / T ) are shown in Figs. 4-6, where the emtpy squares correspond to the molecular simulation results and the solid line corresponds to the results obtained by means of the respective reference equation of state [18-20]. The ARD of the simulation results for the vapor pressure with respect to the values obtained from the corresponding reference equation of state is 5.0% for ethylene, 9.6% for water, and 4.1% for ethanol. The deviation for ethylene (5.0%) turns out to be much larger than the value (2%) reported by Weitz and Potoff [15], using the same potential model for ethylene but a different simulation method. The deviation for water (9.6%) is 33% larger than the value (7.2%) reported by Huang et al. [16], using the same potential model for water but a different simulation method. In contrast, the deviation for ethanol (4.1%) is just 11% larger than the value (3.7%) reported by Schnabel et al. [17], using the same potential model for ethanol but a different simulation method. Simulation results for the vapor pressure and the configurational energies and densities of the coexisting phases at a given temperature were used to compute the heat of vaporization ( H v ) of the pure components from the following thermodynamic relationship:

 1 1 H v  (U T(V )  U T( L ) )  P sat  (V )  ( L )  

  

(2)

16 where U T(V ) and U T( L ) are the total configurational energies of the vapor and liquid phases, respectively, and  (V ) and  ( L ) are their corresponding densities. The resulting graphs of H v vs. T are shown in Figs. 7-9, where the emtpy squares correspond to the molecular simulation results and the solid line corresponds to the results obtained by means of the respective reference equation of state [18-20]. The ARD of the simulation results for the enthalpy of vaporization with respect to the values obtained from the corresponding reference equation of state is 1.7% for water, 4.1% for ethanol, and 12.4% for ethylene. The deviation for water (1.7%) turns out to be 39% smaller than the value (2.8%) reported by Huang et al. [16], using the same potential model for water but a different simulation method. In contrast, the deviation for ethanol (4.1%) turns out to be much larger than the value (0.9%) reported by Schnabel et al. [17], using the same potential model for ethanol but a different simulation method. For ethylene, simulation results exhibit a significant offset (by overestimation) with respect to the equation of state results (see Fig. 7). Since the heat of vaporization was not reported by Weitz and Potoff [15] (proponents of the potential model of ethylene), a direct assessment of the seemingly large deviation (12.4%) obtained in the present work cannot be made. However, similarly large deviations for the heat of vaporization predicted by other LJ plus point-charge potential models have been reported, even when the vapor pressure and the saturated liquid densities are reproduced accurately. For instance, Huang et al. [16] reported an ARD (with respect to the experimental data) of 13.4% for the heat of vaporization predicted by their potential model of ethylene glycol (see Fig. 4 in Ref. [16], in which an offset by overestimation is also shown), even though the saturated liquid densities and vapor pressure of that substance were reproduced with acceptable accuracy by the same potential model (see Figs. 2 and 3 in Ref. [16]). As shown in Figs. 1-9, the statistical uncertainties (error bars) of our GEMC simulation results are small (most of them are within or smaller than the symbol sizes). Therefore, the resulting discrepancy for the average deviations between the present and the previous works is more likely due to the difference between the simulation techniques used. The good agreement of the GEMC simulation results for the vapor pressure and the VLE phase diagrams of the pure components, with respect to calculations carried out by means of the most

17 accurate (reference) equations of state currently available for those components, corroborates the validity of the potential models proposed by Weitz and Potoff [15] for ethylene, Huang et al. [16] for water, and Schnabel et al. [17] for ethanol.

5. Simulation results for the binary mixtures By following the four-stage strategy explained in Section 3.2, binary-mixture VLE molecular simulations were carried out by using the potential models of ethylene, water, and ethanol described in Section 2. Isothermal VLE phase diagrams of the binary subsystems were computed at the following temperatures: 200, 250, and 300°C for ethylene + water, 150, 170, 190, 200, and 220°C for ethylene + ethanol, and 200, 250, 275, and 300°C for ethanol + water. The resulting phase diagrams at 200°C are shown in Figs. 10-13 as graphs of pressure vs. molar composition of the two phases (P-xy diagrams), where the empty squares (with their respective error bars) correspond to the molecular simulation results, the crosses correspond to the available experimental data [21-23], and the solid lines correspond to the results obtained by means of the PRSV2-WS-UNIQUAC thermodynamic model described in our previous study [1]. The resulting phase diagrams at the other temperatures listed above are shown as Figs. S1-S9 in the Supporting Information. 5.1 Results for ethylene + water The isothermal VLE (P-xy) phase diagrams for the binary system ethylene + water are shown in Figs. 10 and 11 at 200°C, and in Figs. S1 and S2 (in the Supporting Information) at 250 and 300°C, respectively. The parameter values reported in Table 2 were used for the calculations with the PRSV2-WS-UNIQUAC thermodynamic model. The solid line and data points at the left side of the diagrams correspond to the liquid phase and the solid line and data points at the right side of the diagrams correspond to the vapor phase. The three sets of data (molecular simulation, thermodynamic model, and experimental measurements) consistently follow the same increasing trend for the mole fraction of ethylene as the pressure increases. At 200 and 250°C, simulation results for the vapor phase exhibit smaller error bars and are closer to both the experimental data and the results obtained from the thermodynamic model, as compared to the simulation results

18 for the liquid phase. At 300°C, as the pressure increases, simulation results for the vapor phase exhibit a monotonic deviation to the left of the corresponding experimental data (see Fig. S2), and the simulation result for the liquid phase at 160 atm is shifted to the right, significantly. A possible interpretation of this behavior is that simulation would be predicting a much narrower phase envelope for ethylene + water at 300°C, which is a surprising result considering the trend of the experimental data. At 200°C, the average absolute deviation (AAD) of the simulation results (ensemble-averages) for the mole fraction of ethylene with respect to the experimental values, is 0.0030 for the liquid phase and 0.0345 for the vapor phase. At 250°C, the AAD is 0.0015 for the liquid phase and 0.0240 for the vapor phase. At 300°C, the AAD is 0.0050 for the liquid phase and 0.0919 for the vapor phase. Additional simulations were carried out in order to assess the effect of including a correction factor  in the Lorentz combining rule for the energy parameter  ij of the LJ potential for unlike interactions, so that  i , j   ( i  j )1 / 2 . By means of a simple random search, several values of the correction factor  were tried for the simulations at 200°C and the best results were obtained with   0.9 , for which the resulting VLE phase diagram is shown in Fig. 11. The AAD of the simulation results for the mole fraction of ethylene at 200°C decreased from 0.0030 to 0.0017 for the liquid phase and from 0.0345 to 0.0189 for the vapor phase, showing an improvement with respect to the use of the default value   1.0 . This correction factor was not applied to the simulations at 250 and 300°C. From an overall examination of the VLE phase diagrams shown in Figs. 10, 11, and S1, it follows that the intermolecular potential models for ethylene and water described in Section 2 are capable of qualitatively predicting the VLE behavior of this binary system at 200 and 250°C.

5.2 Results for ethylene + ethanol

19 The VLE phase diagrams for the binary system ethylene + ethanol are shown in Fig. 12 at 200°C and in Figs. S3-S6 (in the Supporting Information) for temperatures of 150, 170, 190, and 220°C. The parameter values reported in Table 2 were used for the calculations with the PRSV2-WSUNIQUAC thermodynamic model. Once again, the solid line and data points at the left side of the diagrams correspond to the liquid phase and the solid line and data points at the right side of the diagrams correspond to the vapor phase. The three sets of data (molecular simulation, thermodynamic model, and experimental measurements) follow the same trend for the mole fraction of ethylene as the pressure increases. At 150°C and in the middle range of pressures, simulation results for the vapor phase exhibit larger deviations with respect to both the experimental data and the results obtained from the thermodynamic model. At 150°C, the average absolute deviation (AAD) of the simulation results (ensemble-averages) for the mole fraction of ethylene with respect to the experimental values, is 0.0332 for the liquid phase and 0.0781 for the vapor phase. At 170°C, the AAD is 0.0275 for the liquid phase and 0.0476 for the vapor phase. At 190°C, the AAD is 0.0174 for the liquid phase and 0.0196 for the vapor phase. At 200°C, the AAD is 0.0143 for the liquid phase and 0.0189 for the vapor phase. At 220°C, the AAD is 0.0085 for the liquid phase and 0.0156 for the vapor phase. Even though both molecular simulations and the T-P flash calculations made with the help of the thermodynamic model become increasingly difficult in the region near the critical point of the binary mixture, the results associated to these two methods are in very good agreement in that region for some of the phase diagrams, e.g., at temperatures of 170 and 220°C (Figs. S4 and S6). The experimental critical points of this binary system were reported by Tsiklis and Kofman [22] for each of the five values of temperature corresponding to the VLE phase diagrams shown in Figs. 12 and S3-S6. In each of these diagrams the experimental critical point is indicated by the highest pressure cross. At 150 and 170°C (see Figs. S3 and S4), the critical pressures predicted from both molecular simulation and the thermodynamic model are very close to the experimental critical pressures of 127 and 118 atm, respectively [22]. At 190 and 200°C (see Figs. S5 and 12), the critical pressures predicted by extrapolation from the molecular simulation results are very close to the experimental critical pressures of 105 and 98 atm, respectively [22], whereas the

20 critical pressures predicted by means of the thermodynamic model turn out to be slightly higher than the experimental ones. At 220°C (see Fig. S6), both molecular simulation and the thermodynamic model predict a critical pressure about 8% higher than the experimental value of 81 atm [22]. From an overall examination of the VLE phase diagrams shown in Figs. 12 and S3-S6, it follows that the intermolecular potential models for ethylene and ethanol described in Section 2 are capable of qualitatively predicting the VLE behavior of this binary system.

5.3 Results for ethanol + water The VLE phase diagrams for the binary system ethanol + water are shown in Fig. 13 at 200°C and in Figs. S7-S9 (in the Supporting Information) for temperatures of 250, 275, and 300°C. The parameter values reported in Table 2 were used for the calculations with the PRSV2-WSUNIQUAC thermodynamic model. Once again, the solid line and data points at the left side of the diagrams correspond to the liquid phase and the solid line and data points at the right side of the diagrams correspond to the vapor phase. The three sets of data (molecular simulation, thermodynamic model, and experimental measurements) follow the same trend for the mole fraction of ethanol as the pressure increases. Even though simulation results for the vapor phase in the middle range of pressures at 200 and 250°C (Figs. 13 and S7) deviate from both the experimental data and the results obtained from the thermodynamic model, molecular simulation shows to be capable of making an accurate prediction of the azeotropic point of this binary system. At 200°C, the average absolute deviation (AAD) of the simulation results (ensemble-averages) for the mole fraction of ethanol with respect to the experimental values, is 0.0160 for the liquid phase and 0.0557 for the vapor phase. At 250°C, the AAD is 0.0215 for the liquid phase and 0.0254 for the vapor phase. At 275°C, the AAD is 0.0096 for the liquid phase and 0.0122 for the vapor phase. At 300°C, the AAD is 0.0134 for the liquid phase and 0.0180 for the vapor phase.

21 From an overall examination of the VLE phase diagrams shown in Figs. 13 and S7-S9, it follows that the intermolecular potential models for ethanol and water described in Section 2 are capable of qualitatively predicting the VLE behavior of this binary system.

6. Simulation results for the ternary mixture By following the four-stage strategy explained in Section 3.2, ternary-mixture VLE molecular simulations were carried out by using the potential models of ethylene, water, and ethanol described in Section 2. Due to the improved results obtained for the VLE phase diagram of ethylene + water at 200°C by the use of a correction factor   0.9 in the Lorentz combining rule for the LJ energy parameter of unlike interactions (see Section 5.1), such correction factor was also used in the VLE simulations of the ternary mixture. In Section 5, VLE phase diagrams for the binary systems were computed at various temperatures, 200°C being a common value for the three systems with regard to the availability of experimental data [21-23] and simulation results. Due to the availability of experimental data [6] for the ternary system at 200°C as well, the isothermal-isobaric VLE phase diagrams of the ternary system were also computed at 200°C and at the following pressures: 30, 40, 50, 60, 80, and 100 atm. The resulting phase diagrams are shown in Figs. 14-19, where the filled circles correspond to the molecular simulation results, the solid lines correspond to the results obtained from the PRSV2-WS-UNIQUAC thermodynamic model, and the triangles correspond to the few experimental data available [6]. In the ternary diagrams, the upper locus corresponds to the vapor phase at its dew point and the lower locus corresponds to the liquid phase at its bubble-point. The two loci define the saturation envelope for the ternary system, the two phases of which coexist in the region within the envelope. In contrast to the few experimental data points available and included in the phase diagrams, both the simulation and the model predictions extend over larger ranges and follow the same trend for the dew-point and bubble-point lines of the ternary system. At pressures of 40 and 50 atm (Figs. 15 and 16), simulation results accurately match the results obtained by means of the thermodynamic model. At 100 atm (Fig. 19), molecular simulation predicts a slightly narrower two-phase region, indicating that the critical point of the binary mixture ethylene + ethanol at 200°C would occur at a pressure slightly below 100 atm. At that pressure, the bubble-point and

22 dew-point lines of the ternary system would intersect the ethylene-ethanol side of the phase diagram at the composition of that critical point. In contrast, the thermodynamic model predicts that this critical point would occur at a pressure slightly above 100 atm. From the experimental VLE measurements reported by Tsiklis and Kofman [22] for the binary system ethylene + ethanol, the critical point occurs at a pressure of 98 atm at 200°C; therefore, the molecular simulation prediction of a value slightly below 100 atm for the critical pressure at 200°C turns out to be more accurate than the prediction made by means of the thermodynamic model (see Fig. 12).

7. Conclusions From the good agreement of the Gibbs-ensemble Monte Carlo simulation results for the vapor pressure and the VLE phase diagrams of ethylene, water, and ethanol with respect to calculations carried out by means of the most accurate (reference) multiparameter equations of state currently available for those components [18-20], we were able to corroborate the validity of a set of previously published intermolecular potential models for ethylene [15], water [16], and ethanol [17], which had been validated by their proponents from results obtained by means of other simulation methods. These potential models are capable of predicting the available VLE phase diagrams of the binary subsystems ethylene + water [21] (at 200 and 250°C), ethylene + ethanol [22] (at 150, 170, 190, 200, and 220°C) and ethanol + water [23] (at 200, 250, 275, and 300°C). For the binary system ethylene + water, it was shown that the use of a correction factor of 0.9 for the Lorentz combining rule, which is close to the default value of 1.0, improves the calculated phase diagram at 200°C by decreasing by 45% the corresponding deviations for the mole fractions of ethylene in the vapor and liquid phases. Since for most of the simulated binary VLE phase diagrams (11 diagrams out of a total of 12) a correction factor for the Lorentz combining rule was not needed, and the value of the correction factor used for ethylene + water at 200°C is not far from the default value of 1.0, then it can be said that the Lorentz-Berthelot combining rules appear to be appropriate for the three binary systems considered in the present work. This suitability of the

23 Lorentz-Berthelot combining rules for VLE calculations from molecular simulation is in agreement with the findings of previous works (e.g., [31, 32, 45, 46]). Molecular simulation predictions for the VLE phase diagrams of the ternary system at 200°C and pressures of 30, 40, 50, 60, 80, and 100 atm are in very good agreement with predictions that we had previously made [1] by use of a thermodynamic model that combines the Peng-RobinsonStryjek-Vera equation of state, the Wong-Sandler mixing rules, and the UNIQUAC activity coefficient model. The latter agreement is encouraging for the subsequent use of molecular simulation to predict the combined phase and chemical equilibrium of the ternary system and check the validity of predictions that we previously made [1] by means of the thermodynamic model.

Acknowledgement Financial support from the Colombian Administrative Department of Science, Innovation and Technology (COLCIENCIAS), through a research assistanship for doctoral students (Y.M. Muñoz-Muñoz), is gratefully acknowledged.

24

References [1] M. Llano-Restrepo, Y.M. Muñoz-Muñoz, Fluid Phase Equilibr. 307 (2011) 45-57. [2] D.-Y. Peng, D.B. Robinson, Ind. Eng. Chem. Fundam. 15 (1976) 59-64. [3] R. Stryjek, J.H. Vera, Can. J. Chem. Eng. 64 (1986) 323-333. [4] D.S.H. Wong, S.I. Sandler, AIChE J. 38 (1992) 671-680. [5] D.S. Abrams, J.M. Prausnitz, AIChE J. 21 (1975) 116-128. [6] D.S. Tsiklis, A.I. Kulikova, L.I. Shenderei, Khim. Promst. (5) (1960) 401-406. [7] A.Z. Panagiotopoulos, Mol. Phys. 61 (1987) 813-826. [8] A.Z. Panagiotopoulos, N. Quirke, M. Stapleton, D.J. Tildesley, Mol. Phys. 63 (1988) 527-545. [9] A.Z. Panagiotopoulos, M.R. Stapleton, Fluid Phase Equilibr. 53 (1989) 133-141. [10] A.Z. Panagiotopoulos, Int. J. Thermophys. 10 (1989) 447-457. [11] B. Smit, Ph. de Smedt, D. Frenkel, Mol. Phys. 68 (1989) 931-950. [12] B. Smit, D. Frenkel, Mol. Phys. 68 (1989) 951-958. [13] A.Z. Panagiotopoulos, Mol. Simul. 9 (1992) 1-23. [14] A.Z. Panagiotopoulos, J. Phys. Condens. Matter 12 (2000) 25-52. [15] S.L. Weitz, J.J. Potoff, Fluid Phase Equilibr. 234 (2005) 144-150. [16] Y.-L. Huang, T. Merker, M. Heilig, H. Hasse, J. Vrabec, Ind. Eng. Chem. Res. 51 (2012) 74287440. [17] T. Schnabel, J. Vrabec, H. Hasse, Fluid Phase Equilibr. 233 (2005) 134-143. [18] J. Smukala, R. Span, W. Wagner, J. Phys. Chem. Ref. Data 29 (2000) 1053-1121. [19] W. Wagner, A. Pruß, J. Phys. Chem. Ref. Data 31 (2002) 387-535. [20] H.E. Dillon, S.G. Penoncello, Int. J. Thermophys. 25 (2004) 321-335. [21] D.S. Tsiklis, E.V. Mushkina, L.I. Shenderei, Inzh. Fiz. Zh. 1 (8) (1958) 3-7.

25 [22] D.S. Tsiklis, A.N. Kofman, Russ. J. Phys. Chem. 35 (1961) 549-551. [23] F. Barr-David, B.F. Dodge, J. Chem. Eng. Data 4 (1959) 107-121. [24] J.J. Potoff, J.R. Errington, A.Z. Panagiotopoulos, Mol. Phys. 97 (1999) 1073-1083. [25] W.L. Jorgensen, J. Chandrasekhar, J.D. Madura, R.W. Impey, M.L. Klein, J. Chem. Phys. 79 (1983) 926-935. [26] J. Vrabec, H. Hasse, Mol. Phys. 100 (2002) 3375-3383. [27] I. Nezbeda, J. Kolafa, Mol. Simul. 5 (1991) 391-403. [28] J. Vrabec, M. Kettler, H. Hasse, Chem. Phys. Lett. 356 (2002) 431-436. [29] D. Möller, J. Fischer, Mol. Phys. 69 (1990) 463-473. [30] D. Möller, J. Fischer, Fluid Phase Equilibr. 100 (1994) 35-61. [31] J. Carrero-Mantilla, M. Llano-Restrepo, Fluid Phase Equilibr. 208 (2003) 155-169. [32] J. Carrero-Mantilla, M. Llano-Restrepo, Mol. Simul. 29 (2003) 549-554. [33] D. Frenkel, B. Smit, Understanding Molecular Simulation: from Algorithms to Applications, Second Edition, Academic Press, 2002. [34] M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Oxford University Press, 1989. [35] M. Llano-Restrepo, Molecular modeling and Monte Carlo simulation of concentrated aqueous alkali halide solutions at 25°C, Ph.D. Dissertation, Rice University, Houston, 1994. [36] M. Llano-Restrepo, W.G. Chapman, J. Chem. Phys. 100 (1994) 8321-8339. [37] D.C. Rapaport, The Art of Molecular Dynamics Simulation, Cambridge University Press, 1998. [38] J.J. de Pablo, J.M. Prausnitz, Fluid Phase Equilibr. 53 (1989) 177-189. [39] R.P. Fartaria, R.S. Neves, P.C. Rodrigues, F.F. Freitas, F. Silva-Fernandes, Comp. Phys. Comm. 175 (2006) 116-121. [40] S. Mao, Z. Duan, J. Hu, Z. Zhang, L. Shi, Phys. Earth Planet. Inter. 185 (2011) 53-60. [41] L. Shi, S. Mao, Geosci. Front. 3 (2012) 51-58. [42] H. Flyvbjerg, H.G. Petersen, J. Chem. Phys. 91 (1989) 461-466. [43] V. Harismiadis, J. Vorholz, A. Panagiotopoulos, J. Chem. Phys. 105 (1996) 8469. [44]. R. Span, Multiparameter Equations of State: An Accurate Source of Thermodynamic Property Data, Springer Verlag, Berlin, 2000. [45] N. Ferrando, P. Ungerer, Fluid Phase Equilibr. 254 (2007) 211-223. [46] C.G. Pereira, L. Grandjean, S. Betoulle, N. Ferrando, C. Féjean, R. Lugo, J.C. de Hemptinne, P. Mougin, Fluid Phase Equilibr. 382 (2014) 219-234.

26

Vitae Y. Mauricio Muñoz-Muñoz received his undergraduate degree (B.S) in chemical engineering from Universidad Nacional de Colombia, at the Manizales campus, in 2007, and his doctoral degree (in chemical engineering) from Universidad del Valle, Cali, Colombia, in November 2014. He is currently working as a postdoctoral researcher in the Thermodynamics and Energy Technology research group of Professor Jadran Vrabec, at the Faculty of Mechanical Engineering of the University of Paderborn, in Germany.

Mario Llano-Restrepo is Professor of Chemical Engineering at Universidad del Valle in Cali, Colombia. He received his Ph.D. degree (in chemical engineering) from Rice University in 1994. His main research interests are modeling of phase and chemical equilibria, modeling and simulation of separation processes, and molecular simulation. In 2009, he earned a Distinguished Professor recognition from Universidad del Valle for service and excellence in teaching since 1994. Some of the courses he has taught are chemical thermodynamics, statistical thermodynamics, molecular simulation, chemical reaction engineering, mass transfer, separation operations, and separation process modeling and simulation.

27

Figure Captions Figure 1. VLE phase diagram for ethylene. , molecular simulation results of this work; , calculated with the reference multiparameter equation of state of Smukala et al. [18] for ethylene.

Figure 2. VLE phase diagram for water. , molecular simulation results of this work; , calculated with the reference multiparameter equation of state of Wagner and Pruß [19] for water.

Figure 3. VLE phase diagram for ethanol. , molecular simulation results of this work; , calculated with the reference multiparameter equation of state of Dillon and Penoncello [20] for ethanol.

Figure 4. Clausius-Clapeyron plot for the vapor pressure of ethylene as a function of temperature. , molecular simulation results of this work; , calculated with the reference multiparameter equation of state of Smukala et al. [18] for ethylene.

28

Figure 5. Clausius-Clapeyron plot for the vapor pressure of water as a function of temperature. , molecular simulation results of this work; , calculated with the reference multiparameter equation of state of Wagner and Pruß [19] for water.

Figure 6. Clausius-Clapeyron plot for the vapor pressure of ethanol as a function of temperature. , molecular simulation results of this work; , calculated with the reference multiparameter equation of state of Dillon and Penoncello [20] for ethanol.

Figure 7. Heat of vaporization for ethylene as a function of temperature. , molecular simulation results of this work; , calculated with the reference multiparameter equation of state of Smukala et al. [18] for ethylene.

Figure 8. Heat of vaporization for water as a function of temperature. , molecular simulation results of this work; , calculated with the reference multiparameter equation of state of Wagner and Pruß [19] for water.

Figure 9. Heat of vaporization for ethanol as a function of temperature. , molecular simulation results of this work; , calculated with the reference multiparameter equation of state of Dillon and Penoncello [20] for ethanol.

Figure 10. VLE phase diagram of ethylene (1) + water (2) at 200°C. , molecular simulation results of this work; , calculated with the PRSV2-WS-UNIQUAC thermodynamic model; +, experimental data of Tsiklis et al. [21].

Figure 11. VLE phase diagram of ethylene (1) + water (2) at 200°C. , molecular simulation results of this work using a correction factor χ = 0.9 for the Lorentz combining rule; ,

29

calculated with the PRSV2-WS-UNIQUAC thermodynamic model; +, experimental data of Tsiklis et al. [21].

Figure 12. VLE phase diagram of ethylene (1) + ethanol (2) at 200°C. , molecular simulation results of this work; , calculated with the PRSV2-WS-UNIQUAC thermodynamic model; +, experimental data of Tsiklis and Kofman [22].

Figure 13. VLE phase diagram of ethanol (1) + water (2) at 200°C. , molecular simulation results of this work; , calculated with the PRSV2-WS-UNIQUAC thermodynamic model; +, experimental data of Barr-David and Dodge [23].

Figure 14. VLE phase diagram of ethylene + water + ethanol at 200°C and 30 atm. ●, molecular simulation results of this work; , calculated with the PRSV2-WS-UNIQUAC thermodynamic model; ▲, experimental data of Tsiklis et al. [6].

Figure 15. VLE phase diagram of ethylene + water + ethanol at 200°C and 40 atm. ●, molecular simulation results of this work; , calculated with the PRSV2-WS-UNIQUAC thermodynamic model; ▲, experimental data of Tsiklis et al. [6].

Figure 16. VLE phase diagram of ethylene + water + ethanol at 200°C and 50 atm. ●, molecular simulation results of this work; , calculated with the PRSV2-WS-UNIQUAC thermodynamic model; ▲, experimental data of Tsiklis et al. [6].

Figure 17. VLE phase diagram of ethylene + water + ethanol at 200°C and 60 atm. ●, molecular simulation results of this work; , calculated with the PRSV2-WS-UNIQUAC thermodynamic model; ▲, experimental data of Tsiklis et al. [6].

30

Figure 18. VLE phase diagram of ethylene + water + ethanol at 200°C and 80 atm. ●, molecular simulation results of this work; , calculated with the PRSV2-WS-UNIQUAC thermodynamic model; ▲, experimental data of Tsiklis et al. [6].

Figure 19. VLE phase diagram of ethylene + water + ethanol at 200°C and 100 atm. ●, molecular simulation results of this work; , calculated with the PRSV2-WS-UNIQUAC thermodynamic model; ▲, experimental data of Tsiklis et al. [6].

Table 1 Local coordinates (x, y) and charge q of the sites, and size (  ) and energy (  ) parameters of the LJ-site contribution to the planar (z = 0 Å) intermolecular potential models of ethylene [15], water [16], and ethanol [17].

Potential model of ethylene [15] Interaction site

x (Å)

y (Å)

q (e)

 (Å)

 / k B (K)

CH2

0.665

0

0.85

3.72

83.0

CH2

-0.665

0

0.85

3.72

83.0

COM

0

0

-1.70

0

0

Potential model of water [16] O

0

0

0

3.11831

208.08

H

0.7069

0.9133

0.41955

0

0

H

0.7069

-0.9133

0.41955

0

0

NPC

0.20482

0

-0.83910

0

0

Potential model of ethanol [17]

31

CH3

0

0

0

3.6072

120.15

CH2

1.9842

0

0.25560

3.4612

86.291

O

2.0127

1.7156

-0.69711

3.1496

85.053

H

2.9290

1.9683

0.44151

0

0

Table 2

Optimum values of the adjustable parameters for the fit of the PRSV2-WS-UNIQUAC thermodynamic model [1] to the VLE experimental data of Tsiklis et al [21] for ethylene (1) + water (2), Tsiklis and Kofmann [22] for ethylene (1) + ethanol (3), and Barr-David and Dodge [23] for ethanol (3) + water (2).

System (ij)

T (°C)

u ij / R (K)

u ji / R (K)

k ij

12

200

808.94

385.66

1.0181

12

250

-603.75

2007.55

1.1301

12

300

-936.45

1685.95

2.7033

13

150

394.7

94.2

0.115

13

170

306.1

100.0

0.145

13

190

350.0

100.0

0.235

13

200

420.0

100.0

0.260

13

220

155.0

65.0

0.350

32

200

23.67

369.45

0.0085

32

250

51.92

348.74

0.0003

32

275

-22.60

395.02

0.0157

32

300

-293.78

770.99

0.0737

32

280

260

T/K

240

220

200

180

160 0

5

10

ρV , ρL / mol L-1

Figure 1

15

20

700

T/K

600

500

400

300 0

10

20

30

ρV , ρL / mol L-1

Figure 2

40

50

60

550

500

T/K

450

400

350

300

250 0

5

10

ρV , ρL / mol L-1

Figure 3

15

20

4

2

ln( P

sat

/ bar)

3

1

0

-1 3.5

4.0

4.5

5.0

1000 K / T

Figure 4

5.5

6.0

6.5

6

2

ln( P

sat

/ bar)

4

0

-2

-4 1.6

1.8

2.0

2.2

2.4

2.6

1000 K / T

Figure 5

2.8

3.0

3.2

3.4

4

0

ln( P

sat

/ bar)

2

-2

-4

2.0

2.5

3.0

1000 K / T

Figure 6

3.5

4.0

18 16

-1 ∆Hv / kJ mol

14 12 10 8 6 4 2 0 160

180

200

220

T/K

Figure 7

240

260

280

50

-1 ∆Hv / kJ mol

40

30

20

10

0 300

400

500

T/K

Figure 8

600

700

50

-1 ∆Hv / kJ mol

40

30

20

10

0 250

300

350

400

T/K

Figure 9

450

500

160

140

140

120

120

100

100

80

80

60

60

40

40 0.000

0.005

0.010

x1

Figure 10

0.015

0.020 0.4

0.5

0.6

0.7

y1

0.8

0.9

1.0

P / atm

P / atm

160

160

140

140

120

120

100

100

80

80

60

60

40

40 0.000

0.005

x1

Figure 11

0.010 0.4

0.5

0.6

0.7

y1

0.8

0.9

1.0

P / atm

P / atm

160

100

P / atm

80

60

40

20 0.0

0.1

0.2

x1 , y1

Figure 12

0.3

0.4

32 30 28

P / atm

26 24 22 20 18 16 0.0

0.2

0.4

0.6

x1 , y1

Figure 13

0.8

1.0

T = 200 [°C], P = 30 [atm] C2H4 0.0 0.1

1.0 0.9

0.2

0.8

0.3

0.7

0.4

0.6

0.5

0.5

0.6

0.4

0.7

0.3

0.8

0.2

0.9

0.1

1.0

CH3CH2OH

0.0

Figure 14

0.0 0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

H2O

T = 200 [°C], P = 40 [atm] C2H4 0.0 0.1

1.0 0.9

0.2

0.8

0.3

0.7

0.4

0.6

0.5

0.5

0.6

0.4

0.7

0.3

0.8

0.2

0.9

0.1

1.0

CH3CH2OH

0.0

Figure 15

0.0 0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

H2O

T = 200 [°C], P = 50 [atm] C 2H4 0.0 0.1

1.0 0.9

0.2

0.8

0.3

0.7

0.4

0.6

0.5

0.5

0.6

0.4

0.7

0.3

0.8

0.2

0.9

0.1

1.0

CH3CH2OH

0.0

Figure 16

0.0 0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

H2O

T = 200 [°C], P = 60 [atm] C2H4 0.0 0.1

1.0 0.9

0.2

0.8

0.3

0.7

0.4

0.6

0.5

0.5

0.6

0.4

0.7

0.3

0.8

0.2

0.9

0.1

1.0

CH3CH2OH

0.0

Figure 17

0.0 0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

H2O

T = 200 [°C], P = 80 [atm] C2H4 0.0 0.1

1.0 0.9

0.2

0.8

0.3

0.7

0.4

0.6

0.5

0.5

0.6

0.4

0.7

0.3

0.8

0.2

0.9

0.1

1.0

CH3CH2OH

0.0

Figure 18

0.0 0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

H2O

T = 200 [°C], P = 100 [atm] C2H4 0.0 0.1

1.0 0.9

0.2

0.8

0.3

0.7

0.4

0.6

0.5

0.5

0.6

0.4

0.7

0.3

0.8

0.2

0.9

0.1

1.0

CH3CH2OH

0.0

Figure 19

0.0 0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

H2O