2.1. Experimental set-up

Figure 2 illustrates the experimental platform, which consists of a laser-focused unit for generating cavitation bubbles, a schlieren imaging system for capturing shock waves, a dual-angle synchronized imaging system and a three degrees of freedom (3-DoF) mobile platform for precisely adjusting the position of the sample plate. To better visualize bubble collapse patterns and shock waves, a transparent acrylic plate was vertically mounted on the 3-DoF mobile platform, perpendicular to the horizontal plane. Alternatively, the aluminium sample plate was mounted horizontally to study cavitation erosion, as shown in the subfigure in figure 2. These two scenarios are equivalent in terms of bubble dynamics at the corresponding $\gamma$ values since their optical path of the laser focusing is parallel to the plates and the buoyancy effect can be ignored at such small buoyancy parameter $\delta =\sqrt {\rho g R_{max } / p_{\infty }}\sim \mathcal{O}(10^{-2})$ (Reuter et al. Reference Reuter, Deiter and Ohl2022a ). Experiments were conducted under ambient conditions, with temperatures ranging from $20^{\,\circ}$ C to $25^{\,\circ}$ C and a surrounding pressure of $p_\infty =10^5$ Pa.

Figure 2.

The schematic diagram of the experimental platform highlights several key components: a laser-focused cavitation bubble generation unit; a schlieren imaging optical path; a dual-angle synchronized imaging system (plan-view); a bubble-induced material erosion set-up located in the upper right-hand inset (side-view).

A Nd:YAG pulsed laser (Nimma-900, pulse duration 8 ns) with a wavelength of 532 nm and pulse energy ranging from 20 to 40 mJ. The quartz water tank, measuring 100 mm $\times$ 100 mm $\times$ 100 mm, was filled with deionized water to a depth of 90 mm. To achieve a more precise laser focusing, a 30 mm diameter sapphire glass with a thickness of 0.25 mm was installed into a hole on the side of the quartz tank. The laser beam passed through the microscope objective (M Plan Apo L 10 $\times$ , NA = 0.28), the sapphire glass and a convex lens, was focused near the centre of the water tank to generate cavitation bubbles with $R_{max } \sim {O}(1)$ mm. This laser-induced bubble system has been well validated in our previous work, as detailed in Li et al. (Reference Li, Zhao, Zhang and Han2024). Additionally, a number of convex lenses with different focal lengths ( $\phi =25.4$ mm, $f$ = 30 mm, 40 mm and 50 mm) were used in front of the sapphire glass to adjust the focusing angle of the laser, while the laser energy was continuously varied to alter the plasma shape.

A schlieren imaging system, designed to capture the shock waves emanating from bubble collapse, features an optical path depicted in figure 2. A nanosecond high-frequency pulsed laser (Cavilux HF, wavelength 810 nm, pulse duration 10 ns) was employed to enhance the sensitivity of the schlieren imaging instead of the traditional light-emitting diode (LED) point light source. This laser is expanded by a plano-convex lens with a 100 mm focal length to create a parallel beam of approximately 20 mm in diameter that passes through the cavitation region. The wavefront, possessing a density significantly higher than the ambient medium, results in non-uniformities in the refractive index of the surrounding flow field when the bubble collapses, thereby deflecting light rays. The light rays are subsequently focused by a schlieren lens ( $\textrm{F} = 100 $ mm) and intercepted by a sharp knife edge, producing distinctive dark streaks in the imagery. The transient dynamics of bubble nucleation and collapse are captured by a high-speed camera I (Phantom V2012) coupled with a LAOWA macro lens (100 mm, F2.8), at a resolution of $128 \times 64$ pixels, recording at 6 40 000 f.p.s. Simultaneously, high-speed camera II (Phantom V12.1, frame rate 1 80 000 f.p.s.) was used to capture side-view bubble oscillations via the shadow method from another angle of view. Uniform background illumination across the water tank was provided by diffusing two 300 W LED lights through frosted glass. The exposure of camera II was set to 1 $\unicode {x03BC}$ s, and a filter was positioned in front of the lens to shield the sensor.

Figure 3.

Surface characterization of a pure aluminium plate with ultrasmooth treatment: ( $a$ ) the surface roughness before cavitation erosion; ( $b$ ) the plane view of cavitation erosion distribution; ( $c$ ) three-dimensional (3-D) reconstruction of cavitation erosion distribution.

The size of pure aluminium samples employed for cavitation erosion tests is 10 mm $\times$ 10 mm $\times$ 2 mm, and the surface is treated to achieve ultrasmoothness using a precision diamond machining tool (Moore Nanotech UP250L). Figure 3(a) manifests the surface roughness of the processed metal samples, measured by a white light interferometer, indicating a smooth, flat surface with nanometre precision and submicron geometric accuracy. To avoid contamination from water droplets and dust, the metal sample is cleaned with ethanol prior to the erosion tests. During the erosion experiment, side views were captured using camera I with the shadow method ( $128 \times 128$ pixels, 4 30 000 f.p.s.), and the stand-off distance parameters were measured and adjusted.

Surface damage of the metal sample after erosion tests are examined using a super-depth-of-field microscope (Olympus DSX1000). Its measurement accuracy in the $x$ – $y$ plane is $\pm 3\, \%$ , with a repeatability of $\pm 2\, \%$ , while the height measurement repeatability in the $z$ -axis direction is $\sigma _{n-1} \leqslant 1 \ \unicode {x03BC}\text{m}$ . As depicted in figure 3(b), the details of pitting holes are discernible. The microscope can also output high-resolution 3-D images and depth contours through optimized imaging algorithms, as shown in figure 3(c), providing precise quantitative data for subsequent analysis.

2.2. Numerical methodology

The numerical simulation of the dynamic pulsation of cavitation bubbles in proximity to a solid boundary involves solving the compressible and viscous multiphase flow problem. The present study utilized CavBubbleFoam solver, which is a customized and rigorously validated version of CompressibleInterFoam built in OpenFOAM-4.X. Notably, this solver is specifically designed for bubble dynamics, incorporating modifications to the phase fraction field $\alpha$ and ensuring meticulous bubble mass conservation (Koch et al. Reference Koch, Lechner, Reuter, Köhler, Mettin and Lauterborn2016; Zeng et al. Reference Zeng, Gonzalez-Avila, Dijkink, Koukouvinis, Gavaises and Ohl2018, Reference Zeng, An and Ohl2022; Reese et al. Reference Reese, Schädel, Reuter and Ohl2022). Meanwhile, the solver accounts for surface tension but neglects heat and mass transfer. The governing equations comprise the continuity equation and the Navier–Stokes equations,

(2.1) \begin{align} \frac {\partial \rho }{\partial t}+\nabla &\cdot (\rho \boldsymbol{u}) =0 ,\end{align}

(2.2) \begin{align} \frac {\partial \rho \boldsymbol{u}}{\partial t}+\nabla \cdot (\rho \boldsymbol{u} \boldsymbol{u})& =-\nabla p+\nabla \cdot \boldsymbol{\tau }+\boldsymbol{f}_{\sigma },\end{align}

where $\rho$ represents the fluid density, $\boldsymbol{u}$ denotes the velocity, $p$ is the pressure. Here $\boldsymbol{f}_{\sigma }$ represents the surface tension term, modelled using the continuous surface force method (Brackbill, Kothe & Zemach Reference Brackbill, Kothe and Zemach1992), which can be expressed as

(2.3) \begin{equation} \boldsymbol{f}_{\sigma } = \sigma \left [ -\nabla \cdot \left ( \frac {\nabla \alpha _1}{|\nabla \alpha _1|} \right ) \right ] \nabla \alpha _1, \end{equation}

here $\sigma = 0.07\,{\rm N\,m}^{-1}$ is the surface tension coefficient between the gas and water. The viscous stress tensor $\boldsymbol{\tau }$ is given by

(2.4) \begin{equation} \boldsymbol{\tau } = \mu \left (\nabla \boldsymbol{u} + \nabla \boldsymbol{u}^{\mathrm{T}} - \frac {2}{3}(\nabla \cdot \boldsymbol{u}) \boldsymbol{I}\right ), \end{equation}

