3.1. Experimental observations and comparison with simulation results

We present the experimental and simulation results of the particle motion during the growth and collapse of the cavitation bubbles with four typical cases at various stand-off distances and size ratios. Figure 4 presents high-speed snapshots of cavitation bubble interaction with PMMA particles of radius $R_{p}={250}\,\unicode{x03BC}\mathrm{ m}$ with comparable $R_{b,\textit{max}}/R_p\approx 3.0$ and radial offset $r_0/R_p\approx 0.35$ , but varying $\gamma$ . The corresponding axisymmetric simulation results are superimposed. Note that the simulated particle profiles have been aligned horizontally with the experimental snapshots, while their vertical positions are unaltered. Here, red solid lines represent the positions of the solid surfaces, red dashed lines the bubble profiles, and dash-dotted lines the particle profiles in each frame.

Figure 4(a) shows results at $\gamma =2.1$ and $R_{b,\textit{max}}/R_p=3.05$ for both the experiment and the simulation. In the high-speed snapshots, after the bubble is generated, it experiences an explosive growth with a spherical shape. A clearance is observed between the bottom of the particle and the plate before the bubble reaches its maximum radius (approximately ${80.9}\,\unicode{x03BC}\mathrm{ s}$ , see the enlarged view under frame (iv)). This behaviour is consistent with our observations in Ren et al. (Reference Ren, Zuo, Wu and Liu2022). During bubble collapse (see frames (v)–(viii)), the particle exhibits a continuous upward motion. At the late stage of bubble collapse ( ${158}\,\unicode{x03BC}\mathrm{ s}$ ), the bubble is no longer spherical, but it develops a micro-jet towards the boundary. Then the bubble evolves into a vortex ring, and its interactions with the particle have been reported in previous studies (Wu et al. Reference Wu, Li, Zuo and Liu2021; Huang et al. Reference Huang, Ren, Liu and Hu2024). Therefore, we only focus on the dynamics before the first collapse of the cavitation bubble. In the simulation results, figure 4(a), the profiles of the bubble and the particle fit well with the experiments, demonstrating the accuracy of our numerical model.

Figure 4.

High-speed snapshots (grey-level images) and simulation results (red dashed and dash-dotted lines) of interactions between a cavitation bubble and a PMMA particle ( $R_{p}={250}\,\unicode{x03BC}\mathrm{ m}$ ). The experimental conditions for each case are: (a) $\gamma =2.10\pm 0.01$ , $R_{{b,{max}}}/R_{p}=3.05\pm 0.01$ , $r_0/R_p=0.35$ , and (b) $\gamma =2.52\pm 0.08$ , $R_{b,\textit{max}}/R_{p}=3.06\pm 0.01$ , $r_0/R_p=0.37$ . The simulation is axisymmetric with parameters (a) $R_{b,0}={35}\,\unicode{x03BC}\mathrm{ m},\ p_{0}={970}\,\mathrm{MPa}$ , $\gamma =2.10$ , $R_{{b,{max}}}/R_{p}=3.05$ , and (b) $R_{b,0}={35}\,\unicode{x03BC}\mathrm{ m},\ p_{0}={884}\,\mathrm{MPa}$ , $\gamma =2.50$ , $R_{b,\textit{max}}/R_{p}=3.06$ . The horizontal positions of the particle profiles in simulation have been aligned with the experimental snapshots for direct comparison, while the vertical dynamics are preserved. Time zero corresponds to the instant of laser plasma generation. In (a), the frame times for the snapshots are: (i) ${9.52}\,\unicode{x03BC}\mathrm{ s}$ , (ii) ${28.6}\,\unicode{x03BC}\mathrm{ s}$ , (iii) ${57.2}\,\unicode{x03BC}\mathrm{ s}$ , (iv) ${80.9}\,\unicode{x03BC}\mathrm{ s}$ , (v) ${104}\,\unicode{x03BC}\mathrm{ s}$ , (vi) ${124}\,\unicode{x03BC}\mathrm{ s}$ , (vii) ${152}\,\unicode{x03BC}\mathrm{ s}$ , and (viii) ${158}\,\unicode{x03BC}\mathrm{ s}$ . The frame times for the simulation results are: (i) ${7.20}\,\unicode{x03BC}\mathrm{ s}$ , (ii) ${26.4}\,\unicode{x03BC}\mathrm{ s}$ , (iii) ${50.4}\,\unicode{x03BC}\mathrm{ s}$ , (iv) ${79.2}\,\unicode{x03BC}\mathrm{ s}$ , (v) ${103}\,\unicode{x03BC}\mathrm{ s}$ , (vi) ${122}\,\unicode{x03BC}\mathrm{ s}$ , (vii) ${151}\,\unicode{x03BC}\mathrm{ s}$ , and (viii) ${155}\,\unicode{x03BC}\mathrm{ s}$ . In (b), the frame times for the snapshots are: (i) ${4.76}\,\unicode{x03BC}\mathrm{ s}$ , (ii) ${28.6}\,\unicode{x03BC}\mathrm{ s}$ , (iii) ${52.4}\,\unicode{x03BC}\mathrm{ s}$ , (iv) ${76.2}\,\unicode{x03BC}\mathrm{ s}$ , (v) ${100}\,\unicode{x03BC}\mathrm{ s}$ , (vi) ${119}\,\unicode{x03BC}\mathrm{ s}$ , (vii) ${138}\,\unicode{x03BC}\mathrm{ s}$ , and (viii) ${152}\,\unicode{x03BC}\mathrm{ s}$ . The frame times for the simulation results are: (i) ${4.20}\,\unicode{x03BC}\mathrm{ s}$ , (ii) ${27.9}\,\unicode{x03BC}\mathrm{ s}$ , (iii) ${51.9}\,\unicode{x03BC}\mathrm{ s}$ , (iv) ${75.9}\,\unicode{x03BC}\mathrm{ s}$ , (v) ${99.9}\,\unicode{x03BC}\mathrm{ s}$ , (vi) ${119}\,\unicode{x03BC}\mathrm{ s}$ , (vii) ${138}\,\unicode{x03BC}\mathrm{ s}$ , and (viii) ${150}\,\unicode{x03BC}\mathrm{ s}$ . For simulation results, the red solid lines represent the positions of the solid surfaces, the red dashed lines represent the bubble profiles, and the red dash-dotted lines represent the particle profiles in each frame. The regions focused on the bottom of the particle are enlarged and shown under each frame. See supplementary movie 1 for the experimental recordings at https://doi.org/10.1017/jfm.2025.10998.

When we maintain the bubble-to-particle size ratio ( $R_{b,\textit{max}}/R_p=3.06$ ), and increase the dimensionless stand-off distance $\gamma$ to approximately 2.5 (see the snapshots in figure 4 b), the bubble still shows a spherical shape during its explosive growth. The particle exhibits a finite but smaller vertical displacement before the bubble reaches its maximum radius (approximately ${76.2}\,\unicode{x03BC}\mathrm{ s}$ , see the enlarged view under frame (iv)). During bubble collapse ( ${100}\,\unicode{x03BC}\mathrm{ s}$ ) and the late stage ( ${152}\,\unicode{x03BC}\mathrm{ s}$ ), the particle upward displacement is reasonably smaller compared to figure 4(a), indicating that the bubble–particle interactions become weaker as $\gamma$ becomes larger. Still, in this case, the profiles of the bubble and particle in the experimental and simulation results fit well, respectively.

Figures 5 and 6 present cavitation bubble interactions with PMMA particles of radius $R_{p}={50}\,\unicode{x03BC}\mathrm{ m}$ under conditions with a radius ratio $R_{b,\textit{max}}/R_p\gt 10$ . In the high-speed snapshots, the cavitation bubble is partially shown in the upper portion of the images, with the position of the centroid and the radius obtained with the circle-fitting method as shown in figure 2. Using the conversion functions in Appendix A, the corresponding full-boundary parameters for the experimental cases are determined.

Figure 5.

