Fluid Phase Equilibria 144 Ž1998. 343–350 Cavitation processes and negative pressure Tomoyuki Kinjo ) , Mitsuhiro Matsumoto Department of Applied Physics, Nagoya UniÕersity, Nagoya 464-01, Japan Received 13 January 1997; accepted 4 July 1997 Abstract Molecular dynamics simulations of a Lennard–Jones fluid under negative pressure are carried out to investigate the inception of vapor phase in the liquid Ži.e., cavitation. and properties of liquid under negative pressure. The pressure of the system is reduced by uniform system expansion. In the case of larger expansion ratios, the nuclei of bubbles immediately appear, whereas finite waiting time is observed in case of smaller expansion ratio. This means the system reaches a kinetic limit of metastability. The nucleation rate calculated with the waiting time is 8 orders of magnitude larger than that of classical nucleation theory. Isothermal compressibility is also obtained and its density dependence is discussed. q 1998 Elsevier Science B.V. Keywords: Molecular dynamics simulation; Cavitation; Bubble nucleation; Negative pressure 1. Introduction Cavitation is a subject of interest in a number of fields in science and engineering. For example, controlling the cavitation process is relevant to designing a high-performance ship propeller. From a scientific point of view, the cavitation is understood as a typical example of nucleation processes, and is related to thermodynamic properties and the limit of metastable liquid. From a microscopic viewpoint, despite this importance, less attention has been paid to cavitation and properties of liquid under negative pressure. Let us consider a typical equation of state having a van der Waals loop, the isotherm of which is shown in Fig. 1a. Metastable states of liquid, e.g., C, are the states between B and D along the isotherm. The liquid at a stable point A can be expanded beyond the state B which corresponds to boiling point, but can never go over the state D without phase transition. In a state Že.g., E. between D and F, isothermal compressibility is negative Ž k T s yÕy1 Ž E ÕrE P . T - 0., and thus the system within ) Corresponding author. 0378-3812r98r$19.00 q 1998 Elsevier Science B.V. All rights reserved. PII S 0 3 7 8 - 3 8 1 2 Ž 9 7 . 0 0 2 7 8 - 1 T. Kinjo, M. Matsumotor Fluid Phase Equilibria 144 (1998) 343–350 344 Fig. 1. Isotherm at temperature 0.8 Ža. and spinodal curve Žb. of the van der Waals EOS. The units are reduced by the critical point values. B and G obtained by Maxwell’s equal area rule are on the coexistence curve. D and F are on the spinodal line, i.e., on the boundary of metastable state. the region is inherently unstable. These state points are projected to the T–P plane in Fig. 1b where spinodal curve, or the stability limit, QP is also shown. In principle, a liquid quenched to a metastable state undergoes a phase transition. In practice, however, shallow quenching is not sufficient to cause phase transition, because the work required for the creation of surface is larger than the reduction of the free energy of the phase change. For example, it is well known that liquids can withstand modest tensile force, or negative pressure. In the case of liquid-to-vapor phase transition, if the quenching is sufficiently deep, the liquid undergoes explosive phase transition to vapor, which is called vapor explosion for the case of positive pressure, or cavitation for the case of negative pressure w1–3x. The purpose of this study is to investigate the properties of metastable liquid, and to investigate the intrinsic molecular dynamics of phase change. With this objective, we have carried out molecular dynamics Ž MD. simulation of a Lennard–Jones fluid for fairly large systems, and have observed the inception of vapor phase Ž or nuclei of bubbles. in the bulk liquid phase. 2. Simulation methods Microcanonical MD simulations are carried out. The system treated here consists of 10,976 particles. These particles interact each other via the Lennard–Jones potential: s 12 s 6 f LJ Ž r . s 4e y Ž1. r r We perform our calculation with reduced units, i.e., length in s , energy e , and mass m Žmass of a particle.. Thus, the time unit is t s m s 2re , and the pressure unit is esy3. For reference, the unit value for fluid argon are s s 3.4 = 10y10 m, erk B s 120 K, and t s 2.14 ps. The value of the pressure unit is 42.2 MPa. The time step of the simulation is D t s 0.0025, and the cutoff length is 4.0 in these units. Long range corrections are made for potential energy and pressure, assuming that the liquid structure is uniform beyond the cutoff radius. The initial temperature of the system is 0.8. The equilibrated initial configurations are expanded uniformly and instantaneously with various expansion ratios. After this adiabatic expansion, we continue microcanonical MD simulation. ž / ž / ( T. Kinjo, M. Matsumotor Fluid Phase Equilibria 144 (1998) 343–350 345 We do not control the system temperature. In the case of vapor phase nucleation Ž cluster formation., the temperature of the generated clusters is higher than the surrounding vapor due to the latent heat release. In the bubble formation case, however, the temperature of generated bubbles Ž or vapor inside. is not important. The temperature of surrounding liquid is slightly raised because the liquid is contracted during the bubble growth, but this small temperature increase affects little the dynamics of bubble formation. 3. Results 3.1. Limit of metastability We first monitored the time development of the system pressure, and found quite different behavior for fairly small change of the density. Fig. 2 shows three types of pressure change corresponding to metastable state Ž r s 0.708., limit of metastability Ž r s 0.694., and unstable state Ž r s 0.690.. In the case of r s 0.690, pressure immediately increases toward positive value. This increase implies that bubbles appear in the liquid. The case of r s 0.694 is more interesting, in which pressure is maintained in the first stage, and then increases in the second stage. After such waiting time, one bubble Žin this case. appears. The corresponding snapshot is shown in Fig. 3. If a system has higher density than these, we don’t observe bubble formation at least in our simulation time Ž 200t .. Such distinction is sufficiently clear to show the kinetic limit points on phase diagram. We use various initial configurations, and obtain several points of kinetic limit of metastability as shown in Fig. 4. Several empirical equations of state Ž EOS. are known about the Lennard–Jones fluid system, with which we can calculate the spinodal line. In Fig. 4, the line based on EOS of Nicolas et al. w4x is shown; our MD points are near but above the line. More recently obtained EOS Že.g., w5x. gives a similar result at low temperature region. The difference between EOS and our MD simulation may be important, but considering the accuracy of EOS near the spinodal region, we do not discuss it here. Fig. 2. The time developments of the system pressure. T. Kinjo, M. Matsumotor Fluid Phase Equilibria 144 (1998) 343–350 346 Fig. 3. Snapshot of a bubble Žleft. and the corresponding equi-potential surface Žright.. At this kinetic limit, a rapid growth of a bubble after the waiting time is observed, which means the size of the bubble becomes larger than critical one. Hence, one can obtain nucleation rate Jsim from the waiting time tw as follows: 1 Jsim s Ž2. tw V On the other hand, a classical nucleation theory w6x predicts: ( Jclass s r l 2g pm exp y16pg 3 3 k BT Ž P l y Pv . 2 Ž3. With parameters of argon fluid at corresponding temperature, calculated values are Jsim , 4.23 = 10 21 Fig. 4. Kinetic limit of metastability plotted on phase diagram. T. Kinjo, M. Matsumotor Fluid Phase Equilibria 144 (1998) 343–350 347 Fig. 5. The histogram of the local potential energy. Its minimum position y2.5 is used for the definition of a bubble. cmy3 sy1 and Jclass , 1.53 = 10 29 cmy3 sy1. Thus, Jsim is eight orders of magnitude lager than Jclass . Such underestimation of the rate of the classical nucleation theory has been previously pointed out by Zeng and Oxtby w7x, who state that the free energy barrier of bubble formation remains finite near the spinodal region in the approximations of the classical theory w8x. 3.2. Growth of bubble The definition of a bubble is not clear in microscopic scale. We define a ‘bubble’ by calculating the local potential energy as follows. Let us introduce a ghost particle which interact with real particles via the attractive part of Lennard–Jones potential: fatt Ž r . s ½ f LJ Ž r . for r ) 2 1r6 ye for r - 2 1r6 Ž4. We insert the ghost particle in the simulation cell and calculate its potential energy, which we regard Fig. 6. The growth of bubble at a kinetic limit point Ž r s 0.694.. 348 T. Kinjo, M. Matsumotor Fluid Phase Equilibria 144 (1998) 343–350 Fig. 7. Phase diagram and simulated state points. as local potential energy at the position. The histogram of the local potential energy is shown in Fig. 5, from which a threshold of y2.5 seems appropriate to define the bubble region. An equi-potential surface is shown in Fig. 3 Žright. , which captures well the shape of the bubble. The volume of a bubble is plotted against time in Fig. 6, which shows the growth of the bubble. We clearly distinguish two stages, i.e., growth–shrinkage period and explosive growth period. In the classical picture, the bubble in former period seems to be sub-critical and undergo fluctuations. 3.3. Radius distribution function and isothermal compressibility In this section, we briefly investigate the properties of metastable liquid. Four liquid samples Ž three in metastable state and one in normal state. are used, which are shown in Fig. 7. In the deepest quenched state d, a bubble is formed after a finite waiting time; the analysis below is done for the uniform metastable liquid sample. The radial distribution function ŽRDF. is shown in Fig. 8. Although the degree of metastability is quite different, their RDF is almost the same, including that of normal liquid. However, a systematic shift is observed; the overall shape moves to larger r. Fig. 8. Radial distribution function. T. Kinjo, M. Matsumotor Fluid Phase Equilibria 144 (1998) 343–350 349 Fig. 9. Density dependence of isothermal compressibility. To evaluate the difference quantitatively, we introduce a residual function D nŽ r . defined as follows: r H0 D nŽ r . s r g Ž r X . y 1 4p r X 2 d r X Ž5. If we use grand canonical ensemble, D nŽ r . is related to the isothermal compressibility k T as: ` H0 r k BTk T s 1 q r g Ž r X . y 1 4p r X 2 d r X s 1 q D n Ž q` . Ž6. Since our simulation is done in microcanonical ensemble and the system is finite, it is difficult to take the limit of calculated D nŽ r . to evaluate k T . However, D nŽ r . is saturated at large r , 8, we estimate k T with this plateau value. Fig. 9 shows the density dependence of k T . It can be seen that the isothermal compressibility increases with approaching to kinetic limit of metastability; note that the temperature of each state is slightly different. 4. Conclusion In this paper, we found, using MD simulation, the kinetic stability limit of liquid and observed bubble formation at this limit. From the waiting time for a single bubble to form, nucleation rate is calculated, which is much larger than the rate of classical nucleation theory. One of the theoretical defects, i.e., remaining free energy barrier near spinodal region, seems the cause of this underestimation. Bubble growth has two stages, growth and shrinkage, and explosive growth. In the former stage, the system is in fragile equilibrium, which is actually broken up to bubble formation. RDFs of various states are almost the same, but the isothermal compressibility k T has a clear density dependence. Approaching the spinodal, k T increases. We expect that kinetic instability is T. Kinjo, M. Matsumotor Fluid Phase Equilibria 144 (1998) 343–350 350 reached before thermodynamic instability k T ™ `, but accurate calculation from simulation data is difficult at this stage. 5. List of symbols gŽ r. m kB P Pv Pl T V Õ g D nŽ r . r rl Radial distribution function Molecular mass Boltzmann constant Pressure Vapor pressure Pressure on external liquid phase Temperature Volume of the simulation cell Molecular volume Surface tension Residual function, defined in Eq. Ž5. Total number density Number density of liquid phase References w1x w2x w3x w4x w5x w6x w7x w8x C.T. Avedsian, J. Phys. Chem. Ref. Data 14 Ž1985. 695–729. J.P. Heath, G.M. Pound, Condensation and Evaporation–Nucleation and Growth Kinetics, Pergamon, Oxford, 1962. D.H. Trevena, Cavitation and Tension in Liquids, Adam Hilger, Bristol, 1987. J.J. Nicolas, K.E. Gubbins, W.B. Streett, D.J. Tildesley, Mol. Phys. 37 Ž1979. 1429–1454. J.K. Johnson, J.A. Zollweg, K.E. Gubbins, Mol. Phys. 78 Ž1993. 591–618. M. Blander, J.L. Katz, AIChE J. 21 Ž1975. 833–848. X.C. Zeng, D.W. Oxtby, J. Chem. Phys. 94 Ž1991. 4472–4478. J.W. Cahn, J.E. Hilliad, J. Chem. Phys. 28 Ž1968. 258–267.