The reaction ensemble method for the computer

1 downloads 0 Views 164KB Size Report
May 1, 1999 - The reaction ensemble Monte Carlo (REMC) method W. R. Smith and B. Trıska, J. Chem. Phys. .... pointed out some errors in the simulations of Coker and. Watts.14 ..... These entry and exit compositions are denoted by the.
JOURNAL OF CHEMICAL PHYSICS

VOLUME 110, NUMBER 17

1 MAY 1999

The reaction ensemble method for the computer simulation of chemical and phase equilibria. II. The Br21Cl21BrCl system Martin Lı´sal E. Ha´la Laboratory of Thermodynamics, Institute of Chemical Process Fundamentals, Academy of Sciences, 165 02 Prague 6, Czech Republic

Ivo Nezbeda E. Ha´la Laboratory of Thermodynamics, Institute of Chemical Process Fundamentals, Academy of Sciences, 165 02 Prague 6, Czech Republic and Department of Physics, J. E. Purkyneˇ University, ´ stı´ n. Lab., Czech Republic 400 96 U

William R. Smith Department of Mathematics and Statistics, and School of Engineering, College of Physical and Engineering Science, University of Guelph, Guelph, Ontario N1G 2W1, Canada

~Received 5 October 1998; accepted 17 December 1998! The reaction ensemble Monte Carlo ~REMC! method @W. R. Smith and B. Trˇ´ıska, J. Chem. Phys. 100, 3019 ~1994!# is used to study combined reaction and vapor–liquid equilibrium of the Br21Cl21BrCl system. The substances are modeled as nonpolar and dipolar two-site Lennard-Jones molecules with Lorentz–Berthelot mixing rules for unlike atoms. No parameters were fitted to any mixture properties in our calculations. The simulated data are compared with experimental results, and with previous simulation data for the mixture obtained by an indirect semigrand ensemble approach. The REMC method efficiently calculates the complete phase compositions, whereas only a limited subset is available experimentally. The agreement of the simulations with experiment is good. In the course of this work, we used the Gibbs ensemble Monte Carlo method ~which may be regarded as a special case of the REMC method! to calculate the vapor–liquid equilibrium properties of pure BrCl; since this compound is chemically unstable, such data is experimentally inaccessible. © 1999 American Institute of Physics. @S0021-9606~99!50112-8# I. INTRODUCTION

These determine the configurational properties of the mixture, whether or not reaction occurs. One way to use such models to calculate reaction equilibrium is to conduct computer simulation experiments over the complete composition range of the mixture, and use empirical modeling techniques to fit this data to an appropriate EOS or free-energy model. Such a model can then be combined with the ideal-gas reaction terms in the overall system free-energy, and used in the macroscopic approach described above. The reaction ensemble Monte Carlo ~REMC! method6,7 obviates the need to describe the entire compositional behavior of the mixture, and focuses directly on determining the equilibrium composition. The REMC method is ideally suited for the ultimate development of accurate predictive techniques, since the empirical macroscopic-level modeling stage is eliminated. It combines in a single computer simulation the molecular-level model with the ideal-gas driving forces for chemical reaction. It may be applied to any number of simultaneous chemical reactions and simultaneous phases; by directly incorporating the reaction equilibrium conditions ~including those involving phase equilibrium! within the simulation, it avoids any direct calculations of chemical potentials or of fugacities. The non-necessity to calculate these quantities also characterizes the Gibbsensemble method8 for simulation of phase equilibrium ~which may be considered to be a special case of reaction equilibrium5!.

The theoretical prediction of reaction and phase equilibria is of great interest in the chemical industry.1 The combined occurrence of these phenomena underlies the industrially important process of reactive distillation.2 Using a macroscopic-level approach, predictions of reaction and phase equilibrium are customarily carried out using empirical equation-of-state ~EOS! or free-energy models, from which the mixture chemical potentials are determined as a function of the composition; these chemical potentials contain an ideal-gas term and a nonideality term. The former arises from the internal part of the mixture partition function, and the relevant data are available from tables ~e.g., the JANAF Tables3!, or may alternatively be calculated from ideal-gas partition function data.4 These quantities provide the driving forces to move the system composition to its equilibrium value. The nonideal term in the chemical potential may be viewed as arising from the configurational part of the partition function. Given a macroscopic-level chemical potential model for the mixture as a function of composition, the appropriate thermodynamic function is numerically minimized under the overall thermodynamic constraints to determine the equilibrium composition; for example, at specified ~T,P!, the system Gibbs free-energy is minimized.5 Using a microscopic-level approach, a mixture is described by means of a set of intermolecular potential models. 0021-9606/99/110(17)/8597/8/$15.00

