Molecular dynamics simulations in the microcanonical ensemble are performed to study the collapse of a bubble in liquid water using the single-site mW and the four-site TIP4P/2005 water models. To study system size effects, simulations for pure water systems are performed using periodically replicated simulation boxes with linear dimensions, L, ranging from 32 to 512 nm with the largest systems containing 8.7 × 106 and 4.5 × 109 molecules for the TIP4P/2005 and mW water models, respectively. The computationally more efficient mW water model allows us to reach converging behavior when the bubble dynamics results are plotted in reduced units, and the limiting behavior can be obtained through linear extrapolation in L−1. Qualitative differences are observed between simulations with the mW and TIP4P/2005 water models, but they can be explained by the models’ differences in predicted viscosity and surface tension. Although bubble collapse occurs on time scales of only hundreds of picoseconds, the system sizes used here are sufficiently large to obtain bubble dynamics consistent with the Rayleigh–Plesset equation when using the models’ thermophysical properties as input. For the conditions explored here, extreme heating of the interfacial water molecules near the time of collapse is observed for the larger mW water systems (but the model underpredicts the viscosity), whereas heating is less pronounced for the TIP4P/2005 water systems because its larger viscosity contribution slows the collapse dynamics. The presence of nitrogen within the bubble only starts to affect bubble dynamics near the very end of the initial collapse, leading to an incomplete collapse and strong rebound for the mW water model. Although nitrogen is non-condensable at 300 K, it becomes highly compressed and reaches a liquid-like density near the collapse point. We find that the dissolution of nitrogen is much slower than the movement of the collapsing water front, and the re-expansion of the dense nitrogen droplet gives rise to bubble rebound. The incompatibility of the collapse and dissolution time scales should be considered for continuum-scale modeling of bubble dynamics. We also confirm that the diffusion coefficient for dissolved nitrogen is insensitive to pressure as the liquid transitions from a compressed to a stretched state.

The prediction of multiphase or cavitating flows is inherently complex due to the large span in spatial and temporal scales of the distinct physical phenomena involved, ranging from bubble nucleation at molecular scale to sheet or cloud cavitation at the macroscopic scale.1–4 Among them, the collapse of bubbles is of particular interest as it is known to cause erosion in hydroturbines through high-pressure pulses, liquid jets, or local high temperatures generated during collapse.5–7 Experiments have led the investigation of bubble collapse. However, the spatiotemporal resolution is restricted by limitations of the imaging instruments,8 and the measurements of local extreme conditions (e.g., temperature) within the bubble are often indirect.9 In addition to experimental measurements, theoretical analyses and computational studies have contributed to our understanding of bubble collapse. The Rayleigh–Plesset (RP) equation10,11 is an ordinary differential equation describing spherical cavity dynamics in an incompressible infinite medium. It is a simple yet powerful tool for deriving many basic results regarding bubble dynamics.1,12 Meanwhile, computational fluid dynamics (CFD) simulations can provide more details and are capable of probing more complex systems. However, various assumptions are introduced to account for the multiphase nature of the system, such as the source term describing the vapor production rate in a local volume, leading to different number and/or forms of the transport equations being solved.13–15 

For the study of bubble dynamics, molecular dynamics (MD) simulations directly apply classical mechanics to generate trajectories of molecular (particle) systems where the particle interactions are computed using computationally efficient molecular-mechanics force fields. The time and length scales approachable by MD simulations are typically orders of magnitude smaller than those probed by continuum-scale simulations, making MD only suitable for probing molecular scale phenomena, such as nucleation.16–19 However, with increasing computing power, recent MD simulations are capable of probing much larger system sizes.20–23 In this way, system size effects can be carefully estimated and MD simulation results can serve as a validation for the underlying physical models utilized in continuum-scale simulations. Previous MD simulations have looked into bubble collapse in both Lennard-Jones (LJ) fluids and molecular fluids with or without a nearby surface or external field.21–29 To expand the applicability and validate principles seen in MD simulations of aqueous bubble collapse, we study systems with a length scale exceeding a half-micron, compare different water models while considering their thermophysical properties, and examine the effects of non-condensable gas at realistic concentrations.

In this work, we report MD simulation results for the collapse of a bubble in water using both the computationally efficient single-site mW model30 and the very accurate four-site TIP4P/2005 water model.31 System size effects were studied for both water models, but an order of magnitude larger box length (or three orders of magnitude more molecules) can be accessed with the mW model. Results obtained from MD simulations are compared with numerical solutions from the RP equation, which further explains the different collapse dynamics predicted by the two water models. To study the effects of nitrogen on bubble collapse, a single-site mW-compatible nitrogen model was developed, whereas the TraPPE nitrogen model32 without adjustments was used in conjunction with the TIP4P/2005 water model. In addition to system properties, spatially localized properties, such as density and kinetic temperature, were analyzed and visualized with heat maps.

1. Single-site models

The mW water model30 utilizes both two-body and three-body (distance–distance–angle) interaction terms such that
U=ij>iϕ2(rij)+ijik>jϕ3(rij,rik,θijk),
(1)
ϕ2(rij)=AεBσrijmσrijnexpσrijaσ,
(2)
ϕ3(rij,rik,θijk)=λε(cosθijkcosθ0)2expγσrijaσ×expγσrikaσ,
(3)
where rij, rik, and θijk are the distance between particles i and j, the distance between particles i and k, and the angle between particles j, i, and k. The mW model uses the following parameters: A = 7.049 556 277, B = 0.6 022 245 584, m = 4, n = 0, γ = 1.2, a = 1.8, θ0 = 109.47°, λ = 23.15, ϵ/kB = 3114.5 K (where kB is the Boltzmann constant), and σ = 2.3925 Å. Both two-body and three-body interactions are truncated at a distance of (4.3065 Å).
A single-site nitrogen model compatible with the mW water model is developed here. Given the weaker angular dependence of dipole–quadrupole interactions compared to water’s hydrogen-bonding interactions, the nitrogen model does not include a three-body term. To reduce the number of adjustable parameters, the N2–N2 and water–N2 interactions are described by a shifted-force LJ 12-6 potential,
Usf(rij)=U(rij)U(Rc)(rijRc)dU(rij)drij|rij=Rc,
(4)
where
U(rij)=4εijσijrij12σijrij6,
(5)
where ɛij and σij are the well depth and diameter for the unshifted LJ potential. For both N2–N2 and water–N2 interactions, the truncation distance, Rc, is set to 9 Å. Following prior work,33 fluid phase equilibrium calculations are used to determine suitable LJ parameters. Here, ϵN2N2 and σN2N2 are fitted to the critical temperature and density of neat N2. The resulting parameters are ϵN2N2/kB=135 K and σN2N2=3.62 Å, and the model also yields satisfactory saturated vapor pressures (see Fig. S1 in the supplementary material). Given the simplicity of the single-site water and nitrogen models, suitable parameters for the unlike interactions cannot be determined from combining rules. Thus, the ϵwaterN2 and σwaterN2 parameters are fitted to the solubility of nitrogen in liquid water at two state points where we pick the lowest pressure at T = 298 K and the highest pressure at 343 K from the data set in the work of Chapoy et al.34 Using data at elevated pressure reduces the statistical uncertainty for the phase equilibrium simulations (due to the higher mole fraction of nitrogen in the aqueous phase) and should be more relevant as the bubble collapse leads to compression of the nitrogen-rich region. The resulting parameters are ϵwaterN2/kB=123 K and σwaterN2=3.22 Å. Details of the fitting procedure for the N2–N2 and water–N2 interaction parameters are provided in the supplementary material.

2. Multisite models

Both the four-site TIP4P/2005 water31 and the three-site TraPPE nitrogen models32 are nonpolarizable and use a rigid molecular geometry. The pairwise additive interactions are described by a combination of LJ and Coulomb potentials,
U(rij)=4εijσijrij12σijrij6+qiqj4πε0rij,
(6)
where qi and ɛ0 are the partial charge on site i and the relative permittivity of vacuum, respectively. The Lorentz–Berthelot combining rules are used to determine the LJ parameters for the unlike interactions. The Coulomb interactions are computed using the particle–particle particle-mesh Ewald summation technique35 with a relative force error of 10−5. Both the LJ interaction and the real-space contribution from the Ewald summation are truncated at a distance of 14 Å. To test the performance of the combination of the TIP4P/2005 and TraPPE nitrogen models, the solubility of N2 in water is computed at the same two state points as for the single-site models. The TIP4P/2005 + TraPPE combination underestimates the experimental solubility by a factor of about 1.3 (see Table S1), which is sufficiently accurate for the present study to not warrant any adjustment of the unlike interaction parameters.

The MD simulations were performed using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) software package.36 However, minor modifications were required for the largest systems containing more than 5 × 108 molecules to change some variables to 8-byte integer, increase the array dimensions, and initialize velocities from the Maxwell–Boltzmann distribution with zero linear momentum for subregions of the system. All simulations were carried out using cubic boxes with periodic boundary conditions. For the setup of an initially two-phase system, the space is divided into a surrounding liquid region and a spherical bubble region with radius Rw. For the ease of analysis and visualization, the origin of the bubble region is placed at the center of the cubic box. Due to the very low saturated vapor pressures predicted by the mW and TIP4P/2005 water models at 298 K (which underpredict the experimental vapor pressure by factors of 60 and 4, respectively),37 the estimated number of water molecules within the largest bubble would be less than 5. Therefore, for the neat water simulations, all molecules are placed inside the liquid region and the bubble is essentially a void. The setup procedure before starting the bubble collapse simulation can be divided into two steps: (1) First, the water molecules in the liquid region are equilibrated using a simulation in the NVT ensemble at T = 298 K where a spherical LJ wall at Rw prevents the molecules from entering the bubble region. (2) Next, the LJ wall is removed and the system undergoes a very short relaxation for 0.1 ps in the NVT ensemble at T = 298 K. This short period is sufficient to remove a temperature spike upon removal of the LJ wall. The setup for the water/nitrogen mixture adopts a similar procedure, except that the initial configuration is obtained by combining a previously equilibrated pure water system containing molecules only outside of the bubble region with N2 molecules placed randomly inside the LJ wall—that is, a situation where the nitrogen concentrations in the bubble and in the surrounding liquid region are not in equilibrium.

To probe the bubble collapse dynamics unperturbed from thermostating interferences, following the 0.1 ps pre-equilibration, the simulations are carried out in the microcanonical (NVE) ensemble and the timer for the bubble collapse starts at the beginning of the NVE simulations (t = 0). For simulations with the single-site models, the time step used for equilibration and the majority of the collapse period is 5 fs, but it is reduced to 2.5 fs as the collapse point is approached. Otherwise, the large velocities of particles near the bubble wall (see Sec. III C) would lead to problems with energy conservation. The switch to the shorter time step was made at about 0.9 of the bubble collapse time, tc, for the systems started with the larger initial compression and at about 0.97 tc for those with the smaller initial compression. For the simulations with multisite models, the time step is set to 1 fs throughout the entire trajectory. Given the very large system sizes used here, the trajectories need to be broken into smaller parts to fit within the wall-clock limits of the computer systems used. To ensure that the trajectories are not perturbed by these restarts, we used binary restart files, kept the grid structure for the Message Passing Interface (MPI) processes (e.g., a 32 × 32 × 32 MPI grid for the mW-512-0.075 system), and retained LAAMPS settings (e.g., fix, neigh_modify, compute, etc.) in all the input files. It was also checked that no warnings were issued in the log file.