where $\mu$ is the dynamic viscosity and $\boldsymbol{I}$ is the identity tensor. The volume of fluid method is employed to capture and characterize the interface between gas and liquid phases,

(2.5) \begin{align} \frac {\partial \alpha _{1}}{\partial t}+\nabla \cdot (\alpha _{1} \boldsymbol{u})+\nabla \cdot \left (\alpha _{1}\alpha _{2} \boldsymbol{U}_{r}\right )& =\alpha _{1}\alpha _{2}\left (\frac {\psi _{1}}{\rho _{1}}-\frac {\psi _{2}}{\rho _{2}}\right ) \frac {\mathrm{D} p}{\mathrm{D} t}+\alpha _{1} \nabla \cdot \boldsymbol{u} ,\end{align}

(2.6) \begin{align} \alpha _{2}& = 1 - \alpha _{1},\end{align}

where $\alpha _i$ represents the volume fraction of the phase, with subscripts $i=$ 1 and 2 denoting the liquid and gas phases, respectively, and $\boldsymbol{U}_{r}$ is the relative velocity between the two phases. The isothermal compressibility for both fluids, denoted as $\psi _i=\mathrm{D} \rho _i / \mathrm{D} p$ , can be calculated using the Tait equation of state,

(2.7) \begin{equation} p=\left (p_{0}+B\right )\left (\frac {\rho }{\rho _{0}}\right )^{\gamma }-B. \end{equation}

The coefficients of (2.7) are listed in table 1.

Table 1.

The coefficients for Tait state equation of liquid and gas phases.

2.3. Numerical set-ups and validation

Figure 4.

Schematic diagram of the numerical simulation domain. Levels of mesh refinement are indicated by the contours on the symmetry plane.

A schematic diagram shown in figure 4 illustrates the general set-ups for the numerical simulations. The $x$ – $o$ – $y$ plane was designated as a symmetric boundary, allowing for sparing half the computational domain and improving computational efficiency. As a result, the size of the computational domain was set to $94 R_{max } \times 57 R_{max } \times 47 R_{max }$ , with the maximum equivalent radius of the simulated bubble being approximately 530 $\unicode {x03BC}$ m in this study. Similar to Beig, Aboulhasanzadeh & Johnsen (Reference Beig, Aboulhasanzadeh and Johnsen2018), wave-transmissive boundary conditions were implemented on all external boundaries, except the no-slip bottom wall, to emulate an extended fluid domain and minimize the impact of boundary constraints on the simulation outcomes. The simulation utilized a structured mesh with the coarsest cell with $\Delta x=0.73R_{max }$ at the outermost boundary. After eight levels of successive refinement, the finest cell width in the vicinity of the bubble and the adjacent wall was reduced to $\Delta x=0.0057R_{max }$ , ensuring adequate resolution in these critical regions. In line with Reuter & Kaiser (Reference Reuter and Kaiser2019) and Reuter & Ohl (Reference Reuter and Ohl2021), a persistent liquid film, modelled by fixing $\alpha _1=1$ at the wall nodes, prohibits direct contact between the bubble and wall. In the simulation, the initial state of the laser-induced bubble was modelled by introducing a small volume of high-pressure, non-condensable gas into the flow field, with a phase fraction of $\alpha _{2} = 1$ . Initially, we employed the spherical bubble theory (Keller & Kolodner Reference Keller and Kolodner1956; Zhang et al. Reference Zhang, Li, Cui, Li and Liu2023) and iteratively adjusted the initial pressure $P_0$ and radius $R_0$ in an infinite domain. By setting $P_0 = 4800$ bar (within the value range used by Lechner et al. (Reference Lechner, Lauterborn, Koch and Mettin2020), Reese et al. (Reference Reese, Schädel, Reuter and Ohl2022) and Zeng et al. (Reference Zeng, An and Ohl2022)) and $R_0 = 0.055R_{max }$ , we achieved a good agreement between the simulated maximum radius and oscillation period with experimental data. In subsequent simulations, with the bubble’s $R_{max }$ measured at 530 µm, we selected an initial radius $R_0=29$ µm. The effect of phase transition, which contributes to vapour condensation in physical experiments, is accounted for – as demonstrated in the studies by Raman et al. (Reference Raman, Rosselló, Reese and Ohl2022), Fan et al. (Reference Fan, Bußmann, Reuter, Bao, Adami, Gordillo, Adams and Ohl2024 and Ohl, Reese & Ohl (Reference Ohl, Reese and Ohl2024)) – by applying an attenuation coefficient of 0.34 to the bubble’s internal pressure at its maximum volume during the first and second expansions.

Figure 5.

$(a)$ Definitions of non-dimensional stand-off distance $\gamma$ , aspect ratio $\eta$ and the axial asymmetry ratio $\xi$ . $(b)$ Comparison of the time evolution of the bubble radius obtained from the experiment and numerical simulation for the case of $\gamma =0.72$ , $\eta =1.8$ and $\xi =1.9$ . The error bar indicates the uncertainty in the experimental measurements. Non-dimensional times, defined as $t^{*}=t/(R_{max } \sqrt {\rho / P_{\infty }})$ , and the dimensionless radius is expressed as $R^{*}=R/R_{max }$ . $(c)$ Comparison of experimental and numerical bubble profiles, where the red curve represents the numerical bubble isosurface for $\alpha =0.5$ , with experimental time in the bottom right-hand corner and simulation time (in red) in the bottom left-hand corner.

Figure 5 $(a)$ presents the non-dimensional parameters to characterize the cavitation bubble. It has been demonstrated in Liang et al. (Reference Liang, Linz, Freidank, Paltauf and Vogel2022) that the plasma generated by laser usually exhibits a non-spherical shape, resembling a spindle or teardrop shape. In addition, the conventional stand-off distance, the aspect ratio $\eta = (a+c)/2b$ and the axial asymmetry parameter $\xi = {c}/{a}$ are introduced to more precisely reproduce the influence of the initial bubble geometry, which is also the focus of the present study. The initial shape of a laser-induced cavitation bubble can be approximated by two conjoined semiellipsoids sharing an axis of length $2b$ . Specifically, $\eta$ reflects the overall flattening or longitudinal extension of the spindle shape, while $\xi$ quantifies the asymmetry between the left and right semiellipsoids. It should be noted that when dealing with small $\gamma$ values, the maximum bubble radius corresponds to the radius of the bubble’s upper hemisphere at its maximal expansion, as shown in figure 5 $(a)$ , Reuter et al. (Reference Reuter, Deiter and Ohl2022b ) also employed this method, referring to it as a direct image measurement method. In experiments, the initial shape of the bubble is defined at the first frame after the plasma flashing. To validate the accuracy of the numerical simulations, we selected conditions of $\gamma =0.72$ , $\eta =1.8$ and $\xi =1.9$ , and compared the results with experimental data. Figure 5 $(b)$ presents the time evolution of the bubble radius, demonstrating excellent agreement between numerical simulations and experimental data during the first oscillation cycle. Figure 5 $(c)$ further compares the bubble dynamics over the first two oscillation cycles, with the red contours representing the bubble interface from numerical simulations. During the expansion phase, the bubble appeared nearly circular in the top view (frames 1–2) and shrank to its minimum volume at $t^{*}=2.10$ , forming a relatively uniform ring-shaped structure that matched well with the red contour from the experiment. In the second oscillation cycle, the initial asymmetry of the bubble was amplified, leading to increased curvature on both sides of the toroidal bubble (frame 4), which collapsed preferentially (frame 5), and finally converged circumferentially to two points perpendicular to the optical axis (frame 6). It is worth noting that the simulated bubble volume in the second oscillation cycle is slightly larger than the experimental measurements (frame 4), and the collapse speed is slightly lower, possibly due to the simplification of the bubble phase change process in numerical simulations.

3.1. Cavitation erosion patterns

A series of cavitation erosion tests were conducted at 17 different stand-off distances ranging from $\gamma = 0.58$ to $2.20$ , with an interval of approximately 0.1. In terms of the experimental findings presented in this section, the optical path was adjusted to maximize the sphericity of the laser-induced bubbles, aiming for values of both $\xi$ and $\eta$ approaching unity. For each individual test, the laser was activated 900 times (repetition rate 5 Hz) using the same energy level, resulting in an approximate maximum radius of the induced bubbles of 0.97 mm. Based on microscopic observations of all tested specimens, as shown in figure 6, the patterns of cavitation erosion were categorized into five typical types: bipolar; monopolar; annular; solar-halo; central. Detailed description of the erosion patterns are elaborated as follows.

