High-speed video microscopy and numerical modeling of bubble dynamics near a surface of urinary stone

Mechanical action of engineered microbubbles in response to acoustic excitation is increasingly used in medical applications, such as for targeted destruction of biomineralizations.^1–3 To promote preferential action on a target, the bubble shell can be engineered to facilitate accumulation of microbubbles on the surfaces of targets.^1,2 To target and treat urinary stones, Ramaswamy et al.¹ proposed placing engineered microbubbles in the urinary tract or kidney, where the bubbles can accumulate on exposed surfaces of stones (Fig. 1). The accumulated microbubbles then can be energized by one or more energy sources to mechanically erode, pit, and fragment the stone.^1,4,5

With a pulsed laser energy source, a recent in vitro study has shown that adding stone-surface-accumulating (SSA) microbubbles increased the rate of erosion, pitting, and fragmentation of model stones by ∼ 70%.⁴ Other studies have shown that SSA microbubbles can be placed in the urinary tract cystoscopically and energized by an extracorporeal acoustic energy source.⁵ The use of ultrasound with wide beam widths (several centimeters) reduces the burden of precise image guidance (Fig. 1) and, with a cystoscopic delivery of SSA microbubbles, can be performed in diverse clinical settings including the urologist office.

In shock wave lithotripsy, urinary stones are broken with focused shock pulses with pressure amplitudes of 15–150 MPa.⁶ Recent in vitro experiments have shown that urinary stones can be broken by bursts of focused ultrasound with center frequencies in a sub-MHz range and amplitudes of several MPa.⁷ Preliminary experiments with SSA microbubbles suggest that model and urinary stones can be broken by acoustic bursts with pressure amplitudes at the upper levels of diagnostic ultrasound (1.2 ± 0.2 MPa, 0.3–1 MHz).⁵ In these experiments, the stones were implanted in a porcine model and treated by wide beam-width ultrasound with no evidence of urothelium damage and renal parenchymal hemorrhage on histological and gross anatomical examination of post-procedure ureters and kidneys.⁵ 

These encouraging results suggest a need to better understand the dynamics and the mechanical action of SSA microbubbles at these driving conditions. In this work, we used bursts of ultrasound with central frequencies of ∼0.42 MHz and amplitudes of 1.4 ± 0.4 MPa. The dynamics of microbubbles was studied using high-speed video microscopy at a frame rate of up to 10 × 10⁶ fps, capturing complete expansion-collapse cycles of the bubbles. The observed geometry of bubbles at their maximum expansion was used as the initial shape for the numerical modeling of the collapse. For modeling, we solved the multicomponent Euler equations using a two-dimensional (2 D)-axisymmetric code with adaptive time step and mesh refinement algorithms for fine resolution of the gas-liquid interface and shock fronts.^8,9 The use of a 2 D code instead of a one-dimensional (1D) spherical model^10–16 was motivated by the present and previous observations showing that oscillations of microbubbles can be highly nonspherical.^2,17–25

Hsiao and Chahine²⁶ modeled a shelled microbubble subjected to one cycle of 2.5-MHz ultrasound at 1 MPa. The modeling showed nonspherical deformations due to the presence of a nearby rigid boundary and, depending on the standoff distance, one of the following dynamics of the collapse: (1) a single reentrant jet, (2) a ring-type reentrant jet, and (3) a pinching of the bubble.²⁶ Here, we report experimental observations of a circumferential pinching and the numerical modeling showing that the pinching propels two microjets directed away and toward the rigid boundary. When hitting the boundary, the jet was either a single reentrant jet or a ring-type jet, depending on the amount of gas modeled in the bubble.

The present paper is organized as follows. Section II describes the experimental arrangement to record the dynamics of microbubbles at the surface of urinary stones. Section III presents experimental results showing that microbubbles, driven by ∼1.5 MPa, expand up to about 60 μm in diameter, and acquire the shape resembling an oblate spheroid. Section IV describes the numerical modeling of the collapse showing the circumferential pinching and the formation of the two microjets (Sec. IV E). The pressure generated by the collapsing microbubble at the rigid wall (Sec. IV F) exceeded the driving pressure by two-to-three orders of magnitude. Section V shows micro-computed tomography of urinary stones before and after the action of SSA microbubbles illustrating surface erosion and formation of microcracks. The significance and limitations of this work are discussed in Sec. VI. The  Appendix describes numerical details, including an illustration of the adaptive mesh refinement (AMR) algorithm.

SSA microbubbles (Applaud Medical, Inc.) were made of perfluoroalkane gas (C₄F₁₀) encapsulated by lipid shells engineered to accumulate on stone surfaces.¹ The lipid shells incorporate polyethylene glycol structures and synthetic pyrophosphate analogs, collectively functioning to minimize interaction with the urothelium and facilitate accumulation on stone surfaces. The size distribution of the microbubbles was measured with a Coulter counter (Multisizer 4e Coulter Cell Analyzer, Beckman Coulter, Indianapolis, IN) using an aperture size of 30 μm in diameter. The mean diameter of the microbubbles was [mean ± standard deviation (SD)] 1.19 ± 0.04 μm.

These experiments were conducted with surgically retrieved urinary stones. The chemical composition of stones was determined by Fourier-transform infrared spectroscopy and was mostly calcium oxalate monohydrate. Specifically, the surface of the stone shown in the inset of Fig. 2(b) contained 90% calcium oxalate monohydrate, 5% calcium oxalate dihydrate, and 5% apatite; the interior contained 70% calcium oxalate monohydrate, 10% calcium oxalate dihydrate, and 20% apatite.

To position the stone in the focus of a microscope, the stone was glued with a small amount of a Loctite Super Glue to a tip of 0.25-mm thin coverslip (clear vinyl plastic, 18 × 18 × 0.25 mm, VWR International, PA) as shown in the inset of Fig. 2(b). The stone was then kept in water for several weeks for hydration. The hydrated stone was positioned in the water tank (Fig. 2) and aligned at the focus of the microscope using an XYZ-micrometer stage (Thorlabs Inc., NJ).

The test tank was three-dimensionally (3 D) printed from a thermoplastic material (acrylonitrile butadiene styrene) and covered with a waterproof coating (Marine Grade Epoxy 109 Medium, Tap Plastics, CA). For imaging with the inverted microscope, the bottom of the tank had a glass port (75 × 25 × 1 mm microscope slide, VistaVision, VWR International, LLC, Radnor, PA) glued along its edges to the bottom of the tank (Fig. 2). The tank was filled with six liters of water (PURELAB Chorus 1 for Life Science Applications, ELGA, Veolia Water Solutions and Technologies, UK) with an electrical resistivity of 18.2 MOhm-cm and the ultrafiltration to particle size less than 0.05 μm. The water remained in the tank for several days and was in equilibrium with atmospheric gases. During the experiments, the temperature of the water slowly increased from ∼23 °C to ∼28 °C due to the heating by intense light used for the high-speed imaging.