For all the systems studied in this work, the overall water density was fixed at values of 998.0 and 995.0 kg/m3 for the mW and TIP4P/2005 models, respectively. These values are slightly above (by 0.03% and 0.06%, respectively) the predicted values for the saturated liquid density of these models at T = 298 K;37 that is, given sufficient time, these systems would relax into a homogeneous liquid phase. Confining the water molecules to the region outside of the bubble leads to an initial state with a higher density in this region and, consequently, higher (system) pressure. For the mW water model, the effect of initial liquid density was explored by performing two series of simulations with Rw/L = 0.15 and 0.075 (where L is the edge length of the cubic simulation box), yielding bubble volume fractions of 0.0141 and 0.001 77, respectively, at setup. The corresponding initial pressure values for the systems at t = 0 ps are about 800 and 110 atm, respectively. For the TIP4P/2005 water model, simulations were performed only for Rw/L = 0.15, and the initial pressure is about 370 atm. For pure water systems, simulations were performed using various L values to study system size effects. For convenience, pure water systems in this work are denoted by a combination of water model, L in units of nm, and Rw/L. Simulations for the water/nitrogen mixture were performed for mW-128-0.15 and TIP4P/2005-64-0.15 systems and are denoted as mW-128-0.15-N2 and TIP4P/2005-64-0.15-N2, respectively. The nitrogen mole fraction xN2 is about 1.2 × 10−5. This concentration is close to the solubility of N2 at T = 298 K and p = 1 atm, and the resulting initial nitrogen pressure inside the bubble region estimated using the ideal gas equation of state is ∼1 atm. As is true for neat water systems, the equilibrium state for the mixture systems is a homogenous liquid phase with all of the nitrogen molecules dissolved—that is, in the condensed state with small partial molar volume. A summary of the systems studied here is provided in Table I.

TABLE I.

Description of simulated systems. Linear box dimension, L, LJ wall radius normalized by box length, Rw/L, number of water molecules, Nwater, number of nitrogen molecules, NN2, initial bubble radius at the start of the NVE simulation, R0, initial wall velocity, dRdt0, bubble collapse time, tc, and most negative wall velocity captured preceding the collapse point, dRdtBCP. The value of dRdtBCP for the TIP4P/2005-32-0.15 system is not reported due to large noise. A well-defined bubble collapse point is not observed for the TIP4P/2005-64-0.15-N2 system.

LR0dRdt0tcdRdtBCP
Compounds (models)nmRw/LNwaterNN2nmm/spsm/s
Water (mW) 32 0.15 1 093 177 4.75 −31 18.75 −990 
 64 0.15 8 745 414 9.56 −33 33.125 −3120 
 64 0.075 8 745 414 4.81 −7 34.875 −790 
 128 0.15 69 963 308 19.16 −25 63.25 −4060 
 128 0.075 69 963 308 9.64 −11 69.625 −2610 
 256 0.15 559 706 465 38.37 −28 123.625 −8010 
 256 0.075 559 706 465 19.26 −7 148.975 −2840 
 512 0.075 4 477 651 718 38.48 −6 309.00 −4660 
Water (mW)/nitrogen (single-site) 128 0.15 69 963 308 832 19.16 −25 63.625 −1520 
Water (TIP4P/2005) 32 0.15 1 089 891 4.84 −25 69.4 ⋯ 
 48 0.15 3 678 381 7.21 −35 98.5 −320 
 48 0.15 3 678 381 7.21 −34 98.5 −400 
 64 0.15 8 719 125 9.74 −3 121.5 −410 
Water (TIP4P/2005)/nitrogen (TraPPE) 64 0.15 8 719 125 104 9.71 −9 ⋯ ⋯ 
LR0dRdt0tcdRdtBCP
Compounds (models)nmRw/LNwaterNN2nmm/spsm/s
Water (mW) 32 0.15 1 093 177 4.75 −31 18.75 −990 
 64 0.15 8 745 414 9.56 −33 33.125 −3120 
 64 0.075 8 745 414 4.81 −7 34.875 −790 
 128 0.15 69 963 308 19.16 −25 63.25 −4060 
 128 0.075 69 963 308 9.64 −11 69.625 −2610 
 256 0.15 559 706 465 38.37 −28 123.625 −8010 
 256 0.075 559 706 465 19.26 −7 148.975 −2840 
 512 0.075 4 477 651 718 38.48 −6 309.00 −4660 
Water (mW)/nitrogen (single-site) 128 0.15 69 963 308 832 19.16 −25 63.625 −1520 
Water (TIP4P/2005) 32 0.15 1 089 891 4.84 −25 69.4 ⋯ 
 48 0.15 3 678 381 7.21 −35 98.5 −320 
 48 0.15 3 678 381 7.21 −34 98.5 −400 
 64 0.15 8 719 125 9.74 −3 121.5 −410 
Water (TIP4P/2005)/nitrogen (TraPPE) 64 0.15 8 719 125 104 9.71 −9 ⋯ ⋯ 

During the large-scale MD simulations, we compute only the virial pressure of the system from the diagonal components of the system pressure tensor because this does not add to the computational expense. That is, the pressure calculation does not account for spatial variations of the virial in these inhomogeneous systems. Previous works investigating droplets38 and bubbles39 with relatively small systems have calculated the normal and transverse components of contours as a function of radial distance using the Irving–Kirkwood40 and Harasima41 pressure tensors. Here, we have not followed this more complex approach because of (a) the additional computational workload, (b) complications arising from the three-body term used for the mW model, and (c) problems of statistical nature due to rapid evolution of the interface during collapse and the fraction of interfacial molecules being negligible for our large systems. The pressure inside of the bubble is estimated here using the instantaneous number densities of (vapor-like) water and nitrogen molecules within the bubble radius, the average kinetic temperature of these particles, and the ideal gas law.

In addition to the MD simulations probing bubble collapse, which are not run until a given system reaches its equilibrium state due to the prohibitive computational expense for these very large system sizes, we also carried out simulations for smaller systems to investigate nitrogen dissolution and evaporation in bubbly systems and nitrogen diffusion in homogeneous systems.

To assess the time scales associated with nitrogen dissolution from an overcompressed bubble and nitrogen evaporation into a bubble from a supersaturated solution, simulations for the single-site models were conducted under conditions where a bubble is thermodynamically stable for the chosen system size, and the N2 molecules were initialized at locations distributed either wholly inside or outside of the preformed bubble. The system was comprised of 133 155 mW water molecules and 133 mW nitrogen molecules (a N2 mole fraction of 0.001) placed into a cubic box of L = 16 nm (yielding a total mass density of 974 kg/m3). The simulations were carried out in the canonical ensemble with the temperature set to 298 K. Based on prior simulations for neat water,37 a bubble radius of Rw ≈ 2.3 nm (Rw/L = 0.14) is estimated to be stable at this condition, and we assumed that the presence of nitrogen leads to a higher pressure inside the bubble and, hence, a larger bubble. The systems were initialized with a preformed bubble of Rw = 2.4 nm devoid of any water molecules and alternately containing either all 133 N2 molecules for the simulations probing dissolution or no N2 molecules for the simulations probing evaporation. The remaining water and N2 molecules were homogeneously dispersed throughout the remainder of the box volume. A short energy minimization was performed to a force tolerance of 10−6 kcal/mol/Å, after which a 25 ps pre-equilibration was conducted with a time step of 0.5 fs and a spherical LJ wall at Rw in place to prevent water and nitrogen molecules from entering or exiting the bubble region. At this point, the wall was removed and the system was allowed to come to chemical equilibrium over the span of 100 ns. A larger time step of 5 fs was used during this period. The temperature was maintained using the Nosé–Hoover thermostat with a damping time of 100 time steps, and the linear center of mass drift was removed every 100 time steps. Five independent replicas were used for both dissolution and evaporation studies.

MD simulations to determine the diffusion coefficient of nitrogen dissolved in water (i.e., a homogeneous single phase) using both the single-site and multisite models were conducted in the canonical ensemble. The systems contained 10 000 water molecules with either 10 or a singular N2 molecule. The cubic box edge length was adjusted to obtain specific densities, using the total mass of water and nitrogen, ranging from 995 to 1005 kg/m3 for the single-site models and from 992 to 1000 kg/m3 for the multisite models. For the chosen system size with L ≈ 6.7 nm, the lowest densities are much higher than the maximum density allowed for a stable bubble, i.e., the systems are in the homogeneously stretched phase.37 The temperature of 298 K was maintained using the Nosé–Hoover thermostat with a damping time of 100 time steps. To verify that the 10 N2 molecules will not aggregate together at the higher concentration, test simulations were initialized with all N2 molecules in a single cluster. Within less than a nanosecond, this aggregate dispersed entirely, emphasizing that N2 motions are uncorrelated. As will be seen later, statistically indistinguishable diffusion coefficients were obtained for the 10:10 000 and 1:10 000 N2:H2O ratios. Subsequent simulations for the systems with 10 N2 molecules were initialized placing them in dispersed locations. For the single-site models, the initial structure was relaxed for 25 ps with a time step of 0.5 fs, which was followed by a 250 ps equilibration using a 5 fs time step. The diffusion coefficient was computed over a 250 ns production run with a 5 fs time step. For the multisite models, the initial structure was relaxed for 1 ps with a time step of 0.1 fs, followed by a 100 ps equilibration using a 1 fs time step. The production run spanned 20 ns with a 1 fs time step. A total of 12 independent simulation replicas were conducted for the purpose of uncertainty estimation. The average simulation pressure was computed via the molecular virial42 and reported as the average over replicas, with the uncertainty estimated from the standard error of the mean (SEM). Diffusion coefficients were computed via the Einstein relation,42,
6tD=|ri(t)ri(0)|2,
(7)
where the diffusion coefficient, D, is obtained from the mean squared displacement (MSD) (the quantity in angular brackets). The value of D was fitted to MSD data as a function of time, ensuring that the MSD was in the diffusive regime. Because the uncertainty in MSD is heteroscedastic (increasing with time), a weighted linear fit to MSD data from each simulation replica was performed with weights of 1/σ2, where σ2 is the variance in the MSD at each t.43 The variance is the square of standard deviation at a time lag t = τ, which is best represented as the SEM of the observations that contributed to the average value at τ (because the weights 1/σ2 pertain to a normal distribution, rather than the 3D χ-squared statistical distribution associated with an MSD in 3D space). The estimates for D from each simulation replica were then used to compute an average, and the uncertainty in D is reported as the SEM obtained from the 12 replicas.