(i) The bipolar pattern is frequently observed at relatively low stand-off distances ( $\gamma \lt 0.8$ ) within the scope of the present study. It exhibits the most severe local erosion damage symmetrically aligned at the opposite sides of a circular erosion zone, whose connecting line is either parallel to or perpendicular to the axis of the laser beam. Within the circular cavitation zone, there are numerous densely distributed small pitting corrosion pits. However, these indentations have much smaller diameters and depths compared with the erosion in the ‘polar’ regions.

(ii) The monopolar pattern occurs within a large stand-off parameter range ( $\gamma =0.7\sim 1.4$ ), with its primary characteristic being that the most severe erosion appears at one ‘polar’ of the circular erosion zone. The distribution of small pits within the circular erosion area differs between the monopolar and bipolar patterns, with the former exhibiting less overall damage. The erosion distributions of the monopolar pattern at three different $\gamma$ are present in figure 6, indicating significant differences in the extent and range of erosion in the ‘polar’ region. Particularly, there is significant erosion at $\gamma =1.31$ , which is comparable to that observed in the bipolar pattern.

(iii) Annular erosion patterns dominate at $\gamma \approx 1.5$ , with the erosion pits primarily distributed along a ‘ring-like’ circumference. Inside the annular region, there are hardly any significant occurrences of erosion. Compared with the bipolar and monopolar modes, the annular pattern generally possesses an axisymmetric feature, which is more in line with the conventional understanding of cavitation erosion. However, the circumferential asymmetrical patterns, rather than the axisymmetric pattern, manifest the peak erosion damage.

(iv) With the stand-off distance further increased to approximately $\gamma = 1.6\sim 1.8$ , the solar-halo pattern emerges, characterized by the concentrated distribution of pits at the centre and along the circumference of the circular zone. However, no large or deep pits are observed in either of these areas.

(v) For $\gamma$ exceeding 2.0, only small indentations appear at the centre, which is referred to as the centre pattern. At $\gamma$ = 2.2, the surface shows almost no significant damage, with only slight traces visible. It can be anticipated that increasing the stand-off distance further will make the cavitation erosion caused by bubble collapse negligible.

3.2. Quantification of cavitation erosion

In this subsection, quantitative analysis was performed on the cavitation erosion of the tested aluminium specimens. First, single-bubble induced cavitation erosion experiments were conducted, as illustrated in figure 7(a). Since some shallow indentations have depths on the submicron scale, they are almost invisible under bright-field observation and can only be seen under oblique illumination. Therefore, we used a white light interferometer (gbs WLI) for 3-D pit reconstruction and volume measurement. For comparison, the experimental results of Philipp & Lauterborn (Reference Philipp and Lauterborn1998) are reproduced in figure 7(a). It is worth noting that they employed a theoretical estimate of $R_{max }$ in the free field, derived from the Rayleigh equation, to normalize the stand-off distance. Specifically, $R_{max }=T_{\textrm { c}}\sqrt {(p_\infty -p_v)/\rho }/0.915$ , where $T_{\textrm {c}}$ denotes the half-period of the first bubble oscillation cycle. In contrast, this study utilizes directly measured values of $R_{max }$ for normalization. Since the Rayleigh equation tends to overestimate $R_{max }$ for bubbles near a wall, the resulting normalized stand-off distance, defined as $\gamma =d/R_{max }$ , is correspondingly underestimated. To facilitate valid comparison, we also used the Rayleigh equation to calculate $R_{max }$ , and the average value of the data points at the corresponding smallest $\gamma$ to normalize the damage volume. The overall trend of the variation of erosion volume with respect to $\gamma$ observed in our study is consistent with that in Philipp & Lauterborn (Reference Philipp and Lauterborn1998). The maximum erosion volume occurs at the smallest $\gamma$ within the investigated range for both studies. As $\gamma$ increases, the normalized cavitation erosion volume, $V^*$ , initially decreases, then increases and decreases again. Subsequently, another prominent peak appears at $\gamma \approx 1.3$ , followed by a gradual decrease towards zero. However, there are some differences: the peak volume we measured at $\gamma = 0.5$ corresponds to the trough reported in Philipp & Lauterborn (Reference Philipp and Lauterborn1998). Additionally, within the range of $\gamma = 0.5$ to $1.3$ , our results show a greater average value and higher dispersion, especially the peak around $\gamma = 0.9$ .

In addition, we also performed a quantitative measurement of the erosion volume for $N = 900$ as discussed in § 3.1. It should be noted that this study did not directly use the 3-D reconstruction function of the microscope to obtain the damaged fields of all specimens for the $N = 900$ cases, as direct 3-D measurement is highly time-consuming and has a limited field of view. In addition, some specimens exhibited slight overall bending deformation, which could introduce significant errors in direct measurements. Taking the typical bipolar pattern shown in figure 6 as an example, figure 8 illustrates how to obtain the quantitative distribution of cavitation erosion through image-based operations. The detailed method can be found in Appendix A. The specimens were divided into four groups based on different stand-off distances, and different-sized pits were randomly selected as sampling data from each group. The volumes of the pits, $V_{\textit{pit}}$ , were measured using 3-D microscope. According to Dular, Delgosha & Petkovšek (Reference Dular, Delgosha and Petkovšek2013), erosion pits formed on the surfaces of soft metal materials like pure aluminium after plastic deformation are usually approximately semi-ellipsoids. Therefore, the equivalent pit depth can be calculated as $h=3V_{\textit{pit}}/(2\pi R_{\textit{pit}}^2)$ . Figure 9 shows the fitted quadratic relationship between the equivalent radius and depth of the axisymmetric ellipsoid pit based on the sampled data: $h = -0.0004R_{\textit{pit}}^{2} + 0.1244R_{\textit{pit}} + 0.031$ (in micrometres). Regression analysis proves the good predictive ability of this quadratic curve. Finally, a contour map of the cavitation erosion pit depths is provided in figure 8 $(c)$ using the fitted curve. This method offers high detection resolution for small-scale pits.

Figure 7.

Variation of the overall volume of cavitation erosion pits within the range of moderate stand-off distances. Non-dimensionalization is performed using the average damage volume corresponding to the smallest $\gamma$ . $(a)$ In comparison with the experimental results reported by Philipp & Lauterborn (Reference Philipp and Lauterborn1998), the bubble number is $N = 1$ . The $R_{max }$ is calculated using the formula $R_{max }=T_{\textrm { c}}\sqrt {(p_\infty -p_v)/\rho }/0.915$ . $(b)$ A comparison between $N = 1$ and $N = 900$ cases, where $R_{max }$ is directly measured from the images. The solid marks denote the corresponding average values, based on which the corresponding fitted curves are obtained.

Figure 8.

The illustration of the operations to quantify the cavitation erosion. $(a)$ Pixel threshold selection. $(b)$ Binary-based pits image processing with simultaneous edge detection and equivalent circular pits analysis. $(c)$ Contour map of cavitation erosion pit depths.

Based on the quantitative analysis method, figure 7 $(b)$ illustrates the overall variation in erosion damage across moderate stand-off distances, specifically within the range $0.58 \leqslant \gamma \leqslant 2.20$ , and compares with the cases $N = 1$ . To maintain consistency in the methodology throughout the study, the value of $\gamma$ is calculated using direct measurement of $R_{max }$ . Compared with the $N = 1$ case, the variation trend of cavitation volume with $\gamma$ for $N = 900$ is generally consistent, although the peak near $\gamma = 1$ is less pronounced. Moreover, since each specimen underwent hundreds of bubble cycles under the same conditions, random errors were reduced, resulting in significantly lower overall dispersion. At $\gamma = 2.2$ , the extent of cavitation erosion is already considerably small, and it is expected to further diminish monotonically with further increasing distance. Additionally, the fitted curves of the average cavitation erosion volume for both cases exhibit similar shapes. In conjunction with the discussion in Appendix B regarding the impact of bubble repetition count on the cavitation erosion distribution, it is observed that erosion caused by a single bubble collapse is random. However, as the number of bubbles increases, a statistically consistent erosion pattern begins to emerge, indicating that the bubble evolution process and collapse mode remain consistent across occurrences.