We used both continuous and flashlight illumination. The continuous lighting was provided with a fluorescence illumination system EXFO X-cite 120 (XE120, Photonic Solutions Inc., Mississauga, Ontario, CA). This light source uses a 120-W Metal Halide lamp coupled to a liquid lightguide. The end of this lightguide was positioned at about 1 cm above the stone and provided back light to the stone (Fig. 2).

The side lighting was provided with a flashlamp WRF300 (Hadland Imaging, LLC, Santa Cruz, CA). This spark-discharge lamp produced a light pulse with the duration of about 10 μs. The spark light was delivered through a separate liquid lightguide illuminating the side of the stone proximal to the incoming acoustic waves (Fig. 2).

Images were captured with a high-speed camera HPV-X2 (Shimadzu, Kyoto, Japan) operated in one of two modes. In the full-pixel (FP) mode, the camera captures 128 frames of 400- by 250-pixels at a rate up to 5 Mfps. In the half-pixel (HP) mode, the camera captures every other pixel in a checkerboard pattern, interpolating the images into 400- by 250-pixel frames and recording 256 frames at rates up to 10 Mfps. Exposures were 110 ns at 5 Mfps and 50 ns at 10 Mfps.

Magnification was achieved using a microscope (Nikon Eclipse TS100) with a 4× objective (4×/0.13 PhL DL, WD 16.4, Nikon Plan Fluor), a 2.5× projection lens (Nikon CF PL2.5×), and a 34-cm extension tube (Fig. 2). The numerical aperture of the objective was 0.13 giving the diffraction-limited depth of field of 43 ± 12 μm and lateral resolution of 2 ± 0.5 μm. This diffraction-limited lateral resolution of the objective was considered in choosing the magnification of the projection lens and the length of the extension tube. These optical elements were chosen to have a camera resolution of 1 μm/pixel.

In post-processing, recorded images were digitally enlarged (3× by resampling with preserving details in Adobe Photoshop), rotated, and cropped. The acoustic radiation force was directed approximately from right to left in these images (Fig. 4, Mm. 1, and Mm. 2). The buoyant force was directed into the image plane (Fig. 2).

Driving acoustic waves were generated with a custom-made piezoelectric transducer (Sonic Concepts, Inc., Bothell, WA) positioned at ∼9 cm from the proximal surface of the urinary stone (Fig. 2). The active element of the transducer was made of a piezoelectric plate (72.3 × 30.3 × 3.18 mm) divided into eight elements and connected in pairs. In these experiments, each pair was powered by one of four power amplifiers (AP-400B, ENI, USA). The frequency and duration of the acoustic bursts were computer controlled with a specially designed signal generator, allowing us to trigger the HS-camera with TTL pulses while reproducing the frequency modulation used in the clinic⁵ and also to study other driving regimes. In this work, the transducer was driven with a frequency set of 400, 400, 433, and 433 kHz, generating acoustic bursts with a center frequency of 416.5 kHz and ∼30 μs duration of the envelope (Fig. 4).

The spatial characteristics of the acoustic beam (Fig. 3) were measured using a needle hydrophone (HNR-0500, Onda Inc., Sunnyvale, CA) with a diameter of the sensitive element of 0.5 mm. The acoustic field was scanned with a 2-mm step using an XYZ-positioning system assembled with three motorized linear slides (X-LSM150A, Zaber Technologies, Vancouver, BC, Canada). The motion of the positioning system was controlled by a computer through the RS-232 link with a program written in Python, which was also used to acquire hydrophone data recorded by a digital oscilloscope (HDO4024 Teledyne, LeCroy Corp., NY). These measurements were conducted in a 35 × 27 × 20 cm tank with acrylic walls acoustically isolated with sheets of absorptive material 1-cm thick (Aptflex F28, Precision Acoustics, UK) providing echo reduction greater than 25 dB. The tank was filled with 14 liters of water degassed with a pinhole degasser to 1.1 ± 0.5 mg/l (measured with a dissolved oxygen meter DO 6+, Eutech Instruments, Singapore). To prevent the damage of the Onda hydrophone, the spatial characteristics of the acoustic beam were measured at low amplitudes and are shown normalized on the maximum pressure P_max found among the scans in the three orthogonal planes (Fig. 3). At the position of the stone (Z = 90 mm), the acoustic beam had the cross section with the −6-dB width of ∼30 mm and the height of ∼60 mm (solid contour lines, Fig. 3).

Selected measurements were conducted with another needle hydrophone (Y-104, Sonic Concepts, Inc., Bothell, WA). In particular, this hydrophone was used to measure pressure traces at the driving amplitudes (Fig. 4) in the focus of the microscope at the position of the stone (Fig. 2). For these measurements, the stone was replaced with the sensitive tip of the Y-104 hydrophone—the 1.5-mm diameter ceramic crystal enclosed in a 3-mm metal tube. As the diameter of the tip was comparable with the acoustic wavelength (λ ∼ 3.6 mm), the hydrophone's sensitivity depended on its orientation. The angular response of the hydrophone was measured and taken into account, increasing the uncertainty of pressure measurements to ∼30%.

The pressure traces recorded in the HS-camera test tank at the position of the stone were combined with the HS-camera images into multimedia frames using programs written in LabVIEW (National Instruments, Austin, TX) and converted into movies with the H.264-video format using QuickTime Pro 7 (Fig. 4, Mm. 1, and Mm. 2). In these movies, representative hydrophone traces (dark blue) were superimposed on several waveforms (light blue) showing shot-to-shot variability. The vertical arrow shows the timing of the HS-camera with the uncertainty of 0.2 μs. Time t = 0 was positioned at the start of implosion of the largest solitary bubble.

Micro-computed tomography (micro-CT) of urinary stones (Fig. 9) was used to assess the stone damage produced by the SSA microbubbles. The micro-CT images were acquired with a high-resolution X-ray tomography (MicroXCT-200, Xradia, Inc., Pleasanton, CA) before and after the action of SSA microbubbles.

Bubble dynamics at the surface of a urinary stone is shown in movies Mm. 1 and Mm. 2. The movie Mm. 1 was recorded with the HS-camera in FP mode at 5 Mfps and provides full resolution images although at the slower frame rate than Mm. 2, which was recorded in HP mode capturing every other pixel at 10 Mfps.

