Phase equilibria of binary mixtures by molecular simulation and cubic equations of state

PHASE EQUILIBRIA OF BINARY MIXTURES BY MOLECULAR SIMULATION AND CUBIC EQUATIONS OF STATE

V.F.Cabral¹, R.R.C.Pinto^1,2, F.W.Tavares1,^* and M.Castier¹

¹Escola de Química, Universidade Federal do Rio de Janeiro

Cx. P. 68542, 21949-900, Rio de Janeiro - RJ, Brazil

E-mail: [e-mail]

²Centro de Pesquisas e Desenvolvimento Leopoldo Américo M. de Mello,

PETROBRAS, Rio de Janeiro - RJ, Brazil

(Received: March 9, 2000 ; Accepted: May 15, 2001)

Abstract - Molecular simulation data were used to study the performance of equations of state (EoS) and combining rules usually employed in thermodynamic property calculations. The Monte Carlo method and the Gibbs ensemble technique were used for determining composition and densities of vapor and liquid phases in equilibrium for binary mixtures of Lennard-Jones fluids. Simulation results are compared to data in the literature and to those calculated by the t-PR-LJ EoS. The use of adequate combining rules has been shown to be very important for the satisfactory representation of molecular simulation data.

Keywords: excess properties, combining rules, equations of state and molecular simulation

INTRODUCTION

Molecular simulation techniques are increasingly being used to study the thermodynamic properties of systems. Even though these simulations can not be considered to be substitutes for experiments, their results may be particularly useful where experimental measurement is considerably difficult, such as for confined fluids or fluids under extreme temperature and pressure conditions. Molecular simulation data are also used to test the statistical mechanics approximations used in the development of thermodynamic models, as for example, in some of the equations of state (EoS) used in chemical engineering. However, the use of molecular simulation techniques is limited, because of the large amount of computation time required. For this reason, simulated systems typically contain a few hundred particles and use relatively simple semiempirical expressions for the interparticle potential.

Despite evidence that ternary interparticle interactions have an effect on the configurational energy and phase equilibria behavior and the existence of several new theories (Gubbins, 1996), these models require a very large amount of computation time. Therefore, pairwise simplified potentials are often used, and one of these, the Lennard-Jones (LJ) potential is extensively used due to its programming simplicity and the availability of published parameter data for several systems. The LJ potential can represent quite well the qualitative behavior of nonpolar symmetric molecules, even though it provides poor predictions of phase equilibria for systems whose densities are either very low or very high.

There are two main techniques in molecular simulation: the Molecular Dynamics (MD) and Monte Carlo (MC) methods. The former technique is based on the solution of time dependent equations and may be used to estimate transport properties and time-dependent functions. MC methods are based on sampling the phase space and are mainly useful for calculating equilibrium properties. In the field of MC methods, a significant advance was achieved with the Gibbs ensemble technique (Panagiotopoulos, 1987) that allows the direct simulation of phase equilibrium.

Georgoulaki et al. (1994) used the MC method in the Gibbs ensemble to study the ability of conformal solutions theory to predict phase equilibria for binary LJ mixtures. In the Georgoulaki et al. work, simulation data for asymmetric mixtures were compared to theoretical results for cases in which unlike-pair interactions follow the Lorentz-Berhelot rule. Two empirical EoS for the pure LJ fluid were considered in the Georgoulaki et al. study, a Peng-Robinson-type cubic equation and a Benedict-Webb-Rubin-type equation. In the same study, the mixing rule proposed by Harismiadis et al. (1994) was used. Good agreement was found between theory and simulation. In another paper, Liu and Govind (1995) reported results for thermodynamic excess properties obtained from isothermal-isobaric molecular dynamics simulations of binary mixtures of LJ fluids for three parameter ratios, [e₂₂ / e₁₁ = 1 and s₂₂ / s₁₁ = 1.5], [e₂₂ / e₁₁ = 1.5 and s₂₂ / s₁₁ = 1], and [e₂₂ / e₁₁ = 1.5 and s₂₂ / s₁₁ = 1.5].