Figure 9.

The relationship between the equivalent radius and depth of the pits, regarded as axisymmetric ellipsoid shape, on the aluminium specimens after the 900 times of bubble collapse.

3.3. Characteristics of bubble collapse

Figure 10.

The side-view snapshots of laser-induced bubble at different stand-off distances: (a) $\gamma = 0.54$ , $R_{max }$ = 0.96 mm; (b) $\gamma = 1.05$ , $R_{max }$ = 0.96 mm; (c) $\gamma = 1.50$ , $R_{max }$ = 0.93 mm; (d) $\gamma = 1.78$ , $R_{max }$ = 0.95 mm; (e) $\gamma = 1.95$ , $R_{max }$ = 0.90 mm. The red, yellow and blue boxes mark the moments when bubbles collapse to the minimum volume for the first, second and third time, respectively. The green line above each frame represents the position of the wall. The blue arrow denotes the direction of gravity, and the green light path indicates the direction of laser incidence.

The pulsation characteristics of near-wall cavitation bubbles were investigated, focusing on their collapse modes. Figure 10 presents snapshots of bubble evolution at different instants for bubbles with nearly identical $R_{max }$ but varying stand-off distances from the wall. As shown in figure 10( $a$ ), the pulsation behaviour of bubbles nucleating close to the wall resembles the process schematized in figure 1( $b$ ) for bubbles at $\gamma \lt 1$ . However, a notable observation not indicated in the schematic diagram is that the bubble fragmented into a cluster of microbubbles in the later stages of the second collapse. For bubbles with a stand-off distance approximately equal to their maximum radius, as depicted in figure 10( $b$ ), the bubble reached its minimum volume at the end of the first collapse, generally coinciding with the moment when the jet impacted the wall. It is worth noting that, in comparison with the case at $\gamma =0.54$ , the bubble at $\gamma =1.05$ did not make direct contact with the wall during its initial collapse, thereby mitigating the wall’s hindrance to the radial contraction of the bubble. Consequently, this resulted in a more rapid collapse and a smaller final volume of the bubble. Additionally, the bubble was pierced by the jet and torn into two distinct parts, as shown in frame 6 (Philipp & Lauterborn Reference Philipp and Lauterborn1998; Yang, Wang & Keat Reference Yang, Wang and Keat2013). The interaction between these two parts amplifies instability, resulting in a pronounced asymmetry in the bubble pulsation during the second cycle, as shown in figure 11( $b$ ). At $\gamma =1.50$ , as illustrated in figure 10( $c$ ), the bubble pulsation at this greater stand-off distance corresponds to the dynamics shown in figure 1( $a$ ). Whereas, a unique phenomenon of ‘counter jet’ was observed during the second cycle (frames 7 and 8), where a slender jet ejected the central microbubbles away from the wall (Lindau & Lauterborn Reference Lindau and Lauterborn2003; Zhang Reference Zhang2021). Subsequently, the toroidal bubble cluster near the wall collapsed, and even after the second and third collapses, it maintained a relatively regular annular shape, with an increasing radius (figure 11 $c$ ). At $\gamma =1.78$ (figure 10 $d$ ), the bubble’s pulsation during the first cycle resembles the dynamics described in the third case. However, the bubble exhibited faster contraction and jet formation, reaching a smaller minimum volume at the end of the first cycle. During the second collapse, the bubble surface evolved into a ‘cylindrical shape’ and eventually split into two parts in the middle (frame 8). The upper portion then contracted towards the centre of the wall, experiencing central collapse. In the subsequent third cycle, the residual microbubbles in the region away from the wall expanded on the wall, forming a larger ring-shaped microbubble cluster and undergoing a third near-wall collapse. Further increasing $\gamma = 1.95$ (figure 10 $e$ ), the upper portion of the ruptured bubble in the second cycle did not yet contact the wall (frame 9). As the bubble continued to migrate, the microbubble cluster eventually reached the wall in the third cycle (frame 10), but the interaction between the bubbles and the wall weakened.

Figure 11.

The plan-view snapshots of typical regimes of bubble pulsation near a rigid wall: (a) bipolar, $\gamma = 0.59$ , $R_{max }$ = 0.95 mm; (b) mini-annular, $\gamma = 1.06$ , $R_{max }$ = 0.79 mm; (c) annular, $\gamma = 1.52$ , $R_{max }$ = 0.83 mm; (d) solar-halo, $\gamma = 1.77$ , $R_{max }$ = 0.84 mm; (e) central, $\gamma = 1.98$ , $R_{max }$ = 0.85 mm. The red, yellow and blue boxes mark the moments when bubbles collapse to the minimum volume for the first, second and third time, respectively. The blue arrow denotes the direction of gravity, and the green light path indicates the direction of laser incidence.

We further provide the plan view of bubble pulsation under conditions similar to those in figure 10, though not captured simultaneously. From figure 11, we can clearly distinguish the axial asymmetric collapse process of these cavitation bubbles and categorize the bubble collapse modes. We have identified six distinct collapse modes: bipolar; monopolar; annular; mini-annular; solar-halo; central. Specifically, frame 1 of figure 11( $a$ ) shows that the jet tip was not axisymmetric when it contacted the wall, but rather flat, with its long axis parallel to the laser beam. This indicates that flow instability developed rapidly under small distance parameters, yielding significant asymmetry in the initial stage of bubble collapse in the first cycle. Subsequently, the jet head becomes teardrop-shaped due to obstruction by the wall (frames 2 and 3). Although the toroidal bubble at its minimum volume is regular (frame 4), the surface of the rebounded bubble is much rougher compared with the smooth surface in the first cycle (frames 5 and 6). As discussed in Plesset (Reference Plesset1954), this phenomenon was attributed to Rayleigh–Taylor instability. Moreover, the toroidal bubble exhibited a faster expansion rate in the direction perpendicular to the laser beam, then formed a ‘lip-shaped’ bubble during contraction (frame 7). Eventually, the ‘lip-shaped’ bubble was pinched off, and then the separated sections underwent collapse along the circumferential direction during the second cycle. Both sections violently contracted to the middle points on the upper and lower ends (frame 9), thus referred to as the bipolar mode of bubble collapse. It is obvious that the bipolar bubble collapse is strongly related to the bipolar pattern of erosion distribution shown in figure 6. When the collapse intensity at one end of the ring is much higher than the opposite side, or starts at a point and propagates circumferentially, this asymmetric process is known as the monopolar collapse mode (see figure 12 $b$ ). This mode frequently occurs and is associated with the monopolar erosion pattern within the range of $\gamma$ from 0.5 to 0.9. If the laser-induced bubble has good initial sphericity, a ‘lip-shaped’ bubble would not form, and the bubble collapse in the second cycle would remain essentially axisymmetric, namely the annular collapse. The influence of the initial bubble shape on collapse modes will be detailed later in § 4.1. Figures 10( $b$ ) and 11( $b$ ) both show the mini-annular bubble collapse mode at $\gamma$ approximately equal to 1.0. This mode is characterized by the considerably small toroidal bubble radius ( $R_s$ ) formed during the secondary collapse, typically within the range of $R_s \lt 0.5 R_{max }$ . Interestingly, within the range of $\gamma$ dominated by the mini-annular mode, the cavitation erosion pattern is mostly monopolar, and the region where the most severe cavitation occurred is even beyond $R_s$ , similar to the findings reported by Wen et al. (Reference Wen, Yao, Wu, Sun, Yang and Zhong2023). The relevant mechanism will be discussed in § 3.4. When $\gamma$ is around 1.50, as shown in figure 11( $c$ ), the toroidal bubble splits into a cluster of microbubbles during the second contraction cycle, they are generally uniformly distributed along the circumferential direction, corresponding directly to the annular erosion pattern. By comparing figures 11( $d$ ) and 11( $e$ ), similarities can be observed between the solar-halo mode and the central mode during the first and second cycles. However, for the solar-halo mode, damage in the ‘solar’ region was caused by the bubble collapse in the second cycle, while the collapse of the third cycle led to damage in the ‘halo’ region. On the other hand, the central mode exclusively caused damage in the centre as a result of the collapse of the second cycle, and the impact of the third collapse on the wall surface can be disregarded.

