Homogeneous Bubble Nucleation in Stretched Fluids: Cavity

Umbrella sampling Monte Carlo simulations are used to calculate the reversible work of forming cavities of various sizes within the superheated Lennar...
0 downloads 0 Views 96KB Size
Download PDF
Ind. Eng. Chem. Res. 2002, 41, 1113-1121

1113

Homogeneous Bubble Nucleation in Stretched Fluids: Cavity Formation in the Superheated Lennard-Jones Liquid Sudeep Punnathanam and David S. Corti* School of Chemical Engineering, Purdue University, West Lafayette, Indiana 47907-1283

A consideration of various ideas set forth within the scaled particle theory of hard particle fluids, which are also applicable to systems whose particles interact via attractive potentials, suggests that cavity formation plays an important role in the molecular mechanism of the liquid-to-vapor transition. Umbrella sampling Monte Carlo simulations are used to calculate the reversible work of forming cavities of various sizes within the superheated Lennard-Jones liquid maintained at several negative pressures. A critical cavity size is found to occur, beyond which the liquid would phase separate if not for a density constraint that is applied during the simulation. The work of forming this critically sized cavity and its radius is found to decrease as the liquid approaches the spinodal. A critical cavity size is also found for the superheated liquid at positive pressures. The focus on cavity growth in superheated liquids amounts to a new way of studying bubble nucleation and should lead to an improved molecular-based understanding of the kinetics of first-order phase transitions in liquids. 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. Superheated liquids are metastable with respect to the vapor phase and will not remain metastable indefinitely. Spontaneous fluctuations will eventually drive the superheated liquid to evolve to a more stable condition (i.e., the vapor phase) over time. The superheated liquid relaxes toward the more stable vapor phase via the formation of bubbles (i.e., bubble nucleation). Thus, a key quantity of interest is the lifetime of a superheated liquid. Nevertheless, a detailed understanding of the rates and mechanisms of bubble formation in superheated liquids remains incomplete. Several scientific and technological applications depend on accurate knowledge of the lifetimes, or bubble nucleation rates, of superheated liquids. Both in industrial practice and in the laboratory, extensions of a liquid beyond the phase transition point are commonplace.1-3 For example, bubble nucleation is important for understanding industrially hazardous vapor explosions1 and controlling problematic cavitation erosion4,5 and is important in the initial stages of sonoluminescence.6 Bubble nucleation also plays a key role in sonochemistry, in which the formation and subsequent collapse of bubbles is used to enhance the rates of various chemical reactions.7 Knowledge of the kinetic limits of superheated liquids also impacts materials processing, offering insights into the possible improvements of current approaches in the manufacture of new materials.8,9 The formation of a vapor phase from a superheated liquid is an activated process; i.e., a free energy barrier must be surmounted before the phase transition can occur. This is achieved when spontaneous fluctuations in the liquid result in the formation of an embryo of the * School of Chemical Engineering, Purdue Unversity, West Lafayette, IN 47907-1283. E-mail: [email protected].

new phase (microscopic bubbles or small regions where the density is smaller than the bulk). The liquid may contain suspended and dissolved impurities, as well as solid surfaces, that can provide preferential sites for the formation of the embryo of the new phase. This is the case in most practical circumstances, and the process is known as heterogeneous nucleation. In the absence of impurities or solid surfaces, spontaneous fluctuations within the liquid provide the only route toward the formation of the embryo of the gas phase. The process is then called homogeneous nucleation. In both cases, if an embryo is less than some critically sized bubble, the vapor embryo collapses back into the superheated liquid. Once an embryo exceeds this critical size, the bubble spontaneously grows to macroscopic size. Note that homogeneous nucleation always occurs in the superheated liquid. The relative kinetics of heterogeneous or homogeneous nucleation determines which of the two processes is more dominant (not every impurity or solid surface is an effective heterogeneous nucleation site). Homogeneous nucleation is most significant for large superheats or low concentrations of effective heterogeneous nucleation sites. In this paper, we focus on the problem of homogeneous bubble nucleation in liquids. In general, the number of critical nuclei formed per unit time per unit volume, or the nucleation rate, j, can be expressed as1

j ) Bexp[-Wcr/kT]

(1)

where T is the bulk fluid temperature, k is Boltzmann’s constant, Wcr is the reversible work required to form the nucleus of critical size, and B is a kinetic preexponential factor that is believed to be weakly dependent upon temperature. Due to its appearance in the exponential, the nucleation rate is most sensitively dependent upon Wcr. Therefore, previous studies have focused on the determination of Wcr. This paper also concentrates on the estimation of Wcr, which in the case of bubble nucleation corresponds to the reversible work of forming a critical bubble. Once sufficient progress has been made in understanding and accurately estimating

10.1021/ie010554q CCC: $22.00 © 2002 American Chemical Society Published on Web 10/10/2001

1114

Ind. Eng. Chem. Res., Vol. 41, No. 5, 2002

Wcr, our attention will shift to the estimation of the kinetic prefactor B. Within the context of the classical nucleation theory (CNT), the reversible work of forming the critical bubble, assumed to be spherical, is given by1

Wcr )

16πγ3 3(P′ - P)2

(2)

