Zhaoyang Yu,Jing Li*,Xianren Zhang,*

1 State Key Laboratory of Organic-Inorganic Composites,Beijing University of Chemical Technology,Beijing 100029,China

2 College of Biomedical Engineering &the Key Laboratory for Medical Functional Nanomaterials,Jining Medical University,Jining 272067,China

Keywords:Cavitation Dissolved gas Thermodynamics Molecular simulation Surface tension

ABSTRACT Cavitation in water generally takes place at much lower negative pressure than predicted from theories.In this work,we try to stress the discrepancy from the influence of the dissolved gas on cavitation nucleation.By combining molecular dynamics simulation and thermodynamic analysis,we evaluated the lowering of surface tension as a function of density of gas molecules in gas clusters formed in aqueous solution.We found that the obtained surface tension of small gas clusters is much more substantially reduced than expected.The surface tension lowering and the non-ideality of gas molecules in the clusters are then taken into account in determining the nucleation of cavitation,and as a consequence,the required negative pressure for cavitation becomes comparable to experimental values.Thus,we give an alternative explanation for the discrepancy of cavitation pressure between experiment and theory,i.e.,it is the substantially reduced surface tension for small gas nuclei,which have not been taken into account in theory,along with the ideal gas approxiamtion that induce its deviation from the experimental values.

1.Introduction

A liquid can sustain a certain negative pressure in a metastable state due to the nucleation barrier for phase transition.Reducing the external pressure below the threshold associated with the tensile strength of the liquid,a phase transition known as cavitation occurs.Within the framework of nucleation theory for cavitation,small vapor-or gas-filled cavities are needed to be created in the liquid as nuclei,which then initiate the rapid phase transition.

Cavitation in water can be observed abundantly both in macroand micro-scale.Yet,the mechanism for cavitation nucleation remains unclear,partially due to the fact that one is still unable to predict the cavitation-required negative pressure with the nucleation theory.A large discrepancy between the cavitation pressure predicted by classical nucleation theory (CNT) and that from experimental measurement was reported repeatedly [1-4].Different experimental techniques gave a consistent cavitation pressure～-30 MPa for water at room temperature[3,4].But this value is far larger than the theoretical prediction based on classical nucleation theory (CNT) and molecular simulation [3,5].CNT predicted the critical negative pressure for cavitation to be around-140 MPa at room temperature,and recent molecular simulation gave a cavitation pressure lower than -120 MPa.The origin of the large discrepancy is still under debate [1-4].

There exists an explanation that the discrepancy is a consequence of the capillary approximation used in the classical nucleation theory (CNT) [6].Much efforts were devoted to the re123vision of the classical nucleation theory.The representative of this kind of effort is the introduction of the Tolman length correction [7].Lubetkin [8] suggested,alternatively,that the discrepancy was due to a previously unidentified factor: the surface activity of the gases which forms the nuclei.Experiments reported that the measured surface tension was in fact a function of dissolved gas pressure[9-18].This observation leads to the idea that gases can act as surfactants that enables the decreasing of surface tension,although the mechanism is not clear.For cavitation,the possible involvement of dissolved gas in nucleation has been mentioned [8].But,how the dissolved gas affects cavitation nucleation has not been quantitatively evaluated.This is because the extent of surface tension lowering is unknown for small gas nuclei formed during nucleation.

For multi-component (here water and dissolved gas) systems,the bubbles radius dependence of surface tension remains unclear.In fact,most commonly used approximation is γ（r）=γ（r→∞）,with γ（r→∞）the surface tension for the planar interface.From molecular simulation side,most studies on the curvaturedependence of surface tension were in fact performed in a singcomponent,non-polar Lennard-Jones (LJ) system.Rao and Berne[19] used Monte Carlo simulations on an argon-like droplet to determine surface tension.Their study showed that the error induced by the approximating γ（r）with γ（r→∞）will be at most on the order of 10%.The small deviation in surface tension cannot interpret the disagreement of cavitation pressure between experiments and theories.

