Cavitation and Negative Pressure : A Flexible Water Model Molecular Dynamics Simulation

The critical negative pressure for cavitation in water has been theoretically predicted to be in the range of -100 to -200 MPa at room temperature, whereas values around -30 MPa have been obtained by many experiments. The discrepancy has yet to be resolved. Molecular dynamics (MD) is an effective method of observing bubble nucleation, however, most MD simulations use a rigid water model and do not take the effects of intermolecular vibrations into account. In this manuscript we perform MD simulations to study cavitation in water by using a TIP4P/2005f model under volumecontrolled stretching. It is found that the critical negative pressure of water was -168 MPa in the simulation and the critical negative pressure of water containing 50 oxygen molecules was -150 MPa. Hydrogen bonds played a major role in the cavitation process: the breaking of hydrogen bonds promoted bubble generation and growth. The O-H bond could release energy to increase the amount of potential energy in the system, so that cavitation was more likely to occur. When cavitation occurred, the O-H bond could absorb energy to reduce the amount of potential energy in the system, which will promote the growth of bubbles, and stabilise the cavitation bubbles.


Introduction
A liquid can withstand a certain negative pressure in a metastable state until the occurrence of cavitation (Frdric Caupin, 2006), this phenomenon is also known as bubble nucleation.Bubble nucleation can be divided into homogeneous nucleation and heterogeneous nucleation: nucleation in pure liquids is called homogeneous nucleation, nucleation on the surface of impurities or containers is called heterogeneous nucleation (Cochard. H,2013).
Bubble nucleation is the initial stage of the liquid-vapour phase transition, involving the most basic changes in nature, and it is a type of physical phenomenon in that is common as a subject of scientific, and engineering, research.The hazards of cavitation are very serious, and it has been found that cavitation begins with damage to the ship's propeller blades.After the object is exposed to cavitation in liquid motion, the surface will be deformed and the material will be ablated.During the cavitation process, the cavitation rapidly generates, expands, and collapses, forming a shock wave or a high-speed microjet in the liquid.After the metal material is impacted, the surface crystal structure is distorted, chemical instability occurs, and adjacent crystal grains have different potentials, thereby accelerating the electrochemical corrosion process.The mechanical properties of the material in the ablated area deteriorate significantly, resulting in a sharp increase in the amount of cavitation.
Classical nucleation theory (CNT) posits that molecules in a metastable liquid under irregular motion generate density fluctuations.This will result in the formation of areas of lower density within the liquid: the so-called embryo nuclei.When the liquid reaches a certain degree of superheating, the embryo nuclei whose size is smaller than the critical nucleation size disappear, while the embryo nuclei whose size is larger than the critical nucleation size will gradually grow to form macroscopic bubbles.Classical nucleation theory predicts the critical negative pressure for cavitation to be around -140 MPa (Caupin,2015)or -160 MPa (Zheng. Q,1991)around room temperature (300 K).Experimentally, a critical negative pressure of -140 MPa (Alvarenga. A. D. , 1993)was reached at 315 K using the microscopic Berthelot tube method (Azouzi.M. E. M. , 2012); however, in cavitation studies using other experimental methods, such as acoustic cavitation (Arvengas. A. , 2011), shock wave techniques (Quinto-Su, 2013), and the method of artificial trees (Wheeler. T. D. , 2008), the critical negative pressure was about -30 MPa (Frdric Caupin, 2015), which was different from that predicted by classical nucleation theory.
In general, time-scale during bubble nucleation is of the order of nanosecond-microsecond, and the size of vapour cavities is the measured on a scale of nanometre-micron.So far, experiment can observe the phenomenon of micron-scale, but can not directly observe the nano-scale bubble.So new research methods are needed to understand bubble nucleation: molecular dynamics (MD) is the most direct and effective method with which to observe bubble nucleation.Baidakov (Baidakov. V. G. , 2014) studied the homogeneous bubble nucleation process in a Lennard-Jones (L-J) liquid(Ar) under a negative pressure and obtained the relevant parameters in the process of bubble nucleation.The values of J obtained in MD simulation are by 8-20 orders higher than those predicted by CNT.Classical nucleation theory was then modified (MCNT) by incorporating the van der Waals-Cahn-Hilliard gradient theory, explaining the difference between classical nucleation theory and simulated results (Baidakov. V. G. ,2016), but Ar as a non-polar molecule did not truly reflect the nucleation of bubbles in liquid water.Yijin (Mao. Y. ,2013) used the nonequilibrium molecular dynamics (NEMD) method to simulate the process of growth and disruption of nano-bubbles in water (the TIP3P model) and found that the Rayleigh-Plesset (R-P) equation cannot effectively predict the growth of nanobubbles by comparing data with the results of the R-P equation.Wang (Wang. P. , 2017) studied the homogeneous bubble nucleation process in liquid water (TIP4P/2005 model) under negative pressure cavitation, the MD simulation predicts a critical negative pressure for cavitation at -189 MPa under volume-controlled stretching at 300 K, so Wang modified the classic nucleation model, which provided a possible view of the difference between classical nucleation theory and the experiment.The differences between MD results and the tests described above are attributed to several possible facts including the size of the simulated tank, the presence of dissolved gases, or the level of detail inherent to the water model itself.
Okamoto (Okamoto. R. , 2015) studied the formation of nanobubbles in water containing dissolved oxygen molecules, and calculated the bubble free energy.Takenori (Yamamoto. T. , 2011) used nanoscale molecular dynamics (NAMD) techniques to simulate the interfacial properties of helium nanobubbles in water (SPC/FW model) and obtained the relationship between the surface tension and the bubble radius, however, most MD simulations used rigid water molecules, but intermolecular vibration played an important role in bubble nucleation.Takenori used a flexible water model (SPC/FW model), taking vibrations in water molecules into account, but the water model (SPC/FW model) was a simple point charge model, and there was no TIP4P/2005 water model able to describe the surface tension change from the temperature of the three-phase point to the critical temperature with sufficient accuracy (Vega. C. , 2007).
The above research found the difference between the macroscopic bubble nucleation theory and the results of microscopic bubble nucleation, and further perfected the macroscopic bubble nucleation theory according to the microscopic results.Most of them consider the uncoordination of macroscopic theory and did not consider the correctness of micro bubble nucleation results.The object of the above study is simple rigid water molecules or simple single atoms.Because most materials are composed of polyatomic molecules (like water), where a more complicated force field should be accounted for, the bubble nucleation process in liquid water should be studied.Therefore, the above study investigates the complicated process of bubble nucleation by the result of a simple force field.The correctness of the microscopic bubble nucleation process is questionable.This work aims to improve the water model on the basis of the TIP4P/2005 water molecule model using the flexible water molecular model TIP4P/2005f (GonzáLez. M. A. , 2011), to study the liquid water bubble stretching process during nucleation, and on this basis, to study the interaction between water molecules in the bubble nucleation process and the change of energy in the system, and to explain the genesis of micro bubble nucleation from the point of force and energy.