where P′ is the pressure inside the critical bubble, P is the pressure of the surrounding liquid, and γ is the surface tension of the interface (assumed to be independent of the size of the bubble). As the extent of penetration into the metastable region increases, the free energy barrier, or Wcr, becomes smaller. Likewise, the size of the critical bubble also decreases as the degree of superheating is increased. The CNT description of bubble nucleation has several drawbacks, even though CNT predicts reasonably accurate limits of superheating for several liquids at atmospheric or higher pressures.1-3 For example, CNT assumes that the fluid inside the embryos can be described by a uniform bulk phase. This assumption is valid for systems near the coexistence curve but fails for moderate to large extents of superheating. Consequently, CNT fails to predict loss of stability; i.e., the work of forming the critical nucleus does not go to zero at the spinodal (where thereafter, phase separation is a spontaneous process and no activation barrier needs to be surmounted). Hence, CNT predictions of bubble nucleation rates can be in error by several orders of magnitude. Although improvements have been made on its predictive capability,10 the CNT approach relies solely on macroscopic arguments to describe an inherently microscopic process. Important advances in the understanding of the phenomenon of nucleation were obtained through the use of density-functional theory (DFT) by Oxtoby and co-workers.11-14 Within DFT, the barrier to nucleation was obtained by locating the stationary condition of the inhomogeneous grand potential with respect to the density profile across the interface between the bulk metastable phase and the nucleus of the stable phase. DFT represented a large step forward in our understanding of nucleation, being applicable to both vaporto-liquid and liquid-to-vapor nucleation, and even liquidto-solid nucleation. Calculations for bubble nucleation in superheated liquids13-15 revealed that nonclassical effects are much more pronounced in bubble nucleation than in droplet nucleation. Unlike CNT, DFT is able to predict loss of stability. Yet, like CNT, DFT is a continuum theory, although the relevant length scales of this approach are of microscopic size. Also, DFT requires that the shape of the nucleus be specified a priori, like CNT, which is assumed to be spherical (an important limitation of the theory, particularly at high superheatings). The direct molecular simulation of bubble formation in superheated liquids has also revealed important insights into the fundamental understanding of bubble nucleation.16-18 For example, Shen and Debenedetti16 used a biased Monte Carlo (MC) algorithm to determine the reversible work of forming a critical bubble in their simulation cell. Their results correctly predict the loss of stability as the spinodal is approached. Interestingly, as the free energy barrier is mounted, a large nonspherical cavity (ramified and nonspherical region

devoid of particle centers) forms that spans the entire system, casting doubt on previous assumptions that the critical bubble is spherical (at least for large superheatings). The picture of the molecular mechanisms of bubble nucleation that emerges from this study is quite different from that of CNT, and even DFT. The remainder of this paper is structured as follows. Section 2 outlines our new approach to the problem of bubble nucleation in superheated liquids. Recent work on the problem of vapor-to-liquid nucleation and a consideration of various ideas set forth within the scaled particle theory of hard particle fluids suggest that cavity formation plays an important role in the molecular mechanism of the liquid-to-vapor transition. Section 3 discusses the computational methods used to determine the work of forming cavities of various sizes within the superheated Lennard-Jones (LJ) liquid at several negative pressures. Results are presented and discussed in section 4. Conclusions and suggestions for future work are outlined in section 5. 2. Cavity Formation in Superheated Liquids Our approach to the problem of homogeneous bubble nucleation is different from previous studies and follows from earlier work on the molecular theory of vapor-toliquid nucleation.19-21 This work, which has led to significant advances in our understanding of nucleation, focused on the development of a rigorous and consistent molecular theory of vapor-to-liquid nucleation within a dilute supercooled vapor. The statistical mechanical expressions derived in this work are also applicable to the problem of bubble formation in a superheated liquid. One important exception is noted. Although the precise molecular definition of a liquid cluster is arbitrary, there are several generally agreed upon criteria (such as a bond-distance criterion20,21 or an interaction energy criterion19) that are used to determine what configuration of particles may be labeled as a liquidlike cluster. In contrast, what constitutes a bubble on the molecular level is not known, or at least there is no readily apparent physical microscopic definition of the bubbles that participate in the liquid-to-vapor transition. We currently ignore this aspect of the bubble nucleation problem and instead focus on some other relevant issues. For the case of vapor-to-liquid nucleation, the reversible work required to form a cluster of size n in a volume v, W(n/v), is given by19-21

W(n/v) ) -kT ln(Λ-3n3/2V) + Pv - nµ + Uo + F′(n, v, T) (3) where Λ is the thermal de Broglie wavelength, V is the total volume of the system (vapor), P is the pressure of the surrounding vapor phase, µ is the chemical potential of the vapor, Uo is the interaction energy between the cluster and the surroundings, and F′ is the Helmholtz energy of the cluster given some definition of the cluster. The first term on the right side accounts for the translational free energy of the cluster (the cluster can be created anywhere within the volume V), the second term describes the work of forming a cavity (a spherical region devoid of molecular centers) of volume v within the vapor, the third term is the work of removing n particles from the bulk vapor phase, and the fourth term describes the potential energy of interaction between the bulk and the cluster. The separation of translational

Ind. Eng. Chem. Res., Vol. 41, No. 5, 2002 1115