During the first two acoustic cycles in these movies (−3.4 < t < 0.8 μs) the microbubbles expanded to a larger maximum size from the first to the second cycle, concomitant with increasing driving pressure. The second collapse (∼0.8 μs) produced a cluster of bubbles. Bubbles in clusters varied in shape, size, and standoff distance, merging with other bubbles. Here, we focus on the collapse of a solitary bubble starting from its maximum expansion (t ≈ 0 – 0.8 μs, Mm. 1 and Mm. 2).

At maximum expansion [Fig. 4(a) and t = 0, Mm. 1–Mm. 2], the bubbles were nonspherical with a cross-section resembling an ellipse truncated by the stone [Fig. 5(a)]. Average major and minor axes of 16 bubbles were measured to be $2Ry=51±9 μ$m and $2Rx=39±6 μ$m [Fig. 5(b)]. The bubbles had an eccentricity $ϵ=1−(Rx/Ry)2$ of 0.63 ± 0.08 and centers located at h = 12 ± 3 μm from the stone. This standoff distance h was used to model the collapse of the largest bubbles approximating their initial shape as an oblate spheroid¹⁸ with dimensions of 2 Ry = 62 μm and 2Rx = 40 μm.

The collapse of a microbubble was modeled by solving the multicomponent Euler equations with a 2D-axisymmetric code. Specifically, the code models the gas and the liquid as a two-phase compressible flow using the following system of equations:²⁸ 

Here, α₁ and α₂ are the volume fractions of the gas and the liquid; vector u is the velocity of the flow. The gas and the liquid are governed by their equations of state $ek=ek(ρk,pk)$⁠, where ρ_k, e_k and p_k are the density, the internal energy, and the pressure of phase k. The gas in the bubble obeys the ideal-gas equation of state with $γgas=1.1$ (the ratio of specific heats γ for C₄F₁₀ is 1.07).¹⁴ The water obeys the stiffened-gas equation of state with γ_water = 6.12 and $p∞,water=3.43×108$ Pa.²⁹ The gas-liquid mixture has the density $ρ=α1ρ1+α2ρ2$⁠. Likewise, the pressure P in the gas-liquid mixture is the sum of partial pressures $α1p1$ and α₂p₂.

The right-hand side of Eq. (1) describes the relaxation of pressure between the phases, where $δp≡p1−p2$⁠, a coefficient μ determines the speed of relaxation (section 1 of the  Appendix), and interfacial pressure p_I is $(z2p1+z1p2)/(z1+z2)$⁠. Here $zk=ρkck$ is the acoustic impedance of the phase k with c_k being the speed of sound of the corresponding phase.

Due to $p1≠p2$⁠, the total-energy equation of the mixture is replaced by the internal-energy equation for each phase. Nevertheless, the mixture-total-energy equation of the system can be written in the usual form

where $E=(e+1/2)||u||2$ and e are the total and internal energies. Although Eq. (2) may appear redundant (as the internal-energy equations are solved for both phases), it is important to ensure the energy conservation and, hence, to preserve a correct treatment of shocks.

The system of equations (1) was solved with a two-substep approach (section 1 of the  Appendix) using adaptive time steps and an AMR algorithm (section 3 of the  Appendix). In comparison with the numerical modeling in POMA,²⁷ the use of an improved system of equations (1) allows us to take into account expansion and compression of each phase in mixture regions. This improvement and the use of the AMR algorithm (section 3 of the  Appendix) allowed us to increase the accuracy of the results, including the spatial and time resolution of gas-liquid interfaces and shock fronts.

To reduce the computation time and needed resources, the dynamics of the collapsing bubble was modeled with the following simplifications. First, the model was axisymmetric, although the experiments showed some apparent divergency from axial symmetry, in part, due to the irregular surface of urinary stones [inset in Fig. 2(b), Fig. 9]. Second, the wave scattering by the irregular surface of the urinary stone was not modeled in this work and the stone was modeled as a rigid plane. Third, the absolute pressure far from the bubble was assumed to be a constant P_∞.

The assumption of a constant pressure P_∞ was an approximation as the driving acoustic pressure varied during the growth-collapse cycle of the bubbles (Fig. 4, Mm. 1, and Mm. 2). The time for bubbles to reach their maximum size varied from bubble to bubble and increased as the bubbles grew from cycle to cycle (Mm. 1 and Mm. 2). Here, we model the bubble that reached its maximum size at t = 0, which was about 45 degrees into the positive pressure phase of the driving acoustic wave (⁠$−0.3 ≲ t ≲ 0.9 μs$⁠, Fig. 4, Mm. 1 and Mm. 2). The time-average pressure of this “sinusoidal” half-cycle during the time of the collapse (⁠$0 ≲ t ≲ 0.9 μs$⁠) was about 0.7 of the amplitude of the driving acoustic wave. In the experiments, the driving acoustic wave was a superposition of sine waves at two frequencies giving the free-field pressure amplitude of 1.4 ± 0.4 MPa (Fig. 4, Mm. 1, and Mm. 2). The pressure at the urinary stone was higher by ∼ 50%³⁰ due to constructive interference of the incident and scattered waves, increasing the pressure amplitude to ∼2.1 MPa. This amplitude gave the time-average value for the absolute pressure during the collapse (⁠$0 ≲ t ≲ 0.9 μs$⁠) of P_∞ = 1.55 MPa. This time-average pressure was assumed to be spatially uniform within the computational domain (768 × 768 μm, Fig. 10), which was much smaller than the spatial characteristics of the driving acoustic beam (e.g., −6-dB zone of 3 × 6 cm, Fig. 3).

The initial pressure in the bubble P₀ depended on the amount of gas and vapor in the bubble. To bracket the range of possible behaviors, we modeled the collapse with two initial gas pressures: P₀₁ = 10 Pa and P₀₂= 100 kPa. The pressure P₀₁ would occur in the bubble polytropically expanded from a nucleus with an initial radius of ∼1.2 μm to the maximum size observed in the high-speed image (Rx = 20 μm, Ry = 31 μm, and h = 12 μm, Fig. 5) if vaporization, condensation, and gas diffusion were negligible. The pressure P₀₂ would occur in the expanded bubble if vaporization, condensation, and gas diffusion were infinitely fast. The pressure P₀₂ overestimates the pressure and amount of gas in the bubble, potentially cushioning the collapse. In the experiments, the pressure in the expanded bubble (t ≈ 0) was between P₀₁ and P₀₂.