Simulation
The MD simulation using LAMMPS software for the water molecule model using flexible water molecular model TIP4P/2005f was used: in the TIP4P/2005 model, each water molecule has four interacting positions, in addition to the oxygen and hydrogen atoms, and a massless point M on the middle of the bond with the hydrogen and oxygen atoms.Interactions between water molecules include an L-J term and an electrostatic part.The oxygen site carries no charge, but contributes to the L-J term: In the formula, r OO is the distance between two oxygen atoms, σ OO = 3.1644 and ε =0.1852 kcal/mol.The H and M sites are charged (q H = 0.5564 e and q M = -2q H = -1.1128e),but do not contribute to the L-J term.The electrostatic part of the interaction potential between two water molecules is.
where the summation is taken over all pairs of charged sites, r i j is the distance between two charged sites, and k e is the electrostatic constant.
TIP4P/2005f, as a flexible water molecule model, differs from the rigid water molecule model in that there is no fixed bond length and bond angle in the water molecule.The bond length and bond angle are given by: Where r eq is the equilibrium bond length of the OH bond, r eq = 0.9419Å , θ eq is the equilibrium bond angle of the HOH bond angle, θ eq = 107.4Oxygen atoms are not charged, the interaction between them is also the LJ interaction, Ryckaert. J. P. , 1977).In addition, the L-J potential energy parameter between water molecules and oxygen molecules is determined according to the Lorentz-Berthelot formula (Arora. G. , 2005): The total potential of the system is the sum of the L-J and electrostatic interactions between all molecules.The cut-off distance for L-J and electrostatic interactions is set to 8 Å. Electrostatic interactions are calculated using the particleparticle particle-mesh (PPPM) algorithm in LAMMPS (Shi. B. , 2006).The time integration scheme closely follows the time-reversible measure-preserving Verlet and RESPA integrators derived by Tuckerman et al. (Mayo. S. L. , 1990), with a time step of 1 fs.
The simulated water molecule box uses three-dimensional periodic boundary conditions.The temperature is controlled by using a Nose-Hoover thermostat.Each simulation starts with a constant pressure isotropic relaxation (NPT) under atmospheric pressure (p = 1 bar, T = 298 K) to equilibrate the system for 1 ns.Next, the simulation is switched to NVT ensemble to simulate volume-controlled stretching while the pressure is calculated and may become negative.The corresponding temperature and pressure is calculated:

Hydrogen Bonds
The advantage of molecular simulation is that it can reveal detailed micro-structural information.The radial distribution function can effectively describe the microstructure of the liquid.The coordination number (N) is the average number of other atoms in the spherical shell whose radius is r centred on a given atom (or position).The coordination number could be obtained by integrating the radial distribution function: In general, the position of the first trough of the corresponding RDF represented the first coordination layer of the central component.The first valley of the O-O radial distribution function was 3.36 Å, the first peak was 0.75, and the coordination number was 4.56, which means that there were about four oxygen atoms that could be coordinated around each oxygen atom.
In this simulation system, there were four oxygen atoms around each water molecule, and the distance between two water molecules Roo was 2.77 Å, which is similar to the structure of microscopic hydrogen bonds in water.Therefore, it could be said that the liquid water constructed by the TIP4P/2005f model can reflect the microstructure of the hydrogen bonds in the real liquid water.The analysis showed that the hydrogen bond mainly came from the electrostatic interaction.

Pressure Changes in the TIP4P/2005f Model After Stretching
Three stages of relaxation-stretching-constant temperature were simulated for the TIP4P/2005f model, the initial structure contained 5000 water molecules, the initial density was 1.00035 g/cm 3 , the box size was 86.47 Å × 86.47 Å × 20 Å.After relaxation, the density was 0.976 g/cm 3 .After three different stretches, the density was 0.898 g/cm 3 , 0.873 g/cm 3 , and 0.852 g/cm 3 .
As can be seen from Fig. 2, when the liquid water was stretched to ρ = 0.873 g/cm 3 , firstly, the pressure in the liquid was maintained at 850 ps when P = -168 MPa, then, the system pressure began to rise, the final system pressure was maintained at -35 MPa.The simulation results showed that, when the fluid was subjected to a negative pressure of -168 MPa, the liquid began to bubble after having remained in a metastable state for a period of time.When the degree of fluid stretching was slightly increased, the density of the fluid was reduced to 0.8525 g/cm 3 , the time for which the system remained in a metastable state was short (60 ps), then the system pressure increased rapidly.When the degree of fluid stretching was slightly smaller and the fluid was stretched to ρ = 0.898 g/cm 3 , the system was running for 1000 ps at a pressure of -135 MPa, the pressure did not show a significant increase in this situation.Therefore, -168 MPa was deemed to have been the critical negative pressure for liquid water at 298.0 K within 1000 ps in this simulation.
Here, the critical negative pressure obtained by the TIP4P/2005f model was about -168 MPa through MD simulation, compared with the result of Wang(Wang.P. , 2017) who used the TIP4P/2005 water model (-189 MPa): the result obtained in the present research was closer to most of the tests, however, the critical negative pressure obtained from this simulation that used the TIP4P/2005f model were also different from the experimental results, and it was also different from that predicted by classical nucleation theory.For the existence of the above differences, the classical nucleation theory involved a macroscopic liquid and bubbles which was inapplicable to the microscopic scale, so the macroscopic state parameters were not necessarily applicable in microcosm.In addition, comparing with MD simulation, the liquid water in the experiment might have some weak points (such as dissolved gases or bubbles therein), while the liquid water in the simulation was pure water and considered to have had no weak points.
As shown in Fig. 2, there were three stages of system pressure change during the evolution of liquid water cavitation: the first stage was the liquid stretched into a metastable state, the pressure remained unchanged; the second stage was the beginning of the formation of bubbles inside the liquid, the pressure quickly rose to the equilibrium value; in the third stage, the pressure in the liquid remained at its equilibrium value.
Classical nucleation theory posits that the formation of bubbles in the liquid phase needs a certain induction time tw, the pressure during this process remains stable.When gasification cores appear in the liquid phase, the system enters the rapid nucleation stage, then the pressure increases rapidly when the system generates many small gas nuclei.In the late period of nucleation, the small gas nuclei gather to form a larger gas core, the system pressure slowly increases to the equilibrium value.The nucleation rate is the number of critical nuclei produced per unit volume per unit time.The rate of J sim can be obtained from the induction time (Kinjo. T. , 1998): The induction times for cavitation were: 850 ps at ρ = 0.873 g/cm 3 and 60 ps at ρ = 0.8525 g/cm 3 .Thus, at ρ = 0.873 g/cm 3 , J sim = 6.8 × 10 33 m −3 s −1 , while J sim = 9.49 × 10 34 m −3 s −1 at ρ = 0.8525 g/cm 3 .