The calculation of bubble volume follows a similar protocol as our previous work on equilibrium properties of bubbly water systems.37 The analysis uses a mesh-based algorithm with a resolution of 0.2 × 0.2 × 0.2 nm3 to divide space into cells. If a cell contains a liquid-like water molecule (i.e., a water molecule that is surrounded by at least two other water molecules using an O–O distance threshold of 0.33 nm) or its center falls within 0.33 nm of at least two other liquid-like water molecules, then the cell is determined to belong to the liquid region. All non-liquid-like cells are considered as vapor-like. The bubble is defined as the largest contiguous cluster of vapor-like cells, where two cells are considered to be part of the same cluster when they share at least one corner. The bubble volume is calculated from the number of cells belonging to the largest cluster. For the water/nitrogen systems, only the water molecules are considered for sorting cells into liquid-like and vapor-like. To analyze the bubble volume efficiently for systems with large number of molecules, a parallelized neighbor search algorithm with domain decomposition is applied.44 For consistency, the bubble collapse point (BCP) is defined as the time when the bubble volume decays below a cutoff of 0.08 nm3 for the neat water systems (i.e., ten contiguous meshes and roughly equivalent to the volume of three water molecules) or the largest void in the system moves to a new location and as the first volume minimum for the water/nitrogen mixtures. The bubble radius, R, is calculated from the bubble volume assuming a spherical shape, which was found to hold reasonably well even for voids with a radius of about 1 nm.37 The interfacial bubble wall velocity dRdt(t) is obtained numerically from R(t) with a center difference scheme for middle points and one-sided difference scheme for boundary points. These analysis methods are used for both the pure water and the water/nitrogen systems. For the mixtures, the fraction of nitrogen molecules contained within the bubble, αN2, is also of interest. As will be shown later, the water/nitrogen interface becomes somewhat diffuse after the bubble collapse, and a unique definition is not possible. Here, a nitrogen molecule is considered to be inside the bubble when the distance between its center of mass and the bubble center is less than R(t).

The time resolutions near the collapse point for the mW and TIP4P/2005 systems are 0.125 and 0.5 ps, respectively. Given the large system sizes used in the current study, it is not possible to store coordinates during the entire trajectory and much of the analysis is performed on the fly. In cases where the coordinates are needed to obtain spatially resolved images near the collapse point, a short part of the trajectory is re-run using the restart file from a few picoseconds before the collapse point, this time writing the coordinates of all molecules to file.

The initial values (i.e., at the start of the NVE part of the trajectory, t = 0) for the bubble radius, R0, and wall velocity, dRdt0 are listed in Table I. Due to the short thermalization period after the LJ wall is removed, R0 is slightly different from Rw, and dRdt0 assumes a small negative value; that is, the bubble radius is shrinking. The simulations for the mW-L-0.075 systems show consistently less negative values for dRdt0 than the mW-L-0.15 simulations started with a more compressed liquid region.

The Rayleigh–Plesset (RP) equation takes the form
Rd2Rdt2+32dRdt2+4νliqRdRdt+2γρliqR+ΔP(t)ρliq=0,
(8)
where νliq and ρliq represent the kinematic viscosity and density of the surrounding liquid; γ represents the surface tension at the bubble interface; and ΔP(t) is defined as P(t) − Pbub(t), which is the difference between the external pressure infinitely far from the bubble and the pressure within the bubble. Without the viscosity and surface tension terms, the RP equation is named the Rayleigh equation. To map MD simulation conditions into input parameters for solving the RP equation, P is taken to be equal to the system pressure, Psys, while Pbub is assumed to be 0. Although the latter is a satisfactory approximation for the neat water systems, it clearly would not hold for the water/nitrogen mixtures where we find that the nitrogen density inside the bubble becomes quite large as the BCP is approached. For simulations in the microcanonical ensemble, both ρliq and Psys are time-dependent, where ρliq is estimated assuming an infinitely sharp interface. The time evolutions of ρliq, Psys, and their ratio (last term of the RP equation) for the mW-256-0.15, mW-512-0.075, and TIP4P/2005-64-0.15 systems are shown in Fig. S3. The trajectories of ρliq and Psys are elongated by repeating the last data point when solving the RP equation. The initial conditions [R0 and dRdt0] are taken from Table I. The numerical solutions of the RP and Rayleigh equations are obtained using the MATLAB ode45 package.45,46

It was previously noted that some of the thermophysical properties predicted with the mW and TIP4P/2005 water models differ significantly from each other and sometimes deviate substantively from experiment.37 Since the bubble collapse simulations encounter a range of compressed liquid densities, the viscosities for the homogeneous liquid phase represented by the mW and TIP4P/2005 water models were calculated at T = 298 K for the densities spanning 995–1016 kg/m3 (see Fig. S4). We find that the change in νliq with respect to ρliq is not significant and, hence, νliq is assumed to be constant. Values of 3 × 10−7 and 8 × 10−7 m2/s are chosen here for the mW and TIP4P/2005 models, respectively.37 Although the water molecules near the bubble interface heat up significantly during the bubble collapse process (see Sec. III D), the system temperature does not change appreciably and assuming a constant value for νliq appears reasonable. To test the sensitivity of the RP solutions to the viscosity term, the value of νliq is varied. Based on prior work,37, γ decreases by only around 15% when the interface is changed from planar to higher curvature for a roughly spherical bubble with a radius of around 1 nm. Thus, a good first-order approximation is to assume a constant value of γ. To test the sensitivity of the RP predictions to γ, the value predicted for the flat interface at T = 298 K (γ = 0.065 N/m for both the mW and the TIP4P/2005 models)19,37 and a lower value of 0.050 N/m (accounting for curvature and local heating19,37,47) are considered here.

The time evolution of bubble volume, V, bubble radius R, and wall velocity, dRdt for pure water systems with the interactions described by the mW water model are presented in Fig. 1. With increasing system size (L, compare data with same line style but different color) and initial compression of the outside liquid region (compare solid and dashed lines), the most negative wall velocity, dRdtBCP, captured preceding the BCP, becomes more negative, indicating a more violent collapse (for numerical data, see Table I). The most violent collapse is observed for the mW-256-0.15 system, and the corresponding dRdtBCP is about −8 km/s; this corresponds to Mach 3.5 (based on the value of 2.3 km/s for the sound speed of liquid water at T = 298 K and P = 1 atm calculated using the mW model and the method of Lustig48). For this system, a clear rebound is observed with the bubble growing back to R > 3 nm (V > 100 nm3). For the mW-128-0.15 and mW-512-0.075 systems, |dRdtBCP| exceeds 4 km/s, and a minor rebound to R ≈ 0.9 nm (V ≈ 3 nm3) is observed, but the rebound is very short-lived with the bubble vanishing after 12 ps in both cases.

FIG. 1.

Time evolution of the bubble volume (top), bubble radius (middle), and wall velocity (bottom) for pure mW water systems. Different colors denote systems with different initial bubble radii (in nm, see legend). Solid and dashed lines represent data for the mW-L-0.15 (higher initial compression) and mW-L-0.075 (lower initial compression) systems, respectively. For better visualization, data for the wall velocity are only shown up to the most negative value preceding the initial collapse point.

FIG. 1.

Time evolution of the bubble volume (top), bubble radius (middle), and wall velocity (bottom) for pure mW water systems. Different colors denote systems with different initial bubble radii (in nm, see legend). Solid and dashed lines represent data for the mW-L-0.15 (higher initial compression) and mW-L-0.075 (lower initial compression) systems, respectively. For better visualization, data for the wall velocity are only shown up to the most negative value preceding the initial collapse point.

Close modal

To assess system size effects, the time evolution of bubble volume and radius is plotted in reduced units in Fig. 2, where bubble volume and radius are normalized by their initial values at t = 0, and time is normalized by the corresponding bubble collapse time, tc. Note that this type of nondimensionalization process is a common practice in fluid mechanics,49 and bubble collapse dynamics solved on larger scales are often presented in this way as well.50,51 For both the mW-L-0.15 and mW-L-0.075 systems, a converging behavior for the initial bubble collapse process is observed as the system size increases. As mentioned above, the rebound after initial collapse is only observed for a few systems but, where rebounding is observed (i.e., the mW-128-0.15 and mW-256-0.15 systems), there is remarkable consistency between system sizes in reduced units. It should be noted that, under these conditions, the RP equation is not capable of predicting a rebound after the initial collapse point.

FIG. 2.

Reduced bubble volume (top) and bubble radius (bottom) as a function of reduced time for mW-L-0.15 (left) and mW-L-0.075 (right) water systems. Different colors denote systems with different initial bubble radii. The cyan curve shows the extrapolation to infinite system size (see text).

FIG. 2.

Reduced bubble volume (top) and bubble radius (bottom) as a function of reduced time for mW-L-0.15 (left) and mW-L-0.075 (right) water systems. Different colors denote systems with different initial bubble radii. The cyan curve shows the extrapolation to infinite system size (see text).

Close modal

To estimate the limiting behavior (L = ) for the bubble collapse dynamics, an extrapolation based on inverse linear dimension was performed. First, the reduced time range (0 ≤ t/tc ≤ 1) was divided into intervals of length 0.01. For each time slice, a linear fit was performed for V/V0 as a function of 1/L, and the y-intercepts yield the data for L = (see Fig. 3). The correlation coefficients are satisfactory with R2 > 0.995 for 0.1 ≤ t/tc ≤ 0.8. Somewhat larger deviations from linearity are observed at the beginning and near the BCP, which may be explained by the slightly different initial conditions in terms of dRdt0 (see Table I) and the extreme heating and supersonic wall speeds observed for only the larger systems. The estimates of reduced bubble volume and radius for L = are also shown in Fig. 2. For both the mW-L-0.15 and mW-L-0.075 systems, the collapse trajectories for the largest bubble are very close to the limiting behavior, indicating that the collapse behavior of a bubble with initial radius around 40 nm is already a good representative for the collapse of mesoscopic bubbles with radii at least on the order of micrometers, which can be observed experimentally with high-speed cameras.

FIG. 3.

Reduced bubble volume at different reduced times, t/tc, as a function of 1/L obtained for mW-L-0.15 (top) and mW-L-0.075 (middle) systems. Circles and dashed lines represent simulation data and their linear fits. The bottom graph shows the correlation coefficient, R2, of the linear fits as a function of reduced time for the systems with different initial compression of the liquid region.

FIG. 3.

Reduced bubble volume at different reduced times, t/tc, as a function of 1/L obtained for mW-L-0.15 (top) and mW-L-0.075 (middle) systems. Circles and dashed lines represent simulation data and their linear fits. The bottom graph shows the correlation coefficient, R2, of the linear fits as a function of reduced time for the systems with different initial compression of the liquid region.

Close modal