The initial density of gas in the bubble was $ρ01=2.78×10−4 kg/m3$ at P₀₁, and $ρ02=1.2 kg/m3$ (density of air at ∼10–35 °C) for P₀₂. The initial density of water was $ρ∞=1000 kg/m3$⁠. BCs and the computational domain (Fig. 10) are described in section 2 of the  Appendix.

Numerical results are shown in Mm. 3 and Mm. 4.

The content of these movies is illustrated in Fig. 6, showing a representative frame from Mm. 3. In both movies, the top left panels show pressure in color on a logarithmic scale ranging from 1 kPa (blue) to 0.5 GPa (red). Pressures outside of this range are shown either in dark red (>0.5 GPa) or in dark blue (<1 kPa). The color image is overlaid by the volume fraction of gas shown in black and white, with an opacity function to render translucent surfaces. The black shows regions of high gas content. The opaqueness decreases with volume fraction until the gas volume fraction is zero (100% liquid) depicted as 100% transparent.

The bottom left panels in Fig. 6, Mm. 3, and Mm. 4 show numerical schlieren φ of mixture-density gradients

where β is a scaling factor for simultaneous visualization of waves in both fluids (β_gas = 20 and β_water = 200).³¹ 

Thick lines in these images suggest the position of the bubble-water interface. Thin lines typically show moving fronts of pressure waves. The different contrast of shock fronts in these movies show that the bubble with smaller amount of gas (P₀₁, Mm. 3) produced stronger shocks than the bubble with the larger amount of gas (P₀₂, Mm. 4). The difference in contrast of shocks is also seen in Figs. 8(a)–8(b) showing the final stages of the collapse. These results support the notion that the gas in the bubble cushions the collapse, so that bubbles with smaller amount of gas produce stronger collapses.

Both experimental [Fig. 7(a)] and numerical [Fig. 7(b)] images showed a distinctive circumferential narrowing of the bubble (t ≈ 0.4 – 0.7 μs), indicating greater bubble-wall velocities toward the axis of symmetry of the bubble. Figure 7(c) shows bubble-wall velocities v₁ (squares) and v₂ (circles) measured at the moving points 1 (square) and 2 (circle) indicated in the inset. Numerical velocities are shown for both initial gas pressures by thin (P₀₁) and thick (P₀₂) lines. Both bubbles had similar dynamics but the bubble with smaller amount of gas (P₀₁) developed greater bubble-wall velocities and collapsed faster than the bubble with P₀₂. The averaged bubble-wall velocities measured in HS-camera images during the final stage of the collapse (t > 0.6 μs, squares and circles) were in between the numerical curves. This behavior might be expected as in the experiments the gas pressure at t = 0 was between P₀₁ and P₀₂.

Both experimental (dots) and numerical (lines) velocities show that the bubble-wall velocity directed toward the axis of symmetry of the bubble v₁ was greater than the velocity v₂ directed toward the stone (Fig. 7). This disproportion led to the circumferential pinching of the bubble, splitting it into two parts. The distal-to-stone part of the bubble was smaller than the proximal part and collapsed first (Mm. 3 and Mm. 4, Fig. 8). The collapse was intensified by a pressure surge with an amplitude of ∼1.23 GPa at P₀₁ and ∼0.37 GPa at P₀₂. This pressure surge was produced by the converging flow of liquid when it collided at the axis of symmetry of the bubble (0.719 μs in Mm. 3 and 0.763 μs in Mm. 4). The collision formed two axial microjets directed away and toward the stone. The reentrant jet hitting the stone was either a single (P₀₁) or a ring-type (P₀₂) jet (Fig. 8). This jet had the velocity v_jet shown in Fig. 7(c) with the maximum of 4117 m/s at P₀₁ and 1164 m/s at P₀₂.

The impact of the reentrant microjet on the stone produced a hydraulic shock with the water-hammer pressure of ∼0.65 GPa at P₀₁ and ∼0.5 GPa at P₀₂ seen as the first pressure spikes at y = 0 in Fig. 8(c). This water-hammer pressure was followed by a longer somewhat triangular pulse with the stagnation pressure of ∼0.3 GPa at P₀₁ and ∼0.2 GPa at P₀₂ [y = 0, Fig. 8(c)].

The impact of the microjet hitting the stone was intensified by the pressure waves radiated by the collapse of the distal part of the bubble (0.721 μs in Mm. 3 and 0.769 μs in Mm. 4, frames 3 in Fig. 8). These pressure waves also intensified the collapse of the main part of the bubble located proximal to the rigid surface. This proximal bubble toroidally collapsed (frames 4–6, Fig. 8) producing pressure waves at the stone with local maximums initially located at y ≈ 6.2 μm at P₀₁ (0.748 μs, Mm. 3) and y ≈ 9 μm at P₀₂ (0.798 μs, Mm. 4).

The pressure waves radiated by the toroidal collapse converged toward the axis of symmetry of the bubble, producing the overall maximum pressure at the stone at P₀₁ (0.75 μs, Mm. 3). As the toroidal bubble collapsed from the periphery radiating multiple waves (frame 6 in Fig. 8, Mm. 3 and Mm. 4), the location of the maximum pressure was slightly off-axis and not seen in Fig. 8(c). The maximum pressure determined along the entire surface of the stone—rather than at specified locations as in Fig. 8(c)—is shown on the right panel in Fig. 6, Mm. 3 and Mm. 4. The maximum pressure reached 6.88 GPa at P₀₁ (0.75 μs, Mm. 3) and 1.32 GPa at P₀₂ (Mm. 4).

At P₀₂, the greater amount of gas in the bubble cushioned the collapse, so that the pressure waves radiated by the toroidal collapse of the proximal bubble were not as strong as at P₀₁ (frame 6 in Fig. 8, t ≈ 0.798 μs, Mm. 4). These pressure waves collapsed the daughter bubbles that were produced by the distal bubble. In turn, the collapse of these daughter bubbles radiated pressure waves that drove the rebound of the bubble at the stone surface near the axis. The collapse of this rebounding bubble produced the overall maximum pressure at the stone surface at P₀₂ (1.32 GPa at 0.843 μs, Mm. 4).

Micro-CT images of a urinary stone before and after the action of SSA microbubbles showed removal of stone material and formation of microcracks (Fig. 9). Specifically, arrows 1–8 point to places of removed material in the 3 D rendered surface of the stone [Fig. 9(a)], whereas arrows 10 and 11 point to the formation of microcracks visible in the cross sections of the stone [Figs. 9(b) and 9(c)].