Force Parameters
Fig. 4(A) shows that the energy of van der Waals interaction was about 10,050 kcal/mol in the system when ρ was 0.898 g/cm 3 .At ρ = 0.873g/cm 3 , the energy of vdW interaction was maintained at the equilibrium value of about 98,575 kcal/mol in the system.After 850 ps, the energy rapidly increased to 10,350 kcal/mol.As seen in Fig. 4(B), at ρ = 0.898 g/cm 3 , the energy of the electrostatic interaction was maintained at the equilibrium value of about 519,250 kcal/mol in the system: at ρ = 0.873 g/cm 3 , the energy of the electrostatic interaction was maintained at the equilibrium value of about 516,700 kcal/mol.After 850 ps, the energy of the electrostatic interaction rapidly decreased to about 515,500 kcal/mol.It could be found from comparison of electrostatic interaction and van der Waals interaction that the electrostatic interaction in the system was 50 times that of van der Waals interaction, which showed that the electrostatic interactions played a major role in the confined water molecules: in other words, the hydrogen bonds played a major role in the cavitation process.
From all the above phenomena, the following conclusions can be obtained: when the liquid water was stretched from a normal, to a metastable, state, with increased stretching, hydrogen bonds in the water continued to break, and the intensity of the hydrogen bonding was gradually reduced, thereby contributing to the generation of bubbles; when cavitation occurred, the hydrogen bonds in water were further broken, the hydrogen bond strength continued to decrease, thereby promoting the growth of bubbles in liquid water.