The time evolution of the bubble volume, V, bubble radius, R, and wall velocity, dRdt, for pure water systems predicted by the TIP4P/2005 water model are presented in Fig. 4. With the initial liquid pressure value for the TIP4P/2005-L-0.15 systems falling in between those for the mW-L-0.15 and mW-L-0.075 systems (see Fig. S3) and using the same simulation setup, one would expect a similar collapse behavior for the mW and TIP4P/2005 models. However, comparison between Figs. 1 and 4 indicates qualitatively different dynamics. The simulations for the TIP4P/2005 model predict a local minimum in dRdt at the early stage of bubble collapse process. Near the BCP, a local maximum, though less obvious, is also observed in dRdt. These two local extrema correspond to two inflection points in R(t) where the sign of curvature changes, whereas the mW water model predicts a continuously decreasing dRdt (increasing in magnitude) before reaching the BCP. Because of the elongated collapse process, the bubble collapse observed for the TIP4P/2005 water model is overall much less violent.

FIG. 4.

Time evolution of the bubble volume (top), bubble radius (middle), and wall velocity (bottom) for pure TIP4P/2005 water systems. Different colors denote systems with different initial bubble radii (in nm, see legend). For better visualization, data for the wall velocity are only shown up to the most negative value preceding the collapse point.

FIG. 4.

Time evolution of the bubble volume (top), bubble radius (middle), and wall velocity (bottom) for pure TIP4P/2005 water systems. Different colors denote systems with different initial bubble radii (in nm, see legend). For better visualization, data for the wall velocity are only shown up to the most negative value preceding the collapse point.

Close modal

To understand the stark difference between these water models, we turn to the RP equation and recognize model-dependent thermophysical properties.37 By evaluating the signs of different terms on the left-hand side of the RP equation, we observe that it is possible to reach an inflection point in R(t), where d2Rdt2=0, when the negative viscosity term (due to dRdt<0 for a collapse process) balances out the other terms that are positive (due to νliq, γ, and ΔP being positive). As noted earlier, the mW and TIP4P/2005 models yield kinematic viscosity values at ambient conditions of 3 × 10−7 and 8 × 10−7 m2/s, respectively (see also Fig. S4). The corresponding experimental value is 8.9 × 10−7 m2/s.52 Thus, the TIP4P/2005 model underestimates the viscosity by about 10%, whereas the mW model underestimates it by a factor of 3.

Figure 5 presents MD simulation results and numerical solutions of the Rayleigh and RP equation for the largest systems: mW-256-0.15, mW-512-0.075, and TIP4P/2005-64-0.15. Although the Rayleigh equation captures the initial bubble dynamics (i.e., when ΔP is very large) reasonably well for these three systems, it fails to predict the qualitative differences observed for the mW and TIP4P/2005 water models that emerge later in the collapse process. Inclusion of the surface tension term leads to a faster collapse, whereas an opposite effect is observed for the viscosity term. This physically makes sense as reducing the surface area is an additional driving force, while viscosity describes the resistance to flow. Overall, the agreement between MD simulation results (shown in black) and the RP equation utilizing the appropriate kinematic viscosity and surface tension for a given water model (shown in cyan) is very satisfactory, which has also been reported in some previous studies.25,28,53

FIG. 5.

MD simulation results and numerical solutions of the Rayleigh and RP equations for mW-256-0.15 (left), mW-512-0.075 (middle), and TIP4P/2005-64-0.15 (right) systems. The top and bottom rows present data for R(t) and dRdt, respectively. RE and RE + S denote numerical solutions of the Rayleigh equation in its regular form and including the surface tension term, respectively. RP denotes the numerical solutions of the RP equation. The numbers in parenthesis give the values of νliq (in units of 10−7 m2/s) and γ (in units of mN/m).

FIG. 5.

MD simulation results and numerical solutions of the Rayleigh and RP equations for mW-256-0.15 (left), mW-512-0.075 (middle), and TIP4P/2005-64-0.15 (right) systems. The top and bottom rows present data for R(t) and dRdt, respectively. RE and RE + S denote numerical solutions of the Rayleigh equation in its regular form and including the surface tension term, respectively. RP denotes the numerical solutions of the RP equation. The numbers in parenthesis give the values of νliq (in units of 10−7 m2/s) and γ (in units of mN/m).

Close modal

For the mW-256-0.15 system that yields the most violent collapse, the Rayleigh equation with the additional surface tension term yields slightly better agreement with the MD simulations than the Rayleigh or RP equations, but the differences are rather small when the appropriate νliq for the mW model is used. When νliq = 10−6 m2/s (i.e., close to the experimental value) is used, then an inflection point is also observed in the RP prediction for the mW-256-0.15 system, but this qualitatively fails to match the MD trajectory. For the mW-512-0.075 system, the RP equation with the appropriate νliq for the mW model yields the best agreement with the MD trajectory, and the Rayleigh equation with the surface tension term slightly underpredicts the time until the BCP. Again, using νliq close to the experimental value yields an inflection point and extended collapse process that is not in agreement with the MD results. Reducing the value of surface tension to 0.05 N/m while keeping νliq = 3 × 107 m2/s leads to slightly smaller wall velocities and a slight increase in tc.

For the TIP4P/2005-64-0.15 system, only the RP equation with appropriate νliq for the water model predicts collapse dynamics with an inflection point in R(t) that matches the MD trajectory. That is, a larger νliq value is required to yield the inflection point and the stretched collapse dynamics.25 Using a smaller value for γ results in additional stretching of the collapse trajectory. In contrast, the Rayleigh equation and the RP equation with νliq = 2 × 10−7 m2/s yield a much more violent collapse process. However, the RP equation with the appropriate νliq value is not capable of predicting a second inflection point near the BCP that leads to an acceleration of the final segment of collapse. Prior simulations in the canonical ensemble indicate that changes in surface tension19,37 and viscosity37 are relatively small even for (metastable) vapor bubbles with a radius of about 2 nm. We surmise that the acceleration of the final part (R < 2 nm) of the collapse process may be caused by thermodynamic forces as bubbles are not even metastable anymore.18,37

Bubble dynamics obtained for the TIP4P/2005-L-0.15 systems were also examined in reduced units (see Fig. S5) to assess whether the system sizes studied here are sufficient to obtain the limiting behavior (L = ). Linear fits were made for the data at the three system sizes studied for the TIP4P/2005 water model (see Fig. S6). Unlike the findings for the mW systems, where a very satisfactory linear correlation was observed for most of the bubble collapse process, the linear fits for the TIP4P/2005-L-0.15 systems are good (R2 > 0.95) only at intermediate times (0.33 ≤ t/tc ≤ 0.83), and R2 < 0.5 is found for t/tc < 0.27 and for t/tc in the vicinity of 0.92. To test the sensitivity to initial conditions, a second independent run was performed for the TIP4P/2005-48-0.15 system, but both runs yielded very similar data (see Fig. S5). It appears that the two t/tc regions with the worst R2 values are near the two inflection points occurring at somewhat different reduced times for these TIP4P/2005-L-0.15 systems. It is possible that the linear correlation breaks down for bubble collapse processes with inflection points in R(t) or that larger system sizes are needed to reach the regime where a linear correlation holds. With the currently available computing resources, it is not feasible to run simulations for L = 128 and 256 nm for the TIP4P/2005 model, which has more interaction sites per molecule and requires a larger cutoff and Ewald summation together with a shorter time step, while at the same time exhibiting much slower collapse dynamics due to its larger (but closer to experiment) viscosity.

To quantify the temperature and density distributions in time and space, the system is divided into spherical shells containing ∼1000 water molecules. Here, only the translational kinetic temperature is used to allow for a fair comparison of the simulations employing the single-site and multisite models. At the beginning of the collapse process, the temperature is near 300 K for all shells, and the density profile shows strong oscillations reflecting the presence of the wall (see Fig. S7). As the bubble shrinks and the wall velocity increases, the radial temperature profile develops an exponentially decaying behavior with the innermost shell being always found near the highest temperature for any of the shells (see Fig. S7). The change in density profiles is less pronounced (apart, of course, from changes in the location of the bubble wall), but there is a trend toward a slightly increased density (1020 kg/m3) near the bubble wall at t/tc = 0.95. This transient density enhancement may also affect the surface tension.

We focus here on the innermost shell because these molecules are most likely to cause material damage, cavitation erosion, and sonoluminescence. The time evolutions of the average temperature of the innermost shell (T1000) and the threshold temperature for the ten molecules with the highest translational kinetic temperature (T10) for the pure water systems are presented in Fig. 6. In general, the time evolutions of T1000 and T10 are qualitatively similar for the same system as one would expect for a Maxwell–Boltzmann distribution, and T10 ≈ 4 T1000 holds reasonably well. For the mW water systems, a sharp peak is observed in both T1000 and T10 near the collapse point. With increasing system size and initial outside liquid pressure (comparing mW-L′-0.15 and mW-2 L′-0.075 systems starting with the same initial bubble radius), a higher peak temperature in T1000 and T10 is observed. A comparison of Figs. 1 and 6 leads to the conclusion that a more violent collapse [i.e., larger value of dRdtBCP] corresponds to a stronger heating effect for the pure water systems. However, the local heating is more pronounced. As the system size is doubled, we observe that dRdtBCP also approximately doubles (see Table I), whereas T1000 quadruples. For the largest bubble size studied in this work, T1000 peaks at values exceeding 20 000 and 6000 K for the mW-256-0.15 and mW-512-0.075 systems, respectively, which is about 70 and 20 times higher than the initial value of T1000 ≈ 300 K (see Fig. S7). These extreme temperatures well in excess of 1000 K agree with experimental observations.5,54 Our simulation data support the Jarman hypothesis55 that sonoluminescence is thermal in origin, which is based on the empirical premise that sonoluminescent flashes are broad spectrum and would require emitters with temperatures of several thousand Kelvin.56 

FIG. 6.

Time evolution of the average kinetic temperature in the innermost shell containing ∼1000 molecules (top) and of the threshold temperature for the ten hottest molecules (bottom) for the mW-L-0.15 (left), mW-L-0.075 (middle), and TIP4P/2005-L-0.15 (right) systems. Different colors represent systems with different initial bubble radii.

FIG. 6.

Time evolution of the average kinetic temperature in the innermost shell containing ∼1000 molecules (top) and of the threshold temperature for the ten hottest molecules (bottom) for the mW-L-0.15 (left), mW-L-0.075 (middle), and TIP4P/2005-L-0.15 (right) systems. Different colors represent systems with different initial bubble radii.

Close modal

Compared to the mW systems, a qualitatively different heating behavior is observed for the TIP4P/2005 systems studied here. First, the system size does not seem to have an effect on the evolution in reduced time for both T1000 and T10. Second, both T1000 and T10 increase slightly at the beginning of the collapse process, but then their time evolution flattens out at t/tc ≈ 0.2. At the transition point, the corresponding absolute times are around 14, 20, and 24 ps for R0 = 4.8, 7.2, and 9.6 nm, respectively. These times correspond approximately to the initial inflection point in R(t) (see Fig. 4). In other words, after passing the first inflection point when viscosity effects begin to dominate and dRdt starts to increase, the local heating of the innermost shell is no longer observed.