In general, it can be concluded that there is a strong correlation between bubble collapse modes and the associated cavitation erosion patterns, except that monopolar distribution of damage can also be attributed to bubble collapse of the mini-annular mode, in addition to bubble collapse of the monopolar mode.

3.4. Cavitation erosion mechanisms

Figure 12.

Comparison of shockwave front propagation caused by bubble collapse at different $\gamma$ . Each set consists of the plan-view schlieren images on the left and side-view snapshots on the right. The bubble collapse modes and characteristic parameters are: $(a)$ bipolar, $\gamma = 0.63$ , $R_{max } = 1.01$ mm; $(b)$ monopolar, $\gamma = 0.73$ , $R_{max } = 0.97$ mm; $(c)$ mini-annular, $\gamma = 1.08$ , $R_{max } = 0.92$ mm; $(d)$ annular, $\gamma = 1.45$ , $R_{max } = 0.97$ mm; $(e)$ central, $\gamma = 1.85$ , $R_{max } = 0.98$ mm. The blue arrow denotes the direction of gravity, and the green light path indicates the direction of laser incidence. The green dot indicates the laser path is perpendicular to the paper, pointing outward.

Figure 13.

The close-up views of the bipolar pattern of erosion damage. The blue and red circles represent the position of the toroidal bubble at the end of the first and second collapses, respectively. The red arrows indicate the load direction inferred from the plastic deformation of the material surrounding the pit.

Figure 14.

Schematic diagram of the shockwave focusing mechanisms on inducing severe erosion damage in the bipolar and monopolar patterns. The black circular rings indicate the toroidal bubble at the second cycle, while the grey dashed ring of the monopolar II (mini-annular) indicates the separated parts of the annular microbubble cluster.

We demonstrate the link between bubble collapse modes and cavitation erosion patterns, thereby revealing the erosion mechanism at moderate stand-off distances. Initially, we analyse the shock waves produced by four distinct bubble collapse modes through the schlieren images presented in figure 12. The series of five images captures the initial condition of each bubble (frame 1), key moments in the first cycle (frames 2 and 3) and the second cycle (frames 4 and 5).

Figure 12 $(a)$ demonstrates that in the bipolar mode, the axially symmetric collapse of a toroidal bubble during the initial cycle resulted in asynchronous shock waves, forming irregular wavefronts due to complex wave interface. On the other hand, during the second cycle, the asymmetric lip-shaped bubble collapsed along the circumferential direction, generating pronounced shock waves with sharp wavefronts originating from two polar regions. As depicted in figure 13, these two areas precisely correspond to the regions experiencing the most intense cavitation erosion, situated at the upper and lower portions of the red circle. Inspired by Reuter et al. (Reference Reuter, Deiter and Ohl2022a ), the schematic in figure 14 clarifies the shockwave focusing mechanism accountable for this effect. The evenly spaced points along the toroidal bubble represent the emission sources of shock waves induced by bubble collapse. It is assumed that the spherical shock waves propagate at the speed of sound ( $c=1500\,{\rm m\,s}^{-1}$ ) and that the bubble’s circumferential collapse velocity is also $c$ . Initially, the collapse of a bubble shaped like a lip, propagating circumferentially towards the middle at each end. The subsequent shockwave fronts met and superimposed along the direction of the bubble collapse, forming a reinforced wavefront with increased pressure. Furthermore, when the composite wavefronts propagating in opposite directions meet at the middle points, their head-on collision generates a more intense detonation wave focusing, resulting in significantly higher localized pressures at the middle points. As illustrated in figure 13, the stress imposed on the wall due to the focused shock wave far exceeded the yield strength of pure aluminium, resulting in extensive plastic deformation. This deformation manifests as deep pits with diameters greater than 100 $\unicode {x03BC}$ m, along with adjacent wrinkled protrusions. These characteristics indicate that the deformation caused by the impact at high strain rates induced the plastic flow of aluminium, resembling a liquid-like state. Two insets on the right-hand side of figure 13 illustrate that the deformation of small-scale pits outside the blue circle predominantly pointed radially outward. This suggests that the collapse of the toroidal bubble during the first cycle is the primary cause of these pits. However, the deforming direction of small-scale pits inside the blue circle seems more variable. We posit that this variability stems from the complex wave interference inside the toroidal bubble. The animation to further illustrate the cavitation erosion mechanism for the bipolar pattern is provided in the Supplementary movie 1.

The monopolar erosion patterns, as shown in figure 6, are marked by intense cavitation erosion localized to one side of the toroidal bubble. This pattern occurs frequently over a wide range of $\gamma$ values, approximately $0.7 \sim 1.3$ . Unlike the bipolar pattern, which exhibits mirror symmetry, the monopolar pattern breaks this symmetry and may be driven by distinct mechanisms. The first mechanism is similar to the shock wave focusing mechanism observed in the bipolar collapse. As shown in figure 12 $(b)$ , during the second cycle, the toroidal bubble undergoes a circumferentially asymmetric collapse, with a higher curvature on the right-hand side of the ring compared with the left-hand side (frame 4). As a result, the right-hand side of the toroidal bubble was pinched off first, leading to a crescent-shaped bubble, i.e. monopolar collapse mode. As illustrated in the schematic (monopolar I) in figure 14, the shock wavefronts superimposed along the collapse direction of the toroidal bubble and ultimately converged at the left-hand side of the toroidal bubble, inducing the residual bubble to release a strong shock wave (see figure 12 $b$ , frame 5). The second mechanism typically appears in the mini-annular collapse mode. During the contraction phase of the second cycle, the toroidal bubble fragments into multiple discrete small-scale bubbles. These small bubbles then collapsed and emitted shock waves (figure 12 $c$ ). The wavefronts underwent oblique superposition and concentrated in specific regions outside the toroidal bubble, as indicated in figure 14 (monopolar II). This kind of shockwave focusing leads to a lower energy density compared with the first mechanism, thus the overall damage volume in the $\gamma$ range influenced by this mechanism is relatively modest, as depicted in figure 7. At greater distances, the monopolar pattern can also be observed as shown in figure 6 (frame 4). Under these conditions, the toroidal bubble fragmented into several larger sections due to instability (figures 11 $c$ and 12 $d$ ). Upon the bubble collapse, the two aforementioned mechanisms (head-on collision of convergent wavefronts and oblique superposition of wavefronts) may act simultaneously to generate strong shockwave load. As a result, a considerable amount of cavitation erosion volume is apparent in the $\gamma$ range from 1.3 to 1.5.

Overall, the loads that cause severe cavitation erosion at moderate distances are not generated in the first cycle but arise from the circumferentially asymmetric collapse of the toroidal bubble during the second cycle.

3.5. Radial distribution of erosion depth

Figure 15.

Variations in the normalized radius of the annual bubble in the second cycle, $R^*_{s}=R_s/R_{max }$ , with respect to $\gamma$ . The collapse modes are marked, with typical examples shown accordingly. Critical stand-off distance parameters for bubble collapse mode transitions: $\gamma _1 \approx 0.89$ ; $\gamma _2 \approx 1.25$ ; $\gamma _3 \approx 1.56$ ; $\gamma _4 \approx 1.85$ .

