1. Introduction
The collision between particles and dispersed fluid–fluid interfaces (including bubbles and droplets) is commonly encountered in various industrial processes such as spraying (Gong, Li & Kang Reference Gong, Li and Kang2022; Zhang et al. Reference Zhang, Wang, Cheng, Ma, Gao, Brennan and Yan2023), flotation (Hassanzadeh et al. Reference Hassanzadeh, Firouzi, Albijanic and Celik2018; Han et al. Reference Han, Fang, He, Wang, Li, Wu, Cai and Wang2022), catalytic reactions (Chen et al. Reference Chen, Zhang, Xu and Dong2023) and liquid cooling (Chen et al. Reference Chen, Zhang, Jia, Yin, Ni, Xin and Dong2025b ). During collisions, the influence of particles can cause behaviours such as spreading, rebound and breakup of drops. These behaviours significantly affect industrial processes, such as altering heat transfer performance (Fardan, Berndt & Ahmed Reference Fardan, Berndt and Ahmed2021) and the active surface area of catalytic particles (Ge & Fan Reference Ge and Fan2007). Therefore, investigating the dynamic behaviour of collisions is crucial for optimising these processes.
Dispersed fluid–fluid interfaces mainly exist in the form of droplets or bubbles. The particle–droplet collision problem can be categorised into three scenarios based on the initial conditions: droplets impacting stationary particles, particles impacting stationary droplets and the collision of moving particles and moving droplets. There has been extensive research on the impact of droplets on particle surfaces (Pawar et al. Reference Pawar, Henrikson, Finotello, Padding, Deen, Jongsma, Innings and Kuipers2016; Janssen et al. Reference Janssen, Fetzer, Grewing, Burgmann and Janoske2023). After impact, droplets undergo three stages: the initial droplet deformation, an inertia-control stage and a viscous dominance stage (Mitra et al. Reference Mitra, Evans, Doroodchi, Pareek and Joshi2017). The governing factors of these processes can be characterised by two dimensionless parameters, namely the Weber number (We) and Reynolds number (Re) (Yoon & Shin Reference Yoon and Shin2021). Additionally, the surface characteristics of the particles are key factors affecting impact behaviour. Under varying impact conditions, droplets may exhibit different phenomena such as spreading, rebound, breakup and splashing (Wu et al. Reference Wu, Chen, Du, Zhai, Zeng, Yu and Zhou2022). In the case of particles impacting droplets, two behaviours can be observed depending on the impact conditions: partial submersion and complete submersion (penetration) (Mitra et al. Reference Mitra, Doroodchi, Pareek, Joshi and Evans2015). During the impact process, particles are primarily influenced by gravity and capillary forces. If the particle velocity has decreased to zero before complete submersion, this is considered partial submersion. Otherwise, the particle will penetrate the droplet. When suspended particles collide with suspended droplets, four behaviours can occur: agglomeration, separation, penetration and breakage (Singh & Joseph Reference Singh and Joseph2005; Deen, Annaland & Kuipers Reference Deen, Van Annaland and Kuipers2009). Although significant research has been conducted on particle–droplet collision phenomena, many findings do not universally apply to all scaling laws (Cheng, Sun & Gordillo Reference Cheng, Sun and Gordillo2021). Numerous fundamental questions remain to be explored further.
Bubbles in a second-phase fluid typically exist in two forms due to the significant buoyancy forces acting on them: suspended within the fluid and attached to a solid wall (Brabcová et al. Reference Brabcová, Karapantsios, Kostoglou, Basařová and Matis2015; Mohammadi et al. Reference Mohammadi, Nazari, Kayhani and Ahmadi2022). Current research on particle–bubble interaction primarily focuses on flotation processes. The success of particle flotation hinges on two core aspects (Wang et al. Reference Wang, Huang, Yao, Yang, Chen and Wang2025): first, the rupture of the hydration film between the particle and the bubble, which forms a three-phase contact line; second, the external forces acting on the particle must be less than the forces between the particle and the bubble, preventing the particle from detaching from the bubble surface. The flotation process is mainly influenced by factors such as particle size, wettability and the relative velocity between particles and bubbles. It is evident that the flotation process emphasises the movement of particles on the bubble surface, with less attention paid to bubble behaviour.
In boiling processes, the nucleation, growth and detachment of bubbles on heated surfaces have a decisive impact on boiling heat transfer performance (Bureš & Sato Reference Bureš and Sato2022; Chen et al. Reference Chen, Yu, Lu, Wang, Sun, Jiao, Zhang and Tao2024; Park et al. Reference Park, Ahmadi, Foulkes and Boreyko2024). When fluid media contain particles, complex interactions occur between the particles and bubbles (Hsu et al. Reference Hsu, Chiu, Tsai, Chou and Sheng2024). Iwata et al. (Reference Iwata, Takeuchi and Kajishima2010) and Washino et al. (Reference Washino, Chan, Tsujimoto, Tsuji and Tanaka2023) studied the rising behaviour of bubbles in sparse spherical particle clusters and irregular particle clusters, respectively. Yu et al. (Reference Yu, Yang and Wang2024) investigated the rising behaviour of bubbles in dense particle beds. However, these studies did not address the effect of particles on bubble detachment from heated surfaces. In our previous study (Chen et al. Reference Chen, Yin, Zhang, Xin, Ni and Dong2025a ), we investigated the wetting dynamics of a bubble lying between a free particle and a solid horizontal surface. The particle was initially at rest, and the interaction was governed predominantly by capillary effects. Upon contact, capillary waves were generated and propagated along the vapour–liquid interface, driving the motion of the three-phase contact line and eventually inducing bubble detachment from the solid surface. The important role of capillary waves in interfacial transport has also been demonstrated in other multi-phase systems. Huang et al. (Reference Huang, Lori, Yang, Zhuang, Cheng, Chen, Sun, Ming, Cui and Cheng2026) recently demonstrated centimetre-scale droplet jumping on superhydrophobic surfaces induced by capillary waves generated during entrained bubble bursting.
However, in practical particulate boiling systems and vapour–liquid–solid three-phase flows, particles typically collide with bubbles at finite velocities. Under such conditions, particle inertia acts together with capillary forces, giving rise to interfacial dynamics that differs from that arising solely from capillary effects. When particles collide with bubbles on a heated surface, phenomena similar to particles impacting droplets can be observed, where particles may become attached to the bubble surface or submerge within the bubble (Yang et al. Reference Yang, Zhang, Binama, Li and Cai2023). Meanwhile, the bubble may undergo compression, rebound, interfacial-wave propagation or even transitions between different detachment modes, thereby altering the local flow structure and affecting local heat and mass transfer processes (Rahmati et al. Reference Rahmati, Ebrahimi-Dehshali and Hakkaki-Fard2020). Therefore, a systematic investigation of particle-impact-induced bubble detachment is required to elucidate the underlying interfacial dynamics and detachment mechanisms.
To the best of our knowledge, the dynamical mechanisms by which a solid particle impacts a wall-adhered bubble and induces its detachment have not yet been systematically explored. Unlike our previous studies that focused on capillary-driven wetting processes involving a particle initially at rest (Chen et al. Reference Chen, Yin, Zhang, Xin, Ni and Dong2025a ), the present work investigates impact-induced detachment of wall-adhered bubbles under the combined action of particle inertia and capillary effects. Furthermore, the present study considers an axisymmetric particle–bubble system, enabling a more realistic representation of impact-induced bubble deformation and detachment than the two-dimensional planar configuration considered previously. To address this problem, the three-phase flow model based on the coupled volume-of-fluid and level set method, immersed boundary method and discrete element method (VOSET-IBM-DEM) is extended to an axisymmetric framework and employed to numerically investigate the post-impact evolution of a single particle–bubble system, with particular focus on bubble detachment. By incorporating the Kelvin impulse framework together with an energy-based analysis, we identify distinct bubble-detachment regimes and examine how particle properties (including wettability, density ratio and size ratio) affect the detachment behaviour over a range of Weber numbers. The dependence of detachment time and regime transition on the governing physical parameters is elucidated, and the underlying mechanisms are interpreted in terms of interfacial energy release and momentum evolution.
2. Problem statement
We conduct numerical simulations to investigate the dynamics of a particle impacting a bubble attached to a solid wall under an axisymmetric framework. In practical situations, bubbles deform into an ellipsoidal shape due to buoyancy force. To realistically reproduce the bubble morphology, a preliminary simulation of a vapour–liquid two-phase flow without a particle is first performed. The computational configuration for the bubble pre-simulation is shown in figure 1(a). The equivalent bubble diameter is set to D
b
= 1.5 mm, and the contact angle at the bottom wall is fixed at
$\theta _{w}$
= 45°. The initial vertical distance between the bubble centre and the bottom wall, H
2, as well as the effective bubble radius
${R}_{b}'$
, are determined based on the bubble diameter D
b
and the wall contact angle
$\theta _{w}$
. The computational domain has dimensions of L = 1.5D
b
in the horizontal direction and H
1 = 3D
b
in the vertical direction.
Physical properties of water (liquid and vapour).
Table 1. Long description
A table comparing physical properties of water in liquid and vapour states. The table has three rows and four columns. The columns are labeled Medium, Density (kg m3), Dynamic viscosity (Pa s), and Surface tension coefficient (N m1). The rows are labeled Liquid and Vapour. Row 1: Medium, Liquid, Vapour. Row 2: Density (kg m3), 958, 0.6. Row 3: Dynamic viscosity (Pa s), 2.82 x 10^-4, 1.23 x 10^-5. Row 4: Surface tension coefficient (N m1), 0.059.
Schematic of the computational configuration: (a) pre-simulation stage for establishing the equilibrium shape of the surface-attached bubble and (b) particle–bubble impact stage after introducing the particle.
After the bubble reaches a stable configuration, a particle is introduced above the bubble with an initial downward velocity V. The particle diameter is denoted as D p , and the detailed configuration is illustrated in figure 1(b). The bubble diameter D b and the particle impact velocity V are chosen as the characteristic length scale and characteristic velocity, respectively.
Deionised water is selected as the working fluid in this study. The physical properties of the liquid and vapour phases correspond to those at the saturation condition under atmospheric pressure and are listed in table 1. The particle contact angle is denoted by
$\theta _{p}$
, and the density of the particles is referred to as
$\rho _{s}$
. The density ratio between the particles and the liquid phase is defined as
$k_{\rho }=\rho _{s}/\rho _{l}$
, and the size ratio between the particles and the bubble is defined as
$k_{r}= D_{p}/D_{b}$
.
In the boiling of deionised water at low superheating, the entire cycle from nucleation to growth of bubbles spans several tens to hundreds of milliseconds (Siegel & Keshock Reference Siegel and Keshock1964; Kunkelmann Reference Kunkelmann2011; Ling, Li & Tao Reference Ling, Li and Tao2014). The results of this study indicate that the entire process of bubble detachment following particle impact occurs within 20 ms. Therefore, the effect of bubble growth is not considered in this study.
3. Numerical methods
The collision process of a particle into a bubble is simulated using the VOSET-IBM-DEM developed by our group (Chen, Yin & Dong Reference Chen, Yin and Dong2024a
,
Reference Chen, Yin and Dongb
). This method is based on the finite-volume method, which solves for the fluid and particles in a fixed and uniform Cartesian grid. The VOSET method was utilised to capture the vapour–liquid interface (Sun & Tao Reference Sun and Tao2010), while the IBM was applied to resolve the momentum exchange between particles and fluid (Kajishima, Takiguchi & Hamasaki Reference Kajishima, Takiguchi and Hamasaki2001). Additionally, the DEM with soft-sphere collision model was used to calculate the collision forces
$\boldsymbol{f}_{\!\boldsymbol{C}}$
between solids (Tsuji, Kawaguchi & Tanaka Reference Tsuji, Kawaguchi and Tanaka1993). This section provides a brief overview of the hybrid method, with details available in the literature (Chen, Yin & Dong Reference Chen, Yin and Dong2024a
,
Reference Chen, Yin and Dongb
).
The IDEAL algorithm was applied to solve the incompressible mass conservation equation and the Navier–Stokes equations to obtain the velocity field and pressure field (Sun et al. Reference Sun, Qu, He and Tao2008a , Reference Sun, Qu, He and Taob ). The dimensionless governing equations are as follows:
\begin{equation}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}=0,\end{equation}
\begin{equation}\frac{\partial \boldsymbol{u}}{\partial t}+\left(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\right)\boldsymbol{u}=-\boldsymbol{\nabla }\!p+\frac{1}{\textit{Re}}\boldsymbol{\nabla} ^{2}\boldsymbol{u}-\frac{1}{\textit{Fr}}\boldsymbol{j}+\frac{1}{\textit{We}}\boldsymbol{f}_{\!\boldsymbol{\sigma }}+\boldsymbol{f}_{\!\boldsymbol{I}}.\end{equation}
Herein,
$\boldsymbol{u}$
is the velocity vector,
$t$
is time,
$p$
is pressure,
$\boldsymbol{j}=\boldsymbol{g}/| \boldsymbol{g}|$
denotes the unit vector in the direction of the gravity and
$\boldsymbol{g}$
is the gravitational acceleration. The Reynolds, Froude and Weber numbers are defined as
$\textit{Re}=\rho _{l}UL/\mu _{l}$
,
$\textit{Fr}=U^{2}/gL$
and
$\textit{We}=\rho _{l}U^{2}L/\sigma$
, respectively, where
$U$
and
$L$
denote the characteristic velocity and length,
$\mu _{l}$
is the liquid dynamic viscosity and
$\sigma$
is the surface tension coefficient.
The VOSET method is used to capture the vapour–liquid interface. The advection equation for the volume fraction C f was solved to advance the vapour–liquid interface:
\begin{equation}\frac{\partial C_{\!f}}{\partial t}+\boldsymbol{\nabla }\boldsymbol{\cdot }\left(\boldsymbol{u}C_{\!f}\right)=0.\end{equation}
The geometric iterative algorithm was employed to solve the signed distance function
$\phi$
in the vicinity of the interface. This method combines the excellent mass-conservation property of the volume-of-fluid method (Hirt & Nichols Reference Hirt and Nichols1981) with the accurate evaluation of interface geometric properties provided by the level-set method (Sussman et al. Reference Sussman, Almgren, Bell, Colella, Howell and Welcome1996). Meanwhile, the computational cost is reduced by restricting the reconstruction
$\phi$
to only a few grid layers around the interface, rather than updating
$\phi$
throughout the entire computational domain.
The surface tension
$\boldsymbol{f}_{\!\boldsymbol{\sigma }}$
is calculated using our improved continuous surface force model (Popinet Reference Popinet2018; Chen, Yin & Dong Reference Chen, Yin and Dong2024b
):
\begin{equation}\boldsymbol{f}_{\!\boldsymbol{\sigma }}=-\sigma k\boldsymbol{\nabla }H_{\phi }.\end{equation}
The curvature
$k$
of the vapour–liquid interface is solved using the height function technique (Popinet Reference Popinet2009). Function
$H_{\phi }$
is the Heaviside function associated with the signed distance function
$\phi$
.
The interaction force
$\boldsymbol{f}_{\!\boldsymbol{I}}$
between the particles and the fluid is calculated using the immersed solid method proposed by Kajishima et al. (Reference Kajishima, Takiguchi and Hamasaki2001). In this approach, a solid volume fraction field, C
s
, is introduced to represent the local occupancy of the particle within each computational cell. The solid volume fraction is evaluated as (Yuki et al. Reference Yuki, Takeuchi and Kajishima2007)
\begin{equation}C_{s}=\frac{1}{2}\left[1-\tanh \left(\frac{\delta _{p}-R_{p}}{\lambda \xi {\unicode[Arial]{x0394}} x}\right)\right]\!, \end{equation}
\begin{equation}\lambda =\left| n_{p,x}\right| +\left| n_{p,y}\right|\!, \end{equation}
\begin{equation}\xi =0.3+0.05(1-\lambda ^{2}),\end{equation}
where
$\delta _{p}$
is the distance between the grid node and the particle centre,
$R_{p}$
is the particle radius,
${\unicode[Arial]{x0394}} x$
is the grid spacing,
$\lambda$
is an interface-orientation parameter,
$\boldsymbol{n}_{\boldsymbol{p}}=$
(
$n_{p,x},n_{p,y}$
) is the unit normal vector at the solid–liquid interface and
$\xi$
is an adaptive smoothing-width coefficient.
As the particle moves through the computational domain, the solid volume fraction field is recalculated at every time step using the instantaneous particle position and radius. Therefore, cells crossed by the moving particle are naturally treated through the continuous redistribution of C s on the fixed Cartesian grid, without requiring mesh deformation or remeshing.
The local rigid-body velocity
$ \boldsymbol{u}_{\boldsymbol{s}}$
of the particle is defined as
\begin{equation}\boldsymbol{u}_{\boldsymbol{s}}=\boldsymbol{v}_{\!s}+\boldsymbol{\omega }_{\boldsymbol{s}}\times \boldsymbol{r}_{\!\boldsymbol{s}},\end{equation}
where
$\boldsymbol{v}_{\!s}$
and
$\boldsymbol{\omega }_{\boldsymbol{s}}$
denote the translational and angular velocities of the particle, respectively, and
$\boldsymbol{r}_{\!\boldsymbol{s}}$
is the position vector measured from the particle centre.
The interaction force between the particle and the fluid is evaluated as
\begin{equation}\boldsymbol{f}_{\!\boldsymbol{I}}=\frac{C_{s}\big(\boldsymbol{u}_{\boldsymbol{s}}-{\boldsymbol{u}}_{\!\boldsymbol{f}}^{\boldsymbol{F}\boldsymbol{1}}\big)}{{\unicode[Arial]{x0394}} t},\end{equation}
where
${\boldsymbol{u}}_{\!\boldsymbol{f}}^{\boldsymbol{F}\boldsymbol{1}}$
is the intermediate fluid velocity obtained from solving the Navier–Stokes equations without considering the particle effect and
${\unicode[Arial]{x0394}} t$
is the time step.
The corrected velocity field is then given by
\begin{equation}\boldsymbol{u}^{\boldsymbol{n}+\boldsymbol{1}}={\boldsymbol{u}}_{\!\boldsymbol{f}}^{\boldsymbol{F}\boldsymbol{1}}+{\unicode[Arial]{x0394}} t\boldsymbol{f}_{\!\boldsymbol{I}}=C_{s}\boldsymbol{u}_{\boldsymbol{s}}+(1-C_{s}){\boldsymbol{u}}_{\!\boldsymbol{f}}^{\boldsymbol{F}\boldsymbol{1}}.\end{equation}
Herein,
$\boldsymbol{u}^{\boldsymbol{n}+\boldsymbol{1}}$
is the velocity after the immersed boundary correction. Consequently, the immersed boundary forcing constrains the velocity field within the particle region to follow the local rigid-body motion, while the interfacial cells provide a smooth transition between the fluid and solid phases. This treatment effectively enforces the no-slip condition on the moving particle surface.
Assuming the particle is a rigid sphere, its motion in the gas–liquid two-phase flow is governed by the gravitational force, the interaction force
$ \boldsymbol{F}_{\!\boldsymbol{p}}$
from surrounding fluids and the collision forces
$\boldsymbol{F}_{\!\boldsymbol{c}}$
between the particle and the wall. Therefore, the translational and angular velocities of the solid particle can be obtained by applying Newton’s law:
\begin{equation}m_{s}\frac{\text{d}\boldsymbol{v}_{\boldsymbol{s}}}{\text{d}t}=\boldsymbol{F}_{\!\boldsymbol{p}}+\boldsymbol{F}_{\!\boldsymbol{c}}+m_{s}\boldsymbol{g},\end{equation}
\begin{equation}\boldsymbol{I}_{\!\boldsymbol{s}}\frac{\text{d}\boldsymbol{\omega }_{\boldsymbol{s}}}{\text{d}t}=\boldsymbol{T}_{\!\boldsymbol{p}}+\boldsymbol{T}_{\!\boldsymbol{c}},\end{equation}
where
$m_{s}$
is the mass of the particle,
$\boldsymbol{I}_{\!\boldsymbol{s}}$
is the inertia tensor,
$\boldsymbol{T}_{\!\boldsymbol{p}}$
is the interaction torque from surrounding fluids and
$\boldsymbol{T}_{\!\boldsymbol{c}} $
is the torque exerted by the collision force.
The particle–wall collision force
$\boldsymbol{F}_{\!\boldsymbol{c}}$
is computed using the DEM soft-sphere collision model. For particle–wall collisions, a consistent set of contact parameters is employed throughout the simulations. The coefficients of restitution in the normal and tangential directions are both set to 0.9. The normal and tangential spring stiffnesses are 800 and 80 N m−1, respectively, and the friction coefficient is set to 0.3.
The interaction force
$ \boldsymbol{F}_{\!\boldsymbol{p}}$
is determined by the integration of the fluid stress tensor over the particle surface:
\begin{equation}\boldsymbol{F}_{\!\boldsymbol{p}}=\int _{S}\,(-p\boldsymbol{I}+\boldsymbol{\tau })\boldsymbol{\cdot }\boldsymbol{n}_{\boldsymbol{p}}\,\text{d}S+\int _{S}\,(-\boldsymbol{\varPi })\boldsymbol{\cdot }\boldsymbol{n}_{\boldsymbol{p}}\,\text{d}S.\end{equation}
The first term on the right-hand side represents hydrodynamics, which includes pressure and viscous stress, and the second term represents capillary forces. Here,
$\boldsymbol{\tau }$
denotes viscous stress and
$\boldsymbol{\varPi }$
is the capillary pressure tensor. The tensor form of surface tension is introduced, expressed through the capillary pressure tensor (Lafaurie et al. Reference Lafaurie, Nardone, Scardovelli, Zaleski and Zanetti1994), which satisfies
\begin{equation}\boldsymbol{f}_{\!\boldsymbol{\sigma }}=-\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\varPi }.\end{equation}
According to the Gaussian divergence theorem, the integral of the hydrodynamic force terms over the particle surface can be equivalently expressed as the volume integral of these terms over the particle volume; similarly, the capillary force can also be obtained by performing a volume integration of the divergence of the capillary pressure tensor, which yields (Sun & Sakai Reference Sun and Sakai2016)
\begin{equation}\boldsymbol{F}_{\!p}=\int _{P}\,\boldsymbol{\nabla }\boldsymbol{\cdot }(-p\boldsymbol{I}+\boldsymbol{\tau })\,\text{d}V+\int _{P}\,\boldsymbol{\nabla }\boldsymbol{\cdot }(-\boldsymbol{\varPi })\,\text{d}V.\end{equation}
Equation (3.2) in the fluid domain is integrated over the volume occupied by the particle:
\begin{equation}\int _{P}\,\rho _{\!f}\frac{\text{D}\boldsymbol{u}}{\text{D}t}\,\text{d}V=\int _{P}\,\left(-\boldsymbol{\nabla }\!p+\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\tau }+\rho _{\!f}\boldsymbol{g}+\boldsymbol{f}_{\!\boldsymbol{\sigma }}+\rho _{\!f}\boldsymbol{f}_{\!\boldsymbol{I}}\right)\,\text{d}V,\end{equation}
where
$\rho _{\!f}$
is the fluid density.
Substituting the expression for the total force exerted by the fluid on the particle into the above equation and rearranging the volume integral equation:
\begin{equation}\boldsymbol{F}_{\!\boldsymbol{p}}=\int _{P}\,\left(\frac{\text{D}\boldsymbol{u}}{\text{D}t}-\boldsymbol{g}-\boldsymbol{f}_{\!\boldsymbol{I}}\right)\rho _{\!f}\,\text{d}V.\end{equation}
Herein, the first term on the right-hand side is the inertial term caused by the acceleration of the particle motion. In the present study, the contribution of this inertial term is negligible compared with the hydrodynamic and capillary forces, and it is thus omitted. The second term represents the hydrostatic buoyancy force, which is separated from the dynamic interactions. The last term is the fluid–particle interaction force, which includes the hydrodynamic and capillary contributions.
Consequently, the governing equations for particle motion can be written as equations (3.18) and (3.19), where
$ {\unicode[Arial]{x0394}} V$
is the grid volume:
\begin{equation}m_{s}\frac{\text{d}\boldsymbol{v}_{\boldsymbol{s}}}{\text{d}t}=-\sum C_{s}\left(\boldsymbol{f}_{\!\boldsymbol{I}}+\boldsymbol{g}\right)\rho _{\!f}{\unicode[Arial]{x0394}} V+\boldsymbol{F}_{\!\boldsymbol{c}}+m_{s}\boldsymbol{g},\end{equation}
\begin{equation}\boldsymbol{I}_{\!\boldsymbol{s}}\frac{\text{d}\boldsymbol{\omega }_{\boldsymbol{s}}}{\text{d}t}=-\sum C_{s}\boldsymbol{r}_{\!\boldsymbol{s}}\times \left(\boldsymbol{f}_{\!\boldsymbol{I}}+\boldsymbol{g}\right)\rho _{\!f}{\unicode[Arial]{x0394}} V+\boldsymbol{T}_{\!\boldsymbol{c}}.\end{equation}
The curvature at the boundary is corrected based on the height function to satisfy the contact angle boundary condition on the plane (Li & Lin Reference Li and Lin2023). The extension technique proposed by Sussman (Reference Sussman2001) is adopted to satisfy the contact angle boundary condition on the particle surface, which is achieved by solving a specific convection equation:
\begin{equation}\frac{\partial C_{\!f}}{\partial \tau }+\boldsymbol{u}^{\boldsymbol{{extend}}}\boldsymbol{\cdot }\boldsymbol{\nabla }C_{\!f}=0.\end{equation}
Herein,
$\tau =0.5\Delta x$
is the artificial time step. The variable
$\boldsymbol{u}^{\boldsymbol{{extend}}}$
is the extension velocity for interface extrapolation, and the specific solution procedure can be found in the literature (Sussman Reference Sussman2001). In this study, the extension amount in the solid region is set to 4 grid sizes.
The numerical method has been systematically validated in our previous work through a series of benchmark problems, including particle sedimentation and particle impact onto a liquid pool (Chen, Yin & Dong Reference Chen, Yin and Dong2024b ). To further assess the accuracy and robustness of the present axisymmetric framework, additional validation cases are conducted, including particle sedimentation, bubble rising, bubble–particle interaction and particle impact onto a liquid pool. These benchmark cases cover a wide range of interfacial and fluid–particle interaction phenomena and demonstrate the capability of the proposed VOSET-IBM-DEM model to accurately capture particle motion, bubble dynamics, interface deformation and momentum exchange in axisymmetric three-phase flows. Detailed descriptions and quantitative comparisons with reference data are provided in the supplementary material available at https://doi.org/10.1017/jfm.2026.11990.
According to the grid resolution requirements of the validation cases, all simulations presented below are performed with a spatial resolution satisfying
$D_{p}/{\unicode[Arial]{x0394}} x$
= 100, and a fixed time step of
${\unicode[Arial]{x0394}} t$
= 2 × 10−7 s is employed. The capillary time scale
\begin{equation}t_{c}=\sqrt{\rho _{l}{R}_{b}^{3}/\sigma }\end{equation}
is adopted to characterise the bubble dynamics, where
$R_{b}$
is the bubble radius.
4. Results and discussion
4.1. Bubble detachment modes
The dynamic behaviour of a particle impacting a bubble under different Weber numbers (
$\textit{We}=\rho V^{2}D_{b}/\sigma$
) is investigated. Unless otherwise stated, the particle wettability was neutral (
$\theta _{p}$
= 90°), the density ratio was fixed at
$k_{\rho }$
= 5.24 and the size ratio at
$k_{r}$
= 2/3. The numerical results reveal three distinct detachment modes following impact: capillary-wave detachment (CWD), bubble-rebound detachment (BRD) and particle-rebound detachment (PRD). Representative snapshots of these modes are presented in figure 2.
Bubble detachment modes induced by particle impact at different Weber numbers: (a) We = 0.1, CWD; (b) We = 1.0, BRD; (c) We = 10, PRD. (d) Schematic illustration of the key geometric parameters of the bubble.
At the early stage of the interaction, the dynamics bears strong resemblance to the coalescence of two droplets or bubbles (Eggers et al. Reference Eggers, Sprittles and Snoeijer2024). Upon contact, a solid–vapour bridge rapidly forms at the point of impact. The bubble subsequently spreads over the particle surface, in a manner analogous to liquid-bridge growth, accompanied by the release of interfacial energy. A portion of this released surface energy is converted into capillary-wave energy, generating interfacial disturbances that propagate downstream along the vapour–liquid interface (Shi, Cui & Zhang Reference Shi2015). Sustained excitation at the perturbed contact line leads to a transition from an initial single-pulse disturbance to a multi-period wave packet characterised by smaller leading waves followed by larger-amplitude trailing waves. This wave-packet structure is broadly consistent with the self-similar capillary-wave patterns reported by Ding et al. (Reference Ding, Li, Zhang, Sui, Spelt and Thoroddsen2012). The detailed morphology is illustrated in figure 3. Point A denotes the primary crest adjacent to the particle contact line, whose evolution directly reflects local stress concentration. Points B and C correspond to the advancing crest and trough propagating towards the bottom surface, respectively. The wave packet exhibits pronounced spatio-temporal characteristics. In regions not yet reached by the capillary-wave front, the bubble remains quasi-static under the constraint of surface tension. In contrast, in regions traversed by the wave packet, the strong oscillations in interfacial curvature provide the energetic basis for subsequent momentum redistribution. The first column of figure 2 shows the transient configuration prior to disturbance of the bottom contact line. The results indicate that, with increasing We, the time required for the capillary wave to reach the bubble base increases slightly.
Structure of capillary waves.
When the capillary wave packet reaches the bottom surface, the three-phase contact line at the substrate is set into motion under wave-induced forcing. Figure 4 presents the temporal evolution of the dimensionless contact radius,
$k_{\alpha }=\alpha /\alpha _{0}$
. The initial geometric parameters of the bubble–particle system are defined in figure 1, and the characteristic length scales during impact are illustrated in figure 2(d). The contact line first undergoes a retraction driven by the leading small-amplitude wave (second column of figure 2), followed by an initial spreading phase (third column), with the maximum spreading radius increasing monotonically with We. Subsequently, a second retraction occurs under the action of the larger trailing wave (fourth column). At We = 0.1, the bubble detaches directly from the substrate at
$t/t_{c}$
= 0.473. In this regime, detachment is governed primarily by the large trailing wave. This wave-driven mechanism is termed CWD. Based on the contact-line dynamics at the bottom surface, the process can be divided into four stages: a quasi-static stage (I), a first retraction stage (II), a first spreading stage (III) and a second retraction stage (IV). Stages II–IV are all driven by capillary-wave forcing.
Temporal evolution of the dimensionless bottom contact radius.
For We = 1.0 and 10, detachment does not occur directly during stage IV. The extent of the second retraction decreases markedly with increasing We, indicating that particle inertia suppresses the capillary-driven contact-line motion, with the degree of suppression positively correlated with impact velocity. For We = 1.0, the particle velocity decreases to zero at
$t/t_{c}$
= 0.749. Owing to prior compression, the bubble has accumulated potential energy, which is subsequently released to induce a rebound motion, carrying the particle upward synchronously. The bubble eventually detaches at
$t/t_{c}$
= 1.565. This mode, driven by the recovery of bubble deformation, is defined as BRD. Compared with CWD, BRD involves two additional stages: a secondary spreading stage (V), resulting from particle-induced suppression, and a third retraction stage (VI), triggered by bubble rebound.
At We = 10, the substantially larger initial kinetic energy of the particle causes it to strike the substrate directly at
$t/t_{c}$
= 0.441, followed by mechanical rebound. During this collision, the bubble undergoes more severe compression. Thereafter, the release of bubble potential energy acts synergistically with particle rebound, jointly driving the coupled particle–bubble system upward and leading to detachment at
$t/t_{c}$
= 1.030. This mode is defined as PRD. Its evolution also comprises six stages. The first five stages are essentially identical to those of BRD, whereas stage VI corresponds to a third retraction induced by the combined rebound of both the bubble and the particle.
A comparison of detachment times among the three modes shows that CWD, involving the fewest dynamical stages, yields the shortest detachment time (0.473), followed by PRD (1.030), while BRD exhibits the longest detachment time (1.565). The difference between BRD and PRD arises from their distinct rebound-response mechanisms (Chen et al. Reference Chen, Liu, Gao and Ding2019). Figure 5 shows the temporal evolution of the dimensionless particle position,
$k_{{h_{p}}}=h_{p}/H_{p}$
, and velocity,
$k_{{v_{p}}}=v_{p}/V$
. The bubble responds analogously to a soft spring: its compression and rebound are gradual, and the feedback to particle motion is delayed. By contrast, particle impact with the substrate is a typical elastic collision characterised by an impulsive response and extremely high acceleration. Consequently, although the contact area at rebound under We = 10 is larger than that at We = 1.0, the much higher initial rebound velocity shortens the time interval between rebound onset and complete detachment.
Temporal evolution of the dimensionless particle position and velocity.
4.2. Surface energy conversion and geometric suppression
From the perspective of the coupled bubble–particle dynamics, the first four stages of detachment are primarily governed by capillary-wave forcing at the bottom contact line. The underlying energy source can be traced to the surface energy released during bubble spreading over the particle surface. In studies of droplet coalescence dynamics (Eggers et al. Reference Eggers, Sprittles and Snoeijer2024), the evolution of the liquid-bridge radius
$r_{l}$
is commonly employed as a key metric for quantifying the rate of energy conversion. By analogy, in the present configuration the local vapour–liquid interfacial variation is predominantly induced by bubble wetting and spreading over the particle.
To quantify this process, we introduce the spreading angle
$\beta$
to characterise the extent of surface energy release. When the spreading angle reaches
$\beta$
, the corresponding reduction in the local vapour–liquid interfacial area,
${\unicode[Arial]{x0394}} S$
, is equivalent to the area of a spherical cap with central angle
$\beta$
:
\begin{equation}{\unicode[Arial]{x0394}} S=2\pi\kern-1.5pt {R}_{p}^{2}(1-\cos \beta ).\end{equation}
The instantaneous release rate of surface energy is given by
\begin{equation}\dot{E_{s}}=\gamma \frac{\text{d}({\unicode[Arial]{x0394}} S)}{\text{d}t},\end{equation}
where
$\gamma$
is the vapour–liquid surface tension.
Temporal evolution of the spreading angle of the bubble on the particle surface.
Figure 6 shows the evolution of the spreading angle
$\beta$
as a function of the dimensionless time
$t/t_{c}$
for different We. During the early spreading stage,
$\beta$
follows a power-law scaling:
\begin{equation}\beta \sim \left(\frac{t}{t_{c}}\right)^{b}.\end{equation}
For We = 0.1, 1 and 10, the exponent b is 0.46, 0.51 and 0.80, respectively. The increase in b reflects the enhanced wetting dynamics induced by particle impact. For We ≤ 1, the system resides in a capillary-dominated spreading regime, with b
$\approx$
0.5, consistent with the inertial-capillary limit (De Maleprade et al. Reference de Maleprade, Clanet and Quéré2016). When We> 1, mechanical kinetic energy directly participates in interfacial evolution, accelerating the surface energy release rate. As a result, higher We leads to a greater amount of surface energy released per unit time, thereby exciting capillary waves of higher energy density.
In addition to surface energy release, particle impact induces global bubble deformation, thereby modifying its potential energy distribution. To assess this effect, the horizontal (
$R_{r}=w_{b}/W_{b}$
) and vertical (
$R_{z}=h_{b}/H_{b}$
) characteristic dimensions of the bubble are introduced, and their temporal evolutions are shown in figure 7. For We = 0.1, spreading-induced axial elongation dominates, leading to an increase in the vertical extent. This elongation provides a geometrical ‘pre-tension‘ that facilitates subsequent contact-line retraction and promotes capillary-wave-driven detachment. By contrast, at We = 10, strong particle impact leads to dominant compressive deformation. The vertical dimension decreases progressively after impact, which limits the energy density available when capillary waves converge towards the axis. Consequently, despite greater total surface energy release, the conversion into effective upward momentum is hindered. Bubble detachment therefore results from a dynamic competition between surface energy release and mechanically induced compressive suppression. Although large We enhances the energetic source term (surface energy release), the severe geometric deformation and negative momentum injection associated with impact offset the capillary-wave driving effect, rendering capillary-wave-induced detachment less likely.
Temporal evolution of the bubble dimensions in the horizontal and vertical directions.
4.3. Momentum evolution during particle impact
To quantify the momentum accumulation of the bubble during particle impact, the analysis is conducted within the framework of Kelvin impulse theory. According to potential-flow theory, the Kelvin impulse
I
is defined as the surface integral of the velocity potential
$\varphi$
over the bubble surface
$S_{b}$
(Best & Blake Reference Best and Blake1994; Wang et al. Reference Wang, Wu, Zheng, Du and Zhang2022):
\begin{equation}\boldsymbol{I}=\oint _{{S_{b}}}\rho _{l}\varphi \boldsymbol{n}\,\text{d}S,\end{equation}
where
$\boldsymbol{n}$
is the unit normal vector directed into the bubble.
By virtue of the divergence theorem, the surface integral over a closed surface can be converted into a volume integral of its divergence. To extract the vertical component
$I_{z}$
, we construct the auxiliary vector field
$\boldsymbol{A}=\varphi \boldsymbol{k}$
, where
$\boldsymbol{k}$
is the unit vector in the vertical direction. Application of the divergence theorem yields
\begin{equation}\oint _{{S_{b}}}\phi \boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{k}\,\text{d}S=\int _{{V_{b}}}\boldsymbol{\nabla }\boldsymbol{\cdot }(\phi \boldsymbol{k})\,\text{d}V.\end{equation}
The velocity field is defined as
$\boldsymbol{u}=\boldsymbol{\nabla }\varphi$
, and hence
$\boldsymbol{\nabla }\boldsymbol{\cdot }(\varphi \boldsymbol{k})=u_{z}$
, the vertical velocity component. Substitution leads to the volume-integral representation of the Kelvin impulse:
\begin{equation}I_{z}=\rho _{l}\int _{{V_{b}}}u_{z}\,\text{d}V,\end{equation}
which can be interpreted as an equivalent vertical momentum associated with the liquid motion within the bubble volume. Its discrete form reads
\begin{equation}I_{z}=\rho _{l}\sum\limits_{i=1}^{N}{(u_{z}})_{i}F_{i}{\unicode[Arial]{x0394}} V_{i},\end{equation}
where
$F_{i}$
denotes the vapour volume fraction and
${\unicode[Arial]{x0394}} V_{i}$
the cell volume.
Figure 8 shows the temporal evolution of
$I_{z}$
for different We. At the initial stage of impact, increasing We enhances the downward momentum injected by the particle, resulting in a pronounced negative peak of
$I_{z}$
. Subsequently, under the action of capillary waves excited by surface energy release, the negative impulse gradually decays in all cases. For We ≤ 1, the initial downward impulse is relatively small, and the capillary-wave-induced positive momentum accumulation rapidly overcomes inertial suppression, leading to momentum reversal (
$I_{z}$
> 0). The star symbol in the figure marks the transition point of the spreading power law, beyond which contact-line spreading approaches saturation and the spreading rate decreases.
Temporal evolution of the vertical Kelvin impulse.
At We = 0.1, although the spreading transition point is crossed, the impulse evolves smoothly owing to the weak impact. Capillary waves are therefore sufficient to directly induce bubble detachment. At We = 1, capillary waves initially promote partial recovery of the impulse; however, substantial energy is consumed in overcoming substrate adhesion, and continued particle compression imposes geometric suppression. As a result, the total impulse turns negative again after approaching zero. Only when the particle velocity decreases to zero does the stored deformation energy of the bubble dominate the rebound, ultimately triggering detachment. In contrast, at We = 10, the particle remains in the spreading power-law regime prior to substrate impact, and the large impact momentum maintains a persistently negative vertical impulse. Upon collision, mechanical rebound of the particle induces an abrupt and strong momentum reversal, causing
$I_{z}$
to increase sharply with a steep slope and thereby driving rapid detachment.
The temporal rate of change of the Kelvin impulse represents the resultant force acting on the bubble (Best & Blake Reference Best and Blake1994):
\begin{equation}F=\frac{\text{d}I_{z}}{\text{d}t}.\end{equation}
As shown in the figure, compared with the BRD case (We = 1.0), the impulse curve in the PRD regime (We = 10) exhibits a markedly steeper slope. This indicates that under the collision-induced rebound mechanism, the system acquires a stronger instantaneous positive upward impulse, thereby substantially enhancing the efficiency of bubble detachment.
To quantify the effectiveness of energy utilisation under different detachment modes, an energy efficiency ratio
$\eta _{e}$
is introduced:
\begin{equation}\eta _{e}=\frac{E_{k,b}}{E_{\textit{available}}}=\frac{E_{k,b}}{{\unicode[Arial]{x0394}} E_{s}+W_{p,\textit{rebound}}},\end{equation}
where
$E_{k,b}$
denotes the effective detachment kinetic energy defined via the Kelvin impulse:
\begin{equation}E_{k,b}=\frac{{I}_{z,\textit{detach}}^{2}}{2M_{\textit{added}}}=\frac{{I}_{z,\textit{detach}}^{2}}{2\rho _{l}V_{b}}.\end{equation}
Herein,
$I_{z,\textit{detach}}$
denotes the vertical component of the Kelvin impulse at the instant of bubble detachment and
$M_{\textit{added}}$
is added mass. To maintain consistency with the volume-integral form of the Kelvin impulse, the added mass is defined as the mass of liquid displaced by the bubble (
$M_{\textit{added}}=\rho _{l}V_{b}$
). This definition is consistent with potential-flow inertia for bodies moving near boundaries and appropriately quantifies the directed detachment energy transferred to the surrounding liquid.
The available input energy
$E_{\textit{available}}$
consists of the surface energy
${\unicode[Arial]{x0394}} E_{s}$
released during spreading and, in the particle-rebound mode, the mechanical rebound energy
$W_{p,\textit{rebound}}$
of the particle. The capillary-wave energy driving contact-line motion primarily originates from the released surface energy during the first stage:
\begin{equation}{\unicode[Arial]{x0394}} E_{s}=\sigma \big[({A_{lv}}-{A_{sl}}\cos \theta _{p})_{b,\textit{initial}}-\left(A_{lv}-A_{sl}\cos \theta _{p}\right)_{b,\textit{min}}\big],\end{equation}
where
$A_{lv}$
and
$A_{sl}$
denote the vapour–liquid and vapour–solid interfacial areas, respectively. In the present cases,
$\theta _{p}$
= 90°, so that the surface energy variation reduces to (Nam, Kim & Shin Reference Nam, Kim and Shin2013)
\begin{equation}{\unicode[Arial]{x0394}} E_{s}=\sigma [\left(A_{lv}\right)_{b,\textit{initial}}-\left(A_{lv}\right)_{b,\textit{min}}].\end{equation}
The particle rebound energy, included only in the PRD mode, characterises the directional mechanical energy injected into the flow through collision:
\begin{equation}W_{p,\textit{rebound}}=\frac{1}{2}m_{p}{v}_{r,\textit{max}}^{2},\end{equation}
where
$v_{r,\textit{max}}$
is the maximum rebound velocity.
Figure 9 presents the energy components and energy efficiency ratio
$\eta _{e}$
for different We. For We ≤ 10, the released surface energy generally increases with We, consistent with the power-law spreading behaviour discussed previously. This observation further confirms that the surface energy variation in stage I mainly originates from bubble spreading over the particle surface. At We = 25, despite the substantial spreading induced by particle impact, the bubble undergoes pronounced compression. As a result, the net variation of the total interfacial area remains small, leading to a reduced release of surface energy. The kinetic energy at detachment strongly depends on the detachment regime. In the CWD regime, increasing We promotes surface energy release, whereas the detachment kinetic energy decreases, indicating a reduced efficiency of energy conversion. This behaviour can be attributed to the enhanced compression of the bubble, which consumes an increasing fraction of the released surface energy.
Energy and energy utilisation efficiency at different Weber numbers.
After the transition to the BRD regime, both the detachment kinetic energy and
$\eta _{e}$
increase with We. Stronger impacts lead to more significant bubble compression, enabling more potential energy to be stored in the deformed interface and released during the rebound process. Consequently, the upward impulse generated at detachment becomes stronger, leading to a higher energy efficiency ratio
$\eta _{e}$
.
Following the transition to the PRD regime, the rebound of the particle injects substantial directional momentum into the flow, causing the detachment kinetic energy to increase with We. In this regime, the rebound energy of the particle becomes the dominant energy contribution, reaching 10.3 and 102.2 times the released surface energy for We = 10 and 25, respectively. Despite the substantial mechanical energy supplied by particle rebound,
$\eta _{e}$
remains relatively low. A large fraction of the rebound energy is dissipated through viscous effects and particle deceleration during the subsequent upward motion, while only a limited portion contributes to bubble detachment. The further reduction of efficiency from We = 10 to 25 indicates that energy dissipation becomes increasingly important at higher We.
Particle wettability, geometric size and density play decisive roles in governing the post-impact spreading, wetting behaviour and deformation evolution of the bubble. In the following sections, these factors are systematically investigated to elucidate their respective influences on the impact dynamics.
4.4. Synergistic effects of impact We and particle wettability
The synergistic influence of the impact Weber number and particle wettability on system dynamics was studied, as discussed in this section. The parameter ranges considered are 0.1 ≤ We ≤ 25 and 30° ≤
$\theta _{p}$
≤ 150°, with density ratio
$k_{\rho }$
= 5.24 and size ratio
$k_{r}$
= 2/3. The resulting detachment-mode map is presented in figure 10(a). For
$\theta _{p}$
≥ 90° and We ≤ 1, CWD dominates, and its regime expands towards higher We as
$\theta _{p}$
increases. In contrast, for hydrophilic particles (
$\theta _{p}$
= 30°), the system remains in the adhered-sedimentation mode (ASM) for We < 10. At intermediate contact angles and moderate We, BRD prevails. Once We exceeds a critical threshold (approximately 10), the system transitions to PRD. Notably, the onset boundary of PRD exhibits weak sensitivity to wettability, indicating that under strong inertial impact, externally injected mechanical work largely overrides interfacial differences induced by
$\theta _{p}$
.
(a) Detachment regime map and (b) detachment time of the bubble at different Weber numbers and particle wettabilities.
The corresponding evolution of the dimensionless detachment time,
$t_{d}/t_{c}$
, is shown in figure 10(b). Under CWD dominance,
$t_{d}/t_{c}$
remains at a relatively low and nearly constant level, displaying little sensitivity to We. However, as
$\theta _{p}$
decreases from hydrophobic to hydrophilic,
$t_{d}/t_{c}$
increases progressively and becomes increasingly sensitive to wettability. This reflects the highly transient nature of capillary-wave-driven detachment, which is not constrained by prolonged inertial suppression cycles. In the BRD regime,
$t_{d}/t_{c}$
increases significantly and exhibits a trend of increasing with We and decreasing with
$\theta _{p}$
. This behaviour correlates closely with the time required for the particle velocity to decay to zero. As We increases, deeper particle penetration extends the energy storage period, while limited contact-line mobility under less hydrophobic conditions prolongs morphological recovery before detachment can occur. In the PRD regime, enhanced mechanical energy compensation reduces the disparity in
$t_{d}/t_{c}$
across wettabilities. Owing to the impulsive nature of particle–substrate collision, the high-amplitude rebound force generated over a short duration disrupts the quasi-static balance of bubble deformation recovery, leading to overall shorter detachment times than in BRD.
(a) Bubble spreading angle and (b) vertical impulse at
$t/t_{c}$
= 0.127 under different Weber numbers and particle wettabilities.
As established previously, the essence of CWD lies in whether the initial capillary-wave packet can overcome substrate adhesion and geometric suppression induced by particle compression. The energy of this wave packet originates primarily from surface energy released during bubble spreading over the particle (Chen et al. Reference Chen, Yin, Zhang, Xin, Ni and Dong2025a
). According to (4.1) and (4.2), the associated energy-conversion rate scales positively with
${R}_{p}^{2}(1-\cos \beta )$
. Since
$R_{p}$
is fixed in the present study, the value of
$(1-\cos \beta )$
at the end of stage I (
$t/t_{c}$
= 0.127) is used to quantify spreading-induced surface energy variation. Figure 11(a) shows the dependence of
$(1-\cos \beta )$
on
$\theta _{p}$
and We. The results indicate positive correlations with both parameters. Scaling analysis reveals that the normalised spreading factor obeys the power law
\begin{equation}(1-\cos \beta )\propto \theta _{p}\textit{We}^{1/10}.\end{equation}
Thus, the normalised
$(1-\cos \beta )$
increases linearly with
$\theta _{p}$
(slope 1/4) and follows a weak power-law dependence on We (exponent 1/10), indicating stronger sensitivity to wettability than to inertia.
To further quantify capillary-wave-driven detachment capability, the vertical impulse
$I_{z}$
at
$t/t_{c}$
= 0.127 is extracted from figure 11(b). The impulse increases markedly with
$\theta _{p}$
but decreases with increasing We. For We ≤ 1,
$I_{z}$
is primarily governed by contact angle and only weakly affected by inertia. For We ≥ 1, however, the impulse reflects a coupled regulation by both We and wettability. At fixed
$\theta _{p}$
,
$I_{z}$
is more sensitive to We than the spreading factor
$(1-\cos \beta )$
. In particular, when We ≥ 1, the reduction in
$I_{z}$
with increasing We becomes substantial, demonstrating that inertial suppression outweighs spreading enhancement. Under high-We, low-
$\theta _{p}$
conditions,
$I_{z}$
becomes negative, indicating that strong inertial compression fully offsets capillary-wave-driven upward motion.
These observations highlight a competition between ‘spreading gain’ (favouring detachment) and ‘inertial suppression’ (hindering detachment) as We increases at fixed
$\theta _{p}$
. To quantify this balance, a bubble detachment capability factor
$\eta$
is introduced, defined at
$t/t_{c}$
= 0.127 as
\begin{equation}\eta =\mathrm{sgn}\left(I_{z}\right)\frac{E_{k,b}}{{\unicode[Arial]{x0394}} E_{s}},\end{equation}
where
$\mathrm{sgn}(I_{z})$
accounts for impulse direction (+1 for
$I_{z}$
≥ 0,
$-$
1 otherwise),
$E_{k,b}$
is the upward kinetic energy of the bubble and
${\unicode[Arial]{x0394}} E_{s}$
is the released surface energy.
Figure 12 illustrates the variation of
$\eta$
under the combined effects of
$\theta _{p}$
and We. The efficiency increases monotonically with
$\theta _{p}$
, reaching particularly high values in the hydrophobic regime at low We. This confirms, from an energetic perspective, that enhanced hydrophobicity promotes conversion of interfacial potential energy into vertical kinetic energy, thereby facilitating detachment. Conversely,
$\eta$
decreases with increasing We across all wettabilities, and the decline becomes more pronounced at higher We. Under extreme high-We, low-
$\theta _{p}$
conditions,
$\eta$
crosses zero and becomes negative, corresponding to the negative impulse distribution in figure 11(b). In this regime, strong impact inertia not only dissipates released surface energy but also drives downward motion, leading to the ASM state. Taken together with the interfacial evolution and momentum distribution shown in figure 11, it is evident that although increasing We enhances surface energy release via stronger spreading (higher
$(1-\cos \beta )$
), the accompanying downward inertial suppression dominates in magnitude, thereby limiting the potential for capillary-wave-driven detachment.
Bubble detachment capability factor
$\eta$
under different We and
$\theta _{p}$
.
In the PRD regime, the primary driving mechanisms for detachment are the mechanical rebound kinetic energy generated upon particle–substrate collision and the elastic potential energy stored in bubble deformation. To quantify the influence of particle wettability on PRD behaviour, the dimensionless maximum rebound velocity,
$k_{vp\_ r}=v_{p\_ r}/V_{0}$
, and the bubble shape factor at rebound,
$k_{br}=R_{br}/R_{bz}$
, are examined. The results are presented in figure 13, which also includes representative snapshots at We = 10 corresponding to the instant of maximum particle rebound velocity for different contact angles. The maximum rebound velocity increases monotonically with
$\theta _{p}$
. A larger contact angle reduces the spreading resistance of the bubble over the particle surface, thereby facilitating particle penetration through the bubble and diminishing momentum dissipation during descent. When
$\theta _{p}$
= 120°, the particle has already fully penetrated the bubble before attaining its maximum rebound velocity.
Maximum rebound velocity of the particle and bubble shape factor at different
$\theta _{p}$
.
Figure 13. Long description
Panel A: A line graph showing the relationship between particle contact angle (theta_p) and two factors, k_vp_r and k_br, for We = 10 and We = 25. The x-axis represents theta_p in degrees, ranging from 30 to 150. The left y-axis represents k_vp_r, ranging from 0.3 to 0.8. The right y-axis represents k_br, ranging from -2 to 8. The graph includes four data lines: solid blue circles for We = 10, k_vp_r; dashed blue squares for We = 25, k_vp_r; solid orange circles for We = 10, k_br; and dashed orange squares for We = 25, k_br. The blue lines show an increasing trend for k_vp_r with theta_p, while the orange lines show a decreasing trend for k_br with theta_p. Embedded images below the x-axis illustrate the bubble shapes at different theta_p values.
In contrast, the bubble size ratio
$k_{br}$
decreases with increasing
$\theta _{p}$
. A smaller
$k_{br}$
corresponds to stronger bubble compression and hence greater storage of elastic potential energy available for rebound (Liu et al. Reference Liu, Ghigliotti, Feng and Chen2014). These two opposing trends – enhanced particle rebound velocity but reduced bubble compression at higher
$\theta _{p}$
– compensate each other, leading to relatively small variations in detachment time within the PRD regime at high We. It is noteworthy that at
$\theta _{p}$
= 150° and We = 10, a slight local reduction in rebound performance is observed. However, when We increases to 25, the intensified mechanical energy input effectively overcomes this limitation, ensuring efficient bubble detachment under strong-impact conditions.
4.5. Synergistic effects of impact We and density ratio
The coupled influence of the impact Weber number (0.1 ≤ We ≤ 10) and the particle–liquid density ratio (1 ≤
$k_{\rho }$
≤ 10) was systematically investigated, as discussed in this section. The particle contact angle is fixed at 90°, and the particle-to-bubble size ratio is
$k_{r}$
= 2/3. The detachment-mode phase diagram is presented in figure 14(a). An increase in density ratio essentially enhances the mass inertia of the system. At a fixed We, particles with higher
$k_{\rho }$
possess greater downward momentum, intensifying interfacial compression and energy dissipation. Consequently, the density ratio markedly shifts the threshold boundaries between detachment modes. As
$k_{\rho }$
increases from 1.0 to 10.0, the control region of CWD contracts towards lower We, while the ASM gradually emerges at high density ratios. At We = 10, PRD remains dominant. Owing to the stronger inertial effect associated with large
$k_{\rho }$
(= 10), the transition to PRD is triggered earlier, occurring already at We = 2.5.
Figure 14(b) shows the variation of the dimensionless detachment time
$t_{d}/t_{c}$
. In the CWD regime, the detachment time is nearly insensitive to density ratio, indicating that under weak-impact conditions, system dynamics remains governed by capillary-wave energy, and particle mass differences do not yet introduce significant geometric suppression. For high-density particles,
$t_{d}/t_{c}$
exhibits a pronounced nonlinear jump, corresponding to the BRD regime. Deeper penetration by heavier particles prolongs the period of energy storage and morphological recovery, forming a clear dynamical bottleneck. Once the system enters the PRD regime,
$t_{d}/t_{c}$
rapidly decreases and converges across different density ratios. This suggests that for high-mass systems, the strong mechanical work supplied by particle rebound becomes the decisive factor in overcoming interfacial adhesion and shortening the detachment cycle.
(a) Detachment regime map and (b) detachment time of the bubble at different Weber numbers and density ratios.
(a) Bubble spreading angle and (b) vertical impulse at
$t/t_{c}$
= 0.127 at various We and
$k_{\rho }$
.
Figure 15(a) presents the variation of local surface energy change at the end of the initial spreading stage (
$t/t_{c}$
= 0.127), characterised by
$(1-\cos \beta )$
, as functions of
$k_{\rho }$
and We. Over the entire We range, the interfacial energy variation increases monotonically with
$k_{\rho }$
, although the growth rate gradually diminishes and becomes nearly saturated for
$k_{\rho }$
> 7.5. Additionally,
$(1-\cos \beta )$
increases significantly with We. Mechanistically, high-density particles not only carry greater initial kinetic energy but also exhibit stronger resistance to deceleration, maintaining high-speed penetration for a longer duration during bubble compression. This sustained impact forces more pronounced radial expansion of the contact line, promoting surface energy release. In contrast, the vertical impulse
$I_{z}$
displays a distinctly negative evolution trend, as shown in figure 15(b). With increasing
$k_{\rho }$
or We,
$I_{z}$
decreases monotonically. The decline rate becomes less sensitive to density ratio at high
$k_{\rho }$
but intensifies with increasing We. Under combined high-
$k_{\rho }$
and high-We conditions,
$I_{z}$
rapidly crosses zero and enters a deeply negative regime. This indicates that inertia-induced downward compression dominates the vertical momentum balance, and the released capillary-wave driving force is insufficient to counteract the downward traction imposed by heavy particles, thereby promoting transitions between detachment modes.
Figure 16 illustrates the evolution of the bubble detachment efficiency factor
$\eta$
. The efficiency decays approximately exponentially with increasing
$k_{\rho }$
, with the most pronounced attenuation occurring at intermediate We. This suggests that in this range, gravitational (mass-related) and impact inertial effects are comparable in magnitude, rendering the system highly sensitive to density variation. At low We (< 0.25), capillary-wave dynamics dominates momentum response, and
$\eta$
remains relatively high across all density ratios. At higher We (≥2.5), impact inertia prevails. Although the synergistic increase of density and inertia enhances surface energy release, the fraction converted into upward kinetic energy decreases sharply and may even become negative. Particularly when We ≥ 2.5 and
$k_{\rho }$
> 5.0,
$\eta$
approaches zero or becomes negative. In this regime, despite severe interfacial deformation, energy cannot be effectively released in the vertical direction, and the system transitions from CWD to other modes.
Bubble detachment capability factor
$\eta$
under different We and
$k_{\rho }$
.
For the PRD regime at high impact inertia (We = 10.0), figure 17 presents the maximum particle rebound velocity
$k_{vp\_ r}$
and the bubble shape factor at the instant of rebound
$k_{b\_ r}$
for different
$k_{\rho }$
. The results show that
$k_{vp\_ r}$
exhibits a non-monotonic dependence on
$k_{\rho }$
, first increasing and then decreasing. At low
$k_{\rho }$
(
$k_{\rho }$
= 2.5), the particle experiences substantial fluid drag and interfacial reaction forces, leading to rapid deceleration during penetration and a prolonged travel time before reaching the substrate. Consequently, the rebound kinetic energy reserve is limited, explaining the relatively long detachment time of the PRD mode at this density ratio. Meanwhile, the slower compression process allows the bubble more time for morphological recovery, resulting in a larger shape factor
$k_{b\_ r}$
at rebound and a more flattened configuration. As
$k_{\rho }$
increases, the particle more effectively overcomes fluid resistance, significantly shortening the time to reach the substrate and preserving higher impact momentum, which translates into a stronger rebound velocity. However, with further increase in
$k_{\rho }$
,
$k_{vp\_ r}$
declines again. This indicates that under extremely strong impacts associated with very high
$k_{\rho }$
, a larger fraction of momentum is dissipated through inelastic losses during interfacial penetration and substrate collision. For
$k_{\rho }$
> 5.24, both
$k_{vp\_ r}$
and
$k_{b\_ r}$
tend to level off. This behaviour suggests that the mechanical rebound capability approaches a dynamical saturation limit in this parameter range. Owing to the extremely rapid impact process, bubble deformation becomes constrained within a restricted dynamical response window and no longer varies significantly with density ratio.
Maximum rebound velocity of the particle and bubble shape factor at different
$k_{\rho }$
.
(a) Detachment regime map and (b) detachment time of the bubble at different Weber numbers and size ratios.
Figure 18. Long description
Panel A: A scatter plot shows the detachment regime map with the horizontal axis labeled as We (Weber number) and the vertical axis labeled as k_r (size ratio). The plot is divided into regions labeled CWD, BRD, PRD, and ASM, each containing different symbols indicating specific data points. Panel B: A line graph displays the detachment time of the bubble at different Weber numbers and size ratios. The horizontal axis is labeled as We (Weber number) and the vertical axis is labeled as t_d/t_c (detachment time). Different lines represent various size ratios (k_r = 0.5, 0.67, 1.0, 1.6, 2.0), with each line showing how the detachment time changes with the Weber number.
4.6. Synergistic effects of impact We and particle size ratio
This section examines the coupled influence of the impact Weber number (0.1 ≤ We ≤ 10) and the particle-to-bubble size ratio (0.5 ≤
$k_{r}$
≤ 2). The particle contact angle is fixed at 90°, and the particle–fluid density ratio is
$k_{\rho }$
= 5.24. The detachment-mode phase diagram is shown in figure 18(a). In the low-inertia regime (We ≤ 1), the system predominantly exhibits CWD. As
$k_{r}$
increases, the CWD regime expands towards lower We, indicating that larger particles promote capillary-wave-driven detachment under weak-impact conditions. Entering the high-inertia regime (We ≥ 2.5), the detachment mode switches significantly with increasing
$k_{r}$
. For small size ratios (
$k_{r}$
< 1), the system tends to enter the BRD mode. When
$k_{r}$
= 1 at sufficiently large We, the system transitions to PRD. However, for larger size ratios (
$k_{r}$
≥ 1.6), the increased particle mass leads to insufficient bubble rebound energy to overcome inertial suppression, and the system readily falls into the ASM, preventing detachment. Figure 18(b) further presents the variation of the dimensionless detachment time
$t_{d}/t_{c}$
with We and
$k_{r}$
. In the CWD regime, detachment time remains relatively short for all size ratios. Once We exceeds the mode-transition threshold,
$t_{d}/t_{c}$
increases sharply. In the BRD regime, the detachment time increases with both
$k_{r}$
and We, reflecting enhanced penetration depth and prolonged energy recovery. Upon transition to the PRD regime,
$t_{d}/t_{c}$
decreases again due to the strong mechanical rebound contribution.
To clarify the mechanism by which
$k_{r}$
regulates detachment-mode transitions, figure 19 extracts the local interfacial area variation
${\unicode[Arial]{x0394}} S$
caused by bubble spreading over the particle and the vertical impulse
$I_{z}$
at the end of the first stage (
$t/t_{c}$
= 0.127). The results show that
${\unicode[Arial]{x0394}} S$
exhibits a non-monotonic dependence on
$k_{r}$
, first increasing and then decreasing. The peak position is strongly influenced by impact inertia: for We ≤ 1.0, the maximum
${\unicode[Arial]{x0394}} S$
occurs near
$k_{r}$
= 1.6; for We > 1.0, the peak shifts forward to around
$k_{r}$
= 1.0. This peak shift reflects the transition from interfacial-dominated to inertia-dominated dynamics. Under weak inertia, larger particle surface area facilitates greater interfacial extension and potential energy accumulation. However, as We increases, enhanced inertial compression suppresses interfacial deformation, making smaller size ratios more favourable for reaching deformation saturation.
(a) Bubble local interfacial area variation and (b) vertical impulse at
$t/t_{c}$
= 0.127 under different Weber numbers and size ratios.
The vertical impulse
$I_{z}$
decreases monotonically with increasing
$k_{r}$
, and the decline rate intensifies with We. At low We, all size ratios maintain positive impulse, ensuring the spontaneous driving capability of the CWD mode. As
$k_{r}$
increases, however, inertia-induced compression associated with greater particle mass becomes dominant, driving
$I_{z}$
towards negative values. This effect is particularly pronounced at We = 5.0, where the steepest decline of
$I_{z}$
with
$k_{r}$
is observed, highlighting the strong acceleration of downward compression for large particles under high inertia.
Figure 20 further presents the evolution of the bubble detachment capability factor
$\eta$
. The factor decreases continuously with increasing
$k_{r}$
and enters the negative regime at high We. From an energy-partition perspective, this explains why large particles tend to trap the system in the ASM due to insufficient upward momentum despite significant interfacial deformation. Overall, the size ratio regulates detachment-mode transitions by modifying the available deformation space and momentum dissipation pathways. When
$k_{r}$
lies within an intermediate range, the system achieves a more favourable balance between deformation energy storage and maintenance of upward vertical momentum, thereby enhancing detachment robustness.
Bubble detachment capability factor
$\eta$
under different We and
$k_{r}$
.
Maximum rebound velocity of the particle and bubble shape factor at different
$k_{r}$
.
Focusing on the post-transition dynamics in the PRD regime, figure 21 extracts the maximum particle rebound velocity
$k_{vp\_ r}$
and the bubble shape factor at the instant of rebound
$k_{br}$
under the condition We = 10.0, in order to evaluate the size effect on rebound-stage characteristics. The results show that as the size ratio
$k_{r}$
increases, the maximum rebound velocity
$k_{vp\_ r}$
decreases markedly, while the bubble shape factor
$k_{br}$
increases correspondingly. This trend indicates that larger particles experience more severe momentum dissipation during penetration through the bubble and subsequent impact with the substrate. As
$k_{r}$
increases, the reduction in particle velocity prior to substrate collision becomes more pronounced, thereby limiting the rebound velocity achieved during the rebound stage. Meanwhile, the increase in
$k_{br}$
with
$k_{r}$
reflects stronger horizontal expansion of the bubble under the compression of larger particles, implying more intense mechanical suppression. Such a strongly compressed state corresponds to greater storage of deformation potential energy. However, under large-size-ratio conditions, the substantial drop in particle rebound velocity
$k_{vp\_ r}$
means that the accumulated mechanical kinetic energy is insufficient to synergistically combine with the stored deformation energy to overcome adhesive resistance. Consequently, when
$k_{r}$
= 2, the system transitions into the ASM. In contrast, within the range
$k_{r}$
< 1.0, particles possess stronger penetration capability and higher rebound velocities, providing sufficient mechanical compensation to drive bubble detachment. This ensures efficient detachment under high-We conditions.
5. Conclusions
An axisymmetric extension of the VOSET-IBM-DEM was developed to investigate bubble detachment from a solid surface induced by particle impact. Three distinct detachment modes were identified: CWD, BRD and PRD. The vertical Kelvin impulse provides a unified framework for interpreting all detachment modes. Bubble detachment is associated with a reversal of the vertical impulse from negative to positive, establishing a net upward momentum. Based on this mechanism, a bubble detachment capability factor
$\eta$
is introduced to characterise the conversion of released energy into upward detachment motion. Mode transitions are governed by the competition between capillary-wave-driven momentum generation and inertia-induced suppression. The CWD regime is generally confined to We
$\lesssim$
1, whereas the onset of PRD occurs near We
$\approx$
10. The influence of particle properties varies among the three detachment modes. The onset boundary of PRD exhibits only weak sensitivity to particle properties, reflecting the dominance of particle collision, whose associated energy exceeds the released surface energy by one to two orders of magnitude. Within the PRD regime, the rebound velocity exhibits a non-monotonic dependence on density ratio, first increasing and then decreasing, while both the rebound velocity and the bubble deformation factor approach saturation at
$k_{\rho } \gt$
5.24. The size ratio
$k_{r}$
associated with the maximum surface energy release shifts from
$k_{r}\approx$
1.6 at We
$\leq$
1 to
$k_{r}\approx$
1.0 at We
$\gt$
1, indicating a transition from capillary-dominated to inertia-influenced deformation dynamics. These findings establish the detachment capability factor
$\eta$
and the associated regime boundaries as useful criteria for characterising particle-induced bubble detachment.
The present study provides a mechanistic framework for understanding and predicting particle-induced bubble detachment in boiling heat transfer and related gas–liquid–solid multi-phase systems. It should be noted that the present simulations are restricted to axisymmetric configurations and neglect phase change at the vapour–liquid interface. Future work will extend the model to fully three-dimensional boiling systems to investigate the coupled effects of particle motion, bubble growth and phase-change heat transfer under realistic operating conditions.
Supplementary material
Supplementary material is available at https://doi.org/10.1017/jfm.2026.11990.
Acknowledgements
X.C. gratefully acknowledges Professor S. Takagi of the University of Tokyo for valuable discussions and insightful suggestions during the revision of the manuscript.
Declaration of interests
The authors report no conflict of interest.
t/tc
η
θp
θp
t/tc
kρ
η
kρ
kρ
t/tc
η
kr
kr