In this work, MC molecular simulations of pure component and binary mixtures of LJ fluids are used to study systems with the following size and energy parameter ratios: [e₂₂ / e₁₁ = 0.5 and s₂₂ / s₁₁ = 1.0], [e₂₂ / e₁₁ = 0.25 and s₂₂ / s₁₁ = 0.5], and [e₂₂ / e₁₁ = 0.66 and s₂₂ / s₁₁ = 0.35]. Our MC simulation results are validated by comparison with simulation data from the literature and then compared to predictions from the translated Peng-Robinson EoS (Georgoulaki et al., 1994) for LJ fluids. Furthermore, excess properties of liquids determined by MD simulation are compared to predictions of the t-PR-LJ. In both cases, the effect of the combining rules on the EoS predictions is investigated with different combining rules.

THEORY

The Gibbs Ensemble Technique

MC methods are reviewed by Panagiotopoulos(1996) and here we concentrate only on the MC method in the Gibbs ensemble (GMC) (Panagiotopoulos, 1987 and 1988). The basic assumption in the GMC is the coexistence of two phases in equilibrium placed separately in boxes with no physical contact. Some properties of the whole system (box I + box II) are kept constant. The most usual choices are the number of particles of each chemical species, the volume and the temperature (GMC in the constant NVT ensemble) or the number of particles of each chemical species, the pressure and the temperature (GMC in the constant NPT ensemble).

GMC uses three basic trial movements:

(1) a particle displacement inside each box to promote internal thermal equilibrium (like all MC simulations with other ensembles);

(2) a volume rearrangement in each phase to promote mechanical equilibrium (in the NVT version, total volume restriction should be observed);

(3) a particle interchange between the boxes to promote the equalization of chemical potentials.

Particle Displacement

In this step each box is a canonical ensemble (i.e., at constant NVT). In this ensemble, states occur with probability proportional to exp(–E / kT), where E is the configurational energy of the phase. A particle is chosen at random, for example, in box I, and randomly placed in a new trial position in box I. The probability ratio of the new and old states is

where DE^I is the difference in energy for the trial move, k is the Boltzmann constant and T is the absolute temperature of the system. The trial move is accepted with probability given by min(1, P^I_move). An identical procedure is applied in box II.

Volume Rearrangement

In this step, each box is a system from the isothermal-isobaric ensemble. Since the total volume of the system is not constant, changes in volume in the two boxes are independent. The probability for the volume rearrangement step is given by

where b = 1 / kT, P is the system pressure, N^I and N^II are the numbers of particles inside box I and box II, and V^I and V^II are the volumes of box I and box II, respectively. The trial move is accepted with probability given by min (1, P_vol). Equation 2 is identical to the normal constant NPT ensemble simulation acceptance criterion written for two independent regions.

In cases where the total volume of the system (box I + box II) is conserved, a trial move consists of choosing a random change in volume, DV, in box I and a simultaneous change, in volume, -DV, in box II. We assume that the pressure, P, in phase II is the same as that in phase I. In box II, the probability ratio of the new and old states is

Particle Interchange

In this case, each box is a system from the grand canonical ensemble (i.e., at constant mVT). A trial move consists, for example, of an attempted particle creation in box I and an attempted particle destruction in box II. The new position is chosen at random in box I and a particle is chosen at random in box II and destroyed. The probability for the overall interchange step is given by

The trial interchange is accepted with probability given by min(1, P_int).

In the case of mixtures, only the particle interchange step differs from the case of pure components. In equation 4, N^I and N^II are now the numbers of molecules of the species being interchanged.

The main limitation of GMC rests in the particle transfer step for dense fluids because the acceptance probability decreases rapidly as fluid density increases.

Lennard-Jones Fluids

In this work, the LJ potential was used:

where U_ij is the intermolecular potential energy between centers i and j, separated by distance r_ij, and e_ij and s_ij are the energy and size parameters for ij interactions, given by

Equations of State

The calculation of thermodynamic properties of a mixture using an EoS requires mixing rules to relate the EoS parameters of the mixture to those of its components. For a cubic two-parameter EoS, the most commonly used mixing rules are those of the so-called van der Waals one-fluid theory:

where a and b are the attractive and covolume parameters, respectively.

Harismiadis et al. (1991) disagree with the general belief that EoS can only predict the behavior of mixtures of molecules which are similar in size and energy parameters and that they can not predict liquid phase properties correctly. They state that the properties of mixtures of different particles can be predicted effectively with the use of appropriate combining rules. The most commonly employed combining rules for a_ij and b_ij use arithmetic and geometric means as follows:

In this work, the results of calculations using these combining rules were compared to those using the rules proposed by Sadus (1994). There is evidence that the Sadus rules provide qualitatively good results in the prediction of critical diagrams for systems with a large difference in molecular size. The Sadus rules are

where adjustable parameter h was assumed to be equal to 1 in all the calculations.

Harismiadis et al. (1994) also propose an EoS for LJ (t-PR-LJ), based on the translated Peng-Robinson equation presented by Magoulas and Tassios (1990) for hydrocarbons. The new equation (eq. 14) is a modification of the original PR EoS, including a parameter, t, that corrects the liquid phase volume.

The a and b parameters from the Magoulas and Tassios equation were obtained by regression of experimental data for alkanes of different chain lengths. In the Harismiadis et al. equation, parameters a and b are the same as those of Magoulas and Tassios considering the acentric factor of LJ fluids equal to zero due to the high degree of symmetry of the particles. The parameters are given by

Parameter t is given by

where = 0.3074 is the critical compressibility factor given by the original PR equation; , ., and are the critical properties of the LJ fluid; is the reduced temperature ; and a and b are parameters equal to 0.245029 and -25.6632, respectively.

The t-PR-LJ EoS provides good results for pure LJ fluids and has the advantage of computational simplicity. For this reason, it was used in this work, even though other equations such as that presented by Koutras et al. (1992), based on the Carnahan-Starling EoS, and that of the MWBR type with 33 adjustable parameters presented by Nicola et al. (1989), provide very accurate values for LJ fluids.

RESULTS AND DISCUSSION

Critical Properties of Pure Lennard-Jones Fluids

Harismiadis et al. (1994) present a review of the calculation of vapor pressures and critical properties of pure LJ fluids. This work presents simulations using MC methods and the Gibbs ensemble technique. The simulations were carried out with 500 particles and 2x10⁶ configurations were generated, the first 10⁶ of which were discarded to allow the system to reach equilibrium. Table 1 shows the results obtained in our work and several other values from the literature for the critical properties of LJ fluids. The critical temperature, pressure, and density of the pure LJ fluid were estimated by correlating the Gibbs ensemble MC simulation results for the vapor and liquid densities and the vapor pressure with the classical scaling law and the Antoine equation, respectively, following the procedure proposed by Tavares and Sandler (1996). The value between brackets represents the standard deviation of the last digits, i.e., 1.316(6) and 0.128(26) mean 1.316 ± 0.006 and 0.128 ± 0.026, respectively. When absent, it means that value was specified or not present in the original work.

In contrast to the consideration of Harismiadis et al. (1994) that the acentric factor of LJ fluids is equal to zero, simulation at T^* = 0.9212 (T_r = 0.7) led to w = -0.0094, which is within the range of values presented by Reid et al. (1987) for noble gases.

Binary Mixtures

To evaluate the accuracy, stability, and computational requirements of the FORTRAN code used in the simulations of binary mixtures, three systems studied by Georgoulaki et al. (1994) were recalculated in this work. Each simulation system contained 600 particles, and 2x10⁶ to 4x10⁶ configurations were generated, from which 1x10⁶ to 2x10⁶ were discarded to allow the system to reach equilibrium. In each particle change step, 500 insertion trials were done.

Most of the simulations were performed using the NPT version of the GMC method. However, for regions close to the critical point, the NVT version was used, as it proved to be more stable, as also stated by Georgoulaki et al. (1994). Table 2 shows a comparison between data obtained in this work and those of Georgoulaki et al. (1994). Pressures presented with no deviations indicate the use of the NPT version.