In this subsection, we investigate the radial distribution characteristics of cavitation erosion using the same aluminium plate samples discussed in § 3.1. Figure 15 presents the normalized outer radius of the toroidal bubble, $R^*_s$ , and its corresponding collapse mode during the second collapse phase for $0.4 \lt \gamma \lt 2.2$ , with initial bubble shapes having an aspect ratio, $\eta$ , ranging from 1.1 to 1.8. Firstly, the bipolar, monopolar and annular modes coexist for $\gamma$ values between 0.4 and 0.8. Meanwhile, $R^*_s$ initially increases and then decreases, reaching a peak at approximately 0.6. The common feature observed in this range is that the jet impacts the wall before the bubble fully collapses in the first cycle. The interaction between the jet and the wall contributes to the formation of a stable vortex ring, resulting in the toroidal bubble maintaining its intact and continuous shape during the first cycle and enlarging its diameter during the second cycle. In the region dominated by the mini-annular mode of bubble collapse ( $0.89\lt \gamma \lt 1.25$ ), $R^*_s$ is relatively small. This is because the jet impact and collapse to minimum volume occur almost simultaneously in the first cycle before the toroidal bubble is fully formed, leading to a weaker vortex intensity. As the toroidal bubble is already in close proximity to the wall during the second cycle and does not generate strong jets to expand its inner diameter, $R^*_s$ ends up being small. In contrast, during the third range ( $1.25\lt \gamma \lt 1.56$ ), toroidal bubbles that collapse in the mono/annularmodes exhibit higher-intensity jets during the second cycle. These jets, blocked by the wall, result in a radial outward flow that enlarges the inner diameter of the toroidal bubble and contributes to the appearance of a second peak in $R^*_s$ . At $\gamma \gt 1.6$ , most of the potential energy of the bubble is dissipated during the first two collapses, leading to a monotonic decrease in the radius of the bubble ring with increasing $\gamma$ , corresponding to the solar-halo mode. Finally, in the range of $1.8\lt \gamma \lt 2.2$ , a transition occurs to the central mode, characterized by the smallest $R^*_s$ values.

Figure 16.

The radial distribution of average cavitation erosion depth at $\gamma = 0.58$ . A Gaussian multipeak fitting curve is provided, with each Gaussian function for the corresponding peak presented.

Next, we quantitatively analyse the radial distribution characteristics of average erosion depth using the image-based pit reconstruction method proposed in § 3.2, as shown in figure 8. The specific equation used is to obtained average erosion depth,

(3.1) \begin{equation} \overline {h}_i = \frac {\int ^{R_i}_{R_{i-1}}2\pi hR\,{\rm d}r}{\pi \big(R^2_{i}-R^2_{i-1}\big)}, \quad (i=1,2,3, \ldots , n). \end{equation}

Taking the $\gamma = 0.58$ case as an example, figure 16 demonstrates the corresponding radial distribution of cavitation erosion using a multimodal Gaussian function to fit the discrete values of $\overline {h}_i$ $(i=1,2, \ldots , 20)$ . The fitted curve exhibits three peaks, located at $R^*=0$ , 0.35 and 0.65, respectively. The first two smaller peaks primarily result from the impact load of the jet and toroidal bubble collapse in the first cycle, while the largest peak arises from the shock wave focusing resulting from circumferentially asymmetric collapse in the second cycle. Similarly, the fitted curves for the radial distribution of cavitation are plotted for $\gamma =$ 0.58–2.20 in figure 17. To facilitate comparison, the positions of $R^*_s$ and the maximum erosion depth are marked. Overall, for radial positions of $R^*\gt 0.9$ , the erosion depth tends to approach zero, indicating that the cavitation erosion effect mainly occurs within the area where the maximum bubble radius projects onto the wall. There are almost no noticeable erosion pits on the sample surface at $\gamma = 2.2$ , so it is defined as the critical distance for cavitation erosion. When $\gamma \leqslant 0.8$ , the maximum cavity depth is around $R^*=0.6\sim 0.7$ , further substantiating the explanation of the bipolar cavitation erosion mechanism in § 3.3. In addition, the position of the maximum cavity when $\gamma = 0.72\sim 0.8$ is slightly larger than the radius of the toroidal bubble. This phenomenon can be attributed to the combined effects of two shockwave focusing mechanisms: (i) the collision of convergent shockwave front and (ii) the oblique superposition of shock waves. At $\gamma \approx 0.95$ , the erosion level reaches a local maximum in the central region, although its absolute value is relatively small. This is the result of the shockwave focusing within the small toroidal bubble. For $\gamma$ in the range of 1.1–1.2, the toroidal bubble collapses in the mini-annular mode, forming a monopolar cavitation erosion pattern due to the second mechanism of shockwave focusing. The maximum erosion location is significantly greater than $R_s$ , as described in § 3.4. Within the range of $\gamma =$ 1.3–1.5, the erosion depth is significant and the position generally coincides with $R_{s}$ . Within the $\gamma$ range, the monopolar pattern is still primarily caused by the first shockwave focusing mechanism. When $\gamma \gt 1.54$ , we found that the most severely damaged region coincides with the position $R_s$ , indicating that the annular collapse mode is the primary cause of the uniformly distributed pits along the ring. At $\gamma = 1.66$ , the bubble collapse in the solar-halo mode causes relatively severe cavitation erosion damage at a position far larger than $R_{s}$ , indicating that the third collapse of the bubble (figure 11 $d$ , frame 10) begins to influence the distribution characteristics of the pits. For $\gamma \gt 1.8$ , the bubble remains above the wall during the first two pulsation cycles. Its collapse releases spherical shock waves, significantly dissipating energy in the process (figure 12 $e$ ). As a result, cavitation erosion pits predominantly in the central region of the wall, with the peak average depth being less than 0.1 $\unicode {x03BC}$ m.

In conclusion, the radial distribution of cavitation erosion varies significantly at different stand-off distances. However, the regions where severe erosion occurs are mostly located near the toroidal bubble during the second cycle. table 2 summarizes the correspondence between bubble collapse modes and erosion patterns, and lists the salient features associated with each pattern.

Figure 17.

The profiles of radial distribution of cavitation erosion depth from $0.58\leqslant \gamma \leqslant 2.20$ . The dashed line shows the fitted curve of $R^*_{s}$ . The red circles indicate the radial positions with the highest erosion depth at each $\gamma$ .

Table 2.

Correspondence between bubble collapse modes and erosion patterns with typical distribution descriptions.

4.1. The influence of initial non-sphericity on bubble collapse

The results presented earlier, as shown in figure 15, have indicated that the collapse mode of cavitation bubbles is not exclusively influenced by the stand-off distance. Figure 18 compares the collapse modes and shockwave characteristics of bubbles with different initial aspect ratios, $\eta$ , at similar distance parameters in the range of $\gamma$ from 0.6 to 0.7: (i) at $\eta =1.11$ , bubbles with high initial sphericity collapse in the annular mode; (ii) at $\eta =1.41$ , the increased asymmetry in the initial shape leads to the symmetry-breaking instabilities, causing the bubbles to collapse in the monopolar mode; (iii) at $\eta =1.81$ , bubbles with a larger aspect ratio collapse in the bipolar-V mode, with two polar regions located on either side of the laser beam; (iv) at $\eta =2.61$ , bubbles with an extremely elongated shape initially collapse in the bipolar-H mode, with two polar regions distributed axially along the laser beam. The corresponding animations are provided in Supplementary movie 2 for the convenience of readers.

Furthermore, we conducted over 300 experiments within the parameter space of $\gamma =0.4\sim 2.0$ and $\eta =1.0\sim 3.0$ . Based on the experimental observations, we plotted the regime map of bubble collapse modes in figure 19. At $\gamma \lt 0.9$ , the predominant collapse modes are bipolar-H, bipolar-V, monopolar and annular, consistent with the findings in figure 15. In most cases, the monopolar and annular modes exhibit varying degrees of randomness in their occurrence, coexisting within a trapezoidal region of the parameter space. However, when the initial bubble shape is nearly spherical ( $\eta \leqslant 1.1$ ), the collapse exclusively manifests as the annular mode. The bipolar-H is observed only at small stand-off distances ( $\gamma \lt 0.7$ ), with the corresponding critical $\eta$ increasing as $\gamma$ becomes larger. The transition mechanism between the bipolar-H and bipolar-V modes is an intriguing topic that will be further elaborated in § 4.2. For $\gamma$ ranging from 0.9 to 1.2, the mini-annular mode dominates within a parameter region delineated by a yellow triangle, with its dynamic mechanism thoroughly analysed in § 3. As $\gamma$ increases further, mono/annular modes re-emerge in a quadrilateral regime. The bipolar-V modes are extensively distributed within $0.5\lt \gamma \lt 1.5$ , occupying the largest area on the regime map. This indicates that the initial non-sphericity of the bubble has a profound and extensive influence on bubble collapse during the second cycle. The solar-halo and central modes, which dominate at larger stand-off distances, are less sensitive to variations in $\eta$ due to the gradual attenuation of the initial asymmetry during bubble migration towards the boundary.