The combination of HS-video microscopy and numerical modeling employed here showed a pronounced difference between the previously reported dynamics of initially spherical bubbles^22–25,27 and the dynamics of the SSA microbubbles at the surface of urinary stones. The SSA microbubbles expanded nonspherically (Mm. 1 and Mm. 2) and, at their maximum expansion, acquired the shape similar to an oblate spheroid [Figs. 4(a) and 7(a)]. Approximating the cross-section of bubbles as an ellipse truncated by the plane rigid surface [Fig. 5(a)], we modeled the collapse of the microbubbles at the surface of urinary stones. Both the numerical modeling and the experimental observations showed a circumferential constriction pinching the bubbles (Figs. 7–8, Mm. 1–Mm. 4). The pinching was absent in the collapse of the initially spherical bubbles as shown in Fig. 8 in POMA.²⁷ 

Furthermore, whereas for the initially spherical bubbles the axial jet originated from the distal surface of the bubble,^22–25,27 the circumferential pinching of the SSA microbubble split the bubble into two parts producing two axial microjets directed away and toward the stone. The reentrant jet hitting the stone was either a single (P₀₁, Mm. 3) or a ring-type (P₀₂, Mm. 4) jet. The velocity of this jet v_jet reached 4117 m/s at P₀₁ and 1164 m/s at P₀₂ [Fig. 7(c)], supporting the notion that microjets can reach high subsonic or even supersonic speeds.^23,24

The quantitative verification of the predicted velocities of these microjets requires fine spatial (∼1 μm) and temporal (∼1 ns) resolutions (Fig. 8, Mm. 3 and Mm. 4), both of which are experimentally difficult. The numerical frames (Fig. 8) are shown with frame rates up to one billion frames per second. The fastest frame rate of our state-of-the-art camera was 10 × 10⁶ fps (Mm. 2). Although there are cameras that can record a small number (8–15) of frames at rates up to 200 × 10⁶ fps, capturing the exact moment of collapse is difficult even with laser-nucleated bubbles³² and it is practically impossible with the microbubbles at the surface of urinary stones reported in this work. The main reason is that the bubbles are never exactly the same, and tiny variations in bubble's size shifts the moment of collapse, making the capturing of the fine details during the short moment of the collapse technically challenging. Some verification of the numerical model can be provided by comparison with other accepted models and analytical solutions (see section 4 of the  Appendix). The comparison shows a good agreement with the established models and analytical solutions, supporting the trustworthiness of the present numerical results.

It has been suggested that the velocity of microjets is increased by the reflection of the driving wave from the solid surface.²⁴ Here, we showed that microjets produced by SSA microbubbles were also accelerated by the pressure surge when the circumferential flow collides on the axis of symmetry of the bubble (frame 2, Fig. 8). Even though short-lived, this pressure surge was two-to-three orders of magnitude greater than the amplitude of the driving wave (Sec. IV E).

The impact of the jet on the stone produces a hydraulic shock with water-hammer pressure $∝ρcvjet$⁠, where ρ is the density and c is the speed of sound.^21–23 Both the density and the speed of sound grow with pressure,¹⁶ but even at atmospheric pressure (ρ ≈ 1000 kg/m³ and c ≈ 1.5 km/s) the impact of a jet with v_jet ∼1 km/s produces a hydraulic shock of ∼1.5 GPa.

Using geometrical acoustics, Zhong and Chuong³³ developed a theoretical formulation to model the impingement of cavitation microjets and the resultant shock waves in elastic solids, and Zhong et al.³⁴ applied it to consider the impact of the microjets produced by the collapsing bubbles to break urinary stones in shock wave lithotripsy. The impact pressure at the stone boundary and stress-strain at the propagating shock fronts in the stone were calculated for the microjets with a jet diameter of 0.1 mm and velocities in the range of 100–400 m/s. The model predicts that, depending on the contact angles, the compressive stress induced by the jet at the surface of the stone can vary from 0.82 to 4 times that of the water-hammer pressure. The model also predicts that the jet-induced shear stress in the stone can achieve a magnitude of 30%–54% of the water-hammer pressure. Comparison of the model predictions with material failure strengths of urinary stones have suggested that jet impact can lead to stone surface erosion by combined compressive and shear loadings at the jet impacting surface, and spalling failure by tensile forces at the distal surface of the stone.³⁴ 

Numerical simulations of the collapse of a single cavitation bubble near a deformable plane suggest that the material experiences stresses that are much lower than the fluid generated impulsive loads and that the resulting pit characteristics depend not only on the impulsive load amplitude but also on its duration and spatial extent.³⁵ The duration of the hydraulic shock is determined by the ratio of the radius of the jet to the speed of sound in the liquid.²² For example, a microjet with a radius of ∼1.5 μm produces a pressure spike of ∼1 ns. This water-hammer pressure is followed by a longer pulse with the stagnation pressure $∝(1/2)ρvjet2≈0.5$ GPa.²² The duration of this somewhat triangular pulse (∼10–30 ns, Fig. 8) is determined by the ratio of the length to the velocity of the jet. The ratio of the stagnation pressure to the water-hammer pressure is on the order of the Mach number $∼vjet/c$⁠.³³ This ratio was greater at P₀₁ than at P₀₂ due to the greater velocity of the jet at P₀₁ (Fig. 7).

The bubble at P₀₁ collapsed faster and to a smaller minimum size than the bubble with a larger amount of gas, P₀₂. The gas in the bubble cushioned the collapse so that the bubble with P₀₂ produced weaker shocks than the bubble with P₀₁. The difference is apparent by the contrast of shocks in schlieren images and in graphs showing the pressure at the stone (Fig. 8, Mm. 3 and Mm. 4).

As described in Sec. IV F, the maximum pressures (6.88 GPa at P₀₁ and 1.32 GPa at P₀₂) were produced near the axis of symmetry of the bubble. At P₀₁, this maximum was due to constructive interference of pressure waves converging toward the axis from the toroidal collapse of the bubble ring, as previously described.³⁶ The maximum pressure spikes had durations ∼1 ns and quickly decayed off-axis. Longer (∼30 ns) pulses covered larger areas and, regardless of the amount of gas in the bubble, had substantial pressure amplitudes [Fig. 8(c)].