and internal degrees of freedom of the cluster20 is clearly seen in the above expression (i.e., the first term and fifth term on the right side of eq 3). For temperatures below the critical temperature, eq 3 is rigorous except for the second and fourth terms. In general, Uo cannot be separated from F′, and the work of cavity formation is only equal to Pv if the vapor behaves as an ideal gas. For sufficiently dilute vapors, the error introduced by these two terms is quite small (a more rigorous and general expression can be found elsewhere22). Note that eq 3 can be derived without explicitly stating the chosen definition of the liquid cluster. Thus, similar terms would appear in a corresponding expression for the work required to form a bubble of size n in a volume v within a superheated liquid. As the formalism of eq 3 is applied to the bubble nucleation problem, we initially focus on the second term on the right side of eq 3, which now corresponds to the work of forming a cavity within the bulk superheated liquid. We assume, as in eq 3, that the cavity is spherical in shape (a generalization to nonspherical volumes is given elsewhere22). In the case of droplet formation from a low-density vapor, the work of cavity formation is essentially equal to Pv. The relative contribution of this term to the total work of droplet formation is small (since the pressure of the vapor is quite low). In contrast, the work of forming a cavity within a dense liquid can be quite large. As is previously known from the scaled particle theory of hard particle fluids23, the work of forming a cavity whose radius exceeds 1.25d (where d is the diameter of a hard sphere) is greater than about 30kT for hard sphere fluid densities around Fd3 ) 0.85 (where F ) N/V). This quantity is comparable to the work of forming a critically sized bubble within moderately stretched liquids as predicted by CNT via eq 2. The relatively large value of the second term on the right side of eq 3 for liquids suggests that, unlike for a supercooled vapor, cavity formation plays an important, if not dominant, role in the molecular mechanism of the liquid-to-vapor transition. Cavity formation in metastable liquids, and its role in bubble nucleation, has not been, to the authors’ knowledge, studied to any significant extent. Several workers have analyzed cavity formation in stable liquids,24-28 but their motivation dealt with the determination of solvation free energies of various solutes within different solvents (needed, for example, to understand the effects of hydrophobic solutes). In contrast, little research has been done about cavity formation in metastable liquids, particularly liquids under tension, and what role cavities play in bubble nucleation. The statistical geometry of cavities in metastable liquids has received some attention,29-32 but the importance of cavities to phase transitions in liquids was not fully explored. The work of cavity formation within a liquid whose particles interact via some spherically symmetric attractive potential is not known a priori and must be determined via molecular simulation. There are, however, some interesting aspects of cavity formation in metastable liquids that may be deduced from scaled particle theory, some results of which are applicable to systems whose particles interact via attractive potentials. For example, the derivative of the reversible work, W, of forming a cavity with respect to its radius r is given by23,33

∂W ) 4πr2F+kT g 0 ∂r

(4)

where F+ is the local density of surrounding particles in contact with the surface of the cavity. Since F+ by definition is never negative, the work of cavity formation is strictly a monotonic function of the cavity radius. As the cavity approaches macroscopic size, W has the following limiting form23,33

W)

2δ 4πr3 P + 4πr2γ∞ 1 + Ko 3 r

(

)

(5)

where P is the pressure of the bulk fluid, γ∞ is the boundary tension between the bulk fluid and a hard wall (i.e., cavity for which r f ∞), δ is the Tolman length indicating the first-order correction to the boundary tension due to curvature, and Ko is a constant that is required to match the known small cavity limit.23,33 As r f ∞, the work of cavity formation becomes cubic in r, so that the work equals the macroscopic result of Pv. Comparing eqs 4 and 5 reveals that in the large cavity limit, where the cavity becomes equivalent to a hard wall, F+ f P/kT. This limiting value of the contact density in the large cavity limit has some interesting implications for the formation of cavities in dense stable liquids near the coexistence curve.24,27,28,34 For example, in the small cavity limit, the contact density approaches that of the bulk density.23,33 As the cavity increases in size, the surrounding liquid behaves as an elastic medium, resulting in a corresponding increase in the contact density. As the cavity becomes of macroscopic size, however, the contact density must approach the bulk pressure (divided by kT). For dense liquids near the triple point, where bulk pressures are low, the limiting value of the contact density becomes smaller than the bulk density, so that the fluid surrounding the cavity undergoes a drying transition27,28 (a region of low density immediately surrounds the cavity). What is the behavior of the contact density, however, if the liquid is now brought into the metastable region of the phase diagram, specifically to a state in which the liquid is under tension (i.e., P < 0)? In the large cavity limit, eq 5 implies that W f 4πr3P/3 < 0. Thus, W must reach a maximum value. W increases with r in the small cavity limit but must eventually decrease with r beyond some point in order to approach negative values for large r. This conclusion, however, is inconsistent with the required monotonicity of W with r as required by eq 4. Equations 4 and 5 therefore imply that a liquid under tension (negative pressure) cannot be maintained next to a hard wall, or cannot contain a cavity of macroscopic size. If the negative pressure liquid were able to contain a cavity of macroscopic size while remaining homogeneous, eqs 4 and 5 would require that F+ f P/kT < 0, violating the requirement that F+ g 0. We conclude that eq 5 cannot be used to describe the formation of cavities beyond some given size in liquids that are under tension. The apparent breakdown of eq 5 for negative pressure liquids raises the following questions: What is the limit to the size of a cavity that a negative pressure liquid can contain while remaining homogeneous? Must the liquid under tension phase separate when this critical cavity size is reached? The picture for cavity formation in a negative pressure liquid is speculated to be as follows. A superheated liquid at a negative pressure is able to contain small

1116

Ind. Eng. Chem. Res., Vol. 41, No. 5, 2002

cavities without rupturing.29-31,35 Thus, for small cavities, the surrounding liquid will behave as a locally elastic medium so that W will increase with an increase in r. At some intermediate cavity size, however, a drying transition should occur at the cavity surface, where for larger cavity sizes the contact density should begin to decrease with an increase in r. The value of the contact density cannot become negative, but it may become quite small (possibly approaching zero). The sudden decrease in the contact density should cause a change of curvature in the work of cavity formation profile, so that W as a function of r will begin to approach a zero slope (eq 4). A maximum in W is prohibited since a negative slope would imply that F+ < 0. Finally, a state is reached in which the negative pressure liquid cannot accommodate a larger cavity while remaining homogeneous. Any increase in the size of the cavity will cause the metastable liquid to become unstable, i.e., spontaneously phase separate into the corresponding vapor phase. This maximum allowed cavity size, or critical cavity, is analogous to the critical bubble. The above picture of cavity formation in liquids under tension raises the following questions. If the stability of the liquid is governed by a critical cavity, instead of a critical bubble, is the underlying mechanism of bubble nucleation in negative pressure liquids therefore dominated by the formation of cavities? If so, is the free energy barrier to phase separation essentially determined by the work of cavity formation? Considered alone, the growth of cavities leads to the destabilization of the liquid. Does the presence of a few particle centers within the cavity, i.e., a bubble, serve to stabilize the liquid? Or does the cavity instability supersede this? Is similar behavior seen for the growth of cavities in superheated liquids at positive pressures? If cavities are indeed driving the initial stages of the liquid-to-vapor transition, how do these cavities grow and eventually turn into bubbles? This paper discusses some preliminary results gathered from ongoing computational studies of cavity formation in superheated liquids aimed at answering some of the above questions. The focus on cavity growth within liquids amounts to a novel way of looking at bubble nucleation in superheated liquids and should lead to an improved molecular-based understanding of the kinetics of first-order phase transitions in liquids. 3. Simulation Details The reversible work of forming cavities (defined as spherical regions devoid of any particle centers) of various sizes was determined within the superheated Lennard-Jones liquid. The LJ interaction potential between a pair of atoms is given by 12