Overall, the initial shape of the bubble directly influences the distribution of energy and the focusing direction during collapse, particularly for $\gamma \lt 1.0$ . In natural cavitation, bubbles may also exhibit non-spherical shapes upon nucleation. When exposed to vortex flow or non-uniform pressure gradients, bubbles experience compression in certain directions, resulting in deformation into ellipsoidal forms (Dabiri, Sirignano & Joseph Reference Dabiri, Sirignano and Joseph2010; Dular et al. Reference Dular2019; Hsiao, Chahine & Ma Reference Hsiao, Chahine and Ma2020). With increasing $\eta$ , the collapse undergoes a transition from axially symmetric behaviour to circumferentially asymmetric behaviour.

4.2. The mechanism for circumferential asymmetrical bubble collapse

Figure 20.

The evolution of an initially non-spherical bubble collapse near a rigid wall during its first cycle. The instantaneous bubble interface contours sliced from the $x$ – $y$ plane and the $y{-}z$ plane for the bubble ultimately collapsed in ( $a$ ) the bipolar-V mode ( $\gamma =0.72$ , $\eta =1.8$ , $\xi =1.9$ ) and ( $b$ ) the bipolar-H mode ( $\gamma =0.50$ , $\eta =2.2$ , $\xi =2.8$ ). $(c)$ A comparison of the sliced velocity fields for the bipolar-V collapse mode at $t^* = 1.95$ . The white curve showing the numerical bubble isosurface at $\alpha =0.5$ , and the black arrows indicate the velocity vectors.

Three-dimensional numerical simulations were conducted to reveal the mechanism behind circumferentially asymmetric collapse of the bubble, especially the bipolar mode, triggered by initially non-spherical bubbles. Figures 20 $(a)$ and 20 $(b)$ illustrate the time evolution of the contours of the bubbles’ interfaces and flow fields during the first oscillation cycle. The bubbles originated from a ‘teardrop-shaped’ form at nucleation ( $t^*=0$ ) and collapsed in the bipolar-V and bipolar-H modes during the second cycle, respectively.

Figure 21.

Evolution of the vorticity field of the initially non-spherical bubble after jet impact on the wall. ( $a$ ) The bipolar-V mode, $\gamma =0.72$ , $\eta =1.8$ , $\xi =1.9$ ; ( $b$ ) the bipolar-H mode, $\gamma =0.50$ , $\eta =2.2$ , $\xi =2.8$ . The white curve showing the numerical bubble isosurface at $\alpha =0.5$ .

In the bipolar-V mode, during the initial stage of bubble nucleation, the bubble exhibits a larger profile in the $x$ – $y$ plane, requiring the displacement of more surrounding water as it expands. As expected, the greater added mass results in increased inertial resistance, leading to a slightly smaller expansion. When the bubble volume reaches its maximum, its maximum width in the $x$ – $y$ plane is approximately 97 % of that in the $y$ – $z$ plane, leading to a larger curvature at the bubble apex ( $t^{*}=0.94$ ). Based on the Rayleigh equation (Duclaux et al. Reference Duclaux, Caillé, Duez, Ybert, Bocquet and Clanet2007), bubble cross-sections with smaller apex curvature radii generate greater jet momentum during the collapse process. Since the jet maintains a consistent axial velocity across different planes (see figure 20 $c$ ), the jet in the $x$ – $y$ plane, expands more in the transverse direction to increase its volume. In contrast, the jet profile in the $y$ – $z$ plane is relatively narrower. This difference explains the formation mechanism of the elliptical jet tip observed in figure 11 $(a)$ . Additionally, the ‘teardrop shaped’ bubble exhibits asymmetry along its major axis ( $\xi =1.9$ ), causing a deflection angle between the jet axis and the wall normal, as indicated by the blue arrows in figure 20 $(c)$ . A similar phenomenon was observed in frame 4 of the experiment in figure 10 $(a)$ , where the jet typically deflects towards the laser incidence direction. As shown in figure 21 $(a)$ , the elliptical jet tip leads to asymmetric splashing with respect to the wall centre. The upward propagation speed of the annular splashing jet in the $x$ – $y$ plane is slower compared with that in the $y$ – $z$ plane, as shown in frame 1. The impact of the splashing jet on the upper surface of the toroidal bubble results in fragmentation of the vortex ring, where the stronger splashing leads to a relatively larger gas volume being retained within the primary vortex ring, closer to the central axis (see frame 2). Nevertheless, the daughter bubbles coalesce into a larger one as they expand during the second cycle (frame 3). When the bubble reaches its maximum size (frame 4), it is observed that the bubble profile in the $y$ – $z$ plane is larger than that in the $x$ – $y$ plane. Consequently, the toroidal bubble collapses faster and pinches off in the $x$ – $y$ plane. Further, as shown in the 3-D simulation and experimental comparison in figures 22 $(a)$ and 22 $(b)$ , the toroidal bubble exhibits smaller radii of curvature and volume at its left and right ends than at the top and bottom during expansion in the second cycle (see figure 22 $a3$ ). During the subsequent collapse process, the toroidal bubble pinches off into an upper and a lower daughter bubble (see figure 22 $(a4,a5)$ ). From the wall pressure contour maps, it is evident that the water hammer pressure load induced by the jet is only 3.6 MPa, whereas the wall pressure during the first collapse of the toroidal bubble reaches 11 MPa. This indicates that the load during the initial collapse phase of the bubble is relatively small, insufficient to trigger the severe localized erosion phenomena depicted in figure 6. However, as the bubble enters the second cycle, the ‘lip-shaped’ bubble undergoes asymmetric collapse along the circumferential direction, leading to wavefront superposition effects (frame 4). This superposition of wavefronts intensifies the collapse strength at both the upper and lower ends of the bubble, inducing wall pressures up to 63 MPa. This pressure far exceeds the yield strength of pure aluminium, sufficient to cause plastic deformation of the material.

Figure 22.

Comparison between experimental observations and numerical simulation results: bipolar-V, $(a)$ $(b)$ $\gamma = 0.72$ , $\eta = 1.8$ , $\xi = 1.9$ , $R_{max } = 0.53$ mm; bipolar-H $(c)$ $\gamma = 0.61$ , $\eta = 2.61$ , $\xi = 2.3$ , $R_{max } = 0.59$ mm; $(d)$ $\gamma = 0.50$ , $\eta = 2.2$ , $\xi = 2.8$ , $R_{max } = 0.53$ mm. Contour map represents the pressure on the wall surface.

In the bipolar-H mode, as shown in figure 20 $(b)$ , an increase in $\eta$ enhances the curvature difference between the two planes during bubble oscillation, similarly resulting in an elliptical jet tip. The jet in the $y$ – $z$ plane exhibits a significantly wider top radius than in the $x$ – $y$ plane, which limits the volume of the primary vortex ring near the central axis during the collapse phase (figure 21 $b$ , frame 1). When the bubble collapses to its minimum volume, the vortex ring volumes in both planes become nearly identical (frame 2). Note that the position of the cross-section centroid of the bubble in the $y$ – $z$ plane is larger compared with the $x$ – $y$ plane. As a result, as shown in figures 22 $(c)$ and 22 $(d)$ , the top and bottom ends of the toroidal bubble have a smaller oscillation period and collapse earlier in the $y$ – $z$ plane during the second cycle (Li, Prosperetti & van der Meer Reference Li, Prosperetti and van der Meer2020). This is responsible for the pinch-off of the toroidal bubble into a left and a right daughter bubbles. Additionally, as observed in figure 22 $(c\,2)$ , the left end of the toroidal bubble collapses first (possibly due to a larger jet deflection angle), releasing a shockwave, which then converges circumferentially towards the right endpoint. This process ultimately forms a shockwave interference front, similar to the monopolar I mode shown in figure 14. The time phase difference in the collapse of the toroidal bubble on the left- and right-hand sides also carries over to the second cycle. It is worth noting that the wall pressure caused by the bipolar-H mode is only 18 MPa, which is significantly lower than that of the bipolar-V mode.

In summary, the asymmetry of the initial bubble shape results in differences in added mass and inertial resistance during bubble expansion in different planes, thereby forming a jet with an elliptical tip. The variations in the cross-sectional area of the toroidal bubble and the position of the cross-section centroid control the transition between the bipolar-H and bipolar-V collapse modes.