Figure 7 presents the space–time evolution of the specific density and kinetic temperature in the vicinity of the BCP for the mW-128-0.15 and TIP4P/2005-64-0.15 water systems. Within the 2 ps preceding the collapse point, the bubble radius decreases by around 4 nm for the mW-128-0.15 system, while a decrease of less than 1 nm is observed for the TIP4P/2005-64-0.15 system, again reflecting a much more violent collapse for the mW water system with the lower viscosity. For both water models, the moving interfacial region has an intermediate density and a width of around 0.5 nm. As the collapse point is approached, a slight decrease in the interfacial width is observed. For the mW-128-0.15 system, an extremely high density exceeding 1.6 g/cm3 is found at the BCP for shells within 0.6 nm of the center, and a wave of higher density shells moving away from the bubble center is observed for the first 0.5 ps after the BCP. The mW-128-0.15 system exhibits a small rebound (see Fig. 1), but the interface of the rebounding bubble is rather diffuse. The space–time evolution for the kinetic energy shows that the heating is most pronounced at the interface with a well-defined hot spot with temperatures as high as 9000 K at the BCP for the mW-128-0.15 system. Interestingly, the temperature does not show a feature similar to the radiating density wave immediately following the collapse. Cross-sectional heat maps near the BCP are shown for this system in Fig. S8. Before the BCP, the highest temperatures and voxel velocities (all pointing inward) are observed for the interfacial region. After the BCP, however, the region of highest temperature remains at the center, whereas the aftermath of the collapse leads to high density and voxel velocity (all pointing outward) waves radiating away from the center. The collapse process for the TIP4P/2005-64-0.15 system is much less violent, and the space–time evolution shows a more diffuse picture without a large spike in density nor temperature at the BCP (see Fig. 7).

FIG. 7.

Space–time evolution of the specific density in units of g/cm3 (top) and the kinetic temperature in units of Kelvin (bottom) for the mW-128-0.15 (left) and TIP4P/2005-64-0.15 (right) systems. Properties are computed for shells with a thickness of 0.2 nm in the radial dimension.

FIG. 7.

Space–time evolution of the specific density in units of g/cm3 (top) and the kinetic temperature in units of Kelvin (bottom) for the mW-128-0.15 (left) and TIP4P/2005-64-0.15 (right) systems. Properties are computed for shells with a thickness of 0.2 nm in the radial dimension.

Close modal

To investigate the effects of a non-condensable gas on bubble collapse dynamics, we also performed simulations for water/nitrogen mixtures using single-site and multisite models. It is important to emphasize that these simulations were initialized with all of the inert gas molecules placed inside the bubble but that the number of nitrogen molecules is close to the solubility limit and the initial pressure in the bubble is only about 1.2 bars. Figure 8 presents the bubble dynamics results for the water/nitrogen systems and a comparison to the corresponding pure water systems. In the presence of nitrogen, complete bubble collapse is not observed for the mW-128-0.15-N2 and TIP4P/2005-64-0.15-N2 systems, and the minimum bubble radii are found to be 1.5 and 1.2 nm, respectively. Interestingly, the effect of nitrogen on the bubble dynamics is negligible for the overwhelming majority of the trajectory. For the mW-128-0.15-N2 system, the time of 63.6 ps until the first minimum in bubble radius is only 1% longer than the tc for the mW-128-0.15 system. Even for the TIP4P/2005-64-0.15 systems, the bubble radii and wall velocities closely trace each other for about 110 ps (i.e., about 0.9tc for the system without nitrogen). The close agreement between neat water systems and water/nitrogen mixtures for the majority of the collapse trajectory also extends to the density in the liquid region, the system pressure, and their ratio (see Fig. S9) that constitutes the last term on the right-hand side of the RP equation. We also note that the bubbles maintain spherical symmetry throughout the initial collapse process (see Fig. S9).

FIG. 8.

Time evolution of the bubble radius (top) and the wall velocity (bottom) for the mW-128-0.15-N2 and TIP4P/2005-64-0.15-N2 systems and the corresponding pure water systems.

FIG. 8.

Time evolution of the bubble radius (top) and the wall velocity (bottom) for the mW-128-0.15-N2 and TIP4P/2005-64-0.15-N2 systems and the corresponding pure water systems.

Close modal

However, the effect of nitrogen on the final part of the initial collapse process and the subsequent dynamics is very pronounced. For the mW-128-0.15-N2 system, the terminal wall velocity is nearly three times smaller in magnitude than for the corresponding pure water system (see Table I). For the former system, the initial collapse is followed by a much stronger and longer lasting rebound (see Fig. 8). Eventually, the rebounding bubble shrinks again, followed by another rebound, but the second and third collapses are much softer. The R(t) values for the second and third minima are 2.8 and 2.6 nm, respectively. Due to computer time limitations, the trajectory for this large system (containing 7 × 107 water molecules) was stopped after 200 ps because it became obvious that reaching an equilibrium state with no bubble and all nitrogen molecules dissolved and well distributed in the liquid would be unattainable. For the TIP4P/2005-64-0.15-N2, the presence of nitrogen leads to a dramatic dampening of the initial collapse process; that is, there is not a clear minimum in R(t) and also not a second inflection point in dRdt. We observe a trend toward a decreasing bubble size for the last 50 ps of the trajectory for this system but, again, approach toward the equilibrium state with no bubble appears to be extremely slow compared to the initial collapse process (see inset of Fig. 8). Prior large-scale simulations presented in the work of Shekhar et al.21 also investigated the effects of a non-condensable gas, but the initial gas (number) density in the bubble was more than three orders of magnitude larger than in the present work, i.e., the equilibrium state presented in the work of Shekhar et al. study would be a two-phase system.

To provide insight into the pronounced effects on the final phase of the bubble collapse caused by the presence of nitrogen, we calculated a few properties for only the nitrogen molecules contained within the bubble (see Fig. 9). For the mW-128-0.15-N2 system, the fraction of nitrogen molecules within the bubble, αN2, decays only slowly from near unity (831 out of 832 nitrogen molecules contained within the bubble at t0) to a value of 0.9 during the first 50 ps of the trajectory despite that the bubble volume shrinks by a factor of 6. This is followed by a steep decline to a value of 0.43 when the first minimum is reached; that is, there is some degree of mixing during the more violent part of the collapse process. However, the bubble rebound is accompanied by demixing, and αN2 reaches 0.84 before the bubble starts to shrink again. Subsequent bubble size oscillations are also accompanied by mixing as the bubble volume decreases and demixing as the bubble volume increases. Thus, the mixing appears to be mostly limited to an interfacial zone (see below).

FIG. 9.

Time evolution of the fraction of nitrogen molecules contained in the bubble (top), the density of nitrogen in the bubble (top middle), the kinetic temperature of the nitrogen molecules (bottom middle), and the pressure due to nitrogen in the bubble (bottom) for the mW-128-0.15-N2 and TIP4P/2005-64-0.15-N2 systems. The pressure is estimated using the ideal gas law.

FIG. 9.

Time evolution of the fraction of nitrogen molecules contained in the bubble (top), the density of nitrogen in the bubble (top middle), the kinetic temperature of the nitrogen molecules (bottom middle), and the pressure due to nitrogen in the bubble (bottom) for the mW-128-0.15-N2 and TIP4P/2005-64-0.15-N2 systems. The pressure is estimated using the ideal gas law.

Close modal

The density and temperature for nitrogen molecules in the bubble also undergo oscillations driven by initial collapse and subsequent rebounds. For the mW-128-0.15-N2 system, the specific density of nitrogen inside the bubble, ρN2, reaches in excess of 1000 kg/m3, thereby matching the specific density of the surrounding liquid; however, accounting for the difference in molecular mass, the number density of nitrogen is smaller than that of water. During the first rebound, ρN2 falls to about 20 kg/m3, but quite surprisingly it only falls to 200 kg/m3 during the second rebound. The reason is that the first rebound reaches a bubble volume that is about ten times larger than that for the second rebound, whereas αN2 decays by only 10%. The temperature of nitrogen molecules within the bubble, TN2, reaches maxima of 6400, 570, and 450 K during the first, second, and third bubble collapses, respectively, and falls to 300 K during the rebounds. Previous hybrid RP-MD studies also predicted extreme heating exceeding 10 000 K near the collapse point for violent bubble collapse with inert gas content.57,58 Using the density and temperature values, we can also estimate the pressure in the bubble. For convenience, Fig. 9 presents data using the ideal gas equation of state. [We have performed constant-NVT MD simulations for pure nitrogen systems at selected (ρN2,TN2) values observed during the initial collapse for the mW-128-0.15-N2 system and find that the ideal gas equation of state holds well for all but the most extreme state points, but it underestimates the pressure by up to a factor of 4 for TN2>3000 K. Since nitrogen’s Boyle temperature is 325 K,59 the ideal gas equation of state holds rather well for TN2<500 K.] The simulations indicate that, at the BCP of the mW-128-0.15-N2 system, nitrogen is compressed to a density comparable to that observed for a pressure in excess of 1 GPa and pressure equivalents of 370 and 480 bars for the second and third collapses, respectively.

As discussed before, the collapse process for the TIP4P/2005-64-0.15-N2 system is much less violent without a well-defined BCP. This also leads to qualitatively different trends for the properties of nitrogen within the bubble (see Fig. 9). The fraction of nitrogen molecules contained within the bubble decreases at a near constant rate for the first 130 ps of the trajectory and, when the bubble volume becomes steady, αN2 fluctuates around 0.7. Since the majority of nitrogen molecules are still contained in the bubble, but the volume has decreased by a factor of about 600, a very high nitrogen density (ρN2500 kg/m3) and pressure (pN2450 bars) are also reached for the TIP4P/2005-64-0.15-N2 system. Despite the less violent collapse dynamics, it is still much faster than the time required for the dissolution of nitrogen molecules. As also seen for the corresponding pure water system, the temperature of nitrogen molecules only increases until the first inflection point in R(t), and TN2 does not exceed 400 K.

Figure 10 presents the space–time evolution of water density, nitrogen density, and temperature near the initial collapse point for the mixture systems. The V-shaped interface observed in the water density heat map for the mW-128-0.15-N2 system indicates a clear rebound, consistent with the bubble dynamics results presented in Fig. 8. Preceding the collapse, the interfacial width remains at around 0.5 nm, while it grows rapidly during the rebound. Cross-sectional heat maps near the BCP for the mW-128-0.15-N2 system are shown in Fig. S10. As expected from the soft collapse for the TIP4P/2005-64-0.15-N2 system, the radius of the bubble and the interfacial width barely change within the ±2 ps near the first minimum in R at t = 134 ps. For both the single-site and the multisite models, an incomplete mixing is observed near the collapse point, where water molecules do not penetrate deeply into the interior of the bubble and where nitrogen molecules are highly compressed within the bubble. Interestingly, for the mW-128-0.15-N2 system, the radius of the region of high nitrogen density remains nearly constant after the collapse, while the water front is receding. Instead the nitrogen density in this dense core is decreasing. These observations reflect a difference in the time scales between the slow dissolution of nitrogen into the surrounding liquid water and the fast collapsing water front (see Sec. III E). With respect to temperature, the presence of nitrogen leads to slightly more extreme temperatures for the mW-128-0.15-N2 system than the mW-128-0.15 system, while no significant difference is observed for the TIP4P/2005 systems.