6

[(σr) - (σr) ]

uLJ(r) ) 4

(6)

where σ is the interatom separation for which the potential is zero and  is the minimum well depth. The current simulations, however, used the following truncated and shifted LJ potential, uLJT, whereby

uLJT ) uLJ - uLJ(4.0σ) r e 4.0σ ) 0 r g 4.0σ

(7)

The truncated and shifted potential eliminated the need

Figure 1. Pressure-temperature projection of the phase diagram of the Lennard-Jones fluid (eqs 6 and 7). The binodal and the superheated liquid spinodal were obtained using the equation of state of Johnson et al.37 Also shown are the state points at which simulations were performed for Figures 2 and 3.

to apply long-range corrections to the potential. Longrange corrections may be applied if the density of the fluid is uniform beyond the cutoff radius. Since our simulations contain a cavity, necessarily generating an inhomogeneity within the simulation cell, long-range corrections would not be appropriate. Figure 1 displays the P-T projection of the phase diagram of the LJ fluid whose potential is given by eq 7. Figure 1 also contains the superheated liquid spinodal, which was also obtained from the equation of state of Johnson et al.37 The critical parameters of this fluid are estimated to be Tc* ) kT/ ) 1.246, Fc* ) Fσ3 ) 0.308, and Pc* ) Pσ3/ ) 0.118. Also included in Figure 1 are some of the state points used for the present study. The reversible work of forming a cavity was obtained via Monte Carlo simulations within the isothermalisobaric ensemble (NPT). Constant pressure simulations enabled us to maintain the liquid at a fixed negative pressure. Constant pressure simulations also allowed the system volume to vary as the cavity grew in size, ensuring that the local density far from the cavity surface approached the value of the bulk density consistent with the prescribed pressure (in which no cavity was present). The work of cavity formation was determined from the following relation23,33

W(λ) ) -kTln[P(λ)]

(8)

where P(λ) is the probability of successfully inserting a cavity of radius λ at a given point in the fluid. P(λ) can be directly calculated from a standard Monte Carlo simulation, in which the probability of successfully inserting a cavity of size λ is determined. At moderateto-high fluid densities, however, the chance of inserting a cavity in which λ > σ becomes low, and the estimation of P(λ) becomes statistically poor. To overcome this problem, an umbrella sampling scheme38 is employed that favors the formation of large cavities. A cavity of a given size is first introduced into the simulation cell and remains fixed at the center of the cell. In addition to the standard MC displacements of the particles (all particle moves that bring a particle center into the cavity are rejected), another MC trial move is imple-

Ind. Eng. Chem. Res., Vol. 41, No. 5, 2002 1117

mented that attempts to change the radius of the cavity while particle positions are held fixed (an attempted increase in λ is rejected if an external particle center is found within the cavity surface). A biasing potential is added to the standard MC sampling scheme to ensure that large cavity sizes are obtained during the simulation. Since P(λ) is an equilibrium property of the fluid, its ensemble average in the NPT ensemble is given by

∫dV ∫δ(λ - λ*)exp[-β(U + PV)]drN 〈P(λ)〉 ) ∫dV ∫exp[-β(U + PV)]drN

(9)

where drN represents the volume elements of the N particles, P denotes the pressure of the system, V is the volume of the system, U is the potential energy of the system, β ) 1/kT, and δ(λ - λ*) is the Dirac δ function that selects only those configurations for which the cavity size is equal to λ. To sample configurations which contain a large cavity, a biasing potential, ψ or w ) exp(-ψ/kT), is introduced, so that eq 9 can be rewritten as

δ(λ - λ*)

∫dV ∫ w w exp[-β(U + PV)]drN × 〈P(λ)〉 ) ∫dV ∫w exp[-β(U + PV)]drN ∫dV ∫w exp[-β(U + PV)]drN ∫dV ∫w1 w exp[-β(U + PV)]drN

〈

〉

δ(λ - λ*) w ) 1 ww

〈〉

w

(10)

where 〈〉w denotes the average quantities as determined using the biased distribution. Since the resulting ensemble averages are obtained within the biased simulation, the results must be corrected at the end in order to obtain the ensemble averages for the original, unbiased system. From eq 8, good statistics for the biased averages will be obtained if ψ(λ) ) -W(λ). Unfortunately, ψ is not known a priori. Thus, a series of biased simulations are run for several successive windows in which the cavity radius is restricted to remain within some range, ∆λ. The work of insertion is calculated during each window and then curve-fitted to a polynomial in λ. This current estimate of the work of cavity formation is used as an estimate for the biasing potential within the next simulation window. The size of each window was adjusted until a uniform distribution of the cavity radii was obtained, thereby ensuring statistically good estimates of the biased averages. The work curves of each window are finally “linked” together by adjusting the value of the constant in the curve fit to obtain the total W versus λ profile. Each simulation consisted of an equilibration period of 10 000 MC cycles (trial moves per particle) followed by a production run of 100 000 MC cycles. The system size of each simulation was at least 3000 particles. The number of particles within the simulation cell for each window was constantly adjusted to ensure that the density profile far from the cavity surface approached the bulk density near the edge of the simulation cell. Also, since the interparticle potential was truncated at 4.0σ, an initial simulation box size of twice the sum of