In this work,the influence of the dissolved gas on cavitation nucleation has been studied.The lowering of surface tension as a function of density of gas molecules in nuclei is evaluated by combining molecular dynamics simulation and thermodynamic analysis.Unexpectedly,we found that the surface tension of small gas cluster is substantially reduced.The surface tension lowering and non-ideality of the clustered gas molecules are then taken into account for determining cavitation nucleation.We found that the reduced surface tension can give an interpretation on the discrepancy of cavitation pressure between experiments and theories.

2.Molecular Dynamics Simulations

Here we want to stress the essential role of dissolved gas played in cavitation in water,which was most frequently neglected in theoretical analysis.For this purpose,we first perform atomistic molecular dynamics(MD)simulations to demonstrate the clustering of methane in water.Beside,in order to reflect the particularity of water in cavitation [20-22],here we need to employ the allatom model for water,rather than the extensively used LJ model or other coarse-grained models.The simulation results help us to determine the surface tension as a function of size of small clusters mainly composed of gas methane molecules.Then we established a thermodynamic model to analyze the possible nucleation barrier and nucleation rate for the cavitation of methane-saturated aqueous solution.

First,in our MD simulations various gas clusters of different size were achieved by varying the amount of methane molecules dissolved in the water.The single point model [23] for methane was chosen to represent methane molecules,which can reproduce quite well the phase equilibria of methane hydrate [24].Water molecules were simulated with the SPC/E model [25].For the initial simulation configurations,a given number of gas molecules were randomly placed into the simulation box filled by the water molecules.The number of gas molecules was adjusted to obtain different gas clusters.In this work,the number of methane molecules dissolved varied substantially from 70 to 2500 (including 70,100,150,300,400,500,900,1200,1500 and 2500,respectively).

All simulations were performed in the isothermal-isobaric(NPT) ensemble using the program GROMACS 2018 [26] with a 2×10-6ns time step,and neighbor lists were updated every five steps.The Lennard-Jones parameters for the non-bonded interaction between different species were determined with the conventional geometric mean combination rule.Electrostatic interactions were computed using particle-mesh-Ewald algorithm [27],and a cutoff of 1 nm was used for columbic and Lennard-Jones interactions.The molecular geometries were maintained stable using the LINCS algorithm [28] and the water molecules were kept rigid by the SETTLE algorithm [29].Periodic boundary condition was applied along all three directions.V-rescale thermostat [30] was used for fixing temperature at 300 K with a 10-4ns coupling constant,and isotropic Parrinello-Rahman barostat[31]was used for controlling pressure at 0.1 MPa with a coupling constant 2×10-3ns.After a short simulation run for reaching equilibrium,a series of NPT MD simulations of 80 ns were performed at 300 K and 0.1 MPa to investigate the microscopic properties of gas clusters of different sizes.

Our simulation results show that the systems containing 70,100,150,300,400,500,900,1200,1500 and 2500 methane molecules eventually formed the largest gas clusters containing～53,77,113,255,331,441,781,1106,1405 and 2379 methane molecules,respectively(see Fig.1).The typical process for the formation of large gas clusters can be divided into two steps.Methane gas molecules in water first formed two or more small bubbles,and then they merge to form a large bubble (see Fig.1(a) and 1(b)).

After the largest cluster reaches stability,we then determined the radial density distribution of methane molecules in the gas cluster (Fig.2(a)).Depending on the cluster size,the density of methane molecules inside the cluster varies from 8.6 nm-3to 12.4 nm-3.The densities for supercritical methane gas in the clusters are rather large,being even close to the liquid phase density for methane (the liquid density of CH4ranges from 15.8-17.3 nm-3).Based on the methane density obtained in Fig.2a,the methane density in gas clusters ρ as a function the number of methane molecules containednis obtained and shown in Fig.2(b).The figure shows that with the increase of the size of the cluster,the methane molecules density inside the clusters gradually decreases.