FIG. 10.

Space–time evolution of the specific density in units of g/cm3 for water (top) and nitrogen (middle) and the kinetic temperature in units of Kelvin (bottom) for the mW-128-0.15-N2 (left) and TIP4P/2005-64-0.15-N2 (right, using tc = 134 ps) systems. Properties are computed on shells of width 0.2 nm in the radial dimension.

FIG. 10.

Space–time evolution of the specific density in units of g/cm3 for water (top) and nitrogen (middle) and the kinetic temperature in units of Kelvin (bottom) for the mW-128-0.15-N2 (left) and TIP4P/2005-64-0.15-N2 (right, using tc = 134 ps) systems. Properties are computed on shells of width 0.2 nm in the radial dimension.

Close modal

As discussed in Sec. III D, the simulations indicate that the bubble collapse process for the water/nitrogen mixture systems does not lead to destruction of the bubble and dissolution of the nitrogen molecules; instead we observe a long-lived bubble containing a very highly compressed nitrogen gas. The simulations for the bubble collapse used a very large system size (see Table I) that does not allow us to follow the trajectories for more than a few hundred ps. Thus, we conducted simulations for the single-site models using a 83 smaller system (in regard to number of molecules and volume) to assess the time scales that are pertinent to reaching an equilibrium state from initial states representing an overcompressed bubble or an oversaturated solution region. For the system studied here, the initial state with all N2 molecules placed inside the bubble corresponds to a pressure of about 100 bars inside the bubble and the initial state with all N2 molecules dissolved in the liquid region corresponds to a N2 mole fraction of 10−3.

The trajectories for the systems characterized by mostly dissolution or evaporation are shown in Figs. 11 and S11, respectively. During the first 200 ps of the dissolution process, the bubble rapidly grows to Rw = 2.44 nm and loses about 20% of the N2 molecules. That is, the gas bubble with the highly compressed N2 gas prefers a slightly larger volume than a neat water system with the same number of water molecules. This phase of the process reflects the transition from a perfectly spherical bubble with a sharp interface between liquid and gas found for the initial configuration (and enforced by a retaining wall at R = 2.4 nm) to a slightly less spherically symmetric bubble with an interfacial width similar to the molecular size (about 0.3 nm) and all of the N2 molecules that are nominally outside of the bubble still being close to the interface. The radial density profile calculated over the first 1 ns shows some N2 enrichment at the bubble side of the interface (see Fig. 11).

FIG. 11.

Bubble radius (top) and fraction of N2 molecules inside the bubble (bottom) as a function of time for the dissolution process started from an oversaturated bubble. Five independent simulation replicas are shown in orange, red, green, blue, and yellow to illustrate the fluctuations observed in each trajectory. The average over these five trajectories is shown in black. An exponential function of the form A exp(−t/τ) + C is fitted to the average radius and N2 fraction with characteristic time constant τ = 1.0 × 101 ns and coefficients A = 0.11 nm and C = 2.33 nm for the radius and A = 0.46 and C = 0.36 for the N2 fraction. The bottom graph shows the nitrogen radial density profiles obtained over 1 ns time slices starting at the time indicated in the legend. The data for r ≥ 3 nm are scaled up by a factor of 100.

FIG. 11.

Bubble radius (top) and fraction of N2 molecules inside the bubble (bottom) as a function of time for the dissolution process started from an oversaturated bubble. Five independent simulation replicas are shown in orange, red, green, blue, and yellow to illustrate the fluctuations observed in each trajectory. The average over these five trajectories is shown in black. An exponential function of the form A exp(−t/τ) + C is fitted to the average radius and N2 fraction with characteristic time constant τ = 1.0 × 101 ns and coefficients A = 0.11 nm and C = 2.33 nm for the radius and A = 0.46 and C = 0.36 for the N2 fraction. The bottom graph shows the nitrogen radial density profiles obtained over 1 ns time slices starting at the time indicated in the legend. The data for r ≥ 3 nm are scaled up by a factor of 100.

Close modal

Similarly for the evaporation process simulations, the bubble radius rapidly shrinks over the first 200 ps to Rw = 2.22 nm while less than two N2 molecules nominally enter the bubble (see Fig. S11). That is, the oversaturated liquid phase requires a larger volume (i.e., smaller Rw) than the neat water system with the same number of water molecules. The radial density profile calculated over the first 1 ns does not exhibit N2 enrichment at the bubble the interface (see Fig. S11).

After the initial rapid changes, the bubble is seen to shrink or grow more slowly over time as N2 molecules dissolve into or evaporate from the liquid region, respectively. After about 40 ns, the dissolution and evaporation trajectories reach the same equilibrium state with Rw = 2.33 nm and αN2=0.36. This equilibrium state corresponds to a pressure inside the bubble (calculated from the N2 density using the ideal gas law) of about 40 bars and a solution mole fraction of about 6 × 10−4 that are in agreement with the bulk partitioning experiments34 (see also Fig. S2).

The radial density profiles indicate that the slight enhancement of N2 number density at the bubble side of the interface is also present when equilibrium is reached. For the dissolution process, the N2 density in the liquid region shows a decrease with increasing r at 5 and 10 ns, but the liquid region profiles flatten by 20 ns. For the evaporation process, the N2 density profiles in the liquid region are reasonably flat throughout the trajectory. Based on the density profiles in the liquid region not showing a strong gradient well before the bubble reaches equilibrium, we surmise that diffusion is not rate-limiting for the dissolution/evaporation processes.

Figure S12 shows the evolution of system pressure, the pressure inside the bubble, and the surface tension. During the dissolution process of the compressed gas inside the bubble, the system pressure and the pressure in the bubble decrease, and the reverse holds for the evaporation process from the supersaturated solution. Remarkably, the surface tension does not show pronounced changes during the dissolution and evaporation processes, indicating that the surface tension is not very sensitive to the presence of nitrogen over the range of nitrogen densities within the bubble observed during these trajectories, and the limiting value reached at t = 80 ns falls within 1% of the value obtained from equilibrium simulations for vapor bubbles in neat water.37 

As should be expected, the individual trajectories do not follow exactly the same path, but their dissolution and evaporation processes are quite similar. The rapid oscillations that are seen for both Rw and αN2 are largely due to thermal fluctuations, the resolution of the mesh-based bubble detection algorithm, the conversion of mesh volume to Rw assuming a perfect spherical shape, and partitioning of the system into liquid and bubble regions without allowing for an interfacial region. A simple exponential fit of the form A exp(−t/τ) + C to Rw(t) describes the evolution of the average values of Rw and αN2 remarkably well and reveals a characteristic time constant τ of 1.0 × 101 ns for both evaporation and dissolution (see Figs. 11 and S11).

That is, for the systems investigated here, the characteristic time scale for dissolution is two orders of magnitude slower than the time scale for the initial bubble collapse. This is consistent with the observation that there is insufficient time for N2 to equilibrate with the surrounding liquid during the bubble collapse, giving rise to a rebounding effect for the more violent collapse observed for the mW-128-0.15-N2 system, a softening for the collapse of the TIP4P/2005-64-0.15-N2 system, and a long-lived bubble dissolving only on much longer time scales for both systems. The nonequilibrium rates of bubble growth/shrinkage due to dissolution and evaporation from/into a bubble are found here to be comparable, suggesting that traversing the liquid/vapor bubble interface and mass transport away from/toward the interface in the solution phase are rate-limiting. Of course, the extent of separation of the collapse and dissolution time scales may depend on the boundary conditions governing the collapse, but we surmise that it will play a pivotal role during the more violent phase of gas bubble collapse.

As mentioned above, gas diffusion in the solution phase contributes to the separation of time scales between bubble collapse and gas dissolution. Experimental measurements of the diffusion coefficient of dissolved gases in water are challenging because of the low solubility of common gases. In fact, a popular approach to determine gas diffusion coefficients is by observing the collapse dynamics of large bubbles (with diameters larger than 100 μm) usually at ambient pressure.60 The Taylor dispersion method uses a laminar flow tube, and detection of the dispersion curve of the sparingly soluble gases is aided by performing the measurements at elevated pressure.61 Neither of these methods allows for measurements at pressures approaching the saturated vapor pressure of water. This has led to some uncertainty regarding the pressure dependence of the diffusion of dissolved gases, and some researchers studying bubble growth have used the Chapman–Enskog relation62 to extrapolate diffusion coefficients into the low-pressure regime,63 while others have assumed constant diffusion coefficients independent of pressure.64 For MD simulations, there are no conceptual problems with studying the tracer diffusion in a homogeneous solution phase, but the small number of gas molecules in the simulation box leads to a statistical challenge.

Data for N2 diffusion in water at 298 K covering a range of system densities are presented in Fig. 12. The N2 diffusion coefficients obtained for the multisite models exhibit no discernible density dependence while the pressure ranges from −16 to 160 bars, i.e., from the homogeneously stretched to compressed liquid regime. This finding is in agreement with the experimental measurements using the Taylor dispersion method presented in the work of Cadogan et al.61 that yield D = (1.99 ± 0.05) × 10−5 cm2/s for pressures from 120 to 450 bars. Averaging the N2 data from the simulations for the multisite models yields D = (1.97 ± 0.08) × 10−5 cm2/s. Here, it should be noted that the TIP4P/2005 water model also yields an excellent diffusion coefficient at ambient conditions.31,37 Prior experimental measurements of the N2 diffusion coefficient give values of D = (3.1 ± 0.3) × 10−5 cm2/s (interpolated from data at 293 and 303 K) using the bubble collapse method,60, D = (2.30 ± 0.12) × 10−5 cm2/s (interpolated from data at 293 and 303 K) using a stagnant liquid layer method,65 and D = (2.01 ± 0.20) × 10−5 cm2/s using the dispersion technique.66 

FIG. 12.

Diffusion coefficient of N2 in aqueous solution (left axes) and system pressure (right axes) as a function of system density for the single-site (top) and multisite (bottom) models. Diffusion data are depicted as black up and red down triangles for the 10:10 000 and 1:10 000 N2:water molecular ratios, respectively, while the corresponding pressure data are represented as squares and diamonds, respectively.

FIG. 12.

Diffusion coefficient of N2 in aqueous solution (left axes) and system pressure (right axes) as a function of system density for the single-site (top) and multisite (bottom) models. Diffusion data are depicted as black up and red down triangles for the 10:10 000 and 1:10 000 N2:water molecular ratios, respectively, while the corresponding pressure data are represented as squares and diamonds, respectively.

Close modal