Figure 2. Reversible work of formation of a cavity of radius λ in a Lennard-Jones fluid at T* ) 0.80 and at various negative pressures.

the cavity radius and 5.0σ was used. Bulk densities were estimated from the equation of state of Johnson et al.37 An evaluation of the density profile around the cavity for various cavity sizes revealed that our choice of system size within each window was appropriate. The number of windows required to obtain W versus λ increases as λ increases. In addition, the range of λ sampled within each simulation window, ∆λ, decreased with an increase in λ; otherwise poor estimates of P(λ) were obtained. Thus, the number of simulation windows rapidly increased as larger cavity sizes were studied. In conjunction with the relatively large system sizes used in this study, the rapid increase in the number of windows prohibited us from determining the work of forming cavities of radii greater than 3.0σ. 4. Results and Discussion The reversible work of forming a cavity of radius λ within the superheated LJ liquid was determined for different negative pressures along two separate isotherms. Figure 1 shows the state points simulated in this work. At T* ) 0.80, we simulated systems at P* ) -0.3757, -0.4680 and -0.5291, (where the binodal and spinodal pressures are 0.0062 and -0.6, respectively); at T* ) 0.85, we simulated systems at P* ) -0.3679, -0.4058 and -0.4367 (where the binodal and spinodal pressures are 0.00997 and -0.483, respectively). The results are shown in Figures 2 and 3. The work of cavity formation profiles in Figures 2 and 3 clearly exhibit maxima, even though eq 4 prohibits the appearance of such maxima. Section 2 also included a discussion indicating that a liquid under tension will phase separate into the corresponding vapor phase if the liquid includes a cavity beyond some critical size. In fact, for values of λ in which W versus λ began to show a change of curvature (e.g., around λ ) 2.0σ for T* ) 0.80 and P* ) -0.3757 in Figure 2), our simulations revealed that the liquid was no longer able to sample the “liquidlike” regions of configuration space throughout the simulation. Beyond some critical cavity size, the volume of the system was found to increase indefinitely during the simulation (since a negative pressure was imposed during the simulation, the volume of the system never stabilized at a “vaporlike” volume appropriate for a positive pressure system). For

1118

Ind. Eng. Chem. Res., Vol. 41, No. 5, 2002

Figure 3. Reversible work of formation of a cavity of radius λ in a Lennard-Jones fluid at T* ) 0.85 and at various negative pressures.

Figure 4. Effect of the density constraint (used to prevent phase separation) on the reversible work of formation of a cavity for T* ) 0.85 and P* ) -0.4058. The bulk density at this temperature and pressure is Fb* ) 0.69. Fcons* denotes the lower limit on the reduced density imposed during the simulations to prevent the metastable liquid from phase separating.

values of λ smaller than this critical size, the volume of the system remained stable, fluctuating about those states consistent with “liquidlike” densities. To determine the work of forming cavities beyond this critical cavity size, a volume, or density, constraint was imposed during the simulation. Phase separation was prohibited by placing an upper limit on the allowed volume states (or a lower limit to the density) that the liquid could sample. (The use of constraints to obtain the thermophysical properties of metastable systems has been used before.29,30,35) The volume constraint allowed us to generate cavities within the liquid beyond the critical size. The severity of the constraint, however, influences the behavior of the work of cavity formation profile as shown in Figure 4. For severely constrained liquids, the work of cavity formation continues to increase with an increase in the cavity radius. For less constrained systems, W continues to increase with an increase in λ, reaches a maximum value, and then decreases for larger cavity sizes. Note that the appearance of a maximum does not violate eqs 4 and 5, since

the scaled particle theory results are valid for unconstrained systems only. The maxima in Figures 2 and 3 appear during the constrained simulations for the following reason. Prior to reaching the maximum, the liquid samples one set of volume states. Upon reaching the maximum and for cavity sizes beyond that, the liquid samples only those volume states near the imposed upper limit. The differences in relevant volume states, or densities, lead to differences in the probabilities of inserting a cavity, i.e., the work of forming a cavity of fixed size decreases as the volume of the system increases. (Note that eqs 4 and 5 assume that all accessible regions of phase space are fully explored.) The work profiles shown in Figures 2 and 3 are for a volume constraint in which the reduced density, F* ) Fσ3, of the liquid was not allowed to decrease below 0.05 of its bulk value. This choice of the constraint was not too severe to prevent the liquid from sampling its own relevant configuration space, while still allowing the density of particles far from the cavity surface to approach the bulk density consistent with the prescribed pressure. Although the height and location of the maxima are arbitrary, Figure 4 shows that, for not too severely constrained systems, the appearance of a maximum is reproducible. We are, however, currently investigating the possibility that another type of constraint (see examples elsewhere29,30,35) may in fact eliminate the appearance of maxima from the work profiles. Additional simulations utilizing different ways to prevent the liquid from phase separating are under development. Surprisingly, the work of cavity formation profiles shown in Figures 2 and 3 behave in a fashion similar to the work of bubble formation profiles predicted by CNT.1 And like CNT, the value of the work of cavity formation at the maximum and the corresponding cavity radius both decrease with a decrease in the pressure (or increasing penetration into the superheated liquid region). The work of forming the critical bubble was obtained using eq 2, while the size of the critical bubble was obtained via the Laplace relation

P′ - P )

2γ r

(11)

where P′ is the pressure of the critical bubble and r is the radius of the bubble. For a bubble whose vapor is ideal (appropriate for the low temperatures studied), P′ was estimated using the following relation1

P′) Pe(T) exp[(P - Pe)v/kT]

(12)