High-speed snapshots (grey-level images) and simulation results (red solid and dash-dotted lines) of interactions between a cavitation bubble and a PMMA particle ( $R_{p}={50}\,\unicode{x03BC}\mathrm{ m}$ ). The experimental condition is $\gamma =1.3\pm 0.2$ , $R_{{b,{max}}}/R_{p}=15\pm 1$ and $r_0/R_p=2.3$ , where equivalent full-boundary parameters are used. The simulation condition is $R_{b,0}={25}\,\unicode{x03BC}\mathrm{ m},\ p_{0}={3600}\ \mathrm{MPa}$ , $\gamma =1.3$ , $R_{{b,{max}}}/R_{p}=16$ and $r_0/R_p=2.4$ . The horizontal positions of the particle profiles in simulation have been aligned with the experimental snapshots for direct comparison, while the vertical dynamics are preserved. Time zero corresponds to the instant of laser plasma generation. The first frame (no. 1) for the snapshot is ${10}\,\unicode{x03BC}\mathrm{ s}$ , and for the simulation it is ${8}\,\unicode{x03BC}\mathrm{ s}$ . The time interval between frames is ${10}\,\unicode{x03BC}\mathrm{ s}$ , except for frame 16, where the time is ${160}\,\unicode{x03BC}\mathrm{ s}$ for the snapshot and ${159}\,\unicode{x03BC}\mathrm{ s}$ for the simulation. For simulation results, the red solid lines indicate the lower interface of the cavitation bubble, and the red dash-dotted lines represent the particle profiles in each frame. Beneath each snapshot, the difference image (reference frame before time zero minus current frame) highlights particle displacement (bright areas). The white solid line marks the particle’s original position, and the red dotted line shows its profile at each time. In frame 2, the particle’s initial vertical displacement becomes visible, as indicated by the arrow. Scale bars: ${100}\,\unicode{x03BC}\mathrm{ m}$ (main images), ${50}\,\unicode{x03BC}\mathrm{ m}$ (difference images). See supplementary movie 2 for the experimental recording.

Figure 5 presents experimental observations for $\gamma =1.3\pm 0.2$ , $R_{{b,{max}}}/R_{p}=15\pm 1$ and $r_0/R_p=2.3$ . From time zero to frame 8 (approximately ${80}\,\unicode{x03BC}\mathrm{ s}$ ), the bubble’s lower interface descends towards the solid boundary, indicating bubble growth. From frame 8 to frame 15, the lowest position of the interface remains, but the radius of curvature decreases, indicating bubble collapse. At frame 16 ( ${160}\,\unicode{x03BC}\mathrm{ s}$ ), the bubble reaches the late collapse stage, exhibiting a dimple on top that develops into a micro-jet towards the surface. Particle displacement is visualised in the difference images shown beneath each snapshot, generated by subtracting each frame from a reference image before time zero. Bright areas in these images highlight the particle’s displacement. Whilst minimal in frame 1, an upward and radially outward displacement is clearly visible by frame 2, as marked by the arrow. The particle continues to rise throughout the bubble’s growth and collapse, reaching its maximum height (approximately 10 pixels, or ${24}\,\unicode{x03BC}\mathrm{ m}$ ) at frame 16. Quantitative measurements are provided in figure 7. Since the particle has an original radial offset ( $r_0\approx {120}\,\unicode{x03BC}\mathrm{ m}$ ) with respect to the bubble’s axis of symmetry, the particle moves outwards (frame 2–frame 7) and then inwards (frame 8–frame 16), with a radial motion that appears to synchronise with the bubble’s expansion and contraction. The particle’s radial displacement is not the focus of this work. The original radial offset is intentionally minimised (see § 2.1), resulting in a negligible net radial displacement by the end of bubble collapse (frame 16). The simulated bubble profiles for $\gamma =1.3$ and $R_{b,\textit{max}}/R_p=16$ demonstrate close agreement with the experimental observations throughout both the growth and collapse phases. The simulated particle profile ( $r_0/R_p=2.4$ ) also matches the experiment accurately from time zero until frame 15. A discrepancy arises in the late stage of bubble collapse (frame 16), where the simulation predicts a larger vertical displacement of the particle than observed experimentally. While an out-of-plane offset could potentially attenuate the vertical flow, Appendix B indicates that this factor plays only a limited role in influencing the particle vertical motion. The discrepancy is more likely attributed to the simplified force balance model, which may break down in the complex three-dimensional flow at the end of bubble collapse, potentially underestimating dissipative effects (Zeng et al. Reference Zeng, Gonzalez-Avila, Dijkink, Koukouvinis, Gavaises and Ohl2018a , Reference Zeng, An and Ohl2022). This phenomenon requires further investigation. Nevertheless, the simulation successfully captures the primary features of the particle dynamics observed experimentally.

Figure 6.

High-speed snapshots (grey-level images) and simulation results (red solid and dash-dotted lines) of interactions between a cavitation bubble and a PMMA particle ( $R_{p}={50}\,\unicode{x03BC}\mathrm{ m}$ ). The experimental condition is $\gamma =1.8\pm 0.1$ , $R_{b,\textit{max}}/R_{p}=11\pm 1$ and $r_0/R_p=1.0$ , where equivalent full-boundary parameters are used. The simulation condition is $R_{b,0}={25}\,\unicode{x03BC}\mathrm{ m},\ p_{0}={844}\,\mathrm{MPa}$ , $\gamma =1.8$ , $R_{b,\textit{max}}/R_{p}=11$ and $r_0/R_p=1.0$ . The horizontal positions of the particle profiles in simulation have been aligned with the experimental snapshots for direct comparison, while the vertical dynamics are preserved. Time zero corresponds to the instant of laser plasma generation. The first frame (no. 1) for the snapshot is ${10}\,\unicode{x03BC}\mathrm{ s}$ , and for the simulation it is ${8.2}\,\unicode{x03BC}\mathrm{ s}$ . The time interval between frames is ${10}\,\unicode{x03BC}\mathrm{ s}$ , except for frame 10 where the time is ${100}\,\unicode{x03BC}\mathrm{ s}$ for the snapshot and ${96.4}\,\unicode{x03BC}\mathrm{ s}$ for the simulation. For simulation results, the red solid lines indicate the lower interface of the cavitation bubble, and the red dash-dotted lines represent the particle profiles in each frame. Beneath each snapshot, the difference image (reference frame before time zero minus current frame) highlights particle displacement (bright areas). The white solid line marks the particle’s original position, and the red dotted line shows its profile at each time. In frame 4, the particle’s initial vertical displacement becomes visible, as indicated by the arrow. Scale bars: ${100}\,\unicode{x03BC}\mathrm{ m}$ (main images), ${50}\,\unicode{x03BC}\mathrm{ m}$ (difference images). See supplementary movie 2 for the experimental recording.

Figure 6 shows observations for $\gamma =1.8$ , $R_{b,\textit{max}}/R_p=11$ and $r_0/R_p=1.0$ . Similar to figure 5, the bubble lower interface descends towards the solid boundary from time zero to frame 5 (approximately ${50}\,\unicode{x03BC}\mathrm{ s}$ ), indicating bubble growth. The bubble lower interface ascends away from the solid boundary from frame 6 to frame 10, showing an increase in the radius of curvature, indicating bubble collapse. The particle’s vertical displacement becomes visualised as the thin, bright area to the upper left of its original position in the difference image at frame 4. This region is indicated by the arrow below the snapshot in figure 6. From frame 4 to frame 8, the particle rises slightly, as indicated by the expanding bright area in the difference image. Subsequently, at the end of bubble collapse (frames 9–10), this bright area diminishes slightly. Overall, the particle displacement in this case is less than 5 pixels ( $\lesssim{12}\,\unicode{x03BC}\mathrm{ m}$ ), which is notably smaller than that observed in figure 5. Since the particle has a small initial radial offset ( $r_0\approx \,{50}\,\unicode{x03BC}\mathrm{ m}$ ) compared to the maximum bubble radius, the particle shows a minimal displacement in the radial direction until the end of bubble collapse. The simulation adopts $\gamma =1.8$ , $R_{b,\textit{max}}/R_p=11$ and $r_0/R_p=1.0$ . The simulated bubble profiles show good agreement with the experiments during growth and collapse. The simulation also captures the suppressed particle vertical displacement as observed.

Figure 7.