The simulations were extended for each of the three systems by keeping the size parameter ratio constant (1.00; 0.50; 0.35) and varying the energy parameter ratio from 1 to 4, which are values that can be observed in real systems, according to Georgoulaki et al. (1984) and Reid et al. (1987). To allow comparison of the properties obtained in the simulations to those predicted by the t-PR-LJ equation, a FORTRAN program developed by Gandur (1990) was used. The original code provided predictions for the PR and SRK EoS with arithmetic and geometric combining rules (equations 10 and 11). The Sadus combining rules (equations 12 and 13) and the-PR-LJ EoS were incorporated into Gandur`s program. As the purpose of this work is to study the ability of the EoS to predict equilibrium properties and not to fit simulation data, parameter h was always set equal to one in the combining rule defined by equation 12. In Table 2, column I shows the results this work and column II shows the results of Georgoulaki et al. (1994). When the pressure is not shown in column II it means that the value is the same as that established in column I.

Each equilibrium data point from MC simulations of system 1 (s₂₂ /s₁₁ = 1 and e₂₂ /e₁₁ = 0.50) converged in about 2 h in a 133 MHz Pentium. Figure 1 shows the P^* vs. x₁,y₁ and P^* vs. x₁,y₁ and P^* vs. r^* diagrams at T^* = kT / e₁₁ = 1.15, where P^* = Ps³₁₁ / e₁₁ and r = rs³₁₁. As the particles are all of the same size, the use of different combining rules makes no difference in the EoS predictions. As for densities, t-PR-LJ EoS shows good predictions for the vapor phase, but fails to predict the entire liquid phase.

In system 2 (s₂₂ /s₁₁ = 0.35 and e₂₂ /e₁₁ = 0.66), the combining rules were essential for accurate predictions. Arithmetic and geometric rules (equations 10 and 11) estimate a much smaller V-L envelope than that observed in the simulation data. The Sadus rules also underestimate the V-L region, even though they yield a quite good prediction of the critical point and very accurate property predictions for the liquid phase. Figure 2 shows the P^* vs. x₁,y₁ and P^* vs. r^* diagrams at T^* = 1.15.

For system 3 (s₂₂ /s₁₁ = 0.5 and e₂₂ /e₁₁ = 0.25), the mean computational time requirement was 8 h. Due to the large difference in size of the particles, fluctuations were higher in the simulation, specially near the critical point. Phase inversions in the MC simulations were observed several times and the number of equilibrium configurations was much higher than it was for system 1. Figure 3 shows the P^* vs. x₁,y₁ and P^* vs. r^* diagrams at T^*= 1.00, and the combining rules were again very important. The Sadus rules yield a less overestimated vapor-liquid region than the classic rules.

To investigate the predictions of excess properties with EoS, three systems with size and energy parameter ratios equal to [e₂₂ /e₁₁ = 1 and s₂₂ /s₁₁ = 1.5], [e₂₂ /e₁₁ = 1.5 and s₂₂ /s₁₁ = 1], and [e₂₂ /e₁₁ = 1.5 and s₂₂ /s₁₁ = 1.5] at T^*=2 and P^*=2.5 were studied. Figures 4 through 6 show simulation and EoS prediction results.

In the simulation of the system with e₂₂ /e₁₁ = 1.50 and s₂₂ /s₁₁ = 1.00(Fig. 4), where the molecules are all of the same size, g^E / RT and h^E are negative. The values obtained using the t-PR-LJ with the classic rules and the Sadus rules are very similar and represent equally well the values calculated by Liu and Govind (1995).

The excess properties evaluated for the system with equal atom energy (e₂₂ /e₁₁ = 1) using simulation data, g^E / RT, v^E, and h^E, are negative (Fig 5). The corresponding properties determined using the EoS with arithmetic and geometric rules are not in qualitative agreement with the values obtained using MD simulation data (Liu and Govind, 1995). Although the Sadus combining rules can not estimate with precision the behavior of the system, they at least provide results in qualitative agreement with the MD simulation data.

In the system with e₂₂ /e₁₁ = 1.50 and s₂₂ /s₁₁ = 1.50 (Fig. 6), the excess properties calculated by the EoS with the classic rules are shown to be inadequate to represent the system, while the Sadus rules provide a qualitatively correct representation of the system.

CONCLUSIONS

In this work, the ability of t-PR-LJ to predict vapor-liquid equilibrium and excess properties was studied. The effect of the combining rules was found to be very relevant. Arithmetic and geometric rules only give satisfactory predictions for systems with particles of the same size or with small differences in the energy parameter. For other cases, the Sadus rules always give better predictions than the arithmetic and geometric rules, even though these are still inadequate as quantitative predictions. However, it is important to note that all our calculations with EoS in this paper are predictive, i.e., no binary interaction parameters were fitted to the data. Therefore, the better qualitative behavior of the Sadus rules is a very promising aspect for their practical application, since they reduce the need for binary interaction parameters.

ACKNOWLEDGMENTS

The authors acknowledge the financial support provided by FAPERJ, CNPq/Brazil, and PRONEX (Grant no. 124/96).

NOMENCLATURE

a energy parameter in the cubic EoS b volume parameter in the cubic EoS E configurational energy G^E excess Gibbs free energy h^E excess enthalpy k Boltzmann’s constant N number of particles in each phase N_av Avogrado’s number P pressure

P^* reduced pressure, P^* = Ps³₁₁ / e₁₁ critical pressure of component i reduced critical pressure, r_ij interparticle distance R universal gas constant t volume translation parameter of the t-PR-LJ EoS T temperature T^* reduced temperature, critical temperature of component i reduced critical temperature, T_r reduced temperature, T_r = T / T_c U intermolecular potential V volume v molar volume v^E excess molar volume x_i,y_i mole fractions of component i critical compressibility factor given by the original PR equation

Greek Symbols

e_ij LJ energy potential parameter for interaction between components i and j r density r^* reduced density, r^* = rs³₁₁ reduced critical density of component i s_ij LJ size parameter for interaction between components i and j m chemical potential

• Gandur, M. C., Equilíbrio Químico e de Fases com Especificaçăo de Pressăo e Entalpia em Sistemas năo Ideais. MSc Thesis - UFRJ / COPPE, Rio de Janeiro (1990).
• Georgoulaki, A. M., Ntouros, I. V., Tassios, D. P., and Panagiotopoulos, A. Z., Phase Equilibria of Binary Lennard-Jones Mixtures: Simulation and van der Waals 1-fluid Theory. Fluid Phase Equilibria, 100: 153-170 (1994).
• Gubbins, K. E., Molecular Simulation of Fluid Phase Equilibria and Adsorption. Revue de lInstitue Français du Pétrole, 51 (1): 59-71 (1996).
• Harismiadis, V. I., Koutras, N. K., Tassios, D. P., and Panagiotopoulos, A. Z., How Good is Conformal Theory for Phase Equilibrium Predictions, Fluid Phase Equilibria, 65: 1-18 (1991).
• Harismiadis, V. I., Panagiotopoulos, A. Z., and Tassios, D. P., Phase Equilibria of Binary Lennard-Jones Mixtures with Cubic Equations of State, Fluid Phase Equilibria, 94: 1-18 (1994).
• Koutras, N. K., Harismiadis, V. I., and Tassios, D. P., A Simple Equation of State for the Lennard-Jones Fluid. Fluid Phase Equilibria, 77: 13-38 (1992).
• Liu, A. and Govind, R., Determination of Thermodynamic Excess Properties of Mixtures from Computer Simulation. Molecular Simulation, 15: 47-55 (1995).
• Magoulas, K. and Tassios, D. P., Thermophysical Properties of n-alkanes from C1 to C20 and Their Prediction for Higher Ones. Fluid Phase Equilibria, 56: 119-140 (1990).
• Nicolas, J. J., Gubbins, K. E., Street, W. B., and Tildesley, D. J., Equation of State for the Lennard-Jones Fluid. Molecular Physics, 37: 1429-1454 (1979).
• Panagiotopoulos, A. Z., Direct Determination of Phase Coexistence Properties of Fluids by Monte Carlo Simulation in a New Ensemble. Molecular Physics, 61 (4): 813-826 (1987).
• Panagiotopoulos, A. Z., Phase Equilibria in the Gibbs Ensemble: Alternative Derivation, Generalization and Application to Mixture and Membrane Equilibria. Molecular Physics, 63 (4): 527-545 (1988).
• Panagiotopoulos, A. Z., Current Advances in Monte Carlo Methods. Fluid Phase Equilibria, 116: 257-266 1996).
• Reid, R. C., Prausnitz, J. M., and Poling, B. E., The Properties of Gases and Liquids. McGraw-Hill editors, 4^th edition, NY (1987).
• Sadus, R. J., Calculating Critical Transitions of Fluid Mixtures: Theory vs. Experiment. AIChE Journal, 40 (8): 1376-1403 (1994).
• Smit, B., Phase-Diagrams of Lennard-Jones Fluids. J. Chem. Phys. 96: (11) 8639-8640 (1992).
• Tavares, F. W. and Sandler, S. I., Vapour-Liquid Equilibria of Exponential-Six Fluids. Molecular Physics, 87 (6): 1471-1476 (1996).

* To whom correspondence should be addressed

To whom correspondence should be addressed

Publication Dates

• Publication in this collection
  02 Aug 2001
• Date of issue
  June 2001

History

• Received
  09 Mar 2000
• Accepted
  15 May 2001