The present numerical results, however, have the following caveats. First, the greatest pressures were, in part, due to constructive interference of waves under 2D-axisymmetric assumption. Although this assumption significantly reduces the computation time and required resources, the actual bubble dynamics is three dimensional, e.g., due to the non-planar geometry of the stone surface [Figs. 2(b) and 9]. Hence, some highly localized peak pressures simulated here under ideal axisymmetric conditions might not be present at the irregular surface of urinary stones. Second, the stone was modeled as a rigid surface, whereas the finite acoustic impedance and elasticity of urinary stones have been shown to affect the predicted pressure amplitudes.^24,33–36 Third, the pressure during focusing events becomes nearly singular and, to remain finite, depends on actual dissipative mechanisms. The influence of viscosity on pressure generated by the collapsing bubbles was not investigated in this report. Fourth, the maximum pressure converges only slowly and, to some extent, is sensitive to the coarseness of the underlying numerical grid. In the present work, we used an AMR algorithm which showed greater peak pressures than the numerical simulations with a coarser grid.²⁷ The extent to which the coarseness of the grid affects pressure amplitudes was not investigated in this work. For these reasons, we believe that the present modeling results can provide only an order-of-magnitude assessment of bubble-wall velocities and pressures generated during the collapse.

Chahine and Hsiao³⁶ numerically modeled the effect of material compliance on pressure peaks showing that their magnitudes depend on the level of deformation of the material. Softer materials with large deformations experienced impact pressures at lower magnitudes than that at a rigid wall, e.g., the pressure at a rubber plate was one half of the pressure magnitude at the rigid wall.³⁶ Kobayashi et al.³⁰ modeled a shock wave–bubble interaction near rigid and soft boundaries during shock wave lithotripsy. The modeling shows that pressures and velocities of microjets were reduced by approximately one half when the bubble was collapsing near a soft boundary (fat, liver, and gelatin) in comparison with that in front of a urinary stone.

Ohl et al.³⁷ studied the dynamics of a non-equilibrium bubble near bio-materials. The bubble oscillations were initiated “explosively” by increasing the internal pressure in the bubble by orders of magnitude in comparison with the reference pressure in the liquid. The modeling showed that the direction of microjets depended on the ratio of densities of elastic material ρ_m and liquid ρ_l. The bubble (standoff distance ∼1) jetted toward the elastic boundary with $ρm>1.4ρl$ and away from the boundary with $ρm<0.7ρl$⁠. For materials with $ρm≈ρl$⁠, the bubble profiles showed a circumferential narrowing of the bubble and splitting it into smaller bubbles with microjets directed toward and away from the elastic boundary.

Here, we modeled the collapse of a solitary bubble. The bubble dynamics captured by high-speed video microscopy (Mm. 1 and Mm. 2) showed a transition from single to multiple bubbles. Ideally, 3 D models with bubble-bubble interaction will be more appropriate.^36,38 The 3 D modeling shows that interaction between bubbles significantly influences the pressure loading on the material surface and requires extensive studies as both enhancement and negative interference of the interaction on the resulting damage can be seen.³⁶ 

The modeling showed that the solitary bubbles can break into daughter microbubbles (Fig. 8, Mm. 3 and Mm. 4). These modeling results are supported by experimental observations of bubble dynamics in the free field^15,39 and the present observations at the surface of urinary stones showing that the collapse of solitary bubbles can produce daughter microbubbles forming clusters of bubbles.

The daughter microbubbles produced at P₀₂ were larger than at P₀₁ due to larger amount of gas in the parent bubble at P₀₂. In both cases, the daughter microbubbles produced by the distal part of the bubble were seen to collapse during the passage of the pressure waves radiated by the toroidal collapse of the proximal part. These waves (as well as the microjet directed away from the stone) promoted mixing and the flow of daughter microbubbles away from the stone, increasing the cluster size in subsequent acoustic cycles (Mm. 1 and Mm. 2).

It has been suggested that negative pressures (not shown) produced during the collapse can cause secondary cavitation.²⁵ Even though some schlieren images (Mm. 3 and Mm. 4) may suggest nucleation in regions with low pressure, those are visualization artifacts due to different ability of gas and liquid phases to perceive the compression and expansion. Nucleation of nanobubbles was not modeled and would require consideration of surface tension.

Surface tension and the lipid shell were neglected in this model. The surface tension produces the Laplace pressure, which for bubbles greater than several micrometers was about two orders of magnitude smaller than the driving pressure. During the implosion, the Laplace pressure increases but remains small in comparison with the pressures developed in the collapsing bubbles. Including the surface tension in the calculation only alters some details of the collapse and raises slightly the maximum pressure peaks.²⁵ The Laplace pressure is also countered by surface-active molecules in the bubble's shell, which may act to stabilize microbubbles against the dissolution.⁴⁰ Lipid-shell effects on bubble dynamics are generally meaningful only for bubbles oscillating at small amplitudes near their equilibrium radius R_eq, buckling at 0.99 R_eq and rupturing at about 1.05 R_eq.^13,14

Zhao et al.² studied lipid-shell microbubbles adhered to a flexible cellulose boundary and insonified by three cycles of 1.5-MHz ultrasound. The bubbles were seen to collapse toroidally indicative of the formation of a jet and had eccentricity of 0.71 ± 0.06, similar to the eccentricity ϵ = 0.63 ± 0.08 observed in the present study (Fig. 5).

The deviation of the bubble's shape from the spherical geometry (ϵ = 0) is mainly due to the close proximity of the boundary affecting the flow around the bubble. The flow is hindered between the bubble and the boundary, displacing the centroid of the bubble toward the solid wall.²⁵ The displacement may be considered as a result of the secondary Bjerknes forces⁴¹ acting between the growing and collapsing bubble and its mirror image.¹² It has been shown that an isolated hemispherical bubble together with its mirror image behaves like a free spherical bubble.¹² It is interesting, therefore, to approximate the SSA microbubble and its mirror image as a single sphere and estimate the time of collapse τ_c using the Rayleigh expression for a radially symmetric collapse of an empty spherical cavity: $τc≈0.915Rmaxρ/P∞$⁠, where R_max is the initial radius of the sphere at rest and P_∞ is the pressure in the bulk liquid.¹⁰ Taking R_max equal to the major semi-axis of the bubble modeled in this work Ry = 31 μm and P_∞ = 1.55 MPa, one can estimate the time of collapse as τ_c ≈ 721 ns. This τ_c corresponds to the third frame in Fig. 8(a) giving a good estimation for the time of emission of the first shock waves by the collapsing bubble. Although the time of collapse might be assessed with this simple spherical model, the dynamics of the collapse was more nuanced as shown in the present work.

The local curvature of the stone surface and incident wave orientation affect the bubble dynamics and collapse, potentially resulting in a complex translation and deformation of bubbles (Mm. 1 and Mm. 2). The incident wave produces the primary Bjerknes force⁴¹ toward the proximal surface of the stone translating the bubbles along the side and curved surfaces of the stone (not shown). The secondary Bjerknes forces⁴¹ deform and translate the bubbles attracting them into pits, cracks, and crevices.⁴² 