where Pe is the equilibrium pressure at the given temperature and v is the liquid’s molecular volume. Values of the surface tension of the liquid-vapor interface were obtained from the simulation results given by Holcomb et al.36 The simulation and CNT predictions are given in Table 1. The simulation results begin to approach the CNT predictions as the pressure of the liquid increases toward the coexistence pressure. Unlike the CNT results, however, where the work of formation is finite at the spinodal, the work of cavity formation at the maximum appears to vanish at the limit of stability (or at least become very small). The simulation results suggest that the work of cavity formation is able to predict the loss of stability of the superheated liquid. Table 1 also reveals that the work

Ind. Eng. Chem. Res., Vol. 41, No. 5, 2002 1119 Table 1. Comparison of the Size and Reversible Work of Formation of the Critical Cavity and the Critical Bubble Predicted by Classical Nucleation Theorya T*

P*

Fb*

βWCNT

λCNT/σ

βWcavity

λcavity/σ

0.80 0.80 0.80 0.85 0.85 0.85

-0.3757 -0.4680 -0.5291 -0.3679 -0.4058 -0.4367

0.74 0.72 0.70 0.70 0.69 0.68

53.0 34.4 27.0 34.6 28.5 24.7

3.8 3.0 2.7 3.4 3.0 2.8

37 14 6.9 19 12 6.1

2.7 1.7 1.3 2.1 1.6 1.3

WCNT and λCNT are the work of formation and radius of the critical bubble, respectively, as predicted from classical nucleation theory. Wcavity and λcavity are the work of formation and radius of the critical cavity as obtained from the simulations. β ) 1/kT, and Fb* is the bulk density of the metastable liquid. At T* ) 0.80, the spinodal pressure is P* ) -0.6010. At T* ) 0.85, the spinodal pressure is P* ) -0.4828.

Figure 6. Density profile outside a cavity of radius λ in a Lennard-Jones fluid as a function of the distance (R - λ)/σ from the cavity surface. T* ) 0.8, P* ) -0.3757, and Fb* ) 0.74.

Figure 5. Contact density, F+, at the surface of a cavity within a Lennard-Jones fluid as a function of the cavity radius λ for T* ) 0.80 and P* ) -0.3757. The bulk density at this temperature and pressure is Fb* ) 0.74.

of cavity formation is an important component of the total work of forming a critical bubble. In fact, the values of the work of cavity formation are even closer to the work of forming the critical bubble as predicted by DFT.15 The local density of particles surrounding the cavity was also determined for various cavity sizes. The appearance of a drying transition28 is clearly seen in Figures 5 and 6. Figure 5 shows the value of the contact density, F+, versus cavity radius for a liquid maintained at T* ) 0.8 and P* ) -0.3757. The average density within a small bin of 0.01σ from the cavity surface was determined and chosen to represent the value of F+. The depletion of the liquid near the cavity surface is observed for cavity sizes in excess of λ ) 1.8σ, for which F+ rapidly decreases with an increase in the cavity radius. The full density profiles are displayed in Figure 6. The drying transition is again evident by the drastic change in the shape of the density profile for distances close to the cavity surface. Since we were limited to determining the work of cavity formation for radii less than or equal to 3.0σ, we were unable to generate work profiles with maxima for superheated liquids maintained at a positive pressure. To be able to test some of the ideas discussed in section 2 for positive pressure metastable liquids, we instead determined the size of the largest cavity that the liquid could contain without becoming unstable (i.e., the

Figure 7. Maximum cavity radius, λc, that an unconstrained metastable Lennard-Jones fluid could contain at various pressures without phase separating. The arrows indicate the limiting pressures at the spinodal for each of the given temperatures.

critical size beyond which the liquid must be constrained in order to prevent phase separation). Simulations were again carried out in the NPT ensemble, but now a cavity of fixed radius was placed at the center of the simulation box. If the density of the liquid decreased beyond 85% of the bulk density, Fb, within a period of 100 000 MC cycles, the simulation was terminated and the system was considered to be unstable. After some trial and error, the critical cavity radius, λc, was determined to within (0.1σ. Though the value of the critical cavity radius is dependent upon the chosen length of the simulation, the qualitative behavior of the estimated radius as a function of the system pressure did not vary with the length of the MC simulation. Figure 7 shows a plot of the system pressure versus the inverse critical cavity radius λc along five separate isotherms of T* ) 0.75, 0.8, 0.85, 1.1, and 1.13. From this figure, we see that as the pressure decreases, the radius of the maximum allowed cavity decreases as well. This trend is observed not only at negative pressures but also at positive pressures. This suggests that the mechanism of bubble nucleation is similar in both positive and negative pressure liquids. We can also

1120

Ind. Eng. Chem. Res., Vol. 41, No. 5, 2002

conclude that both positive and negative pressure metastable liquids become unstable in the presence of large cavities. Figure 7 also indicates that the size of the critical cavity goes to zero as the pressure of the liquid approaches the spinodal pressure. In fact, any unconstrained simulation performed for conditions right at the spinodal would be considered unstable by the above criteria, predicting a critical cavity size of zero at the spinodal. Thus, we again have an indication that cavity formation is able to predict loss of stability (λc f ∞ at the coexistence curve and λc f 0 at the spinodal), but in a way different from the predictions of DFT. DFT predicts that the work of forming the critical bubble vanishes at the spinodal, while the radius of the critical bubble, unlike the predictions of CNT, diverges at the spinodal.1,11 Our results, on the other hand, indicate that the critical cavity size vanishes at the spinodal. These two results, however, are not at odds with one another. The differences between a cavity (devoid of particle centers) and a bubble (containing particle centers) should lead to different behaviors at the spinodal. DFT predicts a very diffuse interface between the bubble and the surrounding liquid as the spinodal is approached, where fluctuations become infinite in size and the phase transition becomes spontaneous. And since the mechanism of phase separation is no longer an activated process at the spinodal, the liquid should not be able to accommodate a cavity of any size for any appreciable length of time. Finally, we analyzed the results of Figure 7 utilizing a scaling parameter suggested by Shen and Debenedetti15 in their DFT study of homogeneous bubble nucleation in the stretched Lennard-Jones liquid. Shen and Debenedetti showed that various predictions of DFT, including the size of the critical bubble, exhibit scaling behavior; calculations at different temperatures collapsed onto a single curve when the extent of penetration into the coexistence region was normalized by the following ratio