With the density distribution of methane (Fig.2(a)),we can estimate the equilibrium radiusrof the cluster using the hyperbolic tangent fitting ρ=0.5（ρin+ρout）+0.5（ρin-ρout）tanh [2（R-r）/ξ],in which ρinand ρoutare the methane density inside and outside the cluster,respectively,and ξ is the interfacial width.Then the relation between the equilibrium radiusrand the number of methane molecules in the clusternwas determined(see Fig.2(b)).

3.Thermodynamic Analysis

The gas density from MD simulations (Fig.2) enables us to determine the internal pressure of a given cluster,viachoosing an appropriate equation of state.Then,the obtained internal pressure was used to determine the surface tension of the small cluster,with the assumption that the Laplace equation holds.This procedure provides an effective route to surmount the difficulty of molecular simulation in determining the surface tension of highly curved surfaces.

Since the encaged methane molecules in gas clusters are in supercritical state and featured with a high pressure,we employed an equation of state (EOS) [32] that was particularly designed for high pressure supercritical methane gas.Note that if we choose ideal gas law to describe the pressure of gas molecules in the such tiny bubble,as usually done in thermodynamic analysis,a large deviation can be introduced.The EOS developed specially for supercritical gas is in the form of [32].

withkbBoltzmann constant,Tthe temperature,δ the contrast density,δ=ρ/ρc,τ the contrast temperature,τ=Tc/T,andai,di,tibeing the coefficients given in Table 1.Using Eq.(1),along with the obtained gas density ρ from MD simulation,we can obtain internal gas pressure in clusters (see Fig.3(a)).

Table 1 The coefficients used in the EOS [32]

The surface tension γ for a gas cluster was then determined with Laplace equation.

Fig.1.Snapshots of gas cluster formation in the systems containing(a) 300 methane molecules and(b) 1200 methane molecules,respectively.In order to show clearly the gas clusters,the atoms of water molecules are drawn in a much smaller scale.

Fig.2.(a)Radial density distribution for CH4 molecules in various gas clusters that contain a number of methane molecules ranging from 53-2379.(b)The gas density ρ in cluster and the cluster radius r vary with the number of methane molecules contained in gas clusters.

withrthe radius of the gas clusters from MD simulation(Fig.2(b)),Pcthe gas pressure in clsuters(Fig.3(a)),andPextthe external pressure chosen in our MD simulations.Fig.3(b) gives the calculated surface tension according to Eq.(2) as a function of the cluster size,n(the number of methane molecules in the cluster).It is surprising to find that when the number of molecules in the cluster is small,the surface tension can even reduce to 0.35 N·m-1,smaller than half the surface tension for water at ambient pressure,0.72 N·m-1.We must point out that this large decrease of surface tension cannot explained by the Tolman length,which leads to a decreases at most 10%.

Next,with the thermodynamic analysis we discuss how the reduced surface tension affects the nucleation of cavitation.When a given negative pressure is exerted to an initially gas saturated solution,the dissolved gas molecules tend to form gas clusters.For the formation of a gas cluster,the driving force is the chemical potential difference between the metastable state under the exerted negative pressure and the gas saturated state of the aqueous solution before stretching [33].To determine the chemical potential change due to negative pressure,we assumed that for a closed system of methane gas-saturated aqueous solution at an initial pressurePiof 0.1 MPa and a temperature of 300 K,a negative pressurePfis then applied to the closed system.For methane molecules,the change of chemical potential μ due to stretching[33]can be described with the Gibbs-Duhem equation.

in which vm=0.591 nm3is the partial molar volume of a single methane molecule [34].By integrating Eq.(3) to obtain chemical potential change,we have μ=νm（Pf-Pi）.Note that here μ represents the chemical potential change for a dissolved methane molecule due to the exerted negative pressurePf.

Fig.3.(a) Pc and (b)γ as a function of the number of methane molecules in a cluster.

To determine the nucleation barrier for cavitation,we assumed the following path for the formation of a gas cluster: As the first step,a spherical cavity of volumeVcwill generate when the solution is stretched by a negative pressurePf.The creation of the cavity needs additional volume work and surface energy,-（Pc-Pf）Vcand γAcrespectively.HerePcis the internal pressure of the cluster,γ the surface tension andActhe surface area of the cavity.Then,as the second step,the neighboring methane molecules enter the cavity with a change of chemical potential μ.Since the number of neighboring methane molecules entering the cavity are of a rather small fraction when compared to the total number of methane molecules dissolved,the chemical potential for methane molecules in the solution is assumed to keep constant [35].Therefore,the required free energy change for the formation of gas clusters can be wrote as.