Evolution of the cavitation bubble radius and the particle displacement at different $\gamma$ and $R_{b,\textit{max}}/R_p$ . (a) Results of PMMA particles with $R_p={250}\,\unicode{x03BC}\mathrm{ m}$ . Case A (expt): $\gamma =2.10\pm 0.01$ , $R_{{b,{max}}}/R_{p}=3.05\pm 0.01$ , $r_0/R_p=0.35$ . Case B (expt): $\gamma =2.52\pm 0.08$ , $R_{b,\textit{max}}/R_{p}=3.06\pm 0.01$ , $r_0/R_p=0.37$ . Case A (sim): $\gamma =2.1$ , $R_{{b,{max}}}/R_{p}=3.05$ , $r_0/R_p=0$ . Case B (sim): $\gamma =2.5$ , $R_{b,\textit{max}}/R_{p}=3.06$ , $r_0/R_p=0$ . (b) Results of PMMA particles with $R_p={50}\,\unicode{x03BC}\mathrm{ m}$ . For experiments, equivalent full-boundary parameters are used. Case C (expt): $\gamma =1.3\pm 0.2$ , $R_{{b,{max}}}/R_{p}=15\pm 1$ , $r_0/R_p=2.3$ . Case D (expt): $\gamma =1.8\pm 0.1$ , $R_{b,\textit{max}}/R_{p}=11\pm 1$ , $r_0/R_p=1.0$ . Case C (sim): $\gamma =1.3$ , $R_{{b,{max}}}/R_{p}=16$ , $r_0/R_p=2.4$ . Case D (sim): $\gamma =1.8$ , $R_{b,\textit{max}}/R_{p}=11$ , $r_0/R_p=1.0$ . In (a) and (b), panel (i) refers to the evolution of the dimensionless bubble radius $R_{b}/R_{b,\textit{max}}$ with dimensionless time $t/T_c$ with $T_c=R_{b,\textit{max}}\,\sqrt {\rho _{{l}}/\Delta p}$ , and panel (ii) refers to the evolution of the dimensionless vertical particle displacement $\Delta z/R_{p}$ . The grey areas denote the period of bubble explosive growth.

Figure 7 shows the quantitative comparisons of the experimental and simulation results (cases A–D) in terms of the evolution of the cavitation bubble radius and the particle vertical displacement. The time is non-dimensionalised with the scale of bubble lifetime $T_c=R_{b,\textit{max}}\,\sqrt {\rho _{{l}}/\Delta p}$ from the Rayleigh–Plesset equation (Rayleigh Reference Rayleigh1917; Plesset Reference Plesset1949), with $\Delta p=p_{\infty }-p_v$ . Here, $p_{\infty }={1.01\times 10^5}\,\mathrm{Pa}$ is the ambient static pressure, and $p_v={2.3\times 10^3}\,\mathrm{Pa}$ is the vapour pressure at room temperature (Brennen Reference Brennen2014).

For cases A and B, the bubble radii of the simulation results $R_b/R_{b,\textit{max}}$ agree well with the experiments, as shown in figure 7(a)(i). In figure 7(a)(ii), the particle’s dimensionless displacement $\Delta z/R_p$ fits well between the experiments and the simulations. In both cases, the particle detaches from the substrate during the period of the bubble’s explosive growth (before $t/T_c=0.5$ , see the grey area). From the decelerating growth of the bubble radius to the late stage before $t/T_c\approx 1.9$ , the particle experiences an accelerating upward motion, then it decelerates and reaches the maximum height, and finally descends due to the jetting of the cavitation bubble (after collapse).

For cases C and D, the bubble radii of the simulation results $R_b/R_{b,\textit{max}}$ show good agreement with the experimental values equivalent to the full-boundary fitting method; see figure 7(b)(i). The simulation results in figure 7(b)(ii) also align with experiments for case C, where the particle detaches from the substrate before $t/T_c = 0.5$ , and exhibits accelerating upward motion during the bubble’s first cycle with $\Delta z/R_p \lesssim 0.5$ . In contrast, for case D, the particle vertical displacement is strongly suppressed in the experiment, which is also captured by simulation ( $\Delta z/R_p \lesssim 0.1$ ).

3.2. Evolution of the hydrodynamic forces on the particle

To understand the mechanisms behind the particle jump behaviours driven by cavitation bubbles, we show the evolution of the simulated forces on the particle, including the component forces from (2.4)–(2.8) and the resultant hydrodynamic forces, for cases A–D in figure 8. The forces are smoothed using a low-pass Butterworth filter, and non-dimensionalised by the characteristic force scale $F_0=m_pR_p/T_c^2$ , where $m_p$ denotes the mass of the particle. In all cases, the resultant force $F_h$ (thick black dashed line) converts its sign from negative (compressive force) to positive (lift force) during the bubble’s explosive phase (within the grey area). Since the net buoyancy force $F_b$ includes the weight of the particle, the instant for particle detachment is found at $F_h=0$ , as denoted by the arrows. Then the hydrodynamic force serves as a lift force until it decreases by the end of the bubble collapse, which is consistent with the cases for $R_{b,\textit{max}}/R_p\sim 1$ as reported by Ren et al. (Reference Ren, Zuo, Wu and Liu2022).

Figure 8.

Evolution of the dimensionless hydrodynamic forces (smoothed) exerted on the PMMA particle in the vertical direction at different $\gamma$ and $R_{b,\textit{max}}/R_p$ with radii (a,b) $R_p={250}\,\unicode{x03BC}\mathrm{ m}$ and (c,d) $R_p={50}\,\unicode{x03BC}\mathrm{ m}$ . (a) Case A: $\gamma =2.1$ , $R_{{b,{max}}}/R_{p}=3.05$ , $r_0/R_p=0$ . (b) Case B: $\gamma =2.5$ , $R_{b,\textit{max}}/R_{p}=3.06$ , $r_0/R_p=0$ . (c) Case C: $\gamma =1.3$ , $R_{{b,{max}}}/R_{p}=16$ , $r_0/R_p=2.4$ . (d) Case D: $\gamma =1.8$ , $R_{b,\textit{max}}/R_{p}=11$ , $r_0/R_p=1.0$ . The characteristic scale for the force is $F_0=m_pR_p/T_c^2$ . The black dashed lines in (a–d) denote the resultant hydrodynamic force $F_h$ . The arrows denote the instants at $F_h=0$ during bubble explosive growth (period denoted by the grey areas). The insets in (a–c) are the enlarged views of $F_{\textit{drag}}$ and $F_{b}$ .

Dominant forces are identified by comparing the magnitudes of the component forces. For case A, the resultant force $F_h$ is of the same magnitude as the pressure gradient force $F_{\textit{pres}}$ (red solid line) through the first cycle of bubble oscillation; see figure 8(a). The added mass force $F_{add}$ (yellow dash-dotted line) shows the same sign conversion as $F_{\textit{pres}}$ but with an order of magnitude smaller than $F_{\textit{pres}}$ . The lubrication force $F_{\textit{lub}}$ (green dashed line) serves as a negative force after the particle detaches from the substrate (marked by the arrow), and converts its sign to positive at $t/T_c\approx 0.8$ . Following the expression of $F_{\textit{lub}}$ in (2.8), this conversion corresponds to the sign of the term in the square bracket, with increasing gap thickness $\epsilon$ and increasing Reynolds number $\textit{Re}_{v}$ . During bubble growth, the magnitudes of $|F_{\textit{lub}}|$ and $|F_{add}|$ are similar, and their contributions to the resultant hydrodynamic force balance each other. The viscous drag force $F_{\textit{drag}}$ (light blue dash-double-dotted line) remains negative, and its amplitude is only ${\sim}20\,\%$ of $|F_{\textit{lub}}|$ . The net buoyancy force $F_b$ (dark blue dotted line) remains a negative constant, and its absolute value is three orders of magnitude smaller than $|F_{\textit{drag}}|$ . For case B (figure 8 b), the evolution of the component forces is similar to case A, although the magnitudes of the forces are reasonably smaller, about one-half of the magnitudes of the forces for case A.

For case C (figure 8 c), the resultant force $F_h$ is comparable in magnitude to the pressure gradient force $F_{\textit{pres}}$ during the bubble’s explosive growth. After the particle detaches (marked by the arrow), the lubrication force $F_{\textit{lub}}$ serves as a resisting force with an amplitude similar to $F_{\textit{pres}}$ . This allows the added mass force $F_{add}$ to become dominant until $t/T_c\approx 1.5$ , after which $F_{\textit{pres}}$ regains dominance because $F_{add}$ and $F_{\textit{lub}}$ eliminate each other. The viscous drag force $|F_{\textit{drag}}|$ is one order of magnitude smaller than $|F_{\textit{lub}}|$ , while the net buoyancy force $|F_b|$ is three orders of magnitude smaller than $|F_{\textit{drag}}|$ . The significant difference in $F_{\textit{pres}}$ at the end of bubble collapse compared to cases A and B is caused by the distinct flow field and the boundary layer evolution, which will be discussed in § 3.3.

For case D, similar to case C, the added mass force $F_{add}$ dominates during bubble growth, while $F_{\textit{pres}}$ and $F_{\textit{lub}}$ balance each other, as shown in figure 8(d). During bubble collapse, the lubrication force decreases and balances $F_{add}$ , leading to the dominance of $F_{\textit{pres}}$ . Still, the viscous drag force $|F_{\textit{drag}}|$ is one order of magnitude smaller than $|F_{\textit{lub}}|$ , while the net buoyancy force $|F_b|$ is three orders of magnitude smaller than $|F_{\textit{drag}}|$ .