The simulation data for the single-site models yield somewhat lower statistical uncertainties and indicate that the N2 diffusion coefficient may decrease by slightly less than 2% as the pressure increases from −155 to 365 bars, but compared to the data presented in the work of Cadogan et al.61 and multisite model data, the N2 diffusion coefficient of D = (3.59 ± 0.02) × 10−5 cm2/s at |p| ≤ 5 bars is overestimated substantially by a factor of 1.8 (as one may expect from the over prediction of the self-diffusion coefficient in the neat water phase by a factor of 2.9 for the mW water model30,37). Based on the experimental data presented in the work of Cadogan et al.61 and the present simulation study, it is clear that the Chapman–Enskog relation should not be applied to describe the pressure dependence of the diffusion of a dissolved gas. This is because the mean free path length for gas in the solution phase is dictated by the solvent density, which barely changes as a function of pressure for a liquid solvent far from its critical point. The closely scattered data points near zero pressure found here for the single- and multisite models demonstrate that the diffusion coefficient of the dissolved gas does not diverge as the liquid transitions from a compressed to a stretched state. Thus, the observation that the experimentally measured rate of bubble growth upon reduction in pressure can be satisfactorily described under the assumption that the diffusion coefficient is inversely correlated with pressure63 may instead indicate that the growth is not entirely diffusion-controlled.

As also observed in our prior study for neat water,37 the simulations for the mixture indicate a statistically significant change of 1% or 2% for the water diffusion coefficients described by the multi- or single-site models, respectively, over the range of densities investigated here (see Fig. S13). The decrease in water diffusion rate with decreasing pressure can be attributed to the more ice-like behavior exhibited at lower densities.37 

Molecular dynamics simulations are used to probe the collapse dynamics for pure water and water/nitrogen mixtures. For the computationally efficient mW water models, simulations were performed for a range of system sizes (L ranging from 32 to 512 nm) and two values of the initial bubble volume fraction (α = 0.0141 and 0.001 77), yielding initial values for the system pressure of 800 and 100 bars, respectively. We find that, under these conditions, an initial bubble radius of 40 nm is already a good representative for the collapse of micrometer-sized bubbles. Due to the underprediction of viscosity by the mW model, the collapse process is rather violent with a strong rebound for the larger systems and temperatures in excess of 5000 K for the 1000 water molecules belonging to the innermost shell (located at the liquid/vapor interface) of the collapsing bubble. The simulations for the TIP4P/2005 water model, which accurately predicts the water viscosity, yield much elongated collapse dynamics where inflection points corresponding to local extrema in dRdt(t) and only modest heating until the initial inflection point are observed. To understand the difference, MD simulation results are compared with the numerical solutions of the RP equation, and the overall agreement is very satisfactory when all terms in the RP equations are based on the properties of the model systems. It is found that inclusion of the viscosity term in the RP equation is necessary for observing the initial inflection point in R(t), at which the viscosity effect starts to dominate.

The most important finding of this study pertains to the striking differences observed for water/nitrogen mixtures compared to the neat water systems. Even under conditions where all gas molecules would be dissolved in the homogeneous solution phase once equilibrium is reached, we find a clear separation of the time scales for bubble collapse and gas dissolution. Here, the initial collapse is incomplete and results in a gas droplet that is compressed to liquid-like densities. The re-expansion of the dense droplet gives rise to a strong rebound followed by additional incomplete collapse/soft rebound cycles. Based on simulations for smaller systems, we surmise that the time scale for gas dissolution is two orders of magnitude slower than the time scale for the collapse process. We recommend that gas partitioning [i.e., the process of converting a non-condensable gas into a dissolved (condensed) species] and an explicit rate for dissolution/evaporation processes be included in computational fluid dynamics simulations investigating two-phase flow. As an aside, we show that the diffusion coefficient of dissolved gas molecules is insensitive to pressure as the solution phase transitions from a compressed liquid (i.e., large positive pressure) to a homogeneously stretched liquid (i.e., large negative pressure) state.

The supplementary material for this article includes the following: (1) a PDF file containing a description of the development of the single-site nitrogen model, the nitrogen solubility values in aqueous solution using the multisite models, plots of the input data for the Rayleigh–Plesset equation, plots probing the system size effects for the TIP4P/2005 systems, plot of the radial temperature profile at different t/tc values, cross-sectional heat maps near the collapse point for the mW-128-0.15 and mW-128-0.15-N2 systems, plots for the nitrogen evaporation dynamics from oversaturated solution, and the diffusion of water in water/nitrogen mixtures and (2) two MP4 files that show cross-sectional heatmaps for the water and (if present) nitrogen density for systems (a) mW-128-0.15 and mW-128-0.15-N2 and (b) mW-256-0.15.

This work was supported by the United States Office of Naval Research under Grant No. ONR N00014-17-1-2676 with Dr. Ki-Han Kim and Dr. Yin-Lu Young as the program managers. C.K. was supported by the Office of Science, U.S. Department of Energy, under Contract No. DE-AC02-06CH11357. Part of the computer resources was provided by the Minnesota Supercomputing Institute and the Argonne Leadership Computing Facility. We thank Ryan DeFever and Edward Maginn for cross-validation of the shifted-force LJ potential and Joe Katz for alerting us to the controversy regarding the dissolved gas diffusion rate at low pressure.

The submitted manuscript has been created by UChicago Argonne, LLC, Operator of Argonne National Laboratory (“Argonne”). Argonne, a U.S. Department of Energy Office of Science laboratory, is operated under Contract No. DE-AC02-06CH11357. The U.S. Government retains for itself, and others acting on its behalf, a paid-up nonexclusive, irrevocable worldwide license in said article to reproduce, prepare derivative works, distribute copies to the public, and perform publicly and display publicly, by or on behalf of the Government. The Department of Energy will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan. http://energy.gov/downloads/doe-public-access-plan.

The authors have no conflicts to disclose.

Jingyi L. Chen: Conceptualization (equal); Data curation (lead); Formal analysis (lead); Investigation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (supporting). Jesse L. Prelesnik: Conceptualization (supporting); Data curation (supporting); Formal analysis (supporting); Visualization (supporting); Writing – original draft (supporting); Writing – review & editing (supporting). Buyun Liang: Investigation (supporting). Yangzesheng Sun: Investigation (supporting); Methodology (supporting). Mrugank Bhatt: Conceptualization (supporting); Formal analysis (supporting); Writing – review & editing (supporting). Christopher Knight: Software (supporting); Writing – review & editing (supporting). Krishnan Mahesh: Conceptualization (equal); Funding acquisition (equal); Writing – review & editing (supporting). J. Ilja Siepmann: Conceptualization (lead); Funding acquisition (equal); Methodology (equal); Supervision (lead); Visualization (supporting); Writing – review & editing (lead).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