8597

© 1999 American Institute of Physics

8598

Lı´sal, Nezbeda, and Smith

J. Chem. Phys., Vol. 110, No. 17, 1 May 1999

The REMC method has been applied to simple singlephase model systems;6,7,9 a preliminary strategy was presented for using it to model methyl tert-butyl ether ~MTBE! synthesis,10 and it has been used to study model systems at interfaces.11 In this paper, we demonstrate the utility of the method to predict the detailed quantitative thermodynamic behavior of a real system for which a limited set of experimental results is available. In addition to predicting experimentally accessible quantities, we also use the method to predict important, but experimentally inaccessible, quantities. The Br21Cl21BrCl system has been studied experimentally by Karsten12 and by Cheesman and Scott.13 Their results consist of a limited set of vapor–liquid equilibrium data, using the overall Br2 and Cl2 mole fractions as composition variables. Coker and Watts14 simulated this system ~in the liquid phase only! using a modified grand canonical ensemble Monte Carlo method. The combined reaction and phase equilibrium of this system was subsequently calculated using an indirect approach by Kofke and Glandt,15 who employed a semigrand ensemble Monte Carlo method requiring the calculation of a single-species chemical potential and the numerical solution of an equilibrium condition; they also pointed out some errors in the simulations of Coker and Watts.14 Kofke and Glandt15 compared some of their results with the experimental results of Cheesman and Scott.13 Both the experimental and previous theoretical approaches have considered only a limited subset of the reaction and phase equilibrium behavior of the system. The approach of Kofke and Glandt15 makes the calculation of the complete reaction and phase equilibrium diagram a somewhat time-consuming undertaking. ~We remark also that their method cannot be used for cases in which the total number of moles changes due to reaction.! The REMC technique has none of these experimental nor theoretical disadvantages, and can readily calculate complete compositional data for all phases. The purpose of this paper is to employ the REMC method to calculate the equilibrium compositions of the reacting Br21Cl21BrCl system, including complete phase compositional data. In the course of this work, we also calculate the vapor–liquid equilibrium behavior of pure BrCl using Gibbs ensemble Monte Carlo ~GEMC!8 simulations. In addition to being of intrinsic interest, this data is experimentally inaccessible, since BrCl is chemically unstable. In the next section, we describe the intermolecular potential models used. In the subsequent section, we discuss the computational details underlying our approach. The following section contains results and discussion, and the final section contains conclusions. II. INTERMOLECULAR POTENTIAL MODELS

The intermolecular potentials used in our simulations are extensions of the potential models of Singer et al.16 for Br2 and Cl2 to the ternary mixture of these species with BrCl. The interactions among the species are represented by twosite atom–atom potentials. The interaction between atom a and b in different molecules is described by the LennardJones potential

u LJ~ r ab ! 54« ab

FS D S D G s ab r ab

12

2

s ab r ab

6

~2.1!

,

with the potential parameters for Br2 and Cl2 taken from Singer et al.,16 using their constraint of constant bond length l. In Eq. ~2.1!, « ab and s ab are the Lennard-Jones energy and length parameters, respectively, between atoms a and b in different molecules, and r ab is the distance between these atoms. We assumed that the Lennard-Jones potential parameters between unlike atoms satisfy the Lorentz–Berthelot mixing rules « ab 5 A« aa « bb

s ab 5

s aa 1 s bb . 2

~2.2!

Since no potential model for BrCl is currently available, we also assumed it to be described by a two-site LennardJones potential, and as in earlier theoretical studies,14,15 assumed that the Lennard-Jones potential parameters are transferable from the atoms of Br2 and Cl2 to the atoms of BrCl, and that the bond length of BrCl is the arithmetic mean of the Br2 and Cl2 bondlengths. In addition, a dipole–dipole interaction u mm ij 5

H

1 m2 3 3 ~ ei •e j ! 2 2 @~ ei •ri j !~ e j •ri j !# 4pe0 rij rij

J

~2.3!

between BrCl molecules was also included.14 In Eq. ~2.3!, m is the permanent dipole moment of the BrCl molecules, e 0 is the permittivity of free space, ri j is the center-of-mass distance vector between BrCl molecules i and j, r i j 5 u ri j u , and ei and ej are unit vectors along the molecular axes of BrCl molecules labeled i and j, respectively. Values of the potential parameters are given in Table I; we also include for convenience some general properties of Br2, Cl2, and BrCl. In contrast to the potential models used by Coker and Watts,14 and by Kofke and Glandt,15 we did not include an intramolecular contribution due to a flexible bond length. This contribution was not taken into account in the parametrization of the Br2 and Cl2 potentials by Singer et al.,16 and therefore, its influence on the adequacy of the Singer et al. potentials is unknown. Further, we also neglected inductive interactions for BrCl with any other molecule, since its contribution to the intermolecular potential is, on the one hand, less than 1%,15 and on the other hand, its evaluation is very time consuming.17 III. COMPUTATIONAL DETAILS

