Ind. Eng. Chem. Res. 1995,34, 3573-3580
3573
Metastability and Constraints: A Study of the Superheated Lennard-Jones Liquid in the Void-Constrained Ensemble David S. Corti and Pablo G. Debenedetti* Department of Chemical Engineering, Princeton University, Princeton, New Jersey 08544-5263
The rigorous study of metastability requires knowledge of the behavior of systems under various constraints. Such constraints, though impossible to impose analytically for realistic systems, are easily applied in computer simulations. We study the superheated Lennard-Jones liquid in the void-constrained ensemble, in which limits are imposed on the maximum size of voids that are allowed to form. The equation of state is sensitive to the severity of the void constraint, the more so the higher the temperature and the larger the degree of superheating. In contrast, the equation of state is insensitive to bounds placed on the magnitude of allowed density fluctuations in which both voids and clusters are restricted. Constraints that truncate the void size distribution, eliminating the large-void tail, lead to artificially high tensions. Over the range of conditions investigated here, this occurs when the size of the maximum allowed void does not exceed one average interparticle separation. 1. Introduction If a subcritical liquid is held at a pressure lower than its vapor pressure at the given temperature or is heated to a temperature higher than its boiling temperature at the given pressure, it is said to be superheated. Likewise, if a subcritical vapor is compressed to a pressure higher than its condensation pressure at the given temperature or cooled to a temperature lower than its dew point temperature at the given presure, it is said to be supercooled. In the above examples, the homogeneous fluid is metastable. It does not have the lowest Gibbs energy a t the given temperature and pressure. Metastable systems have finite lifetimes. Over times shorter than their lifetimes, they remain homogeneous, and no nuclei of the more stable phase exceeding a critical cut-off size are formed. The system is thus kinetically constrained to sample a restricted region of its phase space. From a statistical mechanical viewpoint, therefore, a metastable system is one that is constrained to explore a limited region of its phase space (Reiss, 1975). Computer simulations are ideally suited to study systems under arbritary constraints. Hence they are useful for the investigation of metastability. Recently, we have studied the superheated Lennard-Jones liquid in the restricted ensemble (Corti and Debenedetti, 1994). In a restricted ensemble simulation, arbitrary bounds are imposed on the magnitude of allowed density fluctuations. In that work, we also presented preliminary results on metastability in the void-constrained ensemble. In the latter ensemble, a limit is placed on the size of voids allowed to form in the superheated liquid. In this paper we investigate the behavior of the superheated Lennard-Jones liquid in the void-constrained ensemble, and we compare these results with those previously obtained in the restricted ensemble. In addition, we propose a criterion for determining the critical severity of a constraint, beyond which the properties of the superheated liquid become unphysical. The long-term goal of this work is to understand what are the appropriate constraints for studying metastable fluids; the approach we have chosen is t o investigate the effects of the type and severity of the constraints
* Author to whom correspondence should be addressed. 0888-5885/95/2634-3573$09.QQI0
on the equation of state of the metastable LennardJonesium in the vapor-liquid coexistence region. Work is also in progress on metastability with respect to ordered phases (i.e., the supercooled Lennard-Jones liquid), but here we restrict our attention to vaporliquid metastability exclusively. A n isolated, homogeneous system is metastable when it does not have the largest entropy consistent with its energy. Spontaneous fluctuations will eventually drive a metastable system toward stable equilibrium via the formation of one or more new phases. In order to measure a property of a metastable system, the observation time must be shorter than the system’s lifetime. When this constraint is satisfied, one finds that metastable systems have experimentally reproducible thermophysical properties (e.g., Skripov et al., 1988). This situation requires, in addition, that the internal relaxation time associated with a given property be much shorter than the observation time. Thus, a thermophysical property of a metastable system is well defined if
where Zlife, robs, and t i n t are the lifetime, observation time, and internal relaxation time, respectively. For example, measuring the isothermal compressibility of a metastable system requires that the observation time be simultaneously much shorter than the system’s lifetime and much longer than the characteristic relaxation time for density fluctuations. Examples of metastability in fluids abound in nature (e.g., clouds, sap ascent in plants, mineral inclusions, supercooling of insects and fish), in technology (e.g., vapor explosions, preservation of labile biochemicals, tension pumps, ultrasonic cleaning), and in the laboratory (e.g., sonochemistry). Despite the obvious importance of metastable fluids, the development of rigorous microscopicallybased theories of metastability remains incomplete (Penrose and Lebowitz, 1971,1987). Since the thermophysical properties of moderately metastable liquids are reproducible (Le., they are true properties of the system), it is sensible to attempt to calculate these properties using the standard techniques of statistical mechanics. However, a fundamental problem arises if one tries to do so. A metastable state is never a 0 1995 American Chemical Society
3674 Ind. Eng. Chem. Res., Vol. 34, No. 10, 1995
condition of maximum entropy of an isolated system. Hence, it does not contribute significantly to the partition function. For example, a superheated liquid at a given density and temperature will, if given enough time, spontaneously separate into a heterogeneous vapor-liquid mixture. If constrained in temperature and pressure, it will spontaneously boil to yield a homogeneous unsaturated vapor. Therefore, in the thermodynamic limit, homogeneous metastable states make negligible contributions to the system’s partition function. In the above example, the most important contributions to the canonical partition function of a fluid whose temperature and bulk density fall inside the vapor-liquid coexistence region are those coming from the inhomogeneous vapor-liquid mixture. In the thermodynamic limit, such contributions are not only the most important ones: they are the only ones that matter. Consequently, the rigorous calculation of the partition function would yield no information whatsoever on metastable states (Hill, 1956). (See, however, Corti and Debenedetti, 1994, for a discussion of an interesting, albeit idealized exception.) In most commonly used theories, such as the van der Waals theory, homogeneous metastable states appear only as result of the mathematical approximations used to “solve” the partition function. For example, in the maximum term method (Hill, 19561, one seeks the homogeneous (single) density that maximizes the generic term in the partition function. Likewise, metastable states are obtained in mean-field treatments, such as the Bragg-Williams or van der Waals theories, in which the density is constrained throughout to be strictly uniform (Hill, 1956). In order to calculate the properties of a metastable system rigorously, we must constrain the evaluation of the partition function by eliminating those configurations that contain large enough nuclei of the more stable phase. For example, we study superheated liquids by eliminating those configurations containing sufficiently large voids. While rigorous, this approach suffers from an unavoidable limitation. In most real systems, partition functions are impossible to calculate exactly. Since configurations cannot be exhaustively enumerated, neither can they be selectively eliminated. Thus, the rigorous application of any type of constraint, defined by specifying which trajectories in phase space are not allowed, is in general impossible for realistic systems. (See Elkoshi et al., 1985, for an illuminating example of the rigorous application of constraints to an idealized, one-dimensional system.) Fortunately, while constraints are impossible to apply analytically, they can be easily imposed in computer simulations. We have studied computationally one such constraint, the restricted ensemble (Corti and Debenedetti, 1994). This ensemble, first defined by Penrose and Lebowitz (1971, 19871, limits the magnitude of allowed density fluctuations within finite regions of space. By dividing the computational cell into a number of smaller subcells of equal size and demanding that the number of particles within each subcell fall within specified limits, the system is constrained to remain homogeneous inside the metastable (and unstable) regions. The restricted ensemble is the computational equivalent of the meanfield theory, except that the resulting equation of state is now exact, having been “derived”by the computer in the course of several simulations. We have also presented preliminary results in the void-constrainedensemble (Corti and Debenedetti, 1994).
In this approach, which is closer in spirit than the restricted ensemble to the rigorous calculation of the partition function via the selective elimination of configurations, we study superheating by preventing the formation of voids (bubbles) larger than a given size. The void-constrained ensemble, unlike the restricted ensemble, is meaningless inside the unstable region. The purpose of this paper is to study the superheated Lennard-Jones fluid in the void-constrained ensemble and to compare the equation of state in the voidconstrained and restricted ensembles. An important question that we seek to answer is whether there exists a “natural” constraint to study superheated liquids. To this end, we use the probability distribution of voids to determine when a constraint is artifically severe. The paper is organized as follows: in section 2, we review the restricted ensemble and revisit some of our previous results. In section 3 we explain the voidconstrained ensemble, and we present new results spanning a wide range of temperatures and densities. We also investigate the severity of void constraints by means of the void distribution function. Section 4 presents our conclusions. 2. Restricted Ensemble In a restricted ensemble simulation, arbitrary bounds are imposed on the magnitude of allowed density fluctuations (Corti and Debenedetti, 1994). These density fluctuations, averaged over a control volume, are bounded by subdividing the system volume into smaller subcells of equal size. If the number of subcells is v and the instantaneous number of particles in each subcell i is Ni, then N , is constrained to satisfy
Ni- 6 N 5 N , INi+ 6N
(2)
where & is the average number of particles in each subcell, NIv (N is the total number of particles), and 6N the maximum allowed deviation from the average value. The choice of 6N is crucial. Too large a value of 6N will result in an unconstrained system (and the system may violate the condition of homogeneity). Too small a value will not allow the system to sample properly its configuration space (and this would in principle affect its calculated equilibrium properties). For an unconstrained single-phase system in the thermodynamic limit, the average magnitude of the root mean square fluctuation, AN, of the number of particles about the mean value, N , in an open control volume is given by the rearrangement of the fluctuation equation in the grand canonical ensemble (McQuarrie, 1976): (3)
where /?= lIkT, k is the Boltzmann constant, P is the pressure, and Q is the number density. In our case, AN at each temperature was chosen to be the value corresponding to the saturated liquid estimated from the equation of state of Johnson et al. (1993). The value of 6N is chosen so that 6N 2. AN, thus letting 6N = AN be an upper bound on the severity of the fluctuation constraint: the most severe constraint thus corresponds t o the natural fluctuation value for the most incompressible condition at any given temperature. Therefore, 6NlAiV is a measure of the severity of the constraint: dNlAN % 1 indicates a severely constrained
Ind. Eng. Chem. Res., Vol. 34,No. 10,1995 3575 .6
I
t
I
-.2
t
--1 1 0 2
1 30
I
\
Y
spinodal ,
,
,
,
,
,
/
35
,
,
,
,
,
.40
I
1
, / , , , , ( , , , , / , , , / ,0.75 ,,,,,,,,,,,,,,,/,,,,,,,, 45 . 5 0 . 5 5 .60 65 . 7 0 .75 80 .85
P‘ Figure 1. Effect of the severity of the uniformity constraint on the equation of state of the superheated Lennard-Jones liquid 0, 6NlAN = 1.0; 0,6NIAN = 3.0. The binodal and spinodal curves were obtained from the equation of state of Johnson et al. (1993). All error bars not shown are contained within the symbols. N = 500, Y = 8. From Corti and Debenedetti, Chem. Eng. Sci. 1994, 49,2117-2734.
system, and larger values of 6NlAN indicate progressively weaker constraints. We have implemented this constraint within the canonical ensemble (constant N , V, and T,where V is the volume and T is the absolute temperature) using the standard Metropolis Monte Carlo algorithm (Allen and Tildesley, 1992). The Monte Carlo algorithm follows the standard Metropolis method guidelines except that the violation of the density constraint takes precedence over all other acceptance criteria. Monte Carlo simulations were performed at four subcritical temperatures: T* = kTk = 0.75,0.9,1.1,and 1.28where E is the Lennard-Jones well depth. The best estimate of the critical parameters are 1.316f 0.006 and e,*= eca3 = 0.304 f 0.006,where u is the separation corresponding t o zero potential (Smit, 1992). The LennardJones potential was truncated at 2.50, and a long-range correction (Allen and Tildesley, 1992) was added to obtain the configurational energy and pressure (U*= UIN and P* = Pa31c). Each simulation was performed a t a specific density, e* = ea3, number of boxes, Y , and severity of constraint, 6NlAN. Every simulation was preceded by an initial equilibration period of 1000 Monte Carlo steps per particle (MCS). Initial particle positions were chosen at random (Allen a_nd Tildesley, 1992) inside each subcell, with a t least Ni - 6N particles in each subcell. No two particles were allowed t o be separated by a distance less than 0.7~.ARer equilibration, ensemble averages (pressure and energy) were calculated from 3000 MCS, corresponding to 1.5 x lo6 configurations for N = 500. Ensemble averages and standard deviations were obtained from block averages accumulated during the simulations. Pressure uncertainties in the superheated liquid region were within &lo%, except for densities near saturation in which fluctuations sometimes exceeded this value. Figure 1 (Corti and Debenedetti, 1994) shows the effect of the severity of the constraint, dNlAN, on the equation of state of the superheated Lennard-Jones liquid. As long as the constraint is severe enough t o prevent incipient phase separation, the value of 6NlAN
r=
has a negligible effect on the equation of state. This result is not surprising for mildly superheated liquid states, since the most severely constrained system (6Nl AN = 1.0)is still allowed to sample phase space properly (remembering that AN is chosen to be the natural value for the saturated liquid). However, this result is quite surprising for highly superheated liquid states where density fluctuations are quite severe. For these states, the allowed particle fluctuation, 6N, is smaller than the natural fluctuation for the saturated liquid. Yet in spite of this, we see an insensitivity to the severity of the constraint a t all superheated liquid densities. The main finding of our studies in the restricted ensemble is that the equation of state is insensitive to the severity of the restriction on the allowed density fluctuations as long as the size of a subcell exceeds three average interparticle separations (see Corti and Debenedetti, 1994, for a detailed discussion). We now compare these results with those obtained by preventing the formation of voids that exceed a specified size.
3. Void-Constrained Ensemble The restricted ensemble imposes bounds on the magnitude of density fluctuations within a chosen sample volume. This constraint is a zeroth-order approximation: we are not interested in the actual configuration of the particles within our control volume. On the other hand, one can simulate a superheated liquid by limiting the maximum size of a void (or cavity) allowed to form. This constraint is now localized (as opposed to averaged over a subcell) and biased (restrict voids, not clusters). In this case, metastability is maintained by suppressing cavitation (a method quite analogous to the experimental techniques used to study superheated liquids). The success of this procedure hinges upon the development of an efficient void-counting algorithm. To define void sizes, we performed a VoronoiDelaunay tessellation (Tanemura et al., 1983). The definition of the Voronoi and Delaunay dual tessellation is as follows: a point x (where x is a vector) belongs to the Voronoi cell of atom i located a t position xi if it is closer to xi than to any other point xj of the system (Riedinger et al., 1988)
where Vi denotes the Voronoi polyhedron which surrounds atom i. The dual Delaunay construction is a tiling of space by simplices (d-dimensional triangle where d is the system’s dimension) whose vertices are the original atom points xi while the centers of the spheres circumscribing these simplices are the Voronoi vertices. The diameter of the circumsphere minus the Lennard-Jones diameter, u, is taken as the effective size of a void about a Lennard-Jones particle. For convenience, these voids were visualized as spherical particles with diameters equal to their effective lengths. Figure 2 illustrates the dual tessellation and the definition of a void size. In the void-constrained ensemble, the relevant parameter is the diameter of the maximum allowed void, dga = d,,Iu. The superheated liquid will be prevented from sampling those configurations which contain voids greater than dga. If we reference dgaxto the average interparticle separation, ( ~ * ) - l ’ then ~, we obtain a dimensionless parameter, b [ = C ~ Z ~ ( Q * ) ~ ’ ~ I ,
3576 Ind. Eng. Chem. Res., Vol. 34,No. 10,1995 .8
.4 -
.e
0 -
*&
-.I
-
-.8
-
-.2
-
-.a -1.0
-
A' 0.75
I
.
-1.2
,30 .35 .4Q . 4 5 .50 ,55 .60 .65 .70 .75 .80 .85
P' Figure 3. Effect of the severity of the void-size constraint on the equation of state of the Superheated LennardJones liquid A, b = 1.0: 0,b = 1.5; 0,b = 2.0. N = 256. The solid line passes
through each T' = 1.1isotherm, the dashed line through each P = 0.9 isotherm, the dot-dashed line through each l" = 0.75 isotherm. The binodal and spinodal curves were obtained from the equation of state of Johnson et al. (1993). All error bars not shown are contained within the symbols.
Figure 2. (top) Two-dimensional illustration of the Vomnoil Delaunay dual construction. The central atom i is surrounded by atoms j. The solid lines form the Vomnoi polygon about atom i (which contains all p i n t s in the plane which are closer to atom i than to atom j ) . The dashed lines form the Delaunay triangles whose circumspheres are centered a t the corresponding vertices of the Vomnoi polygon. (bottom) The determination of the size (diameter) of "void particles" about a LennardJones particle in the superheated liquid. The LennardJones particles lie on the vertices of the Delaunay triangles (smaller shaded circles). The void particles are centered on the vertices of the Vomnoi polygon. For clarity, only one void particle is shown (larger shaded circle). Its diameter equals that of the sphere circumscribingthe Delaunay triangle and centered on the vertex of the Voronoi polygon, minus the LennardJanes diameter, v
which quantifies the severity of the constraint. The value of dzm by itself is not an effective indicator of the severity of the void constraint. Rather, the size of the maximum void allowed in relation to the average size of voids naturally occurring in the metastable liquid is a more relevant indicator of the severity of the constraint. Hence, we use the parameter b, which scales dEm to the average interparticle separation, to quantify the severity of the constraint. The smaller the value of b, the more severely constrained the superheated liquid. Note that b is a purely geometric constraint and is not a function of the system temperature. We are currently studying density- and temperature-dependent constraints (e.g., dEm < dLt, where dLt is the critical bubble diameter estimated from nucleation theory); in this paper we restrict our attention to the purely geometric constraint defined above. Simulations of the superheated LennardJones fluid in the void-constrained ensemble were performed using the canonical Monte Carlo algorithm (constant N,V,and
T). Normal Monte Carlo procedures were followed except that a particle move was rejected if that move created a void larger than some specified value. All simulations were started from a n initial random configuration of particles which had no voids greater than An initial equilibration period the chosen limit, d:* of 1000 MCS was performed. Ensemble averages were calculated over 4000 MCS. For all simulations, the system size was N = 256. We determined the equation of state of the superheated liquid a t three subcritical temperatures: T+= 0.75,0.90,and 1.1. In each case, we chose void constraints of b = 1.0,1.5,and 2.0. In so doing, we were interested in determining the effect of the severity of the constraint on the resulting equilibrium properties, specifically pressure and configurational energy, of the superheated liquid. Uncertainties in the pressure (again determined from fluctuations in block averages) were within i5% except for densities near or greater than saturation in which they sometimes exceeded *lo%. In all cases, uncertainties in the configurational energy were less than f2%. Work is underway on system sizes of N 2 500. Our initial choice of N = 256 was determined by the long time required for the VoronoiDelaunay tessellation. For the quantities reported here, however, system size effects are not expected because the equation of state in the restricted ensemble is system-size independent (Corti and Debenedetti, 1994) and calculated pressures in the voidconstrained ensemble for systems which are not overconstrained are consistent with restricted ensemble and equation of state calculations (Johnson et al., 1993),as will be shown below. Figure 3 shows the equation of state of the superheated Lennard-Jones liquid in the void-constrained ensemble for various values of b. Unlike the restricted ensemble, the void-constrained ensemble is quite sensitive to the severity of the constraint. In addition, the constraint has quite different effects on the pressure for different temperatures and extents of superheating. At low temperatures and for densities close to the satura-
Ind. Eng. Chem. Res., Vol. 34,No. 10,1995 3577 -3.6
4 .4 .2
-4.4 -
z *\
D
-4.6
0
-
-.2
-4.0 -
a -.4 -.6
-5.2 -5.0
-5.4
-
-5.6
-
-5.0
-.6
I
-6 0 " ' " ' " 45 50
-1.0
-
'\ spinodal
~
'
'
55
'
'
"
"
'
60
'
~
'
"
~
65
"
'
1
1
70
1
'
1
"
1
75
'
60
85
65
P' Figure 4. Effect of the severity of the void-size constraint on the configurational energy of the superheated Lennard-Jones liquid: A, b = 1.0; 0 , b = 1.5; 0,to b = 2.0. Arrows indicate the spinodal densities which were obtained from the equation of state of Johnson et al. (1993). All error bars not shown are contained within the symbols. N = 256.
tion curve, b has little effect on the equation of state. At the highest temperature simulated, on the other hand, b = 2.0 was not a severe enough constraint to prevent the system from partially phase separating to whatever extent is possible in a finite system (note the positive pressure). For all the temperatures studied, b = 1.0 yielded significantly lower pressures than the other less severe constraints at densities close to the spinodal. Such a severe pressure decrease occurs up to densities progressively closer to the binodal as the temperature increases. In fact, b = 1.0 leads to distinctively lower pressures well into the one-phase liquid region for F = 1.1. Similar trends are observed for the configurational energy per particle, U*/N,they are shown in Figure 4. For densities close to the saturation curve, the severity of the constraint has no effect on the configurational energy. However, sufficiently close to the spinodal, b does affect the configurational energy. For severe superheating, a decrease in b (increase in the severity of the constraint) causes an increase in the configurational energy. Again, this sensitivity to b persists closer to the binodal as the temperature is increased. Such trends were not observed in the restricted ensemble where U*/N was essentially independent of the severity of the constraint. We see that the superheated liquid responds quite differently to space-averaged unbiased constraints (the restricted ensemble) than it does to localized biased constraints (suppression of voids, but not of aggregation). However, as shown in Figure 5, there is a range of constraints in the void-constrained ensemble which yields pressures consistent with those of the restricted ensemble. These results, together with Figure 3,suggest that the pressures obtained in the void-constrained ensemble for b = 1.0 at all temperatures studied are artificially low due to the seventy of the constraint. The question naturally arises as to whether there is an objective way of determining when a constraint is too severe. To answer this question, we analyze the probability density distribution of void sizes in the superheated
Figure 6. Comparison of the equations of state of the superheated Lennard-Jones liquid obtained from the restricted and voidconstrained ensembles: 0,simulations in the restricted ensemble where N = 500, 6 N I m = 1.0, and v = 8; 0, simulations in the void-constrained ensemble where N = 256 and b = 2.0, except for r" = 1.1 where b = 1.5. The binodal and spinodal curves were obtained from the equation of state of Johnson et al. (1993). All error bars not shown are contained within the symbols. The dashed line joining void-constrained data closer to the spinodal denotes partial phase separation. I
i -2
1.0 -
e0
'0
Q
2h 1.2 d
h
Q
0
Q
-
LI Q
.t 0
'-
I'
1
0
.2
.4
.6
.0
1.0 1.2 1.4 1.6 1.0 2.0 2 2
Figure 6. Void-size probability density distributions of the superheated Lennard-Jones liquid for various constraints in the void-constrained ensemble at !P = 0.75 and e* = 0.7: A, b = 1.0; 0,b = 1.5; 0,b = 2.0. N = 256. Solid and dashed lines are drawn through data points for b = 2.0 and b = 1.0, respectively, as a guide to the eye.
liquid for various constraints. Figure 6 displays the void-size distribution for T* = 0.75 and e* = 0.70 for each of the constraints simulated. The distribution is a probability density, obtained by dividing the probability distribution by the histogram sampling width of 0.10. Void sizes for the entire system were determined every 5 MCS (or 800 times during a simulation of 4000 MCS). For b = 2.0, voids as large as 2.00 contribute to the superheated liquid's equilibrium properties. For b = 1.5,voids up to 1.690 are sampled. In contrast, we see
3678 Ind. Eng. Chem. Res., Vol. 34,No. 10, 1995 2.4 I,
2.2
I
-
1.8 2 w 1.6 2.0
~
c
4
1.4-
2
1.2
2
1.0
(d
L, 0 L
a
-
'8 ,6
-
.4
-
.2 O m , m l
0
2
"
4
"
6
'
I
8
'
1 0
1 2
1 4
0
1 6
0
2
4
6
m, I , 8 1 0 1 2 1 4 1 6 1 8 2 0 2 2 2 4
d/. Figure 7. Void-size probability density distributions of the superheated Lennard-Jones liquid for various constraints in the void-constrained ensemble at Tc = 0.75 and e* = 0.8: A, b = 1.0; 0, b = 1.5; 0,b = 2.0. N = 256. Solid and dashed lines are drawn through data points for b = 2.0'and b = 1.0, respectively, as a guide to the eye.
that the statistical void distribution for b = 1.0 is artificially truncated beyond 1.20: the distribution is now qualitatively different, there being no large-void tail. As long as the system is able to sample a void distribution including the large-void tail, the calculated pressures are statistically indistinguishable (see Figure 3 for b = 2.0 and b = 1.5). For example, in Figure 6, the void-size distributions for b = 1.5 and b = 2.0 are statistically identical (except for void diameters greater than 1 . 6 9 where ~ the distribution corresponding t o b = 1.5 is truncated; note, however, the very low values of the probability density); hence, the corresponding system pressures are statistically indistinguishable. For b = 1.0, on the other hand, the artificial truncation of the void distribution causes a drastic reduction in the system's pressure, a result entirely consistent with the fact that the mechanism for relieving tension in a liquid is cavitation. Figures 7-9 show analogous calculations for (F= 0.75, e* = 0.8),(F= 1.1,e* = 0.71, and (F = 1.1,e* = 0.6).Again, we see that forb = 1.0 the largevoid tail of the distribution is absent, resulting in pressures lower than expected. Deviations between b = 1.0 and the other void-size distributions are greatest for F = 1.1 and e* = 0.6. The distributions corresponding to b = 1.5 and b = 2.0 are indistinguishable in every case. Comparison with Figure 3 again shows that the pressures corresponding to b = 1.5 and b = 2.0, both of which allow for a large-void tail in the distribution, are statistically indistinguishable. We have calculated the radial distribution function of the superheated liquid under various constraints. Figures 10 and 11show the radial distribution function of the superheated liquid a t T* = 1.1,e* = 0.6 and r" = 1.1,e* = 0.7,respectively, for b = 1.0, 1.5, and 2.0. As a consistency check, void-constrained pressures obtained from the virial were found to match, within error bars, the pressures obtained by integrating the radial distribution function (McQuanie, 1976). The sensitivity of the pair correlation function to the void constraint depends on the state being simulated. While the curves forb = 1.5 and b = 2.0 are indistinguishable regardless of the degree of superheating (which is t o be
u u Figure 8. Void-size probability density distributions of the superheated Lennard-Jones liquid for various constraints in the void-constrained ensemble at Tc = 1.1and e* = 0.7: A, b = 1.0, 0, b = 1.5; 0, b = 2.0. N = 256. Solid and dashed lines are drawn through data points for b = 2.0 and b = 1.0, respectively, as a guide to the eye. 2.2 I
1.8
h 1.6
5
VI
aJ
1.4
I
1 -
ah 1 . 2 -
5
5(d
1.0
-
e
0 C
a
-
.6 .4
-
2 0
0
.2
.4
.6
.8
1.0 1 . 2 1 . 4 1 . 6 1 . 6 2 . 0 2 . 2 2 4
d/o Figure 9. Void-size probability density distributions of the superheated Lennard-Jones liquid for various constraints in the void-constrained ensemble at Tc = 1.1and e* = 0.6: A, b = 1.0; 0,b = 1.5; 0, b = 2.0. N = 256. Solid and dashed lines are drawn through data points for b = 2.0 and b = 1.0, respectively, as a guide to the eye.
expected, since the equation of state is insensitive to the severity of the constraint for b z 1.5), differences with respect t o the overconstrained curve, b = 1.0, depend sensitively on the proximity to the spinodal. In Figure 10, where the density is close to the superheated liquid spinodal, the b = 1.0 curve is clearly different from the b = 1.5 and b = 2.0 functions. On the other hand, in Figure 11,where the density corresponds t o the stable liquid, the radial distribution function for b = 1.0 is only slightly different from the b = 1.5 or b = 2.0 curves. These small differences, however, translate into appreciable differences in pressure (see Figure 3). The radial distribution function, then, does not appear to be an effective indicator of artificial (overconstrained) systems. Furthermore, it appears to provide less in-
Ind.Eng. Chem. Res., Vol. 34, No. 10, 1995 3679
0.0 I 0.0
I
I , 1 .o
2.0
3.0
4.0
TI0
Figure 10. Radial distribution function of the superheated Lennard-Jones liquid for various constraints in the void-constrained ensemble at Tc = 1.1 and e* = 0.6: -, b = 2.0;- - -, b = 1.5;*.*, b = 1.0. 2.5
I
0.51 i 0.0 L 0.0
I
I ,
1 .o
2.0
3.0
4.0
do
Figure 11. Radial distribution function of the superheated Lennard-Jones liquid for various constraints in the void-constrained ensemble a t F = 1.1 and e* = 0.7: -, b = 2.0;- - -, b = 1.5; b = 1.0. .*a,
formation about what constitutes an unphysical system than the void-size probability distribution. Even for those densities in which the radial distribution function is quite sensitive t o the choice of the severity of the constraint, there is little information that one can infer about changes in the short- and long-range structure in the superheated liquid. 4. Conclusions
The rigorous statistical mechanics of metastability is the statistical mechanics of constrained ensembles. As long as a metastable system is studied over time scales shorter than its lifetime and longer than its characteristic internal relaxation times, the corresponding thermophysical properties are reproducible. The rigorous calculation of the equilibrium properties of metastable systems requires that constraints be imposed that serve to prevent the sampling of configurations which would lead the system away from homogeneity (i.e., phase separate). For nontrivial systems, configurationscannot be completely enumerated and, hence, cannot be exhaustively eliminated. For real systems, the systematic elimination of unwanted trajectories in phase space is
an intractable problem. However, in contrast to analytical calculations, the removal of unwanted configurations, or a constraint, is easily imposed in computer simulations. Here, we have studied one such constraint: the voidconstrained ensemble. In this ensemble, limits are placed on the maximum size of voids allowed to form in the liquid. The equation of state of the superheated, voidconstrained LennardJones liquid is sensitive to the severity of the constraint. This sensitivity increases with temperature and with the extent of superheating. Comparison with restricted ensemble simulations, in which bounds are imposed on the magnitude of allowed density fluctuations within a control volume, show that the system responds quite differently to space-averaged unbiased constraints (restricted ensemble) than it does to localized biased constraints (void-constrained ensemble). If the maximum allowed voids are equal to, or larger than, 1.5 times the average interparticle distance, then the pressure in the void-constrained ensemble is consistent with that obtained in the restricted ensemble. The statistical distribution of voids was used t o determine when a constraint leads to unphysical results. As long as the void distribution has a well-developed large-void tail, the pressure is independent of the constraint. Constraints that result in abruptly interrupted void distributions lead to significantly lower pressures, a finding that is consistent with the fact that the mechanism in which tension is relieved in a stretched liquid is cavitation. Over the range of conditions examined so far, artificially low pressures (or high tensions) result when the maximum allowed void does not exceed one average interparticle separation. We also looked at the radial distribution function of the superheated liquid under various constraints. We found that these curves yielded no further information than could be obtained from the simulation pressures, failing to provide an additional means by which we could determine whether a system was overconstrained. In ongoing work, we are studying the connectivity of voids under various constraints. We plan to investigate the relationship between measures of void connectivity, severity of the constraint, and extent of superheating. Such studies, combined with further comparisons between different constrained ensembles along the lines presented here, should greatly improve our understanding of the thermophysical properties of constrained fluids.
Acknowledgment The financial support of the US. Department of Energy, Office of Basic Energy Sciences, Division of Chemical Sciences, Office of Energy Research, through Grant DE-FG02-87ER13714, is gratefully acknowledged.
Literature Cited Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, 1992. Corti, D. S.; Debenedetti, P. G. A Computational Study of Metastability in Vapor-Liquid Equilibrium. Chem. Eng. Sci. 1994,49,2717-2734. Elkoshi, Z.; Reiss, H.; Hammerich, A. D. One-Dimensional Rigorous Hole Theory of Fluids: Internally Constrained Ensembles. J . Stat. Phys. 1985, 41, 685-708.
3580 Ind. Eng. Chem.Res.,Vol. 34, No. 10,1995 George, M. F.; Becwar, M. R.; Burke, M. J. Freezing Avoidance By Deep Supercooling of Tissue Water in Winter-Hardy Plants. Cryobiol. 1982,19,628. Hill, T. L. Statistical Mechanics. Principles and Selected Applications; McGraw-Hill: New York, 1956. Johnson, J. K.; Zollweg, J. A.; Gubbins, K. E. The Lennard-Jones Equation of State Revisited. Mol. Phys. 1993,78, 591-618. McQuarrie, D.A. Statistical Physics; Harper & Row: New York, 1976;pp 61-62. Penrose, 0.; Lebowitz, J. L. Rigorous Treatment of Metastable States in van der Waals-Maxwell Theory. J . Stat. Phys. 1971, 3,211-236. Penrose, 0.;Lebowitz, J. L. Towards a Rigorous Molecular Theory of Metastability. In Fluctuation Phenomena; Montroll, E., Lebowitz, J. L.; Eds.; North Holland: Amsterdam, 1987; Chapter 5. Reiss, H. The Existence of the Spinodal, An Incompletely Solved Problem in the Thermodynamics of Solids. Ber. Bunsen-Ges. Phys. Chem. 1975,943-957. Riedinger, R.; Habar, M.; Oelhafen, P.; Guntherodt, H. J.
About the Delaunay-Voronoi Tessellation. J . Comput. Phys. 1988,74, 61-72. Skripov, V. P.; Sinitsyn, E. N.; Pavlov, P. A.; Ermakov, G. V.; Muratov, G. N.;Bulanov, N. V.; Baidakov, V. G. Thermophysical Properties of Liquids in the Metastable (Superheated) State; Gordon and Breach New York, 1988. Smit, B. Phase Diagrams of Lennard-Jones Fluids. J . Chem. Phys. 1992,96,8639-8640. Tanemura, M.; Ogawa, T.; Ogita, N. A New Algorithm for ThreeDimensional Voronoi Tessellation. J.Comput. Phys. 1983,51, 191-207.
Received for review January 24,1995 Accepted April 28,1995*
IE950071K
* Abstract published in Advance ACS Abstracts, August 15, 1995.