1.
M. S.
Plesset
and
A.
Prosperetti
, “
Bubble dynamics and cavitation
,”
Annu. Rev. Fluid Mech.
9
,
145
185
(
1977
).
2.
G.
Wang
,
I.
Senocak
,
W.
Shyy
,
T.
Ikohagi
, and
S.
Cao
, “
Dynamics of attached turbulent cavitating flows
,”
Prog. Aerosp. Sci.
37
,
551
581
(
2001
).
3.
A.
Singhal
,
M. M.
Athavale
,
H. Y.
Li
, and
Y.
Jiang
, “
Mathematical basis and validation of the full cavitation model
,”
J. Fluids Eng.
124
,
617
624
(
2002
).
4.
C. F.
Delale
,
K.
Okita
, and
Y.
Matsumoto
, “
Steady-state cavitating nozzle flows with nucleation
,”
J. Fluids Eng.
127
,
770
777
(
2005
).
5.
E. B.
Flint
and
K. S.
Suslick
, “
The temperature of cavitation
,”
Science
253
,
1397
1399
(
1991
).
6.
P.
Kumar
and
R. P.
Saini
, “
Study of cavitation in hydro turbines—A review
,”
Renewable Sustainable Energy Rev.
14
,
374
383
(
2010
).
7.
B. K.
Sreedhar
,
S. K.
Albert
, and
A. B.
Pandit
, “
Cavitation damage: Theory and measurements—A review
,”
Wear
372–373
,
177
196
(
2017
).
8.
W.
Lauterborn
and
T.
Kurz
, “
The bubble challenge for high-speed photography
,” in
The Micro-World Observed by Ultra High-Speed Cameras
(
Springer
,
2018
), pp.
19
47
.
9.
M. P.
Brenner
,
S.
Hilgenfeldt
, and
D.
Lohse
, “
Single-bubble sonoluminescence
,”
Rev. Mod. Phys.
74
,
425
(
2002
).
10.
Lord
Rayleigh
, “
On the pressure developed in a liquid during the collapse of a spherical cavity
,”
London, Edinburgh Dublin Philos. Mag. J. Sci.
34
,
94
98
(
1917
).
11.
M. S.
Plesset
, “
The dynamics of cavitation bubbles
,”
J. Appl. Mech.
16
,
277
282
(
1949
).
12.
J.
Franc
,
Fluid Dynamics of Cavitation and Cavitating Turbopumps
(
Springer
,
2007
), pp.
1
41
.
13.
A.
Gnanaskandan
and
K.
Mahesh
, “
A numerical method to simulate turbulent cavitating flows
,”
Int. J. Multiphase Flow
70
,
22
34
(
2015
).
14.
R. K.
Shukla
,
C.
Pantano
, and
J. B.
Freund
, “
An interface capturing method for the simulation of multi-phase compressible flows
,”
J. Comput. Phys.
229
,
7411
7439
(
2010
).
15.
M.
Pelanti
and
K.
Shyue
, “
A mixture-energy-consistent six-equation two-phase numerical model for fluids with interfaces, cavitation and evaporation waves
,”
J. Comput. Phys.
259
,
331
357
(
2014
).
16.
J. L. F.
Abascal
,
M. A.
Gonzalez
,
J. L.
Aragones
, and
C.
Valeriani
, “
Homogeneous bubble nucleation in water at negative pressure: A Voronoi polyhedra analysis
,”
J. Chem. Phys.
138
,
084508
(
2013
).
17.
G.
Menzl
,
M. A.
Gonzalez
,
P.
Geiger
,
F.
Caupin
,
J. L. F.
Abascal
,
C.
Valeriani
, and
C.
Dellago
, “
Molecular mechanism for cavitation in water under tension
,”
Proc. Natl. Acad. Sci. U. S. A.
113
,
13582
13587
(
2016
).
18.
S. H.
Min
and
M. L.
Berkowitz
, “
Bubbles in water under stretch-induced cavitation
,”
J. Chem. Phys.
150
,
054501
(
2019
).
19.
I.
Sanchez-Burgos
,
M. C.
Muniz
,
J. R.
Espinosa
, and
A. Z.
Panagiotopoulos
, “
A deep potential model for liquid–vapor equilibrium and cavitation rates of water
,”
J. Chem. Phys.
158
,
184504
(
2023
).
20.
J.
Diemand
,
R.
Angélil
,
K. K.
Tanaka
, and
H.
Tanaka
, “
Large scale molecular dynamics simulations of homogeneous nucleation
,”
J. Chem. Phys.
139
,
074309
(
2013
).
21.
A.
Shekhar
,
K.
Nomura
,
R. K.
Kalia
,
A.
Nakano
, and
P.
Vashishta
, “
Nanobubble collapse on a silica surface in water: Billion-atom reactive molecular dynamics simulations
,”
Phys. Rev. Lett.
111
,
184503
(
2013
).
22.
S.
Zhan
,
H.
Duan
,
L.
Pan
,
J.
Tu
,
D.
Jia
,
T.
Yang
, and
J.
Li
, “
Molecular dynamics simulation of shock-induced microscopic bubble collapse
,”
Phys. Chem. Chem. Phys.
23
,
8446
8455
(
2021
).
23.
R.
Cabriolu
,
B. G.
Pollet
, and
P.
Ballone
, “
Effect of organic ions on the formation and collapse of nanometric bubbles in ionic liquid/water solutions: A molecular dynamics study
,”
J. Phys. Chem. B
127
,
1628
1644
(
2023
).
24.
C.
Xiao
,
D. M.
Heyes
, and
J. G.
Powles
, “
The collapsing bubble in a liquid by molecular dynamics simulations
,”
Mol. Phys.
100
,
3451
3468
(
2002
).
25.
R.
Hołyst
,
M.
Litniewski
, and
P.
Garstecki
, “
Large-scale molecular dynamics verification of the Rayleigh-Plesset approximation for collapse of nanobubbles
,”
Phys. Rev. E
82
,
066309
(
2010
).
26.
F.
Lugli
,
S.
Höfinger
, and
F.
Zerbetto
, “
The collapse of nanobubbles in water
,”
J. Am. Chem. Soc.
127
,
8020
8021
(
2005
).
27.
H.
Fu
,
J.
Comer
,
W.
Cai
, and
C.
Chipot
, “
Sonoporation at small and large length scales: Effect of cavitation bubble collapse on membranes
,”
J. Phys. Chem. Lett.
6
,
413
418
(
2015
).
28.
V. H.
Man
,
M. S.
Li
,
P.
Derreumaux
, and
P. H.
Nguyen
, “
Rayleigh-Plesset equation of the bubble stable cavitation in water: A nonequilibrium all-atom molecular dynamics simulation study
,”
J. Chem. Phys.
148
,
094505
(
2018
).
29.
S. H.
Min
,
S.
Wijesinghe
,
E. Y.
Lau
, and
M. L.
Berkowitz
, “
Damage to polystyrene polymer film by shock wave induced bubble collapse
,”
J. Phys. Chem. B
124
,
7494
7499
(
2020
).
30.
V.
Molinero
and
E. B.
Moore
, “
Water modeled as an intermediate element between carbon and silicon
,”
J. Phys. Chem. B
113
,
4008
4016
(
2009
).
31.
J. L. F.
Abascal
and
C.
Vega
, “
A general purpose model for the condensed phases of water: TIP4P/2005
,”
J. Chem. Phys.
123
,
234505
(
2005
).
32.
J. J.
Potoff
and
J. I.
Siepmann
, “
Vapor–liquid equilibria of mixtures containing alkanes, carbon dioxide, and nitrogen
,”
AIChE J.
47
,
1676
1682
(
2001
).
33.
J. I.
Siepmann
,
S.
Karaborni
, and
B.
Smit
, “
Simulating the critical behaviour of complex fluids
,”
Nature
365
,
330
332
(
1993
).
34.
A.
Chapoy
,
A. H.
Mohammadi
,
B.
Tohidi
, and
D.
Richon
, “
Gas solubility measurement and modeling for the nitrogen + water system from 274.18 K to 363.02 K
,”
J. Chem. Eng. Data
49
,
1110
1115
(
2004
).
35.
R. W.
Hockney
and
J. W.
Eastwood
,
Computer Simulation Using Particles
(
CRC Press
,
1988
).
36.
S.
Plimpton
, “
Fast parallel algorithms for short-range molecular dynamics
,”
J. Comput. Phys.
117
,
1
19
(
1995
).
37.
J. L.
Chen
,
B.
Xue
,
K.
Mahesh
, and
J. I.
Siepmann
, “
Molecular simulations probing the thermophysical properties of homogeneously stretched and bubbly water systems
,”
J. Chem. Eng. Data
64
,
3755
3771
(
2019
).
38.
S. M.
Thompson
,
K. E.
Gubbins
,
J. P. R. B.
Walton
,
R. A. R.
Chantry
, and
J. S.
Rowlinson
, “
A molecular dynamics study of liquid drops
,”
J. Chem. Phys.
81
,
530
542
(
1984
).
39.
H.
Rezaei Nejad
,
M.
Ghassemi
,
S. M.
Mirnouri Langroudi
, and
A.
Shahabi
, “
A molecular dynamics study of nano-bubble surface tension
,”
Mol. Simul.
37
,
23
30
(
2011
).
40.
J. H.
Irving
and
J. G.
Kirkwood
, “
The statistical mechanical theory of transport processes. IV. The equations of hydrodynamics
,”
J. Chem. Phys.
18
,
817
829
(
1950
).
41.
A.
Harasima
, “
Molecular theory of surface tension
,” in , edited by P. Debye Prigogine (Wiley,
1957
), Chap. 7.
42.
M. P.
Allen
and
D. J.
Tildesley
,
Computer Simulation of Liquids
(
Oxford University Press
,
2017
).
43.
X.
Michalet
, “
Mean square displacement analysis of single-particle trajectories with localization error: Brownian motion in an isotropic medium
,”
Phys. Rev. E
82
,
041914
(
2010
).
44.
D.
Brown
,
J. H. R.
Clarke
,
M.
Okuda
, and
T.
Yamazaki
, “
A domain decomposition parallelization strategy for molecular dynamics simulations on distributed memory machines
,”
Comput. Phys. Commun.
74
,
67
80
(
1993
).
45.
J. R.
Dormand
and
P. J.
Prince
, “
A family of embedded Runge-Kutta formulae
,”
J. Comput. Appl. Math.
6
,
19
26
(
1980
).
46.
L. F.
Shampine
and
M. W.
Reichelt
, “
The MATLAB ODE suite
,”
SIAM J. Sci. Comput.
18
,
1
22
(
1997
).
47.
F.
Magaletti
,
M.
Gallo
, and
C. M.
Casciola
, “
Water cavitation from ambient to high temperatures
,”
Sci. Rep.
11
,
20801
(
2021
).
48.
R.
Lustig
, “
Direct molecular NVT simulation of the isobaric heat capacity, speed of sound and Joule–Thomson coefficient
,”
Mol. Simul.
37
,
457
465
(
2011
).
49.
B. R.
Munson
,
T. H.
Okiishi
,
W. W.
Huebsch
, and
A. P.
Rothmayer
,
Fluid Mechanics
(
Wiley
,
Singapore
,
2013
).
50.
D.
Obreschkow
,
M.
Bruderer
, and
M.
Farhat
, “
Analytical approximations for the collapse of an empty spherical bubble
,”
Phys. Rev. E
85
,
066303
(
2012
).
51.
E.
Johnsen
and
T.
Colonius
, “
Numerical simulations of non-spherical bubble collapse
,”
J. Fluid Mech.
629
,
231
262
(
2009
).
52.
W. M.
Haynes
,
CRC Handbook of Chemistry and Physics
(
CRC Press
,
2014
).
53.
H.
Okumura
and
N.
Ito
, “
Nonequilibrium molecular dynamics simulations of a bubble
,”
Phys. Rev. E
67
,
045301
(
2003
).
54.
H.
Xu
and
K. S.
Suslick
, “
Molecular emission and temperature measurements from single-bubble sonoluminescence
,”
Phys. Rev. Lett.
104
,
244301
(
2010
).
55.
P.
Jarman
, “
Sonoluminescence: A discussion
,”
J. Acoust. Soc. Am.
32
,
1459
1462
(
1960
).
56.
W. B.
McNamara
III
,
Y. T.
Didenko
, and
K. S.
Suslick
, “
Sonoluminescence temperatures during multi-bubble cavitation
,”
Nature
401
,
772
775
(
1999
).
57.
D.
Schanz
,
B.
Metten
,
T.
Kurz
, and
W.
Lauterborn
, “
Molecular dynamics simulations of cavitation bubble collapse and sonoluminescence
,”
New J. Phys.
14
,
113019
(
2012
).
58.
A.
Bass
,
S. J.
Ruuth
,
C.
Camara
,
B.
Merriman
, and
S.
Putterman
, “
Molecular dynamics of extreme mass segregation in a rapidly collapsing bubble
,”
Phys. Rev. Lett.
101
,
234301
(
2008
).
59.
G.
Coccia
,
G.
Di Nicola
,
S.
Tomassetti
,
M.
Pierantozzi
, and
G.
Passerini
, “
Determination of the Boyle temperature of pure gases using artificial neural networks
,”
Fluid Phase Equilib.
493
,
36
42
(
2019
).
60.
D. L.
Wise
and
G.
Houghton
, “
The diffusion coefficients of ten slightly soluble gases in water at 10–60°C
,”
Chem. Eng. Sci.
21
,
999
1010
(
1966
).
61.
S. P.
Cadogan
,
G. C.
Maitland
, and
J. P. M.
Trusler
, “
Diffusion coefficients of CO2 and N2 in water at temperatures between 298.15 K and 423.15 K at pressures up to 45 MPa
,”
J. Chem. Eng. Data
59
,
519
525
(
2014
).
62.
E. N.
Fuller
,
P. D.
Schettler
, and
J. C.
Giddings
, “
New method for prediction of binary gas-phase diffusion coefficients
,”
Ind. Eng. Chem.
58
,
18
27
(
1966
).
63.
J.
Li
,
H.
Chen
,
W.
Zhou
,
B.
Wu
,
S. D.
Stoyanov
, and
E. G.
Pelan
, “
Growth of bubbles on a solid surface in response to a pressure reduction
,”
Langmuir
30
,
4223
4228
(
2014
).
64.
I.
Akhatov
,
O.
Lindau
,
A.
Topolnikov
,
R.
Mettin
,
N.
Vakhitova
, and
W.
Lauterborn
, “
Collapse and rebound of a laser-induced cavitation bubble
,”
Phys. Fluids
13
,
2805
2819
(
2001
).
65.
P. T. H. M.
Verhallen
,
L. J. P.
Oomen
,
A. J. J. M. v. d.
Elsen
,
J.
Kruger
, and
J. M. H.
Fortuin
, “
The diffusion coefficients of helium, hydrogen, oxygen and nitrogen in water determined from the permeability of a stagnant liquid layer in the quasi-s
,”
Chem. Eng. Sci.
39
,
1535
1541
(
1984
).
66.
R. T.
Ferrell
and
D. M.
Himmelblau
, “
Diffusion coefficients of nitrogen and oxygen in water
,”
J. Chem. Eng. Data
12
,
111
115
(
1967
).

Supplementary Material