For the pure BrCl GEMC simulations, we used N 5512 particles in cubic boxes, the minimum image convention, periodic boundary conditions, and cutoff radius equal to half the box length. The Lennard-Jones long-range correction for the configurational energy was included,18 assuming that the radial distribution functions are unity beyond the cutoff radius. The dipolar long-range correction of the configurational energy was treated by the reaction-field method with the dielectric constant e RF set to infinity.19 The GEMC simulations were organized in cycles as follows. Each cycle consisted of three steps: n D translation and rotational moves, n V volume moves, and n T particle transfers. The three types

Lı´sal, Nezbeda, and Smith

J. Chem. Phys., Vol. 110, No. 17, 1 May 1999

8599

TABLE I. Potential parameters for the atoms Br and Cl,16 and general properties of Br2, Cl2, and BrCl. e and s are the Lennard-Jones energy and length potential parameters, respectively, k B is Boltzmann’s constant, l is the bond length, m is the permanent dipole moment, M is the molecular weight, T c is the critical temperature, P c is the critical pressure, and r c is the critical molar density. For Br2, the critical parameters were taken from Lide, ~Ref. 27! and for Cl2 from Angus et al. ~Ref. 28! For BrCl, the critical parameters were estimated from the Gibbs ensemble simulation data of Table II.

Br Cl Br2 Cl2 BrCl

e /k B (K)

s ~Å!

257.2 178.3 l ~Å!

3.538 3.332 m ~C•m!a

M 3103 (kg•mol21)

T c (K)

r c (mol•m23)

P c ~kPa!

¯ ¯ 1.87•10230

159.808 70.9054 115.3567

588 416.956 518.6

7874.0 8134.5 7252.3

10340.0 7991.4 9319.2

2.229 2.099 2.164

a

The dipole of BrCl is located a distance 0.3073l from the Br atom along the interatomic axis.

of moves were selected at random with fixed probabilities, chosen so that the appropriate ratios of each type of move were obtained. The ratio n D :n V :n T in the cycle was set to N:1:1000. Acceptance ratios for translation and rotational moves, and for volume changes, were adjusted to approximately 30%. After an initial equilibration period of 122 3104 cycles, we generated ~depending on the thermodynamic conditions! between 0.5213105 cycles to accumulate averages of the desired quantities. The precisions of the simulated data were obtained using block averages, with 500 cycles per block. In additional to ensemble averages of the quantities of direct interest, we also monitored convergence profiles of the thermodynamic quantities, in order to keep the development of the system under control.20 For the Br21Cl21BrCl system, reaction and phase equilibrium at fixed ~T,P! requires equilibrium for the reactions Br2~ l ! 5Br2~ g ! ,

~3.1!

Cl2~ l ! 5Cl2~ g ! ,

~3.2!

BrCl~ l ! 5BrCl~ g ! ,

~3.3!

Br21Cl252BrCl,

~3.4!

where l and g denote liquid and vapor phases, respectively ~in the last reaction, each species may be in either phase!. We remark that we have used the convention that a chemical species is characterized by both its chemical formula and its phase, enabling the treatment of phase equilibria on exactly the same basis as reaction equilibria.5 There are thus four independent stoichiometric reactions for this system, and the above form a linearly independent set. The vapor and liquid phases are represented in the simulations by two separate boxes, labelled I and II. An arbitrary simulation box in the following will be denoted by the symbol a. In the REMC method, the required conditions of vapor–liquid and chemical-reaction equilibrium are ensured by performing a combination of four steps: particle displacements, volume changes, interphase particle transfers, and reaction moves. The transition probability k→l for a particle displacement in box a is18 a PD kl 5min@ 1,exp~ 2 b DU kl !# ,

~3.5!

and the transition probability k→l for a volume change of box a is18 P Vkl 5min

H F 1,exp

2 b DU akl 2 b P ~ V al 2V ak ! 1N a

ln

V al V ka

GJ

. ~3.6!