4.3. The impact pressure induced by circumferential asymmetrical bubble collapse

To determine whether the localized high-pressure phenomenon discussed in § 4.2 occurs only when the initial bubble is non-spherical, we conducted five comparative numerical simulations with different initial bubble shapes under the same stand-off distance. The results indicate that variations in the initial shape have little impact on the bubble oscillation period, and detailed analysis is provided in Appendix C. We then compared the time evolution of the pressure field during the first and second collapse stages, as shown in figure 23.

For the first bubble collapse, the pressure curve exhibits a three-peak feature, corresponding to three significant physical phenomena: jet impact on the wall; liquid ‘splashing’; the toroidal bubble collapse. First, as indicated by Peak 1, the jet penetrates the bubble’s lower surface and impacts the wall at a velocity of approximately 70 m s⁻¹, forming a localized high-pressure region. Notably, a thin liquid film of approximately 10.6 $\unicode {x03BC}$ m exists between the bubble’s lower surface and the wall. According to Reuter & Kaiser (Reference Reuter and Kaiser2019), using total internal reflection shadowgraphy, it was observed that for $\gamma \lt 1.0$ , the liquid film thickness during bubble collapse ranges from 5 to 40 $\unicode {x03BC}$ m, suggesting that the jet impact does not directly reach the wall, thereby mitigating the wall pressure to some extent. Through numerical simulations, we found that the maximum jet impact pressure is attenuated by approximately 50 % when transmitted to the centre of the wall through a water film with a thickness of 10.6 $\unicode {x03BC}$ m. A comparison of the curves for different $\eta$ values shows that as the bubble deviates from a spherical shape, the wall pressure induced by the jet decreases as $\eta$ increases. After the jet penetrates the bubble, it spreads radially along the wall, forming an annular liquid film in a direction different from the jet. This phenomenon, often referred to as the ‘splashing effect’ or ‘secondary jet’ by researchers (Tong et al. Reference Tong, Schiffers, Shaw, Blake and Emmony1999; Brujan et al. Reference Brujan, Keen, Vogel and Blake2002), generates significant pressure (Peak 2) when the rapidly expanding splashing edge collides with the upper surface of the collapsing bubble. However, since this phenomenon occurs far above the wall, it induces high local pressure but significantly attenuates before reaching the wall. When the bubble collapses to its minimum volume, Peak 3 is formed, with its peak value approximately ten times that of Peak 1. However, since the collapse of the toroidal bubble is not simultaneous at all locations, multiple pressure peaks are generated at different positions. Notably, when $\eta = 1.8$ , the ellipsoidal initial shape ( $\xi = 1.0$ ) induces a higher flow field pressure compared with the case with $\xi = 1.9$ .

During the second collapse of the bubbles, the peak pressure of the flow field induced by bubble collapse was higher than that of the first collapse, consistent with the analysis of the cavitation erosion mechanism in § 3.4. Furthermore, except for $\eta = 1.0$ , the flow field pressure during the second collapse exhibited a distinct double-peak feature. Figure 24 presents three typical frames from the second collapse cycle: the top row shows the morphological characteristics of the toroidal bubbles during the contraction process (including front and side views); the middle row depicts the tearing state of the bubbles in the late collapse stage; and the bottom row illustrates the moment when the bubbles collapse to their minimum volume.

(i) At $\eta =1.0$ , the toroidal bubble exhibited a uniform thickness and contracted at similar speeds across its surface. The curvature within the ring showed minimal variation, with only slightly elevated local curvatures at the four corners. The bubble tore into four sections and collapsed simultaneously at $t^{*} = 3.49$ , exhibiting a single-peak characteristic. This process generated a significant pressure in the flow field, but the wall pressure reached only 21 MPa.

Figure 23.

Comparison of the temporal evolution of the pressure induced by bubbles ( $R_{max } \approx 0.53$ mm, $\gamma \approx 0.72$ ) with varying initial shapes: $(a)$ the first collapse; $(b)$ the second collapse. A counterpart animation is provided in the supplementary movie 2.

Figure 24.

The 3-D numerical results ( $R_{max } \approx 0.53$ mm, $\gamma \approx 0.72$ ) demonstrate the influence of the initial bubble shape on the second collapse mode. The characteristic parameters are: (1) $\eta = 1.0$ , $\xi = 1.0$ ; (2) $\eta = 1.2$ , $\xi = 1.3$ ; (3) $\eta = 1.8$ , $\xi = 1.9$ ; (4) $\eta = 1.8$ , $\xi = 1.0$ ; (5) $\eta = 2.2$ , $\xi = 2.8$ .

(ii) For $\eta =1.2$ , the bubble ring also maintained a relatively uniform thickness, with slightly increased curvature at the four corners. The left- and right-hand sides collapsed earlier than the top and bottom sides. Assuming the shock wave propagates at the speed of sound (1500 m s⁻¹) towards the top and bottom residual bubbles, over a distance of approximately 410 $\unicode {x03BC}$ m, the estimated time is 0.27 $\unicode {x03BC}$ s. The time interval between Peak 4 and Peak 5 is approximately 0.4 $\unicode {x03BC}$ s. This indicates that the shock waves from the left- and right-hand sides enhanced the collapse of the top and bottom bubbles, resulting in a maximum pressure of 192 MPa in the flow field. However, the collapse occurred at a distance of approximately 75 $\unicode {x03BC}$ m from the wall (figure 23), leading to a wall pressure peak of only 28 MPa, significantly lower than the flow field pressure peak.

(iii) For $\eta =1.8$ and $\xi =1.9$ , the curvature on the left- and right-hand sides of the bubble increased significantly during the secondary oscillation. The side view showed that the central thickness of the bubble ring was approximately 62 % of the thickness at the top and bottom. The thinner and narrower regions collapsed and disappeared first, producing a pressure peak (Peak 4) in the flow field, as shown in the subplot for $t^{*} = 3.47$ in figure 23 $(b)$ . The final lip-shaped bubble collapse generated a localized wall pressure peak reaching 63 MPa (Peak 5).

(iv) Notably, for $\eta =1.8$ and $\xi =1.0$ , the overall collapse pattern of the bubble was similar to the case of $\xi =1.9$ , indicating that axial asymmetry had a limited influence on the overall collapse behaviour. However, the wall pressure peak caused by $\xi =1.0$ reached 92 MPa, nearly equal to the maximum pressure in the flow field and approximately 1.46 times that for $\xi =1.9$ .

(v) For $\eta =4.0$ , the asymmetry between the top/bottom and left/right ends of the bubble became more pronounced, with the width approximately twice that of the sides and the thickness approximately 1.7 times greater. The bipolar-V collapse mode was observed again. However, the wall pressure peak was significantly reduced, reaching only 32 MPa. This is also due to the attenuation effect of water on the shockwave load.

Overall, at moderate stand-off distances, the flow field pressure induced by the secondary collapse of the bubble is higher than that of the initial collapse. As $\eta$ increases, the circumferential asymmetry of bubble collapse becomes more pronounced. Additionally, under the same aspect ratio, elliptical initial bubbles generate higher wall pressure compared with teardrop-shaped bubbles, reaching up to three times that of spherical bubbles.

Finally, we replicated an extreme case studied by Rodriguez Jr et al. (Reference Rodriguez, Beig, Barbier and Johnsen2022) (figure 2 $a$ therein), with $\gamma = 1.25$ , and the peak pressure is approximately 400 MPa (the comparison results are not shown here). Our numerical results exhibit strong similarities in bubble collapse patterns and jetting behaviours. While the timing of the peak pressures between the two models varies by only 1 %, the amplitude of the peak pressure shows a deviation of approximately 40 %. This discrepancy may arise from the lower accuracy of the present finite volume method in capturing the shockwave front compared with higher-order methods. Additionally, thermal effects were not considered in our model (Beig et al. Reference Beig, Aboulhasanzadeh and Johnsen2018; Zhou & Prosperetti Reference Zhou and Prosperetti2020; White, Beig & Johnsen Reference White, Beig and Johnsen2023), which could also contribute to the observed differences. Future work will focus on investigating the influence of thermodynamic effects and ambient temperature on cavitation erosion.