Various directions of the microjet and the migration of the bubble centroid have been illustrated in a recent 3 D modeling of a bubble at various standoff distances from a rigid flat wall and subjected to ultrasound propagating parallel to the wall.³⁸ In this situation, the bubble was acted on by the primary Bjerknes force parallel to the wall and the secondary Bjerknes force toward the wall. The modeling showed that the bubble centroid was moving toward the wall and along the ultrasound direction. A decrease of the standoff distance was associated with an increase of the centroid migration to the wall and a decrease of its translation along the direction of ultrasound. Furthermore, a decrease of the standoff distance was associated with a progressive change of the direction of the microjet toward the wall. The authors noted that the bubble system absorbs the energy from the ultrasound and transforms the uniform momentum of the ultrasound parallel to the wall to the highly concentrated momentum of a high-speed liquid jet pointing to the wall.³⁸ This notion suggests the potential significance of microbubbles for the breakage of side surfaces of urinary stones.

Urinary stones are typically curved and have various flaws (e.g., microcracks and crevices) both at the surface and within the stone [Figs. 2(b) and 9].⁴³ Some flaws can be seen with X-ray computer tomography (Fig. 9) and their presence correlates with the efficiency of stone breakage in shock wave lithotripsy.⁴³ We speculate that the mechanical action of microbubbles collapsing in cracks and crevices may cause stone erosion to a greater extent than the bubbles at the planar surfaces due to a greater stress concentration within cracks and crevices.

Borkent et al.²⁰ modeled the collapse of a bubble with the shape of a spherical segment stemmed from an air-filled cylindrical pit with the radius five times smaller than the initial radius of the bubble. Similar to our observations, their modeling shows rushing of liquid toward the axis and pinching-off the bubble with the formation of fast needlelike water jets directed toward and away from the substrate surface. The pinching produced a small distal bubble which was neglected from further modeling.²⁰ Here, we showed that the collapse of the distal bubble radiated pressure waves that can intensify the microjet toward the stone as well as the collapse of the proximal bubble. Furthermore, the distal bubble produced daughter microbubbles, which rebounded driven by the collapse of the proximal bubble. As described in Sec. IV D, the rebound of these daughter bubbles produced pressure waves that, in turn, drove a rebound of a microbubble at the stone. At P₀₂, this rebounding bubble produced the overall maximum pressure at the stone.

The present numerical simulations predicted substantial pressures in the close proximity to the collapsing bubbles (Fig. 8). The direct verification of these predictions would require measurements of highly localized short pressure spikes with temporal resolution on the order of 1 ns and spatial resolution on the order of a few micrometers. Such measurements are technically difficult due to limited frequency response and finite dimensions of pressure sensors. In this work, indirect evidence of the substantial pressure produced by the microbubbles were provided by the micro-computed tomography of urinary stones before and after the action of the SSA microbubbles. The examination of the 3 D surface [Fig. 9(a)] showed removal of stone material, creation of pits and crevices, whereas the cross sections of the stone [Figs. 9(b) and 9(c)] showed the formation of microcracks.

In principle, cavitation bubbles are capable of breaking even the strongest materials, such as stainless steel, within a fraction of a second.¹⁵ Such intense cavitation, however, has been observed with bubbles driven at lower frequencies (∼20 kHz) and higher acoustic amplitudes (tens of MPa).^15,32 In the present work, the microbubbles were driven with ultrasound in the sub-MHz range and pressure amplitudes on the order of 1.5 MPa. Even with these driving amplitudes, the numerical modeling (Fig. 8, Mm. 3 and Mm. 4) suggests that the SSA microbubbles can produce local pressures exceeding the amplitudes of shock pulses used in shock wave lithotripsy to break urinary stones.

In shock wave lithotripsy, shock pulses consist of the positive-pressure phase of ∼0.5–2 μs with peak pressures of ∼15–150 MPa and the negative-pressure phase of ∼2–5 μs with pressures of 5–20 MPa.^6,44–46 Lithotripter pulses are delivered at ∼0.5–2 Hz; in contrast, the driving considered here produces thousands of growth-collapse cycles per second.⁴⁷ Furthermore, in shock wave lithotripsy, cavitation bubbles collapse under static pressure of 0.1 MPa.^10,11,48 In contrast, the collapses investigated here are intensified by the ∼1.5-MPa pressure of the driving wave. Whereas the ∼1.5-MPa driving wave is insufficient to cause stone breakage by itself,⁷ the SSA microbubbles can produce local pressures sufficient to facilitate breakage of urinary stones.⁵ 

It has been shown that a collapsing bubble can concentrate the relatively weak and disperse energy density of the driving acoustic wave by 12 orders of magnitude.⁴⁹ Let us consider the energy conversion from the incident acoustic wave to bubble expansion and to jet impact during one acoustic cycle (⁠$−1.5 ≲ t ≲ 0.9 μ$s, Fig. 4, Mm. 1 and Mm. 2). The negative pressure phase of the incident acoustic wave (⁠$−1.5 ≲ t ≲−0.3 μ$s) expanded the bubble providing the momentum to the surrounding liquid. The inertia of the liquid continued to expand the bubble into the positive pressure phase (⁠$−0.3 ≲ t ≲ 0.9 μ$s) driving the bubble to reach its maximum volume V_b under a positive pressure P_L—resulting in a potential energy $PLVb$ in the surrounding liquid. The inertial growth was followed by the inertial collapse, concentrating the majority of the energy onto a small region during a short interval at the final stage of the collapse. This concentration of energy by the growing and collapsing bubble produced energy flux densities orders of magnitude greater than that of the incident acoustic wave.

Specifically, the modeling suggests that the microjet impinging the stone had the velocity $vjet∼1$ km/s (Fig. 7), thereby having the kinetic energy flux density $ρvjet3/2∼108$ W/cm². In comparison, the incident acoustic wave with pressure amplitude $Pa∼1.5$ MPa (Fig. 4, Mm. 1 and Mm. 2) had the peak instantaneous intensity $Pa2/(ρc)∼150$ W/cm²—six orders of magnitude smaller than the energy flux density produced by the microjet.

The kinetic energy of the impinging jet was orders of magnitude greater than the energy delivered to the impact area during the whole period T = 2.4 μs by the incident acoustic wave. For example, the microjet impacting the area A_j with the radius of ∼1 μm (Fig. 8) for $τjet∼7$ ns (e.g., 724–731 ns, Mm. 3) had a kinetic energy on the order of $τjet Aj ρvjet3/2∼10$ nJ. In comparison, the energy transmitted with the incident acoustic wave during a period T through the same cross-sectional area A_j was $AjTPa2/(2ρc)∼0.01$ nJ—three orders of magnitude smaller than the kinetic energy of the microjet.