µliq - µsat(T) ∆µ ) ∆µspin µspin(T) - µsat(T)

(13)

where µliq is the chemical potential of the liquid, µspin(T) is the chemical potential at the liquid spinodal at the given T, and µsat(T) is the chemical potential of the saturated liquid. ∆µ/∆µspin goes to zero at the coexistence curve and approaches unity at the spinodal. ∆µ/∆µspin was found to be a natural parameter in the correlation of data for bubble nucleation. We rescaled the size of the critical cavity using the above ratio of chemical potentials for the lowest three temperatures studied (T* ) 0.75, 0.8, and 0.85). The chemical potentials were estimated using the equation of state of Johnson et al.37 The chemical potentials obtained from this equation of state are not accurate for temperatures close to the critical temperature but are reasonably accurate for systems far from the critical point. The results are shown in Figure 8. Interestingly, the critical cavity sizes plotted versus ∆µ/∆µspin for the three lowest temperatures do appear to be, within the error bars, consistent with the collapse of the data onto one curve. The use of the natural parameter ∆µ/∆µspin followed from an analysis of the formation of bubbles within the superheated liquid. The scaling of the size of the critical parameter with the extent of penetration into the

Figure 8. Maximum cavity radius, λc, that an unconstrained metastable Lennard-Jones fluid could contain without phase separating as a function of the ratio ∆µ/∆µspin. The values of T* are 0.75 (star), 0.80 (circles), and 0.85 (triangles).

superheated region again suggests that cavity formation is an important component of the process of bubble nucleation. Additional data and more accurate chemical potential estimations, however, are needed to verify the trend exhibited by Figure 8. 5. Conclusions Bubble nucleation was studied from the viewpoint of understanding cavity formation within superheated liquids. The reversible work of forming cavities of various sizes was determined for liquids under tension at various temperatures. A critical cavity size was found to occur, beyond which the liquid would phase separate if not for the imposition of a volume constraint during the simulations. The radius of this critical cavity was found to decrease with an increase in the extent of penetration into the metastable region of the phase diagram, approaching a zero radius at the spinodal. Unlike stable liquids, the presence of large cavities destabilizes metastable liquids, leading to further spontaneous growth of the cavity, i.e., phase separation. The size of this critical cavity obtained at several temperatures also scaled with the ratio ∆µ/∆µspin, suggested by Shen and Debenedetti15 to be a natural parameter to correlate data dealing with the problem of bubble nucleation. To determine the work of forming cavities beyond the critical cavity size, a volume constraint was imposed during the simulations. For at least negative pressure liquids, the work of cavity formation eventually displayed a maximum, decreasing for larger cavity radii. The dependence of the work of cavity formation and radius at the maximum on the temperature and pressure of the metastable liquid is qualitatively similar to that of a critical bubble. The simulations results become comparable to the CNT predictions as the liquid approaches the coexistence curve, indicating that the work of cavity formation is an important contribution to the work of forming the critical bubble. The effect of imposing different constraints, however, on the work of cavity formation in superheated liquids needs to be further investigated and is the focus of ongoing work. Overall, our simulations suggest that cavity formation is an important part of the process of bubble nucleation

Ind. Eng. Chem. Res., Vol. 41, No. 5, 2002 1121

in metastable liquids (at either positive or negative pressure). Whether cavity formation plays a dominant role in the mechanism of bubble nucleation remains to be seen. A bubble, which is characterized by a region of very low density, is similar to a cavity, a region completely empty of particle centers. Yet, the size of a bubble needed to contain even a few particles, such that the density of the bubble is near that of the corresponding vapor, is large. Thus, for highly superheated liquids, bubble nucleation may proceed via the formation of a small cavity that first serves to destabilize the liquid. As the cavity grows in size, particles may diffuse into the cavity, thereby forming a bubble. Afterward, bubble formation and growth would be driving the phase transition. This scenario may not be relevant for moderately to slightly superheated liquids, where the estimated critical cavity size is quite large. In this case, a bubble may form that could destabilize the liquid before the critical cavity size is reached, or even serve to stabilize the critical cavity. Nevertheless, the work of cavity formation would still be a large contribution to the overall work of forming the critical bubble. Additional computational studies are needed to investigate the importance of cavity formation and growth in metastable liquids. In the end, the focus on cavity growth should lead to an improved understanding of the molecular mechanism of bubble nucleation in superheated liquids. Acknowledgment This work was supported by a grant from the Purdue Research Foundation. Literature Cited (1) Debenedetti, P. G. Metastable Liquids; Princeton University Press: Princeton, New Jersey, 1996. (2) Skripov, V. P. Metastable Liquids; John Wiley & Sons: New York, 1974. (3) 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 Science Publishers: New York, 1988. (4) Gardellin, D. Minimize Valve Cavitation by Proper System Design. Chem. Eng. Prog. 1996, 92, 52. (5) Chen, Y. L.; Israelachvili, J. New Mechanism of Cavitation Damage. Science 1991, 252, 1157. (6) Putterman, S. J. Sonoluminescence: Sound into Light. Sci. Am. 1995, 272, 46. (7) Pandit, A. B.; Moholkar, V. S. Harness Cavitation to Improve Processing. Chem. Eng. Prog. 1996, 92, 57. (8) Abe, T.; Messner, W. C.; Reed, M. L. Effects of Elevated Temperature Treatments in Microstructure Release Procedures. J. Microelectromech. Sys. 1995, 4, 66. (9) Song, K. H.; Xu, X. Explosive Phase Transformation in Excimer Laser Ablation. Appl. Surf. Sci. 1998, 127, 111. (10) Kiselev, S. B. Kinetic Boundary of Metastable States in Superheated and Stretched Liquids. Phys. A 1999, 269, 252. (11) Oxtoby, D. W. Homogeneous Nucleation: Theory and Experiment. J. Phys.: Condens. Matter 1992, 4, 7627. (12) Oxtoby, D. W.; Evans, R. Nonclassical Nucleation Theory for the Gas-Liquid Transition. J. Chem. Phys. 1988, 89, 7521. (13) Zeng, X. C.; Oxtoby, D. W. Gas-Liquid Nucleation in Lennard-Jones Fluids. J. Chem. Phys. 1991, 94, 4472. (14) Talanquer, V.; Oxtoby, D. W. Nucleation of Bubbles in Binary Fluids. J. Chem. Phys. 1995, 102, 2156. (15) Shen, V. K.; Debenedetti, P. G. Density-Functional Study of Homogeneous Bubble Nucleation in the Stretched LennardJones Fluid. J. Chem. Phys. 2001, 114, 4149.