a 5U al 2U ak is the change in In Eqs. ~3.5! and ~3.6!, DU kl configurational energy in box a, b 51/(k BT), k B is Boltzmann’s constant, V a is the volume of box a, and N a is the total number of molecules in box a. Particle transfers between phase boxes were implemented using the algorithm of Green et al.,21 which is computationally efficient when the compositions of the coexisting phases are quite different.21,22 It involves choosing the donor and recipient boxes at random, then randomly choosing the particle which is to be transferred from the donor box regardless of its type, and then transferring it to a random position in the recipient box. The transition probability k →l for the transfer of a particle from box I into box II is

H F

I II P I→II kl 5min 1,exp 2 b DU kl 2 b DU kl 1ln

N IV II ~ N II11 ! V I

GJ

. ~3.7!

The ideal-gas driving term for a reaction j is given by6

F

G j [exp 2

DG 0j ~ T, P 0 ! RT

G

,

~3.8!

where R is the universal gas constant, and DG 0j (T, P 0 ) is the standard Gibbs energy change for reaction j at the system temperature T and the standard-state pressure P 0 , given by n

DG 0j ~ T ! [

( n i j DG 0f i~ T ! ,

i51

~3.9!

where the dependence of DG 0 on P 0 has been suppressed, n is the number of species in the system, n i j is the stoichiometric coefficient of species i in reaction j, and DG 0f i (T) is the molar standard ideal-gas Gibbs energy of formation of species i at the system temperature T and standard-state pressure P 0 . G is zero for a phase equilibrium reaction ~reactions ~3.1!–~3.3! above!, and nonzero otherwise. For reaction ~3.4!,

8600

Lı´sal, Nezbeda, and Smith

J. Chem. Phys., Vol. 110, No. 17, 1 May 1999

TABLE II. Vapor–liquid equilibrium data for pure BrCl from the Gibbs ensemble computer simulations of this paper. T is the temperature, U is the internal energy, r is the molar density, and P s is the vapor pressure. Superscripts g and l denote the vapor and liquid phases, respectively. The simulation uncertainties are given in the last digits as subscripts. T ~K!

U g (kJ•mol21)

U l (kJ•mol21)

r g (mol•m23)

r l (mol•m23)

P s ~kPa!

488.68 462.96 385.8 334.36 283.0 257.2

23.7461 22.3537 20.7113 20.259 20.0810 20.0813

213.83 215.32 219.02 221.11 223.11 224.21

2806.54604 1675.72736 451.4640 151.8381 38.6117 16.284

12 620.93410 13 925.71472 16 824.11162 18 338.9983 19 673.7864 20 333.7668

6543.84595 4519.12391 1271.71536 399.2949 88.4264 34.1178

DG 0 ~ T ! 52DG 0f ,BrCl~ T ! 2DG 0f ,Br ~ T ! 2DG 0f ,Cl ~ T ! . 2 2

~3.10!

Reaction moves for reaction ~3.4! need be carried out in only one phase box for the system at hand, since Eqs. ~3.1!– ~3.4! form a maximal linearly independent set of reactions. For convenience, we performed these moves in the liquid box, where more particles are presented; however, other strategies are possible which will affect the convergence ~but not the final results! of the method. Reaction moves were carried out in forward and reverse directions according to a preset probability of 0.5. In the forward direction, a Br2 and a Cl2 molecule ~reactants! are chosen at random, and an attempt is made to change their identity into two molecules of BrCl ~products!. The attempt is accepted with transition probability k→l given by6 P FR kl 5G

a a N Cl N Br 2

2

a a 11 !~ N BrCl 12 ! ~ N BrCl

exp~ 2 b DU akl ! .

~3.11!

Similarly, in the reverse direction, two molecules of BrCl ~products! are chosen at random, and an attempt is made to change their identity into Br2 and Cl2 molecules ~reactants!. This attempt is accepted with transition probability k→l given by6 P RR kl 5

a a 21 ! N BrCl 1 ~ N BrCl exp~ 2 b DU akl ! . a a G ~ N Br 11 !~ N Cl 11 ! 2 2

~3.12!

a a a In Eqs. ~3.11! and ~3.12!, N Br , N Cl , and N BrCl are the 2 2 numbers of molecules of Br2, Cl2, and BrCl in box a, respectively. As for the pure BrCl GEMC simulations, the REMC simulations were similarly organized in cycles as follows. Each cycle consisted of four steps: n D translation and rotational moves, n V volume moves, n T particle transfers, and n j forward and reverse reaction moves. The four types of moves were selected at random with fixed probabilities, chosen so that the appropriate ratios of each type of move were obtained. The ratio n D :n V :n T :n j in the cycle was set to N:2:N:5000. The other simulation details were the same as in the case of the GEMC simulations. The initial equilibration period required 0.5213105 cycles and final production runs were carried out using between 1223105 cycles.

IV. RESULTS AND DISCUSSION

The vapor–liquid equilibria ~VLE! of pure LennardJones diatomics corresponding to Br2 and Cl2 have been

studied by Romano and Singer23 and for Cl2 by Galassi and Tildesley.24 Both sets of workers deemed the intermolecular potentials to provide a good representation of the experimental data for these compounds. Our new GEMC simulation data for pure BrCl are given in Table II. We estimated the critical temperature T c and density r c of BrCl from a leastsquares fit of the rectilinear diameter law25

r g1 r l 5 r c 1C 1 ~ T2T c ! , 2

~4.1!

and the critical scaling relation25

r l 2 r g 5C 2 ~ T c 2T ! 1/3.

~4.2!

The estimated critical temperature and density are given in Table I. We have fitted the vapor-pressure curves for all pure components to the Antoine equation,26 log10 P s 5A2

B . T1C

~4.3!

The Antoine constants for all pure fluids are given in Table III. The critical pressure, P c , of BrCl was obtained from the Antoine equation and is also given in Table I. The vapor– liquid coexistence curves for all three pure fluids are shown in Fig. 1. Karsten,12 and Cheesman and Scott13 experimentally examined VLE in Br21Cl2 mixtures by means of total vaporpressure measurements. They studied the VLE behavior as a function of the overall atomic compositions, z Br2 , and z Cl2 , defined for each phase by

TABLE III. Constants A, B, and C for the Antoine vapor-pressure equation @Eq. ~4.3!# for Br2, Cl2, and BrCl. Pressure is in kPa and temperature is in K. For Br2, we used the data of Romano and Singer ~Ref. 23!, for Cl2 that of Galassi and Tildesley ~Ref. 24!, and for BrCl, our GEMC simulation results of Table II.

Br2 Cl2 BrCl

A

B

C

6.387 19 6.721 73 6.780 62

1347.206 93 1268.379 69 1582.753 54

223.324 05 32.018 79 44.416 00

Lı´sal, Nezbeda, and Smith

J. Chem. Phys., Vol. 110, No. 17, 1 May 1999

FIG. 1. Vapor–liquid coexistence curves for pure Br2, Cl2, and BrCl. Open circles are the simulation results of Romano and Singer ~Ref. 23! for Br2 and filled circles are the Br2 experimental data ~Ref. 23!. Open triangles denote GEMC data of Galassi and Tildesley ~Ref. 24! for Cl2, and filled triangles denote the Cl2 experimental data ~Ref. 28!. Diamonds represent GEMC results of this paper for BrCl. ~a! Temperature-density coexistence envelopes; the lines are drawn through simulation data only as a guide. ~b! Vapor-pressure curves; the curves through simulation data represent fits of the data to the Antoine equation @Eq. ~4.3!#.

z Br25x Br21 21 x BrCl ,

~4.4!

z Cl25x Cl21 21 x BrCl512z Br2 ,

~4.5!

where x denotes the usual species mole fraction. They studied z Br2 and z Cl2 in each phase as a function of the overall T T and z Cl , defined by atomic compositions in the system, z Br 2

T z Br 5z Br2~ g ! f 1z Br2~ l 2

!~ 12 f ! ,

T T z Cl 512z Br , 2

2

~4.6! ~4.7!

2

where f is the molar fraction of the system in the vapor phase. We applied the REMC method6 for the system Br21Cl21BrCl at P51 atm at the temperatures $263.15, 283.15, 293.15, and 303.15 K%. We used G values interpo-

FIG. 2. Typical convergence profiles at T5263.15 K in the production period of ~a! coexistence densities, and ~b! coexistence mole fractions of Cl2.

lated from molar Gibbs free-energy data available in standard stables.3 Our simulation values of the thermodynamic quantities after the equilibration period did not show any systematic trends away from the average values. Convergence profiles such as the ones shown in Fig. 2 were monitored to verify that convergence was satisfactory for each production period. The final results are summarized in Table IV. In a nonreacting binary system, the atomic compositions coincide with the mole fractions of each component. In a ternary nonreacting system, there are two independent composition variables, and phase equilibrium is described by a 4-dimensional ( P,T,x 1 ,x 2 ) surface, where x 1 and x 2 are any two composition variables in a phase; at fixed P or fixed T, phase equilibrium is described by three-dimensional dewand bubble-point surfaces. The presence of an equilibrium chemical reaction in a ternary system reduces the number of degrees of freedom by one, resulting in a corresponding dimensional reduction in the dependence of the phase equilibrium behavior, producing dew- and bubble-point curves on the aforementioned surfaces. The overall atomic composi-

TABLE IV. Vapor–liquid equilibrium data of the reacting system Br21Cl21BrCl from reaction ensemble Monte Carlo simulations at P51 atm. T is the temperature, G is given by Eqs. ~3.8! and ~3.10!, using thermochemical data interpolated from the JANAF tables ~Ref 3!, U is the internal energy, r is the molar density, x Br2 , x Cl2 , and x BrCl are the Br2, Cl2, and BrCl mole fractions, respectively, z Br2 is the overall atomic composition of Br2, f is the molar fraction of the system in a phase, and g and l denote the vapor and liquid phases, respectively. The simulation uncertainties are given in the last digits as subscripts. U (kJ•mol21)

r (mol•m23)

8601

x Br2

x Cl2

x BrCl

z Br2

f

T5263.15 K, G58.40 47.029 g 20.0758 l 222.476235 20 475.688

0.007 929389 0.114 1 109

0.8121200 0.3172197

0.1800197 0.568797

0.097 93138 0.398 5 157

0.4476319

T5283.15 K, G57.97 43.733 g 20.08113 l 224.272245 20 232.488

0.039 44136 0.326 8 219

0.5475226 0.1144118

0.4131259 0.5588109

0.2460265 0.6062273

0.4187246

T5293.15 K, G57.79 42.916 g 20.08038 l 225.085186 19 868.7106

0.095 2587 0.470 1 155

0.3983128 0.058 1268

0.5064132 0.4718141

0.3485153 0.7060226

0.3538172

T5303.15 K, G57.62 40.521 g 20.08170 l 225.812173 19 546.3104

0.166 7 105 0.618 0 72

0.253787 0.025 90 28

0.5796126 0.356176

0.4565168 0.7960110

0.4038267

8602

J. Chem. Phys., Vol. 110, No. 17, 1 May 1999

FIG. 3. Temperature-atomic composition diagram of the reactive bromine– chlorine mixture. The dew- and bubble-point lines drawn through the experimental data of Cheesman and Scott ~Ref. 13! are represented by dashed and solid curves, respectively. Filled circles denote the experimental data of Karsten ~Ref. 12!, squares denote the simulation results of pure Br2 and Cl2 ~Refs. 23, 24!, triangles denote the Kofke and Glandt simulation results ~Ref. 15!, and diamonds denote the simulation results of this paper.

tions of each phase are transformations of the molar compositions on these curves, but they do not in themselves provide a complete description of the phase behavior. In Fig. 3, we show the vapor–liquid phase envelope at P51 atm, which shows the coexistence values of z Br2(l) and z Br2(g) as a function of temperature. Although this data does not describe the complete compositional dependence of the coexisting phases, it is the most convenient data obtainable experimentally. Our results are compared with the experimental data of Cheesman and Scott ~CS!13 and with the simulation results of Kofke and Glandt ~KG!.15 The experimental results of Karsten12 ~obtained in 1907! are in acceptable agreement with the CS results on the liquid side, but ~as concluded by Cheesman and Scott13! they are in error on the vapor side; we have included them only for completeness. On the vapor side, the overall agreement among the CS, KG, and our results is only fair. At T5263.15 K, our dew-point composition is in excellent agreement with that of CS. At higher temperatures, our dew-point compositions underestimate those of CS by about 10%. Further, our simulations did not reproduce the experimentally observed slight indentation of the bubble-point curve near the equimolar composition. On the liquid side, the overall agreement among CS, KG, and our results is excellent. Our bubble-point compositions underestimate those of CS by about 3%. In Fig. 4, we show the projection in the ternary composition triangle of the reaction and phase equilibrium diagram. Shown are the vapor–liquid coexistence curve, several equilibrium tie-lines, and a typical complete reaction equilibrium curve ~the latter at T5283.15 K!. z Br2 is a monotonically increasing parameter along the length of this latter curve. For a given temperature between the saturation temperatures of

Lı´sal, Nezbeda, and Smith

FIG. 4. Projection in the ternary composition triangle of the complete reaction and phase equilibrium diagram of the reactive bromine–chlorine mixture. Filled symbols denote the experimental vapor-phase compositions of Cheesman and Scott ~Ref. 13!. Symbols with a cross hair denote the liquid and vapor saturation simulation data of Kofke and Glandt ~Ref. 15!. Open symbols connected by an equilibrium tie-line denote the results of this paper. The vapor–liquid coexistence curve drawn through the simulation data of this paper is represented by the solid curve. The dotted curves are the single-phase portions of the reaction equilibrium curve at T5283.15 K, calculated in this paper.

pure Cl2 and Br2, the reaction equilibrium curve has singlephase regions, which intersect the nonreacting ternary phase equilibrium diagram in a point; the locus of such points as temperature varies forms the phase envelope shown. ~Outside this T range, reaction equilibrium occurs only in a single phase.! The curve at T5283.15 K begins in the gas phase at low z Br2 values, enters the saturation region, and then exits the saturation region in the liquid phase at a higher z Br2 value. These entry and exit compositions are denoted by the open circles in the figure. We also show in Fig. 4 the saturation liquid and vapor-phase compositions of Kofke and Glandt15 and the corresponding vapor-phase compositions of Cheesman and Scott.13 The liquid-phase results of the former authors were given in their paper, and we numerically calculated the vapor-phase results of both sets of workers assuming the ideal-gas EOS and using a single-phase chemical equilibrium calculation in both cases, using the stated z Br2(g) values. Although Kofke and Glandt15 used the virial EOS with the second virial coefficient, we found that the ideal-gas EOS results typically agreed with our REMC vapor-phase compositions to within 2%. The overall agreement among the CS, KG, and our results is very good. Our vapor-phase compositions are in better agreement with the CS results than those of KG. Our liquid-phase compositions agree with the KG liquid-phase compositions to within about 3%. In Fig. 5, we show the detailed dependence of the complete liquid and vapor compositions on the overall system composition at T5283.15 and 303.15 K. Also shown are the vapor-phase curves ~again obtained from equilibrium calculations assuming the ideal-gas EOS and using the stated z Br2(g) values!, the CS saturation vapor-phase compositions,

Lı´sal, Nezbeda, and Smith

J. Chem. Phys., Vol. 110, No. 17, 1 May 1999

FIG. 5. Dependence of the liquid and vapor coexistence compositions of the reactive bromine–chlorine mixture on the overall system composition at ~a! T5283.15 K and ~b! T5303.15 K. Filled circles denote the experimental vapor-phase compositions of Cheesman and Scott ~Ref. 13!. Diamonds denote vapor-phase compositions from the simulations of Kofke and Glandt ~Ref. 15!. Circles denote the results of this paper; corresponding phaseboundary points are connected by equilibrium tie-lines. Curves on the vapor side correspond to single-phase chemical equilibrium calculations assuming the ideal-gas EOS at the given overall composition; curves on the liquid side are drawn through the simulation results of Kofke and Glandt ~Ref. 15!.

and the KG liquid-phase compositions. On the vapor side, our simulation results agree within their statistical uncertainties with those obtained from the equilibrium calculations assuming the ideal-gas EOS. On the liquid side, the agreement between our results and the KG results is excellent, except in the vicinity of the liquid boundary at T 5283.15 K, where the agreement is slightly worse. Overall, our calculations differ slightly from the KG results. We believe that this is due primarily to the use of different values of the thermochemical data parameters G by KG. Our G values ~Table IV! were interpolated from data available in the JANAF thermochemical tables,3 a widely accepted source of accurate data; the KG values: $ (G,T) % 5 $ (10.69,263.15),(10.01,283.15),(9.47,303.15) % , which they calculated from partition function data, are approximately 20% larger than these values. Further evidence that the KG values are in error is the fact that the approximate ~constant! value of G quoted by Cheesman and Scott13 (1./0.382 56.93) is closer to the JANAF values than to the KG values. Finally, to estimate the effect of G values on our calculations, at T5263.15 K, we have found numerically that a 5% increase in G results in approximately a 1% increase in the atomic compositions of each phase. This is consistent with the observed differences between our results and the KG results. We conjecture that similar observations hold at the other temperatures, and that the differences in G account for the bulk of the differences between the KG results and ours.

V. CONCLUSIONS

We have used the reaction ensemble Monte Carlo ~REMC! method to directly simulate combined reaction and phase equilibrium for the Br21Cl21BrCl system over a range of temperatures and overall system compositions. We compared these calculations with the available experimental

8603

results12,13 and with those previously calculated by Kofke and Glandt15 using an indirect semigrand ensemble approach. The agreement of our results with the experimental results is good. We attribute differences between ours and the previously calculated theoretical results15 largely to inaccurate values used by the latter for the standard reaction freeenergy change data. We have shown that the REMC method can be used to calculate data that is not available experimentally, including the complete reaction and phase equilibrium diagram for the system. The Br21Cl21BrCl system studied here is fairly simple at the molecular level, and it is likely that results of comparable accuracy to those obtained here by the REMC method could also be obtained using an appropriate engineeringbased equation of state ~EOS!. However, the REMC method can be applied, in principle, to much more complex systems ~including molecularly confined systems!, for which a reasonable EOS is not available. In addition, the REMC method permits the direct testing of the adequacy of a given intermolecular potential model to predict observed thermodynamic properties, and obviates the need to construct an explicit EOS or free-energy model of the system.

ACKNOWLEDGMENTS

This research was supported by the Grant Agency of the Czech Republic under Grant No. 203/98/1446, and by the Natural Sciences and Engineering Research Council of Canada under Grant. No. OGP1041.

W. D. Seider and S. Widagdo, Fluid Phase Equilibria 123, 283 ~1996!. D. Barbosa and M. F. Doherty, Chem. Eng. Sci. 43, 529 ~1988!; S. Ung and M. F. Doherty, ibid. 50, 3201 ~1995!. 3 M. W. Chase, Jr., NIST-JANAF Thermochemical Tables, 4th ed., J. Phys. Chem. Ref. Data Monogr. 9, ~1985!. 4 D. A. McQuarrie, Statistical Mechanics ~Harper & Row, New York, 1976!. 5 W. R. Smith and R. W. Missen, Chemical Reaction Equilibrium Analysis: Theory and Algorithms ~Wiley-Interscience, New York, 1982; reprinted, with corrections, Krieger, Malabar, FL, 1991!. 6 W. R. Smith and B. Trˇ´ıska, J. Chem. Phys. 100, 3019 ~1994!. 7 J. K. Johnson, A. Z. Panagiotopoulos, and K. E. Gubbins, Mol. Phys. 81, 717 ~1994!. 8 A. Z. Panagiotopoulos, Mol. Simul. 9, 1 ~1992!. 9 W. R. Smith, I. Nezbeda, M. Strnad, B. Trˇ´ıska, S. Labı´k, and A. Malijevsky´, J. Chem. Phys. 109, 1052 ~1998!. 10 W. R. Smith, I. Nezbeda, and M. Strnad, paper 13a, AIChE Annual Meeting, Chicago, IL, 1996. 11 M. Borowko, A. Patrikiejew, S. Sokolowski, R. Zagorski, O. Pizio, and I. N. Abramenko, Czech. J. Phys. 48, 371 ~1998!. 12 B. J. Karsten, Z. Anorg. Chem. 53, 365 ~1907!. 13 G. H. Cheesman and D. L. Scott, Aust. J. Chem. 21, 287 ~1968!. 14 D. F. Coker and R. O. Watts, Chem. Phys. Lett. 78, 333 ~1981!; Mol. Phys. 44, 1303 ~1981!. 15 D. A. Kofke and E. D. Glandt, Mol. Phys. 64, 1105 ~1988!. 16 K. Singer, A. Taylor, and J. V. L. Singer, Mol. Phys. 33, 1757 ~1977!. 17 B. J. Costa Cabral, J. L. Rivail, and B. Bigot, J. Chem. Phys. 86, 1467 ~1987!. 18 M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids ~Clarendon, Oxford, 1987!. 19 B. Saager, J. Fischer, and M. Neumann, Mol. Simul. 6, 27 ~1991!. 20 I. Nezbeda and J. Kolafa, Mol. Simul. 14, 153 ~1995!. 1 2

8604 21

Lı´sal, Nezbeda, and Smith

J. Chem. Phys., Vol. 110, No. 17, 1 May 1999

D. G. Green, G. Jackson, E. de Miguel, and L. F. Rull, J. Chem. Phys. 101, 3190 ~1994!. 22 L. F. Rull, G. Jackson, and B. Smit, Mol. Phys. 85, 435 ~1995!. 23 S. Romano and K. Singer, Mol. Phys. 37, 1765 ~1979!. 24 G. Galassi and D. J. Tildesley, Mol. Simul. 13, 11 ~1994!. 25 J. V. Sengers and J. M. H. Levelt Sengers, Progress in Liquid Physics, edited by C. A. Croxton ~Wiley, New York, 1978!, pp. 103–174.

26

R. C. Reid, J. M. Prausnitz, and B. E. Poling, The Properties of Gases & Liquids, 4th ed. ~McGraw-Hill, New York, 1987!. 27 D. R. Lide, editor, CRC Handbook of Chemistry and Physics ~CRC, Boca Raton, FL, 1992!. 28 Chlorine, International Thermodynamic Tables of the Fluid State, edited by S. Angus, B. Armstrong, and K. M. De Reuch, Chemical Data Series, Vol. 8, No. 31 ~Pergamon, New York, 1984!.