Energy Parameters
Fig. 5(A) indicated that when flexible water model was stretched to 0.898 g/cm 3 , the potential energy inside the system was maintained at the equilibrium value of -47,850 kcal/mol, and when the flexible water model was stretched to 0.873 g/cm 3 , the potential energy inside the system was maintained at its equilibrium value of -47,520 kcal/mol for some time, after 850 ps, the system potential energy decreased rapidly, and finally, the potential energy of the system was balanced at -48,100 kcal/mol.) shows that, when the flexible water model was stretched to 0.898 g/cm 3 and 0.873 g/cm 3 respectively, the kinetic energy in both systems was maintained at the equilibrium value of 13,320 kcal/mol.Fig. 6(A) shows that the energy of O-H bond could be maintained at the equilibrium value of about 6050 kcal/mol in the system when liquid water was stretched to 0.898 g/cm 3 .When liquid water was stretched to 0.873 g/cm 3 , the energy of O-H bond could be maintained at about 6000 kcal/mol.After 850 ps, the bond energy in the system increased rapidly.Finally, the bond energy was maintained at about 6100 kcal/mol.
As shown in Fig. 6(B), the energy of the H-O-H angle of the system could be maintained at the equilibrium value of 1920 kcal/mol when liquid water was stretched to 0.898 g/cm 3 .When ρ = 0.873 g/cm 3 , the energy of H-O-H angle was always about 1920 kcal/mol, even though cavitation occurred at 850 ps, the angle energy did not change thereafter.
In summary, there was a difference in correlation energy parameters of the different metastable liquids before bubble nucleation.The angle energy and kinetic energy had the same equilibrium value in two metastable liquids when ρ = 0.898 g/cm 3 and ρ = 0.873 g/cm 3 , while the potential energy and the bond energy all had different equilibrium values.With the increase of the degree of stretching, the lower density metastable liquid had higher potential energy and lower O-H bond energy, the rates of increase, and decrease, were 0.69% and 0.83%, respectively.
With increased stretching in a metastable liquid, the distance across each water molecule gradually increased, while the van der Waals interaction of system was weakened, the O-H bond in the water molecule released energy, thus, the activity of each water molecule increased, and the whole water molecule system became more unstable, so the potential energy of the water molecule system increased.
After bubble nucleation, for a metastable liquid at 0.873 g/cm 3 , the angle energy and kinetic energy of system remained constant, while the potential energy decreased by 1.22% and the bond energy of the system increased by 1.67%.So it could be speculated that when cavitation occurred, with increasing bubble volume, the distance between water molecules in the liquid phase was closer than that when unstretched, so the van der Waals interaction increased: the O-H bond also absorbs energy, and the water molecules in the system were more stable, so the potential energy decreased.
In summary, hydrogen bonds played a major role in cavitation process, the reason why the bubble generation and growth was that the hydrogen bonds were constantly being broken.Furthermore, the O-H bonds of water molecules also played on the XY-plane at t = 1000 ps an indispensable role in the cavitation process, when the liquid water was stretched from the normal state to the metastable state, the energy released in the O-H bond increased the potential energy of the system, making cavitation more likely to occur, but when cavitation occurred, the O-H bonds absorb energy to reduce the system potential energy and promote the growth of bubbles, finally stabilising the bubbles.

Liquid Water Systems Containing Oxygen Molecules
The simulation of a liquid water system containing oxygen molecules consists of three stages: a relaxation-stretchingconstant temperature stage, whose initial structure contained 4912 water molecules and 50 oxygen molecules.Here the initial density was 1.00035 g/cm 3 , the box size was 86.47 Å × 86.47 Å × 20 Å , and the density was 0.978 g/cm 3 after relaxation, while the density of the liquid water (flexible water molecules) system with oxygen molecules was 0.888 g/cm 3 , and 0.876 g/cm 3 respectively after stretching.It should be noted that only the vdW interaction of the oxygen molecules and water molecules are considered, while the electrostatic interaction of the oxygen molecules and water molecules to form hydrogen bonds are not considered in the simulation.
As shown in Fig. 7, when ρ = 0.873 g/cm 3 the system pressure remained at around -168 MPa, and after a waiting time (850 ps), cavitation began to occur and the system pressure increased rapidly to a final pressure of -35 MPa.For the dissolved oxygen model, the system pressure was maintained at around -150 MPa at ρ = 0.888 g/cm 3 and 0.876 g/cm 3 and then the pressure quickly rose to the equilibrium value of -35 MPa after cavitation.The critical negative pressure was -150 MPa.
Cavitation occurred at different times under different tensions: cavitation occurred at 440 ps when ρ was 0.888 g/cm 3 , while it occurred at 70 ps when ρ was 0.876 g/cm 3 .The nucleation rate for bubbles was 8.37 × 10 34 m −3 s −1 at ρ = 0.876 g/cm 3 and 1.34 × 10 35 m −3 s −1 at ρ = 0.888 g/cm 3 .So at the same initial density, the initial relaxation pressure, the temperature and the same tension, liquid water containing 50 oxygen molecules had a higher critical negative pressure and a larger bubble nucleation rate than that in liquid water without oxygen molecules.During bubble nucleation, the oxygen atoms dissolved in the water reduced the critical negative pressure and increased the bubble nucleation rate.
Table 2 lists the energy parameter changes during cavitation when the density of liquid water was 0.873 g/cm 3 (with no oxygen molecules present therein): Table 3 lists the energy parameter changes during cavitation when the density of liquid water was 0.876 g/cm 3 (with 50 oxygen molecules present therein).2 and 3, compared with those systems without oxygen molecules, the system that had 50 oxygen molecules dissolved therein had greater potential energy, a lower bond energy, and a lower hydrogen bond strength during cavitation.The cause of this phenomenon may have been that, during cavitation, on the one hand, oxygen molecules weakened the hydrogen bonds in liquid water; on the other hand, because of the presence of covalent bonds in oxygen molecules, the bonding in the system released more energy after the tensile liquid entered a metastable state, which imparted a higher potential energy to the metastable fluid, thus the system with dissolved oxygen molecules had a lower critical negative pressure and a greater bubble nucleation rate than a system without oxygen molecules.