In summary, the dominant force governing particle jump shifts with the bubble–particle size ratio. At small ratios ( $R_{b,\textit{max}}/R_p\sim 3$ ), the pressure gradient force $F_{\textit{pres}}$ dominates, while the added mass force $F_{add}$ and the lubrication force $F_{\textit{lub}}$ balance each other. At large ratios ( $R_{b,\textit{max}}/R_p\gtrsim 10$ ), an enhanced lubrication force $F_{\textit{lub}}$ counteracts $F_{\textit{pres}}$ , allowing the added mass force $F_{add}$ to dominate during bubble growth. This dynamic competition among $F_{\textit{pres}}$ , $F_{add}$ and $F_{\textit{lub}}$ is the distinguishing feature of large size ratio cases. In all cases, the effects of $F_{\textit{drag}}$ and $F_{b}$ are negligible.

3.3. Evolution of the bubble-induced boundary layer flows

To further understand how the evolution of the hydrodynamic forces on the particle is determined by the bubble-induced boundary layer flow, we look into the evolution of the flow velocity. Since we use the one-way fluid–particle coupling simulation, we show the simulation results of the undisturbed velocity distributions at $r=r_0$ between the bubble bottom and the substrate for different $\gamma$ and $R_{b,\textit{max}}/R_p$ (cases A–D); see figure 9. Panels (i) show $U_z$ versus $z/R_p$ , where the different line colours denote the dimensionless time $t/T_c$ for each case.

Figure 9.

Simulation results of the undisturbed velocity distribution at $r=r_0$ with different $\gamma$ and $R_{b,\textit{max}}/R_p$ for PMMA particles with radii (a,b) $R_p={250}\,\unicode{x03BC}\mathrm{ m}$ and (c,d) $R_p={50}\,\unicode{x03BC}\mathrm{ m}$ . (a) Case A: $\gamma =2.1$ , $R_{{b,{max}}}/R_{p}=3.05$ , $r_0/R_p=0$ . (b) Case B: $\gamma =2.5$ , $R_{b,\textit{max}}/R_{p}=3.06$ , $r_0/R_p=0$ . (c) Case C: $\gamma =1.3$ , $R_{{b,{max}}}/R_{p}=16$ , $r_0/R_p=2.4$ . (d) Case D: $\gamma =1.8$ , $R_{b,\textit{max}}/R_{p}=11$ , $r_0/R_p=1.0$ . In (a–d), panel (i) refers to the profile of the wall-normal velocity $U_z$ , and panel (ii) refers to the vertical velocity gradient $\partial \,U_z/\partial \,z$ . The colour bar on the right denotes the dimensionless time $t/T_c$ for each case. The grey circles in panel (ii) indicate the boundary layer thickness $\delta /R_p$ at the different times. In (a), the times $t/T_c$ are 0.037, 0.28, 0.52, 0.76, 1.00, 1.24, 1.48 and 1.72; in (b), the times $t/T_c$ are 0.025, 0.27, 0.51, 0.75, 0.99, 1.23, 1.47 and 1.71; in (c), the times $t/T_c$ are 0.013, 0.26, 0.51, 0.75, 1.00, 1.25, 1.49 and 1.74; in (d), the times $t/T_c$ are 0.045, 0.29, 0.53, 0.77, 1.02, 1.26, 1.50 and 1.74.

Figure 9(a)(i) shows the typical evolution of the undisturbed velocity profiles between the bubble bottom and the solid wall for case A ( $\gamma =2.1$ ). At $t/T_c=0.037$ (dark blue line), the negative velocity indicates a front stagnation flow towards the substrate at $\approx {1}\ \mathrm{m}\ \text{s}^{-1}$ at level $z/R_p=2$ . As the bubble grows explosively, the fluid velocity (at the same level) increases ( $t/T_c=0.28$ ), indicating the accelerating growth of the bubble. Then the bubble growth decelerates ( $t/T_c=0.28$ –1.00) until the bubble grows to its maximum size when the velocity becomes approximately zero. During the bubble’s collapse, the velocity becomes positive and increases again, indicating a rear stagnation flow ( $t/T_c=1.00$ –1.72). Compared with case A, the velocity profiles are similar for case B ( $\gamma =2.5$ , see figure 9 bi) and case D ( $\gamma =1.8$ , see figure 9 di), indicating that the evolution of the flow field between the bubble bottom and the rigid wall is similar. For case C ( $\gamma =1.3$ , see figure 9 ci), which is a particular representative for a bubble with flat bottom during growth and collapse, the bubble’s explosive growth phase can be obviously observed before $t/T_c=0.26$ . Yet before $t/T_c=1.49$ (orange line), the velocity remains negative, indicating that the front stagnation flow lasts until the late stage of the bubble collapse.

Then we analyse the evolution of the boundary layer thickness for the different cases based on the velocity profiles. From the solution to the axisymmetric stagnation (Homann) flow (Homann Reference Homann1936), we know that the potential flow solution out of the boundary layer is $U_z=-2Bz$ and $U_r=Br$ , where $B$ is the strain rate of the flow, $z$ is the vertical distance from the wall, and $r$ is the radial distance from the axis. Within the boundary layer, the component velocities read $U_z=-2B\,f(z)$ and $U_r=Br\,{\text{d}f}/{\text{d}z}$ , where $f(z)$ is a function determining the velocity profile. Outside the boundary layer, $z\gg \,\delta$ , the function $f(z)$ fulfils $f(z)\,\to \,z$ and ${\text{d}f}/{\text{d}z}\,\to \,1$ . The boundary layer thickness $\delta$ is defined as the normal distance from the wall where the radial velocity component is equal to 99 % of the corresponding free-stream value, i.e. $({\text{d}f}/{\text{d}z})|_{(z=\delta )}=0.99$ . We prove in § 3.4 that with $r,z\ll L,R_{b,\textit{max}}$ the analytical solution to the bubble-induced stagnation flow has the same form as the Homann solution. Based on this, we identify the boundary layer thickness $\delta$ by finding the position where the vertical velocity gradient $\partial \,U_z/\partial \,z$ reaches 99 % of a constant value. As shown in panels (ii) of figure 9, the boundary layer thickness $\delta /R_p$ for each case is denoted by the grey circles at the selected times.

The evolution of the unsteady boundary layer thickness with time ( $\delta /R_p$ versus $t/T_c$ ) is plotted in figure 10, which is superimposed with the central position of the particle, $z_p/R_p=1+\Delta z/R_p$ , for comparisons. For case A (figure 10 a), the boundary layer thickness $\delta /R_p$ increases slowly from zero to approximately 0.25 at $t/T_c\approx 1.0$ . Then the boundary layer grows fast to a peak value $\delta /R_p\approx 0.9$ (close to the particle centre) at $t/T_c\approx 1.1$ as the bubble reaches its maximum size (see figure 7 ai). As the bubble collapses, the boundary layer thickness decreases and approaches a constant $\delta /R_p\approx 0.5$ until $t/T_c=1.95$ . The boundary layer thickness remains less than 50 % of the particle radius during most time of the bubble’s oscillation, indicating that the particle is mainly subjected to the force induced by the inviscid flow field, corresponding to the leading role of the pressure gradient force $F_{\textit{pres}}$ .

Figure 10.

Simulation results of evolution of the dimensionless boundary layer thickness $\delta /R_p$ at $r=r_0$ with different $\gamma$ and $R_{b,\textit{max}}/R_p$ for PMMA particles with radii (a,b) $R_p={250}\,\unicode{x03BC}\mathrm{ m}$ and (c,d) $R_p={50}\,\unicode{x03BC}\mathrm{ m}$ . (a) Case A: $\gamma =2.1$ , $R_{{b,{max}}}/R_{p}=3.05$ , $r_0/R_p=0$ . (b) Case B: $\gamma =2.5$ , $R_{b,\textit{max}}/R_{p}=3.06$ , $r_0/R_p=0$ . (c) Case C: $\gamma =1.3$ , $R_{{b,{max}}}/R_{p}=16$ , $r_0/R_p=2.4$ . (d) Case D: $\gamma =1.8$ , $R_{b,\textit{max}}/R_{p}=11$ , $r_0/R_p=1.0$ . The boundary layer thickness $\delta /R_p$ is compared with the position of the particle centre $z_p/R_p$ in each case. For the experimental data, bubble parameters obtained from the partial-boundary fitting method are converted into their full-boundary equivalents. The grey areas in (a–d) denote the period of bubble explosive growth.