(16) Shen, V. K.; Debenedetti, P. G. A Computational Study of Homogeneous Liquid-Vapor Nucleation in the Lennard-Jones Fluid. J. Chem. Phys. 1999, 111, 3581. (17) Kusaka, I.; Oxtoby, D. W. Identifying Physical Clusters in Bubble Nucleation. J. Chem. Phys. 1999, 111, 1104. (18) Kinjo, T.; Matsumoto, M. Cavitation Processes and Negative Pressure. Fluid Phase Equil. 1998, 144, 343. (19) Schaaf, P.; Senger, B.; Reiss, H. Defining Physical Clusters in Nucleation Theory From the N-Particle Distribution Function. J. Phys. Chem. B 1997, 101, 8740. (20) Senger, B.; Schaaf, P.; Corti, D. S.; Bowles, R.; Voegel, J.C.; Reiss, H. A Molecular Theory of the Homogeneous Nucleation Rate. I. Formulation and Fundamental Issues. J. Chem. Phys. 1999, 110, 6421. (21) Senger, B.; Schaaf, P.; Corti, D. S.; Bowles, R.; Pointu, D.; Voegel, J.-C.; Reiss, H. A Molecular Theory of the Homogeneous Nucleation Rate. II. Application to Argon Vapor. J. Chem. Phys. 1999, 110, 6438. (22) Schaaf, P.; Senger, B.; Voegel, J.-C.; Reiss, H. Extended (n/v)-Stillinger Cluster for Use in the Theory of Homogeneous Nucleation. Phys. Rev. E 1999, 60, 771. (23) Reiss, H.; Frisch, H. L.; Lebowitz, J. L. Statistical Mechanics of Rigid Spheres. J. Chem. Phys. 1959, 31, 369. (24) Stecki, J.; Toxvaerd, S. Hard Sphere Cavity in a LennardJones Liquid. J. Chem. Phys. 1990, 93, 7342. (25) Irisa, M.; Nagayama, K.; Hirata, F. Extended Scaled Particle Theory for Dilute Solutions of Arbitrary Shaped Solutes. An Application to Solvation Free Energies of Hydrocarbons. Chem. Phys. Lett. 1993, 207, 430. (26) Floris, F. M.; Selmi, M.; Tani, A.; Tomasi, J. Free Energy and Entropy for Inserting Cavities in Water: Comparison of Monte Carlo Simulation and Scaled Particle Theory Results. J. Chem. Phys. 1997, 107, 6353. (27) Lum, K.; Chandler, D.; Weeks, J. D. Hydrophobicity at Small and Large Length Scales. J. Phys. Chem. B 1999, 103, 4570. (28) Huang, D. M.; Chandler, D. Cavity Formation and the Drying Transition in the Lennard-Jones Fluid. Phys. Rev. E 2000, 61, 1501. (29) Corti, D. S.; Debenedetti, P. G. Metastability and Constraints: A Study of the Superheated Lennard-Jones Liquid in the Void-Constrained Ensemble. Ind. Eng. Chem. Res. 1995, 34, 3573. (30) Corti, D. S.; Debenedetti, P. G.; Sastry, S.; Stillinger, F. H. Constraints, Metastability, and Inherent Structures in Liquids. Phys. Rev. E 1997, 55, 5522. (31) Sastry, S.; Debenedetti, P. G.; Stillinger, F. H. Statistical Geometry of Particle Packings. II. “Weak Spots” in Liquids. Phys. Rev. E 1997, 56, 5533. (32) Vishnyakov, A.; Debenedetti, P. G.; Neimark, A. V. Statistical Geometry of Cavities in a Metastable Confined Fluid. Phys. Rev. E 2000, 62, 538. (33) Reiss, H. Scaled Particle Theory of Hard Sphere Fluids to 1976. In Statistical Mechanics and Statistical Methods in Theory and Application; Landman, U., Ed.; Plenum: New York, 1977. (34) Stillinger, F. H. Structure in Aqueous Solutions of Nonpolar Solutes from the Standpoint of Scaled-Particle Theory. J. Solution Chem. 1973, 2, 141. (35) Corti, D. S.; Debenedetti, P. G. A Computational Study of Metastability in Vapor-Liquid Equilibrium. Chem. Eng. Sci. 1994, 49, 2717. (36) Holcomb, C. D.; Clancy, P.; Thompson, S. M.; Zollweg J. A. A critical study of simulations of the Lennard-Jones liquidvapor interface. Fluid Phase Equil. 1992, 75, 185. (37) Johnson, J. K.; Zollweg, J. A.; Gubbins, K. E. The LennardJones equation of state revisited. Mol. Phys. 1993, 78, 591. (38) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford Science Publications: Oxford, 1987.

Received for review June 28, 2001 Revised manuscript received August 13, 2001 Accepted August 13, 2001 IE010554Q