Discussion
MD simulations were conducted to study cavitation of water: the TIP4P/2005f model, under volume-controlled stretching, was used.The liquid water with, and without, oxygen molecules were subjected to different degrees of stretching before analysing the changes in energy parameters during cavitation.A few key findings are summarised as follows: Compared with the critical negative pressure of -189 MPa calculated by Wang(Wang.P. , 2017) using the TIP4P/2005 rigid water molecular model simulation, the critical negative pressure calculated using the TIP4P/2005f flexible water molecular model in simulated bubble nucleation was -168 MPa.The critical negative pressure calculated for a flexible water molecule containing 50 oxygen molecules was -150 MPa.It was found that the calculated critical negative pressure was closer to the experimental value when using the flexible water molecule model, therefore, it was deemed more appropriate to use flexible water molecules in the study of bubble nucleation.Although the critical negative pressure obtained was much less than the experimental value, it can be seen from the calculation of the critical negative pressure obtained by use of the flexible water molecule containing 50 oxygen molecules that the improved initial configuration of the simulation, such as the addition of dissolved gas molecules, etc. was worthwhile.The liquid water model was thus improved and gave a simulated negative pressure closer to that found experimentally.
Hydrogen bonding plays a major role in the cavitation process: due to the continuous breaking of hydrogen bonds, bubble generation and growth are promoted.However, comparing rigid water molecules with flexible water molecules, it can be found that the OH bonds in the water molecules contributed to the cavitation process.The process also played an indispensable role in liquid water stretching from a normal state to a metastable state.The OH bond energy was released which made cavitation easier by virtue of the associated potential energy released.When cavitation occurred, the O-H bond energy absorption capacity reduced the potential of the system to promote the growth of bubbles, and finally cavitation bubbles were stabilised.
The different numbers of dissolved oxygen molecules in the system exerted some influence on the initial cavitation of flexible water molecules.The effect of oxygen covalent bonding significantly reduced the maximum critical negative pressure and increased the bubble nucleation rate.On the one hand oxygen molecules weakened the hydrogen bonds in the liquid water; on the other hand, due to the presence of covalent bonds in the oxygen molecules, the bonds in the system could release more energy after the stretched liquid reached a metastable state Metastable fluids have higher potential energy, thus, a system containing dissolved oxygen molecules has a lower maximum critical under-pressure and a larger bubble nucleation rate than a molecule not containing oxygen.

Figure 2 .
Figure 2. Pressure changes after stretching

Figure 5 .
Figure 5. Potential energy and kinetic energy of the system under different metastable states

Figure 7 .
Figure 7. Variation of pressure after liquid water stretching • .The type of the OH bond is MORSE and both D r and β are OH bond parameters, where D r = 103.3436kcal/mol, β = 2.287Å −1 , the type of HOH bond angle is harmonic, and k a is the HOH bond angle parameter, where k a = 43.9349Kcal/mol•rad 2 .

Table 2 .
The energy parameters changes in the system during cavitation (0-0.876)