With Eq.(4),we can determine how the reduced surface tension affects the nucleation barrier for cavitation,and how it improves the agreement of cavitation pressure between experiments and theories.Since experiments reported cavitation pressure ranges generally from -20 MPa to -30 MPa,here we chose the negative pressure for cavitation in methane-saturated aqueous solution to-24 MPa,namelyPf=-24 MPa.The required free energy change for gas cluster formation,ΔGas given in Eq.(4),shows a local maximum which is denoted as the nucleation barrier ΔG*(Fig.4).The figure shows that the energy barrier of cavitation nucleation ΔG*under the experimental pressure is rather low,mainly because of the significant decrease of surface tension for gas clusters and also because of the appropriate EOS used.

Fig.4.ΔG as a function of the number of methane molecules in a cluster.

Then,the rate of cluster formationJcan be calculated with.

in whichJ0can be approximately calculated [16] according toJ0=ρwxe0.5with ρwthe number density of the solvent and × the molar concentration of the solute.At the given pressure of-24 MPa,the nucleation rate is aboutJ=5.5871×1023number/s,which means that nucleation occurs instantaneously at the given negative pressure of -24 MPa.

Then we explored the negative pressurePfthat is needed for the cavitation nucleation,by considering the influence of both molecule numbernand negative pressurePf.The free energy change,ΔG,for forming a gas cluster of sizenat the given external pressurePfcan be determined with Eq.(4).Fig.5 gives the contour of ΔGas a function ofnandPf,in which the critical nuclei at a given negative pressure was identified when,ΔGcorresponds to a local maximal asnchanges.For the determined critical nuclei as a function of the given negative pressure are represented by the dished line in Fig.5.It demonstrates that whenPf＜-15 MPa,the nucleation barrier decreases sharply with the decrease of pressure,so that in this range the cavitation becomes commonly observed.If we assumeJ＞1 as the condition for observable cavitation nucleation within the experimental reachable time,the upper bound of cavitation pressure can be located at -11.6 MPa.

Fig.5.The contour of free energy change ΔG as a function of the exerted negative pressure Pf and the number of gas molecules in gas clusters n.In this figure,the red dashed line gives the size of critical nuclei n* at given negative pressure.Along the line,the nucleation barrier,ΔG*,increases with n*.The star symbol represents the highest value (upper bound) of negative pressure Pf for experimentally observable cavitation (with J=1).

4.Conclusions

In this work we stress the essential role of dissolved gas played on cavitation nucleation,which was most frequently neglected in theoretical analysis.For this purpose,we first perform all-atom molecular dynamics (MD) simulations to show the clustering of methane molecules in water.By combining molecular dynamics simulation and thermodynamic analysis,we evaluated the lowering of surface tension of gas clusters formed in aqueous solution as a function of gas density in the clusters.It is surprising to find that the surface tension of small gas clusters is more significantly reduced than we expected.Furthermore,the effect of lowering surface tension is enhanced by the increase of density for gas molecules in the nucleus,or by the decrease of the cluster size.

Cavitation in water is generally takes place at much lower negative pressure than predicted from the theory.To interpret the discrepancy,the surface tension lowering found here is then taken into account in determining the nucleation pressure.Further,we also included the non-ideality of encaged gas molecules in clusters by choosing an equation of state that is particularly designed for high-pressure supercritical gases.Our results show that the determined cavitation pressure becomes comparable to experimental values.Thus,we suggest that the non-ideality of gas molecules confined in small gas nuclei and the substantially reduced surface tension,both of which have not been taken into account in theory,induce the deviation of predicted cavitation pressure from the experimental value.

Data availability

Data will be made available on request.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

This research was supported by the National Natural Science Foundation of China (21978007).