For case B (figure 10 b), the evolution of the boundary layer thickness during bubble growth is similar to case A, with a smaller peak value ( $\delta /R_p\approx 0.4$ ) at the bubble’s maximum size ( $t/T_c\approx 1.0$ ). During bubble collapse, the boundary layer thickness remains a constant $\delta /R_p\approx 0.2$ , less than that of case A. The evolution of $\delta /R_p$ also indicates that the particle’s motion is dominated by the inviscid flow field.

For case C (figure 10 c), the boundary layer thickness increases from zero to a peak value $\delta /R_p\approx 1.8$ at $t/T_c\approx 1.8$ , exceeding the position of the particle centre. Then it drops to $\delta /R_p\approx 1$ by the end of bubble collapse. The peak of the boundary layer thickness is delayed because the front stagnation flow is largely retarded due to the lasting downward motion of the bubble bottom, as already seen in figure 5. During this time period, the particle motion is largely influenced by the boundary layer flow, leading to a large lubrication force $F_{\textit{lub}}$ . After $t/T_c\approx 1.8$ , the particle centre is suddenly exposed to the outer flow again, causing an abrupt increase in $F_{\textit{pres}}$ as observed in figure 8(c).

For case D (figure 10 d), the boundary layer thickness $\delta /R_p$ grows fast from zero to a peak value of approximately 1.0 at $t/T_c\approx 1.1$ , then it maintains at the level 0.9 until the end of bubble collapse, almost immersing the particle centre. During bubble growth, compared with cases A and B, the front stagnation flow in case D is retarded near the particle centre, allowing the fast increase of the lubrication force $F_{\textit{lub}}$ shown in figure 8(d), which suppresses the particle vertical motion.

3.4. Scaling analysis of the particle displacement at the end of bubble collapse

To elucidate the scaling law for the particle’s displacement at the end of bubble collapse, we perform a model analysis of the dynamics of the cavitation bubble–particle interaction. Figure 11 presents a schematic of the model for the flow fields induced by a cavitation bubble without particles. In figure 11(a), we assume an axisymmetric stagnation point flow between the bubble’s bottom interface and the rigid wall. This flow includes a viscous boundary layer with a time-dependent thickness $\delta$ , reflecting the unsteady nature of the cavitation bubble dynamics. The velocity components within the boundary layer ( $z\lt \delta$ ) are denoted by $(u_r,u_z)$ , while the velocity components outside the boundary layer are $(U_r,U_z)$ . For simplification, the velocity field outside the boundary layer is assumed to be approximately equivalent to the solution obtained from the potential flow model, as shown in figure 11(b). Following the potential flow theory framework established in Best & Blake (Reference Best and Blake1994), the flow field induced by the cavitation bubble is superimposed with its mirror image symmetric about a horizontal rigid wall. This approach satisfies the boundary condition of zero normal velocity at the rigid wall.

Figure 11.

Schematic of the theoretical model for cavitation bubble-induced flow near a rigid boundary. (a) The flow field near the stagnation point under a growing cavitation bubble. (b) The flow field outside the boundary layer established by a spherical cavitation bubble near a solid wall with the method of images.

In cylindrical coordinates, the flow field is symmetric about the $z$ -axis passing through the centre of the cavitation bubble. Consistent with prior works (Best & Blake Reference Best and Blake1994; Wu et al. Reference Wu, Li, Zuo and Liu2021; Ren et al. Reference Ren, Zuo, Wu and Liu2022), the cavitation bubble is modelled as a sphere with a time-dependent radius $R_{b}(t)$ . A point source of strength $Q(t)=4\pi R_{b}^2\dot {R}_{{b}}$ is placed at the bubble’s centre, where $\dot {R}_{b}$ represents the velocity at the bubble interface. A second, identical point source is positioned at the location of the mirror image cavitation bubble, symmetric about the plane $z=0$ . The flow velocity $\boldsymbol{U}=(U_r,U_z)$ at any point $(r,z)$ with an arbitrary rotation angle about the $z$ -axis is the superposition of the velocities induced by both point sources:

(3.1) \begin{align} U_r=R_{b}^2\dot {R}_{{b}}\,\left \{\frac {r}{\big [r^2+(L-z)^2\big ]^{{3}/{2}}}+\frac {r}{\big [r^2+(L+z)^2\big ]^{{3}/{2}}}\right \} \!, \\[-28pt] \nonumber \end{align}

(3.2) \begin{align} U_z=R_{b}^2\dot {R}_{{b}}\,\left \{\frac {-(L-z)}{\big [r^2+(L-z)^2\big ]^{{3}/{2}}}+\frac {L+z}{\big [r^2+(L+z)^2\big ]^{{3}/{2}}}\right \}\!. \\[-2pt] \nonumber \end{align}

The time evolution of the cavitation bubble radius is governed by the modified Rayleigh equation near the solid wall (Best & Blake Reference Best and Blake1994):

(3.3) \begin{align} R_{b}\ddot {R}_{b}+\frac {3}{2}\dot {R}_{b}^2+\frac {R_{b}}{2L}\big (R_{b}\ddot {R}_{b}+2\dot {R}_{b}^2\big )=-\frac {\Delta p}{\rho _{{l}}}, \end{align}

where $\ddot {R}_{b}$ denotes the acceleration of the bubble interface, with the driving pressure difference $\Delta p=p_{\infty }-p_{{v}}$ , where $p_{\infty }={1.01\times 10^5}\,\mathrm{Pa}$ is the pressure in the far field, and $p_{v}={2.3\times 10^3}\,\mathrm{Pa}$ is the vapour saturation pressure. The initial conditions for (3.3) are $R_{b}(0)=R_{b,\textit{max}}$ and $\dot {R}_{b}(0)=0$ . Solving (3.3) yields $R_{b}(t)$ , which describes the collapse of a cavitation bubble from its maximum radius. To obtain the complete $R_{b}(t)$ curve encompassing both growth and collapse phases, the solution curve is reflected about the time point where $R_{b}=R_{b,\textit{max}}$ and the time zero is shifted to the instant of bubble generation (Obreschkow, Bruderer & Farhat Reference Obreschkow, Bruderer and Farhat2012; Wu et al. Reference Wu, Zuo, Stone and Liu2017).

Within the region close to the bottom of the cavitation bubble ( $r,z\ll L,R_{b,\textit{max}}$ ), the potential flow field described by (3.1) and (3.2) can be further simplified:

(3.4) \begin{align} \left \{ \begin{array}{l} \displaystyle U_r=\frac {2R_{{b}}^2\dot {R}_{{b}}}{L^3}\,r,\\ U_z=-\dfrac {4R_{{b}}^2\dot {R}_{{b}}}{L^3}\,z. \end{array}\right . \end{align}

Here, we define $B={2R_{b}^2\dot {R}_{b}}/{L^3}$ with units s $^{-1}$ . With this substitution, (3.4) becomes

(3.5) \begin{align} \left \{ \begin{array}{l} \displaystyle U_r=Br,\\ U_z=-2Bz. \end{array}\right . \end{align}

The simplified inviscid flow field solution in (3.5) resembles the outer flow field of the steady axisymmetric stagnation flow solution presented by Homann (Reference Homann1936). However, our case involves a highly unsteady flow where $B$ varies with time due to the cavitation bubble’s growth and collapse.

Following a method similar to that of Homann (Reference Homann1936), the Navier–Stokes equations governing viscous flow are established. Substituting the inviscid solution (3.5) into the Euler equations in cylindrical coordinates:

(3.6) \begin{align} \left \{\begin{array}{l} \displaystyle \frac {\partial U_r}{\partial t}+U_r\,\frac {\partial U_r}{\partial r}+U_z\,\frac {\partial U_r}{\partial z}=-\frac {1}{\rho _{l}}\frac {\partial p}{\partial r},\\ \displaystyle \frac {\partial U_z}{\partial t}+U_r\,\frac {\partial U_z}{\partial r}+U_z\,\frac {\partial U_z}{\partial z}=-\frac {1}{\rho _{l}}\frac {\partial p}{\partial z}. \end{array}\right . \end{align}

After simplification, the equations governing the inviscid pressure field are obtained:

(3.7) \begin{align} \left \{\begin{array}{l} \displaystyle \frac {\partial p}{\partial r}=-\rho _{l}\big (\dot {B}+B^2\big )r, \\[6pt] \displaystyle \dfrac {\partial p}{\partial z}=-\rho _{l}\big (-2\dot {B}+4B^2\big )z, \end{array}\right . \end{align}

with $\dot {B}={\text{d}B}/{\text{d}t}$ . Integrating (3.7) and introducing the stagnation pressure $p_{\textit{st}}$ yields

(3.8) \begin{align} p=p_{\textit{st}}-\frac {1}{2}\rho _{l}\big [\big(\dot {B}+B^2 \big)r^2+\big(-2\dot {B}+4B^2 \big)z^2\big ]. \end{align}

The viscous flow solution is assumed to be

(3.9) \begin{align} \left \{\begin{array}{l} \displaystyle u_r=Br\,\frac {\mathrm{d}f}{\mathrm{d}z},\\ u_z=-2B\,f(z). \end{array}\right . \end{align}

Here, $f(z)$ is an undetermined function of spatial coordinate $z$ . The pressure field of the viscous flow solution is also assumed to be

(3.10) \begin{align} p=p_{\textit{st}}-\frac {1}{2}\rho _{l}\big [\big(\dot {B}+B^2 \big)r^2+\big(-2\dot {B}+4B^2 \big)\,g(z)\big ], \end{align}

where $g(z)$ is an undetermined function. We then substitute the viscous flow solution from (3.9) into the continuity equation and Navier–Stokes equations in cylindrical coordinates:

(3.11) \begin{align} \left \{\begin{array}{l} \displaystyle \frac {\partial u_r}{\partial r}+\frac {u_r}{r}+\frac {\partial u_z}{\partial z}=0,\\\displaystyle \frac {\partial u_r}{\partial t}+u_r\,\frac {\partial u_r}{\partial r}+u_z\,\frac {\partial u_r}{\partial z}=-\frac {1}{\rho _{l}}\frac {\partial p}{\partial r}+\nu _{l}\left (\frac {\partial ^2 u_r}{\partial r^2}+\frac {1}{r}\,\frac {u_r}{r}-\frac {u_r}{r^2}+\frac {\partial ^2 u_r}{\partial z^2}\right )\!,\\\displaystyle \frac {\partial u_z}{\partial t}+u_r\,\frac {\partial u_z}{\partial r}+u_z\,\frac {\partial u_z}{\partial z}=-\frac {1}{\rho _{l}}\frac {\partial p}{\partial z}+\nu _{l}\left (\frac {\partial ^2 u_z}{\partial r^2}+\frac {1}{r}\,\frac {u_z}{r}+\frac {\partial ^2 u_z}{\partial z^2}\right )\!, \end{array}\right . \end{align}

where $\nu _{l}$ is the kinematic viscosity of water. The continuity equation is automatically satisfied. Below, we simplify the Navier–Stokes equation in the $r$ -direction.

In the following derivation, ${\mathrm{d} f}/{\mathrm{d} z}$ is abbreviated as $f^{\prime }$ , ${\mathrm{d}^2 f}/{\mathrm{d} z^2}$ as $f''$ , and ${\mathrm{d}^3 f}/{\mathrm{d} z^3}$ as $f'''$ . The simplified Navier–Stokes equation in the $r$ -direction becomes

(3.12) \begin{align} f'''+\frac {2B}{\nu _{l}}ff''+\frac {B}{\nu _{l}}\big (1-f'^2\big )+\frac {\dot {B}}{\nu _{l}B}\big (1-f'\big )=0. \end{align}

This equation governs the axisymmetric viscous stagnation flow induced by the cavitation bubble.

The boundary conditions for (3.12) at time $t$ are given by

(3.13) \begin{align} \left \{\begin{array}{l} \displaystyle z=0,\quad f'(0)=f(0)=0,\\ z\to \infty,\quad f'(z)\to \,1.\\ \end{array}\right . \end{align}

The no-slip boundary condition ( $f = f' = 0$ ) applies at the solid wall ( $z = 0$ ). The far-field boundary condition indicates that the velocity profile approaches the inviscid flow ( $f' = 1$ ) as the distance from the wall $z$ increases far beyond the boundary layer thickness.

The time-dependent nature of parameter $B$ in (3.12) presents a significant challenge for obtaining a numerical solution. As noted by Kolomenskiy & Moffatt (Reference Kolomenskiy and Moffatt2012), solutions to similar unsteady viscous stagnation flow equations are highly sensitive to the parameter $\dot {B}/B^2$ . Depending on the range of $\dot {B}/B^2$ , the equation may exhibit various solution behaviours, including single-valued, bifurcated, multi-valued, or even no solution at all. These complexities likely reflect instabilities within the boundary layer of the actual flow.

Therefore, a detailed analysis of solving (3.12) and the characteristics of its solutions is deferred for future work. Instead, we leverage a simplified form of (3.12) to identify key parameters and derive scaling laws for particle motion.

To link the flow field to the displacement of the particle, we leverage the momentum theorem. The particle’s velocity $v_{p,\textit{end}}$ at the end of bubble collapse $T_{\textit{end}}$ (before jet penetration) reads

(3.14) \begin{align} v_{p,\textit{end}}=\frac {\int _0^{T_{\textit{end}}}F_{z}\,\text{d}t}{m_{p}}, \end{align}

where $F_z\ (=F_h+F_N)$ represents the vertical force acting on the particle. We denote the time of particle detachment from the solid surface as $t_0$ . Since the vertical force acting on the particle remains zero ( $F_z=0$ ) before $t_0$ , the integral in (3.14) reduces to

(3.15) \begin{align} \int _0^{T_{\textit{end}}} F_z\, \text{d}t = \int _{t_0}^{T_{\textit{end}}} F_z\, \text{d}t. \end{align}

After particle detachment ( $\Delta z\gt 0$ ), from figure 8 we find that the resultant force on the particle equals $F_h$ , which is dominated by the sum of the pressure gradient force $F_{\textit{pres}}$ , the added mass force $F_{add}$ , and the lubrication force $F_{\textit{lub}}$ . Here, we estimate the lower and upper bounds of the integral (3.15) as

(3.16) \begin{align} \int _{t_0}^{T_{\textit{end}}} F_{add}\, \text{d}t\lesssim \int _{t_0}^{T_{\textit{end}}} F_z\, \text{d}t \lesssim \int _{t_0}^{T_{\textit{end}}} F_{\textit{pres}}\, \text{d}t. \end{align}

From the Navier–Stokes equation, the pressure gradient force (2.5) is related to the fluid unsteady inertia as

(3.17)

where $\boldsymbol{u}$ is the undisturbed fluid velocity in the stagnation flow field. In (3.17), the viscous terms are omitted for the limited contribution of the viscous drag force. This expression is equivalent to the undisturbed force in Annamalai & Balachandar (Reference Annamalai and Balachandar2017). By simplifying the right-hand side of (3.17) with the undisturbed fluid velocity at the particle centre ( $z=z_p$ ), we obtain

(3.18) \begin{align} \begin{aligned} F_{\textit{pres}}\approx \frac {4}{3}\pi R_{p}^3\rho _{{l}}\frac {\text{D}{u_z}}{\text{D}t}\Bigg |_{(z=z_p)}. \end{aligned} \end{align}

The added mass force is linked to the pressure gradient force as $F_{add}\approx ({1}/{2})F_{\textit{pres}}-({1}/{2})( {\rho _l}/{\rho _p})F_z$ .

From (3.14)–(3.18), we obtain

(3.19) \begin{align} \begin{aligned} v_{p,\textit{end}}&\lesssim \frac {\rho _{{l}}}{\rho _{p}}\left [u_z{\big(z=z_p,t=T_{\textit{end}} \big)}-u_z{ \big(z=z_p,t=t_0 \big)}\right ]\!,\\v_{p,\textit{end}}&\gtrsim \frac {\rho _{{l}}}{2\rho _{p}+\rho _{{l}}}\left [u_z{\big(z=z_p,t=T_{\textit{end}} \big)}-u_z{\big(z=z_p,t=t_0 \big)}\right ]\!. \end{aligned} \end{align}

Recalling (3.9), we recognise that at the end of the bubble collapse ( $R_{b}\to 0$ ), the parameter $B(T_{\textit{end}})$ approaches zero. Thus we can infer $u_z{(z=z_p,t=T_{\textit{end}})}\to 0$ . From the evolution of the particle’s displacement with time (figure 7), we approximate the particle’s velocity at the end of bubble collapse $v_{p,\textit{end}}$ as twice the average velocity ${\Delta z_{\textit{end}}}/{\Delta T}$ , with $\Delta T=T_{\textit{end}}-t_0$ .

To solve for the right-hand side of (3.19), we now need to derive the flow velocity $u_z(z=z_p,t=t_0)$ at the centre of the particle at time $t=t_0$ when the particle is about to leave the substrate ( $z_p=R_p$ ). By comparing the position of the particle centre and the boundary layer thickness, we consider two regimes: (i) when the particle centre is far beyond the boundary layer at $t_0$ , i.e. $\delta (t_0)/R_p\ll 1$ , the particle detachment is dominated by the inviscid flow; and (ii) when the particle centre position is comparable to or lower than the boundary layer thickness at $t_0$ , i.e. $\delta (t_0)/R_p\gtrsim 1$ , the particle detachment is dominated by the viscous flow. In regime (i), we recall (3.5) and have

(3.20) \begin{align} u_z \big(z=R_p,t=t_0 \big)=-2B(t_0)\,R_p\quad \text{for}\ \delta (t_0)/R_p\ll 1. \end{align}

In regime (ii), we recall (3.9) and have

(3.21) \begin{align} u_z\big(z=R_p,t=t_0 \big)=-2B(t_0)\,f(R_p)\quad \text{for}\ \delta (t_0)/R_p\gtrsim 1. \end{align}

To estimate the vertical liquid velocity $u_z(z=R_p,t=t_0)$ for both regimes, we need to approximate the value of $B(t_0)$ .

For regime (i) with $\delta (t_0)/R_p\ll 1$ , from figures 8(a) and 8(b) we find that at $t=t_0$ , the hydrodynamic force $F_h$ is dominated by the pressure gradient force, thus with $F_{\textit{pres}}=0$ . Substituting (3.18) into $F_{\textit{pres}}(t_0)=0$ yields solving $\text{D}u_z/\text{D}t=0$ at $t=t_0$ .

For regime (ii) with $\delta (t_0)/R_p\gtrsim 1$ , from figures 8(c) and 8(d) we find that at $t=t_0$ , the hydrodynamic force $F_h$ is dominated by the added mass force, thus with $F_{add}=0$ . Substituting (2.6) and $a_p(t_0)=\text{d}v_p/\text{d}t=0$ into $F_{add}(t_0)=0$ also yields solving $\text{D}u_z/\text{D}t=0$ at $t=t_0$ .

We now focus on regime (ii). Substituting (3.9) into $\text{D}u_z/\text{D}t=0$ at $z=R_p$ and $t=t_0$ , we obtain

(3.22) \begin{align} \big (\dot {B}f-2B^2ff'\big )\Big |_{(z=R_{p},\ t=t_0)}=0. \end{align}

Since both $f(R_p)$ and $B(t_0)$ are non-zero, solving (3.22) for $t_0$ reduces to solving ${\dot {B}}/({2B^2})=f'(R_p)$ .

When the particle centre position is comparable to the boundary layer thickness at $t_0$ , i.e. $\delta (t_0)/R_p\sim 1$ , we assume $f'(R_p)\approx 1$ . This leads to

(3.23) \begin{align} \left (R_{{b}}^2\ddot {R}_{{b}}+2R_{{b}}\dot {R}_{{b}}^2-\frac {4}{L^3}R_b^4\dot {R}_b^2\right )\Bigg |_{(t=t_0)}=0. \end{align}

Eliminating $\ddot {R}_{{b}}$ in (3.3) by using (3.23), we express the bubble wall velocity at $t=t_0$ as

(3.24) \begin{align} \dot {R}_{{b}}|_{(t=t_0)}=\sqrt {\frac {\Delta p}{\rho _{l}}\left [-\frac {2}{\gamma ^4}\left (\frac {R_b}{R_{b,\textit{max}}}\right )^4-\frac {4}{\gamma ^3}\left (\frac {R_b}{R_{b,\textit{max}}}\right )^3+\frac {1}{2}\right ]^{-1}}\Bigg |_{(t=t_0)}. \end{align}

To further eliminate $\dot {R}_b$ , we can leverage the solution provided by Best & Blake (Reference Best and Blake1994) for the modified Rayleigh equation for a growing cavitation bubble near a solid wall:

(3.25) \begin{align} \dot {R}_{{b}}\,\sqrt {\frac {\rho _{{l}}}{\Delta p}}=\sqrt {\frac {2}{3}\left (1+\frac {1}{2\gamma }\frac {R_{{b}}}{R_{{b,{max}}}}\right )^{-1}\left [\left (\frac {R_{{b,{max}}}}{R_{{b}}}\right )^3-1\right ]}. \end{align}

Combining (3.24) and (3.25), we obtain

(3.26) \begin{align} & \left [\frac {4}{\gamma ^4}\left (\frac {R_{{b}}}{R_{{b,{max}}}}\right )^7+\frac {8}{\gamma ^3}\left (\frac {R_{{b}}}{R_{{b,{max}}}}\right )^6-\left (\frac {3}{2\gamma }+\frac {4}{\gamma ^4}\right )\left (\frac {R_{{b}}}{R_{{b,{max}}}}\right )^4\right. \nonumber \\& \quad \left. -\left (4+\frac {8}{\gamma ^3}\right )\,\left (\frac {R_{{b}}}{R_{{b,{max}}}}\right )^3+1\right ]\Bigg |_{(t=t_0)}=0. \end{align}

The solution to (3.26) depends on the dimensionless stand-off distance $\gamma$ . Numerical calculations show that for the range of $\gamma$ considered in this work (1.0–5.0), the solution $ ({R_{b}}/{R_{b,\textit{max}}} )|_{(t=t_0)}$ falls within the range $({0.419},0.618)$ , and the non-dimensional time $t_0/T_c$ ranges from 0.0608 to 0.159, both increasing with increasing $\gamma$ . Then the value of $B(t_0)$ can be calculated through

(3.27) \begin{align} B(t_0)=\frac {2}{\gamma ^3}\,G(\gamma )\,\sqrt {\frac {\Delta p}{\rho _{{l}}}}\,R_{{b,{max}}}^{-1}, \end{align}

with a function

(3.28) \begin{align} G(\gamma )=\frac {\big(R_b/R_{b,\textit{max}}\big)|_{(t=t_0)}^2}{\sqrt {-\dfrac {2}{\gamma ^4}\left (\dfrac {R_b}{R_{b,\textit{max}}}\right )\Bigg |_{(t=t_0)}^4-\dfrac {4}{\gamma ^3}\left (\dfrac {R_b}{R_{b,\textit{max}}}\right )\Bigg |_{(t=t_0)}^3+\dfrac {1}{2}}}. \end{align}

The function $G(\gamma )$ increases with increasing $\gamma$ , and its value ranges between 0.463 and 0.544 with $\gamma \in ({1.0},5.0)$ . Therefore, within this range, we can approximate $G(\gamma )$ as a constant, i.e. $G(\gamma )\approx 0.5$ . Consequently, the value of $B(t_0)$ is approximated as

(3.29) \begin{align} B(t_0)\,\approx \,\frac {1}{\gamma ^3}\sqrt {\frac {\Delta p}{\rho _{{l}}}}\,R_{{b,{max}}}^{-1}=\gamma ^{-3}\,T_c^{-1}. \end{align}

When the particle centre is far below the edge of the boundary layer at $t_0$ , i.e. $\delta (t_0)/R_p\gg 1$ , we assume $f'(R_p)\to 0$ . This leads to

(3.30) \begin{align} \big (2R_{{b}}\dot {R}_{{b}}^2+R_{{b}}^2\ddot {R}_{{b}}\big )\Big |_{(t=t_0)}=0. \end{align}

Following a similar procedure from (3.24) to (3.29) (for details see Appendix D), we arrive at

(3.31) \begin{align} B(t_0)\approx \gamma ^{-3}T_c^{-1}. \end{align}

The result above coincides with (3.29) for $f'(R_p)\to 1$ .

To estimate the undisturbed liquid velocity $u_z(z=R_p,t=t_0)$ for $\delta (t_0)/R_p\gtrsim 1$ , we can simplify the solution to (3.12) by approximating $f'''(z)$ as a linear function, i.e. $f'''(z)\approx -f'''(0)(z/\delta _0-1)$ , where $\delta (t_0)$ is denoted as $\delta _0$ for simplification. By integration while combining the boundary conditions ( $f'(0)=f(0)=0$ ) and $f'(\delta _0)=0.99$ , we obtain

(3.32) \begin{align} f(z)&\approx 0.2475{\delta _0}\left (\frac {z}{\delta _0}-1\right )^4+0.99z-0.2475{\delta _0}, \nonumber\\f'(z)&\approx 0.99\left (\frac {z}{\delta _0}-1\right )^3+0.99,\nonumber\\f''(z)&\approx \frac {2.97}{\delta _0}\left (\frac {z}{\delta _0}-1\right )^2 \!,\nonumber\\f'''(z)&\approx \frac {5.94}{\delta _0^2}\left (\frac {z}{\delta _0}-1\right )\!, \end{align}

with $f'''(0)\approx -5.94/\delta _0^2$ under the condition $0\leqslant z\leqslant \delta _0$ . Therefore, $f(R_p)$ is a function of $R_p$ and $\delta _0$ . To determine $\delta _0$ , from (3.12) and (3.22), we notice $f'''(0)=-( {B(t_0)}/{\nu _{l}})-[\dot{B}(t_0)/(\nu_{l} B(t_0))] =-({B(t_0)}/{\nu _{l}}) [1+2\,f'(R_p) ]$ . Combining this relation with $f'(R_p)\approx 0.99(R_p/\delta _0-1)^3+0.99$ and the expression for $B(t_0)$ (both (3.29) and (3.31) give the same result), we derive an equation for $\delta _0/R_p$ :

(3.33) \begin{align} \left (\frac {\delta _0}{R_p}\right )^3+5.94\left (\frac {\delta _0}{R_p}\right )^2-5.94\left [1+\gamma ^3\left (\frac {R_{b,\textit{max}}}{R_p}\right )^2\textit{Re}_b^{-1}\right ]\left (\frac {\delta _0}{R_p}\right )+1.98=0, \end{align}

with the bubble Reynolds number based on the maximum bubble radius

(3.34) \begin{align} \textit{Re}_{b}=\frac {R_{b,\textit{max}}\sqrt {\Delta p/\rho _{l}}}{\nu _{l}}. \end{align}

The solution to (3.33) can be evaluated by numerical calculation.

For $\delta _0/R_p\ll 1$ , we have $f'(R_p)=1$ , and then $f'''(0)=- {3\,B(t_0)}/{\nu _{l}}\approx - {5.94}/{\delta _0^2}$ , leading to $\delta _0\approx \sqrt { {1.98\nu _{l}}/{B(t_0)}}$ , or

(3.35) \begin{align} \frac {\delta _0}{R_p}\approx 1.41\gamma ^{{3}/{2}}\left (\frac {R_{b,\textit{max}}}{R_p}\right )\textit{Re}_b^{-{1}/{2}}. \end{align}

The solutions shown in figure 12 indicate that the linear relationship (3.35) holds approximately for the range of $\delta _0/R_p$ from 0 to 1.5. In figure 12, $\delta _0/R_p$ from the simulation results (figure 10) shows a similar trend to the lines, validating the above theoretical derivations.

Figure 12.

Boundary layer thickness at particle detachment ( $\delta _0/R_p$ ) as a function of the stand-off distance $\gamma$ , the bubble–particle radius ratio $R_{b,\textit{max}}/R_p$ , and the bubble Reynolds number $\textit{Re}_b$ . The solid grey line denotes the combination of the linear solution (3.35) and the numerical solution to (3.33). The dashed black line denotes the linear relationship only. The data points denote $\delta _0/R_p$ in simulations (cases A–D) read from figure 10.

Finally, to simplify (3.19), we write the non-dimensionalised forms of (3.20) and (3.21), and combine them as

(3.36) \begin{align} \frac {u_z(z=R_p,t=t_0)}{\sqrt {{\Delta p}/\rho _{{l}}}}\approx -2\gamma ^{-3}\left (\frac {R_{b,\textit{max}}}{R_p}\right )^{-1}\mathcal{F}, \end{align}

with a function

(3.37) \begin{align} \mathcal{F}=\begin{cases} \displaystyle 1 & \text{for}\ \delta _0/R_p\ll 1,\\ \dfrac {f(R_p)}{R_p} & \text{for}\ \delta _0/R_p\gtrsim 1. \end{cases} \end{align}

From figure 12, we identify $\mathcal{F}$ as a function of $\gamma ^{{3}/{2}} ( {R_{b,\textit{max}}}/{R_p} )\,\textit{Re}_b^{-{1}/{2}}$ . Therefore, we arrive at

(3.38) \begin{align} \begin{aligned} \frac {v_{p,\textit{end}}}{\sqrt {{\Delta p}/\rho _{{l}}}}&\lesssim 2\frac {\rho _l}{\rho _p}\gamma ^{-3}\left (\frac {R_{b,\textit{max}}}{R_p}\right )^{-1}\mathcal{F},\\ \frac {v_{p,\textit{end}}}{\sqrt {{\Delta p}/\rho _{{l}}}}&\gtrsim \frac {2\rho _l}{2\rho _p+\rho _l}\gamma ^{-3}\left (\frac {R_{b,\textit{max}}}{R_p}\right )^{-1}\mathcal{F}. \end{aligned} \end{align}

Then we simplify the relationship $v_{p,\textit{end}}\sim \,2\,\Delta z_{\textit{end}}/\Delta T$ by approximating $T_{\textit{end}}/T_c$ as $1.83(1+0.205/\gamma )$ (Rattray Reference Rattray1951; Reuter et al. Reference Reuter, Zeng and Ohl2022). For the range of $\gamma$ considered in our work (1.0–5.0), and recalling the range of $t_0/T_c$ , the ratio $\Delta T/T_c=(T_{\textit{end}}-t_0)/T_c$ is approximately equal to 2. Then we have

(3.39) \begin{align} \frac {2\rho _l}{2\rho _p+\rho _l}\gamma ^{-3}\mathcal{F}\lesssim \frac {\Delta z_{\textit{end}}}{R_p}\lesssim 2\frac {\rho _l}{\rho _p}\gamma ^{-3}\,\mathcal{F}. \end{align}

Based on this, figure 13 shows the relationship between $({\Delta z_{\textit{end}}}/{R_p}) \,({\rho _p}/{\rho _l}) \,\gamma ^3$ and $\gamma ^{{3}/{2}} ( {R_{b,\textit{max}}}/{R_p} )\,\textit{Re}_b^{-{1}/{2}}$ , with the results of experiments and simulations. When $\gamma ^{{3}/{2}} ( {R_{b,\textit{max}}}/{R_p} )\,{\textit{Re}}_b^{-{1}/{2}}$ is below $\sim 0.2$ , the inviscid flow dominates the particle motion (regime (i), corresponding to $\delta _0/R_p\ll 1$ ), leading to $ {\Delta z_{{end}}}/{R_p}({\rho _p}/{\rho _l})\gamma ^3\sim 2$ (dash-dotted black line). This regime is supported by the experimental cases of PMMA particles with $R_p\,={250}\,\unicode{x03BC}\mathrm{ m}$ and the simulation cases A and B. When $\gamma ^{{3}/{2}} ( {R_{b,\textit{max}}}/{R_p} )\,{\textit{Re}}_b^{- {1}/{2}}$ exceeds $\sim 0.2$ (corresponding to a larger boundary layer thickness comparable to $R_p$ ), the dominating mechanism for particle motion transfers to regime (ii), leading to the decrease of $({\Delta z_{\textit{end}}}/{R_p})\,({\rho _p}/{\rho _l}) \,\gamma ^3$ with the increase of $\gamma ^{{3}/{2}} ( {R_{b,\textit{max}}}/{R_p} )\,\textit{Re}_b^{-{1}/{2}}$ , as supported by the experimental cases of PMMA particles with $R_p={50}\,\unicode{x03BC}\mathrm{ m}$ , and PS particles with $R_p={25}\,\unicode{x03BC}\mathrm{ m}$ , as well as the simulation cases C and D. The two lower bounds of the shaded region (regime (ii)) arise from the different density ratios $ {\rho _l}/{\rho _p}$ , but the two curves almost overlap each other for the small density differences between PMMA and PS.

Figure 13.

Relationship between $({\Delta z_{\textit{end}}}/{R_p}) \,({\rho _p}/{\rho _l})\,\gamma ^3$ and $\gamma ^{{3}/{2}}( {R_{b,\textit{max}}}/{R_p})\,\textit{Re}_b^{-{1}/{2}}$ . Data points for experiments are differentiated by markers representing the type and radius of the particles, and the simulation results (cases A–D) are presented, as shown in the legend. The shaded region in regime (ii) represents (3.39) with the upper and lower bounds for different particle densities. For the experimental data, bubble parameters obtained from the partial-boundary fitting method are converted into their full-boundary equivalents.

It is important to clarify the limitations of the theoretical analysis in this work. The assumptions made here, namely that the particle has negligible influence on the flow field around the cavitation bubble, and that the flow field coordinates $(r,z)$ at the particle’s location satisfy $r,z\ll L,R_{{b,{max}}}$ , may not hold for the cases with ${R_{{b,{max}}}}/{R_{p}}\lesssim 1$ or bubble–particle contact. In such scenarios, the pressure distribution on the particle surface should be considered, and the particle’s motion may be dominated by the net pressure difference caused by the vapour pressure inside the bubble during contact, which has been described by Ren et al. (Reference Ren, Zuo, Wu and Liu2022).