The kinetic energy of the jet was only a fraction of the potential energy of the bubble $PLVb$⁠. Specifically, the bubble modeled in Mm. 3 and Mm. 4 had the initial volume $Vb≈72·10−9$ cm³, while the absolute pressure in the liquid was P_L = 1.55 MPa, so that the potential energy of the bubble was $PLVb≈112$ nJ. A significant portion of this energy was concentrated and released during the toroidal collapse of the main (proximal to the stone) portion of the bubble, radiating intense pressure waves right at the surface of the stone (Fig. 8, Mm. 3 and Mm. 4). Due to energy dissipation (e.g., on heating of gas in the bubble), the intensity of these pressure waves decreased with the amount of gas in the bubble and was smaller at P₀₂ than at P₀₁ (Fig. 8, Mm. 3 and Mm. 4).

In summary, the present results suggest that, for a broad range in the amount of internal gas, the collapsing bubbles produced peak pressures exceeding the driving pressure by two-to-three orders of magnitude [Fig. 8(c), plots in Mm. 3 and Mm. 4]. This focusing of pressure is the basis for the development of treatment modalities with driving pressures substantially lower than those required without accumulating microbubbles. The significance of these results is likely relevant to many applications dealing with growth and collapse of microbubbles in close proximity to rigid boundaries,^1–7,47 and future work will explore such applications.

Physical experimentation and numerical modeling have been combined to gain a detailed understanding of the dynamics of microbubbles engineered to accumulate on the surface of stones. The nonspherical form assumed by SSA microbubbles has been precisely characterized, with an observed form of truncated spheroid broadly consistent with previously presented observations of bubbles in close proximity to solid and elastic boundaries.^1,2,18,19 For the observed standoff distances and microbubble geometry, the numerical model showed a distinctive circumferential pinching forcing a microjet to the stone. SSA microbubbles can concentrate the relatively weak and disperse energy of driving acoustic waves, focusing this energy into the stone with pressure spikes some two-to-three orders of magnitude greater than the amplitude of the driving wave. These results provide insights on the use of SSA microbubbles in combination with extracorporeal acoustic energy sources with peak pressures on the order of 1 MPa as a potential approach to treat urinary stones.

We thank Dr. Marshall Stoller for valuable discussion, suggesting and arranging micro-CT, and for providing urinary stones, Dr. Sunita Ho for micro-CT of urinary stones, Dr. Ryoji Shiraki for chemical analysis of stone composition, Dr. Kyle Morrison (Sonic Concepts) for piezoelectric transducers, Dr. Matt Hopcroft for the development of the driving acoustic system and software for scanning of acoustic beams, David Bell for help with programming of the signal generator, and Todd Rumbaugh (Hadland Imaging) for advice on high-speed imaging. T.K. is a founding member of Applaud. Y.P., W.B.P., and D.L. are employees/investigators for Applaud. Experiments where performed at Applaud, the numerical modeling in Caltech. T.C. and K.S. acknowledge support from the Office of Naval Research (ONR, N0014-18-1-2625). T.C. and K.M. also acknowledge support from the ONR (N00014-17-1-2676) and the NIH (P01-DK043881).

The numerical approach consists of two steps. First, the system of equations (1) is solved without the relaxation terms (⁠$μ→0$⁠). The resultant (pressure-disequilibrium) hyperbolic equations are solved using a Godunov-type method; the Riemann problem is solved with a HLLC solver.²⁸ Then, in the pressure-relaxation step, $μ→+∞$ and the numerical method converges to solving the following five equations of the mechanical-equilibrium model of Kapila et al.:⁵⁰ 

Here, $K∇·u$ accounts for the compression and expansion of both phases in mixture regions, where K is the ratio

Further details about the pressure-relaxation algorithm and the re-initialization procedure ensuring the conservation of total energy can be found in Saurel et al.²⁸ 

The 2 D-axisymmetric assumption allowed using a computational domain in the form of an xy-plane with the x axis directed along the axis of symmetry of the bubble (Fig. 10). The left boundary of the domain was at the stone surface modeled as an acoustically-rigid plane. The top and the right boundaries assumed non-reflective BCs mimicking an infinite domain.

Spatial and time resolutions of shock fronts and gas-liquid interfaces were improved by using the AMR algorithm.⁹ This shock- and interface-capturing algorithm was combined with a MUSCL-like scheme, which provides second-order accuracy in both space and time. In addition to the AMR and high-order methods, the THINC interface-sharpening method⁵¹ was used to maintain the interface relatively sharp.

Figure 10 (inset) shows the computational grid at t = 0 illustrating the AMR algorithm. At every time step, the AMR re-meshes the computational domain by progressively subdividing the numerical grid into smaller steps reaching the maximum grid-line density in specified regions. These regions are defined by significant transition of calculated parameters (such as pressure, density, or volume fraction) and include contact discontinuities, shock fronts and other fast changing pressure waves, as well as gas-liquid interfaces. If the region with fine grid-line density moves to another location (e.g., following the motion of the bubble wall or shock waves), the mesh at the previous location is coarsened. In this work, the maximum density of grid lines was set to 400 cells per unit length giving a spatial resolution at the finest refinement level of 0.05 μm per grid line.

The present model and numerical method have been validated for 1 D single phase and multi-phase Riemann problems, and several selected multidimensional flows: liquid/gas shock tube,^8,9,52 surface-tension test,^9,53 and cavitation in air/water mixture.⁵² Regarding bubble dynamics, we compared the present numerical modeling results for spherical bubble collapse (on a 3 D grid) with analytical solutions in the weakly compressible limit.⁵² The comparison showed a good agreement with discrepancies less than or around one percent. We also compared (data not shown) maximum pressure at the wall for different standoff distances using the present model and other established models and methods.^54–56 The other methods used reduced models with thermodynamic equilibrium assumption (equilibrium of velocity, pressure and temperature between phases in mixture regions) and the Tait equation of state. In the present study, we use a model assuming mechanical equilibrium (velocity and pressure) but no thermal equilibrium, and the stiffened-gas equation of state. All modeling results were performed with the same grid resolution (100 cells per diameter of the bubble) and yielded the same order of magnitude maximum pressures at the rigid wall and similar trends in function of the standoff distance, supporting the trustworthiness of the present modeling results.