1. Introduction
The transport and dynamics of rigid or deformable particles and their interaction with surfaces are of considerable importance in various scientific and industrial fields due to their potential for efficient, and sustainable applications. Understanding these interactions is crucial for optimising processes ranging from wastewater treatment to mineral beneficiation and food processing. In many engineering applications, a deformable surface, or interface separating two fluids, affects the motion of rigid or deformable bodies. Furthermore, in oceanography, particles and marine snow settle in a background fluid that consists of large density gradients caused by variations in salinity and temperature levels, which may greatly affect their fate. Bubbles rising in the ocean may also interact with sedimenting particles, which may play a role in their transport (Masry et al. Reference Masry, Rossignol, Roussel, Bourgogne, Bussière, R'mili and Wong-Wah-Chung2021). Comprehensive reviews have been carried out on the motion in viscosity and/or density stratified fluids by Govindarajan & Sahu (Reference Govindarajan and Sahu2014), Ardekani, Doostmohammadi & Desai (Reference Ardekani, Doostmohammadi and Desai2017), Magnaudet & Mercier (Reference Magnaudet and Mercier2020) and More & Ardekani (Reference More and Ardekani2023).
At the core of these interactions is the dynamic behaviour of bubbles. When a bubble impacts and interacts with a horizontal (Klaseboer et al. Reference Klaseboer, Manica, Hendrix, Ohl and Chan2014; Manica, Klaseboer & Chan Reference Manica, Klaseboer and Chan2015, Reference Manica, Klaseboer and Chan2016) or curved (Basařová, Zawala & Zedníková Reference Basařová, Zawala and Zedníková2019; Esmaili et al. Reference Esmaili, Shukla, Eifert and Jung2019) solid surface, a hydrodynamic dimple (Ascoli, Dandy & Leal Reference Ascoli, Dandy and Leal1990) might form on the surface of the bubble closest to the solid boundary. If the approach velocity of the interaction is sufficiently high, the bubble may rebound, or bounce, after colliding with the solid surface (Zawala et al. Reference Zawala, Krasowska, Dabros and Malysa2007; Chan, Klaseboer & Manica Reference Chan, Klaseboer and Manica2011; Zawala & Dabros Reference Zawala and Dabros2013). If the approach velocity is low, a trapped liquid film in the gap region between the bubble and the solid surface will start to drain. Klaseboer et al. (Reference Klaseboer, Chevaillier, Gourdon and Masbernat2000) and Chan et al. (Reference Chan, Klaseboer and Manica2011) have studied the film drainage and coalescence of bubbles and drops, which present further complexities. The deformation of these bodies can be characterised across six orders of magnitude in length scales, which presents difficulties when studying the dynamic nature of the interaction both experimentally and numerically. Moreover, as a bubble impacts a solid wall, it exerts forces that can dislodge contaminants or biofilms (Gómez-Suárez, Busscher & van der Mei Reference Gómez-Suárez, Busscher and van der Mei2001; Sharma et al. Reference Sharma, Gibcus, van der Mei and Busscher2005a,Reference Sharma, Gibcus, Van Der Mei and Busscherb; Parini & Pitt Reference Parini and Pitt2006; Esmaili et al. Reference Esmaili, Shukla, Eifert and Jung2019). Menesses et al. (Reference Menesses, Belden, Dickenson and Bird2017) have shown that a continuous stream of single bubbles can effectively prevent biofouling by generating consistent shear stress. This principle extends beyond cleaning, influencing the behaviour of bubbles in different media and applications. The adhesion behaviour of bubbles or droplets with solids is influenced by the properties of the deformable interface and the surrounding fluid (Danov et al. Reference Danov, Stanimirova, Kralchevsky, Marinova, Stoyanov, Blijdenstein, Cox and Pelan2016), the presence of contaminants (Basarova, Soušková & Zawala Reference Basarova, Soušková and Zawala2018; Legawiec et al. Reference Legawiec, Kruszelnicki, Zawadzka, Basařová, Zawala and Polowczyk2023) and the surface characteristics of the solid, such as its geometry or wettability (Krasowska, Zawala & Malysa Reference Krasowska, Zawala and Malysa2009). As a bubble adheres to a solid surface, it forms a three-phase contact line, which is the interface where the bubble, solid and liquid meet. The dynamics and stability of the contact line are crucial in understanding how bubbles interact with solid surfaces.
In the context of mineral beneficiation, froth flotation relies on bubbles to separate valuable minerals from waste material. The process of flotation depends heavily on the characteristics of the separated particles, and the bubbles. Specifically, the surface wettability of the solid particles is crucial for allowing the bubbles to selectively separate hydrophobic from hydrophilic particles (Xie et al. Reference Xie2021), and the bubbles must be appropriately sized to enhance the probability of collision between the bubbles and the particles (Yoon & Luttrell Reference Yoon and Luttrell1989). The main causes behind the detachment of particles from bubbles during flotation separation systems were summarised by Klassen & Mokrousov (Reference Klassen and Mokrousov1963) and Wang et al. (Reference Wang, Wang, Yu, Yao, Zhao, Liu, Chen and Wang2023). These include the sliding of the particle on the bubble surface, the overall hydrodynamic behaviour of the bubble, such as its lateral migration or deformation, or the collision of the aggregate with a solid wall. When a bubble–particle aggregate collides with a solid horizontal wall, Wang et al. (Reference Wang, Wang, Yu, Yao, Zhao, Liu, Chen and Wang2023) have found that the detachment of the particle is dependent on the inclination angle of the major axis of the bubble with the solid wall during the collision. Comprehensive reviews of recent advancements in understanding the hydrodynamics and surface interactions in froth flotation systems can be found in Wang & Liu (Reference Wang and Liu2021) and Xie et al. (Reference Xie2021).
In liquid–solid suspensions, the behaviour of bubbles is equally important. Bubbles interact with the solid particles, affecting their rise velocity and overall dynamics. Hooshyar et al. (Reference Hooshyar, Van Ommen, Hamersma, Sundaresan and Mudde2013) have considered the interaction between a rising bubble and a neutrally buoyant particle suspension. The particles were characterised based on their diameter, and when the size of the particles is sufficiently small, the bubble rises without collision. The work has shown that as the rising bubble interacts with larger particles, the bubble's rise velocity decreases, its surface deforms, and it imparts momentum on the surrounding particles. The bubble will then continue to rise until it interacts with another particle, repeating the cycle of the interaction.
Many experimental (Clift, Grace & Weber Reference Clift, Grace and Weber1978; Bhaga & Weber Reference Bhaga and Weber1981) and computational (Sussman & Puckett Reference Sussman and Puckett2000; Dijkhuizen, van Sint Annaland & Kuipers Reference Dijkhuizen, van Sint Annaland and Kuipers2010) studies have considered the hydrodynamic behaviour of particles or bubbles settling or rising in a quiescent background fluid. At small Reynolds number, 
$Re=D_p U_t\rho _l/\mu _l$, the particle terminal speed is 
$U_t=\frac {2}{9}(\rho _p-\rho _l) gR^2/\mu _l$ where 
$\rho _p$ is the particle density, 
$R=D_p/2$ is the particle (bubble) radius, 
$g$ is the gravitational acceleration, and 
$\rho _l$ and 
$\mu _l$ are the fluid density and viscosity, respectively. At larger 
$Re$, (Mordant & Pinton Reference Mordant and Pinton2000; ten Cate et al. Reference ten Cate, Nieuwstad, Derksen and Van den Akker2002) considered the inertial regime, where 
$1.5< Re<7700$, such that 
$U_t$ is unknown a priori. An alternative way of characterising the particle motion in a homogeneous fluid is through the Archimedes number, 
$Ar=\rho _l(\zeta _p g)^{1/2}R^{3/2}/\mu _l$, where 
$\zeta _p=\rho _p-\rho _l/\rho _l$, which corresponds to a particle Reynolds number based on a gravitational velocity scale (Magnaudet & Mercier Reference Magnaudet and Mercier2020). Tripathi, Sahu & Govindarajan (Reference Tripathi, Sahu and Govindarajan2015) recently conducted three-dimensional direct numerical simulations of a rising bubble in a quiescent, Newtonian fluid. The study found that the rising bubble has five distinct regimes, governed by two non-dimensional parameters, which are the bubble Galilei, 
$Ga=\rho _l g^{1/2} R^{3/2} /\mu _l$, and Bond, 
$Bo=\rho _lgR^2/\gamma$, numbers, where 
$\gamma$ is the surface tension. The computational study was in wide agreement with the experimental work of Bhaga & Weber (Reference Bhaga and Weber1981), and the five distinct regimes are (i) axisymmetric bubble geometry at small 
$Ga$ and small 
$Bo$, (ii) skirted spherical cap bubble at small 
$Ga$ and large 
$Bo$, (iii) a spiralling bubble at large 
$Ga$ and small 
$Bo$, and (iv) peripheral or (v) central breakup at large 
$Ga$ and 
$Bo$. The bubble dynamics can be further classified based on surface mobility; mobile bubbles with low contamination, by the presence of surfactants for example, achieve higher terminal velocities, while immobile, contaminated bubbles rise more slowly due to the absence of tangential velocity (Esmaili et al. Reference Esmaili, Shukla, Eifert and Jung2019).
The aim of this paper is to investigate the interaction between a freely rising, deformable bubble and a freely settling spherical particle in a quiescent, Newtonian fluid. This study examines the impact of several parameters on the interaction behaviour, including the Galilei number, the Bond number, the density contrast between the solid and the fluid 
$(\zeta _p)$, the ratio of the particle radius to the bubble radius 
$(\varOmega )$ and the surface wettability characteristics of the particle. The paper implements a three-dimensional particle-resolved direct numerical simulation using the level contour reconstruction method, which is a hybrid level-set/front-tracking method to accurately capture the motion and the interaction of the bubble and the particle. The outline of the paper is as follows: § 2 contains the methodology and the numerical technique; § 3 provides a discussion of the results; and § 4 contains the conclusions and future directions.
2. Methodology
This section contains the problem formulation, governing equations and the numerical method used in this paper. We perform highly resolved numerical simulations of the incompressible Navier–Stokes equations in a three-dimensional cubical Cartesian domain with 
$L_x=L_y=L_z=16R$, as shown in figure 1. The computational domain is divided into parallel subdomains, which are each uniformly discretised by a 
$64^3$ finite-difference mesh. The equations governing Newtonian and incompressible two-phase flows correspond to mass and momentum conservation, respectively, expressed in the present work using a single-field formulation,
\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}=0, \end{gather}
\begin{gather}\rho\left(\frac{\partial \boldsymbol{u}}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}\right)={-} \boldsymbol{\nabla}{p}+\rho\boldsymbol{g}+\boldsymbol{\nabla}\boldsymbol{\cdot}\mu (\boldsymbol{\nabla}\boldsymbol{u}+\boldsymbol{\nabla}\boldsymbol{u}^T)+ \boldsymbol{F}+\boldsymbol{F}_{FSI}, \end{gather}
where 
$\boldsymbol {u}$, 
$\rho$, 
$p$, 
$\mu$, 
$\boldsymbol {g}$, 
$\boldsymbol {F}$ and 
$\boldsymbol {F}_{FSI}$ denote the velocity, density, pressure, viscosity, gravitational acceleration, the local surface tension force at the interface, and the direct forcing term for fluid–structure interaction, respectively. Here, 
$\boldsymbol {F}_{FSI}$ can be computed based on the work of Fadlun et al. (Reference Fadlun, Verzicco, Orlandi and Mohd-Yusof2000) as
\begin{equation} \boldsymbol{F}_{FSI}=\rho \left( \frac{\boldsymbol{u}_s-\boldsymbol{u}}{\Delta t} \right), \end{equation}
where 
$\boldsymbol {u}_s$ is the velocity of the immersed solid boundary. Direct forcing can be considered as imposing the particle velocity directly on the boundary which is equivalent to applying 
$\boldsymbol {F}_{FSI}$ inside the solid structure (Pathak & Raessi Reference Pathak and Raessi2016). (Feedback forcing (Goldstein, Handler & Sirovich Reference Goldstein, Handler and Sirovich1993) to obtain the necessary force to satisfy the desired velocity of the particle is an alternative approach, however, direct forcing is found to be faster and more efficient (Yang & Stern Reference Yang and Stern2012; Pathak & Raessi Reference Pathak and Raessi2016).)
Figure 1. (a) The three-dimensional Cartesian cubic domain with a size of 
$16R \times 16R \times 16R$, showing the subdomain decomposition and the initial problem set-up; (b) the initial separation distance between the particle and the bubble is shown.
We employ a hybrid formulation for the surface tension force, 
$\boldsymbol {F}$ (Shin et al. Reference Shin, Abdel-Khalik, Daru and Juric2005; Shin, Chergui & Juric Reference Shin, Chergui and Juric2017):
\begin{equation} \boldsymbol{F}=\gamma\kappa_{H_{f}}\mathcal{H}_f.\end{equation}
Here 
$\gamma$ is the surface tension, considered constant, 
$\mathcal {H}_f$ is the Heaviside function, and 
$\kappa _{H_{f}}$ is twice the mean interface curvature calculated on the Eulerian grid using the following equations:
\begin{gather} \kappa_{H_{f}}= \frac{\boldsymbol{F}_L\boldsymbol{\cdot}\boldsymbol{N}}{\boldsymbol{N} \boldsymbol{\cdot}\boldsymbol{N}}, \end{gather}
\begin{gather}\boldsymbol{F}_L=\int_{\varGamma(t)}\kappa_f\boldsymbol{n}_f \delta_f(\boldsymbol{x}-\boldsymbol{x}_f)\, {\rm d} S, \end{gather}
\begin{gather}\boldsymbol{N}=\int_{\varGamma(t)}\boldsymbol{n}_f\delta_f (\boldsymbol{x}-\boldsymbol{x}_f)\,{\rm d} S. \end{gather}
In (2.6) and (2.7), 
$\boldsymbol {x}_f$ is a parameterisation of the interface, 
$\varGamma (t)$, and 
$\delta _f(\boldsymbol {x}-\boldsymbol {x}_f)$ is a three-dimensional Dirac delta function that is non-zero only at the interface, 
$\boldsymbol {x}=\boldsymbol {x}_f$. Here 
$\boldsymbol {n}_f$ is the unit normal vector to the interface, where the interface has an area element, 
${\rm d} S$; 
$\kappa _f$ is again twice the mean interface curvature, however, it is now obtained from the Lagrangian interface structure (Shin et al. Reference Shin, Chergui and Juric2017; Kahouadji et al. Reference Kahouadji, Nowak, Kovalchuk, Chergui, Juric, Shin, Simmons, Craster and Matar2018). The geometric information in 
$\boldsymbol {N}$, such as the unit normal 
$\boldsymbol {n}_f$ and the interface infinitesimal area 
${\rm d} S$, are directly computed from the Lagrangian interface, then distributed onto the Eulerian grid using a discrete delta function and the immersed boundary method developed by Peskin (Reference Peskin1977). The Lagrangian elements of the interface are advected by integrating the following equation using a second-order Runge–Kutta method:
\begin{equation} \frac{{\rm d} \boldsymbol{x}_f}{{\rm d} t}=\boldsymbol{V};\end{equation}
here, the interface velocity 
$\boldsymbol {V}$ is interpolated from the Eulerian velocity. Further details on calculating the force, constructing the function field 
$\boldsymbol {N}$ and the Heaviside function, and advecting the interface can be found in Shin et al. (Reference Shin, Abdel-Khalik, Daru and Juric2005), Shin & Juric (Reference Shin and Juric2007), Shin (Reference Shin2007), Shin & Juric (Reference Shin and Juric2009a), Shin & Juric (Reference Shin and Juric2009b) and Shin, Yoon & Juric (Reference Shin, Yoon and Juric2011).
The particle's motion and interaction with the fluid are modelled using a fictitious domain method, such that the Navier–Stokes equations are applied to the entire flow field. Tracking the solid body's motion is achieved by introducing a distance function 
$\phi _s$. At the start of the simulation, 
$\phi ^0_{s}$ is calculated at the cell centre of the Eulerian grid 
$(x^0_{i},y^0_{i},z^0_{i})$ to create the solid volume, which is then updated to calculate 
$\phi ^n_{s}$ at 
$(x^n_{i},y^n_{i},z^n_{i})$ by tracking the solid's centre position and rotation over time. The solid object's movement is updated by its momentum-averaged translational and rotational velocities, which are calculated in the following:
\begin{gather} m_p \boldsymbol{u}_p = \int_{\mathcal{V}} \rho \boldsymbol{u} \, {\rm d} \mathcal{V}, \end{gather}
\begin{gather} \boldsymbol{\mathsf{I}}_p {\boldsymbol{\omega}_p}=\int_{\mathcal{V}} \boldsymbol{r}\times \rho \boldsymbol{u} \, {\rm d} \mathcal{V}, \end{gather}
where 
$m_p$ is the mass of the solid, 
$\boldsymbol {u}_p$ is its translational velocity, 
$\mathcal {V}$ is the solid volume, 
$\boldsymbol {r}$ is the radial vector from it's centre, 
$\boldsymbol{\mathsf{I}}_p$ is the moment of inertia and 
$\boldsymbol {\omega _p}$ is the rotational velocity. The moment of inertia can be calculated using
\begin{equation} \boldsymbol{\mathsf{I}}_p =\int_{\mathcal{V}} \rho_p |(\boldsymbol{r}\boldsymbol{\cdot} \boldsymbol{r})\boldsymbol{\mathsf{I}}-\boldsymbol{r} \otimes \boldsymbol{r}| \, {\rm d} \mathcal{V}, \end{equation}and the position of the solid can be tracked by the following velocity field:
\begin{equation} \boldsymbol{u}_s (x,y,z)=\boldsymbol{u}_p + \boldsymbol{\omega}_p \times \boldsymbol{r}. \end{equation}
Utilising equation (2.12), the updated solid distance function can be calculated at the current position 
$(x^n_{i},y^n_{i},z^n_{i})$ by tracing back to the original position 
$(x^0_{i},y^0_{i},z^0_{i})$ using
\begin{equation} \phi_s(x^n,y^n,z^n)=\phi_s(x^n-{u}_s\Delta t,y^n-{v}_s\Delta t,z^n-{w}_s\Delta t)=\phi_s(x^0,y^0,z^0). \end{equation}Rigid body constraints are applied to the solid volume for the momentum-averaged translational and rotational velocities, and an additional high viscosity, i.e. 500 times that of water, coefficient is utilised for the solid volume.
Within the single-fluid formulation, the material properties in the domain are calculated as follows:
\begin{gather} \rho(\boldsymbol{x},t)=[\rho_b+(\rho_l-\rho_b)\mathcal{H}_f(\boldsymbol{x},t)] \mathcal{H}_s(\boldsymbol{x},t)+\rho_p(1-\mathcal{H}_s(\boldsymbol{x},t)), \end{gather}
\begin{gather}\mu(\boldsymbol{x},t)=[\mu_b+(\mu_l-\mu_b)\mathcal{H}_f(\boldsymbol{x},t)] \mathcal{H}_s(\boldsymbol{x},t)+\mu_p(1-\mathcal{H}_s(\boldsymbol{x},t)), \end{gather}
where the subscripts 
$(l,b,p)$ denote the background liquid, the bubble and the solid phase, respectively. A smoothed three-dimensional Heaviside function 
$\mathcal {H}_{f}(\boldsymbol {x},t)$ for the fluid phases is defined as zero in the bubble phase, and unity in the background liquid phase. A similar approach has been used for defining the properties of the solid particle, 
$\mathcal {H}_s(\boldsymbol {x},t)$ which is considered as zero in the solid phase and unity for both fluid phases (the bubble and the surrounding liquid). The sharp transition between the two phases is resolved numerically with a steep, yet smooth transition across three to four grid cells (Shin et al. Reference Shin, Chergui and Juric2017; Kahouadji et al. Reference Kahouadji, Nowak, Kovalchuk, Chergui, Juric, Shin, Simmons, Craster and Matar2018).
To non-dimensionalise the governing equations, we introduce the following scaling:
\begin{equation} \left.\begin{array}{c}
\displaystyle (x,y,z)= R(\tilde{x},\tilde{y},\tilde{z}),
\quad \boldsymbol{u}=\sqrt{gR} \tilde{\boldsymbol{u}},\quad
t=\left(\dfrac{R}{\sqrt{gR}}\right)\tilde{t},\\ p= (\rho_l
gR) \tilde{p},\quad \mu=\mu_l\tilde{\mu},\quad
\rho=\rho_l\tilde{\rho}, \end{array}\right\}
\end{equation}
where the tildes, introduced temporarily, designate dimensionless variables. After dropping the tildes, the non-dimensional governing equations read
\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}=0, \end{gather}
\begin{gather}\rho\left(\frac{\partial \boldsymbol{u}}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}\right)={-}\boldsymbol{\nabla}{p}-\rho\boldsymbol{e}_z+\frac{1}{Ga} \boldsymbol{\nabla}\boldsymbol{\cdot}[\mu(\boldsymbol{\nabla}\boldsymbol{u}+\boldsymbol{\nabla}\boldsymbol{u}^T)] +\frac{\boldsymbol{F}}{Bo}+\boldsymbol{F}_{FSI}, \end{gather}
where 
$\boldsymbol {e}_z$ denotes the unit vector in the 
$z$-direction (see figure 1), 
$Ga$ is the Galilei number and 
$Bo$ is the Bond number. The non-dimensional material properties then read
\begin{gather} \rho(\boldsymbol{x},t)=\left[\frac{\rho_b}{\rho_l}+ \left(1-\frac{\rho_b}{\rho_l}\right)\mathcal{H}_f (\boldsymbol{x},t)\right]\mathcal{H}_s(\boldsymbol{x},t)+ \frac{\rho_p}{\rho_l}(1-\mathcal{H}_s(\boldsymbol{x},t)), \end{gather}
\begin{gather}\mu(\boldsymbol{x},t)=\left[\frac{\mu_b}{\mu_l}+\left(1- \frac{\mu_b}{\mu_l}\right)\mathcal{H}_f(\boldsymbol{x},t)\right] \mathcal{H}_s(\boldsymbol{x},t)+\frac{\mu_p}{\mu_l}(1- \mathcal{H}_s(\boldsymbol{x},t)). \end{gather} To mitigate the stress singularity at the contact line (Sui, Ding & Spelt Reference Sui, Ding and Spelt2014), the generalised Navier boundary condition is employed in the framework of the level contour reconstruction method to model the contact line motion (Yamamoto et al. Reference Yamamoto, Ito, Wakimoto and Katoh2013; Shin, Chergui & Juric Reference Shin, Chergui and Juric2018), which may be crucial for the interaction between bubbles and particles over a certain range of parameters as will be discussed in the following section. Derived from molecular dynamics simulations (Qian, Wang & Sheng Reference Qian, Wang and Sheng2003), the boundary condition indicates that the interaction force in a thin layer adjacent to the solid wall includes contributions from both viscous shear stress and an uncompensated Young's stress. The boundary condition thereby combines Navier slip, which relates the contact line displacement to the shear stress at the wall, with the uncompensated Young's stress. This approach models microscale slip information at the macroscale, addressing the infinite shear stress problem near the contact line observed under no-slip conditions. The non-dimensional slip velocity 
$u_{cl}$ is defined as follows:
\begin{equation} u_{cl}= \lambda_s \left( \left.\frac{\partial u}{\partial n}\right|_{wall} + \frac{\cos \theta_{{cont}} - \cos \theta_{{ext}} } {Ca \Delta x}\right), \end{equation}
where 
$\lambda _{s}$ is the slip length, set at a quarter of the size of a grid cell in this work, 
$Ca$ is the capillary number defined as 
$Ca=\mu _l \sqrt {gR}/ \gamma$, 
$\theta _{{cont}}$ is the apparent contact angle governed by the tangential component of the interface element 
$\boldsymbol {t}_{{cont}}$, and 
$\theta _{{ext}}$ is the extended contact angle (see Appendix B) for an extended interface inside the wall modelled by tracing a new vector 
$\boldsymbol {t}_{{ext}}$ in the extended surface. The extended interface contact angle 
$\theta _{{ext}}$ is used to account for the contact angle hysteresis, such that the hysteresis model is as follows:
\begin{equation} \theta_{{ext}} = \begin{cases} \theta_{{ext}} = \theta_A & \text{if } \theta_{{cont}} > \theta_A, \\ \theta_{{ext}} = \theta_R & \text{if } \theta_{{cont}} < \theta_R, \\ \theta_R < \theta_{{ext}} < \theta_A & \text{otherwise}, \end{cases} \end{equation}
where 
$\theta _A$ and 
$\theta _R$ are the advancing and receding contact angles, respectively.
The numerical method involves temporally discretising equation (2.2) into the following form:
\begin{equation} \frac{\boldsymbol{u}^{n+1}-\boldsymbol{u}^n}{\Delta t}={-}(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u})^n+ \frac{1}{\rho^n}[\boldsymbol{G}^n-\boldsymbol{\nabla} p], \end{equation}
where 
$\boldsymbol {G}$ contains the gravitational, surface tension, fluid–structure interaction and viscous forces. The momentum solver computes the velocity and pressure variables on a fixed and uniform Eulerian mesh by means of Chorin's projection method (Chorin Reference Chorin1968). For spatial discretisation, the well-known staggered mesh, MAC (marker and cell) method, is used (Harlow & Welch Reference Harlow and Welch1965). The nonlinear term is spatially discretised using a second-order essentially non-oscillatory scheme (Shu & Osher Reference Shu and Osher1989; Sussman et al. Reference Sussman, Fatemi, Smereka and Osher1998), whereas the other terms are spatially discretised using standard second-order centred differences.
The time integration of (2.23) is split into two substeps. An intermediate unprojected velocity, 
$\boldsymbol {u}^*$, is first calculated neglecting the pressure gradient,
\begin{equation} \frac{\boldsymbol{u}^*-\boldsymbol{u}^n}{\Delta t}={-}(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u})^n+ \frac{\boldsymbol{G}^n}{\rho^n},\end{equation}
followed by the calculation of the final velocity, 
$\boldsymbol {u}^{n+1}$,
\begin{equation} \frac{\boldsymbol{u}^{n+1}-\boldsymbol{u}^*}{\Delta t} ={-}\frac{1}{\rho^n}\boldsymbol{\nabla} p,\end{equation}
where 
$\Delta t$ is the time step. By enforcing the divergence-free condition on 
$\boldsymbol {u}^{n+1}$, the elliptic pressure Poisson equation,
\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \left(\frac{1}{\rho^n} \boldsymbol{\nabla} p\right)=\frac{\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}^*}{\Delta t}, \end{equation}
is solved using a multigrid iterative method whence 
$\boldsymbol {u}^{n+1}$ is obtained:
\begin{equation} \boldsymbol{u}^{n+1}=\boldsymbol{u}^*-\frac{\Delta t}{\rho^n}\boldsymbol{\nabla} p.\end{equation} The temporal integration scheme for all the simulations performed is based on a second-order Gear method (Tucker Reference Tucker2013), with implicit solution of the viscous terms of the velocity components (Kahouadji et al. Reference Kahouadji, Nowak, Kovalchuk, Chergui, Juric, Shin, Simmons, Craster and Matar2018). The time step, 
$\Delta t$, at each temporal iteration is set to satisfy the following criterion:
\begin{equation} \Delta t=\min \{\Delta t_{cap},\Delta t_{CFL},\Delta t_{int},\Delta t_{vis}\},\end{equation}
where 
$\Delta t_{cap},\Delta t_{CFL},\Delta t_{int},\Delta t_{vis}$ are the capillary, Courant–Friedrichs–Lewy (CFL), interfacial and viscous time steps, respectively. These time steps are defined as follows:
\begin{equation} \left.\begin{array}{c} \displaystyle{\Delta t_{cap} \equiv \dfrac{1}{2}\sqrt{\dfrac{(\rho_{b} + \rho_{l}) \Delta x_{min}^3}{{\rm \pi} \sigma}},} \quad \displaystyle{\Delta t_{CFL} \equiv \min_j \left( \min_{domain} \left(\dfrac{\Delta x_j}{u_j}\right) \right),} \\ \displaystyle{\Delta t_{int} \equiv \min_j \left( \min_{\varGamma(t)} \left( \dfrac{\Delta x_j}{\|{\boldsymbol V}\|}\right) \right),} \quad \displaystyle{\Delta t_{ vis} \equiv \min\left(\dfrac{\rho_{b}}{\mu_{b}},\dfrac{\rho_{l}}{\mu_{l}}\right)\dfrac{\Delta x_{min}^2}{6},} \end{array}\right\} \end{equation}
where 
$\Delta x_{min}=\min _j (\Delta x_j)$ is the minimum size 
$x$ for cell 
$j$. The mesh size used in this study is 
$512^3$, which corresponds to a resolution of 
$R/\Delta x=32$. No-penetration boundary conditions are imposed on the lateral domain walls, and no-slip, no-penetration boundary conditions are imposed on the top and bottom boundaries.
In this paper, the viscosity ratio 
$\lambda =\mu _l/\mu _b=100$, and the density ratio 
$\beta =\rho _l/\rho _b=1000$ between the bubble and the background fluid are kept constant. The bubble and the particle radii are kept constant 
$\varOmega =R_b/R_p=1$ unless otherwise stated, with an initial distance of 
$6R$ from the centre of the bubble to the centre of the particle. The non-dimensional parameters governing the problem are 
$Ar$, 
$Bo$ and 
$Ga$. The parameter 
$Ar$, governed by 
$\zeta _p$, is used instead of 
$Re$ as the particle terminal velocity is unknown a priori. Unless otherwise stated, the particle Archimedes number 
$(Ar=10)$ is kept constant by setting 
$\zeta _p=1$. This study explores the bubble parameter range of 
$Bo=[0.5\unicode{x2013}5]$ and 
$Ga=[10, 20]$; 
$Ga$ is kept relatively low to explore the regime wherein the bubble follows a rectilinear path and interacts with the settling particle.
The numerical method utilised in this study has been extensively validated with experimental works found in the literature. In figure 2, two different validation results are provided. Figure 2(a,b) show a bubble rising from rest in a quiescent Newtonian fluid, approaching a terminal velocity and interacting with a solid, immobile and impermeable wall. The results are compared with the experimental work of Kosior et al. (Reference Kosior, Zawala and Malysa2014). The bubble terminal velocity found using this numerical method is approximately 
$350\ \textrm {mm}\ \textrm {s}^{-1}$, which is in excellent agreement with the terminal velocity found in the experiments. Figure 2(a) shows the temporal evolution of the bubble horizontal over vertical deformation, superposed with the results found in the experiments; figure 2(b) shows a comparison between the bubble shapes obtained from the experiments (top-half) and our numerical method (bottom-half); figure 2(c) shows the results of the numerical framework of a rising deformable bubble in a quiescent, Newtonian fluid compared with the experimental results of Bhaga & Weber (Reference Bhaga and Weber1981). When comparing the terminal shape of the bubble predicted by the simulations with that observed in the experiments, excellent agreement is found. Furthermore, the numerical framework was validated against the experimental works on settling solid particles (Mordant & Pinton Reference Mordant and Pinton2000; ten Cate et al. Reference ten Cate, Nieuwstad, Derksen and Van den Akker2002), as seen in figure 3. The particle velocity profiles obtained from our numerical method are in agreement with the results found in the experiments, which ultimately inspires confidence in the accuracy and reliability of our computations. Additional validation case studies are provided in Appendix B.
Figure 2. Validation of the numerical technique: (a,b) show a comparison of our predictions with the experimental data of Kosior, Zawala & Malysa (Reference Kosior, Zawala and Malysa2014) for a freely rising bubble interacting with a solid wall; (a) compares the bubble deformation as a ratio of the horizontal over the vertical bubble axis, and (b) depicts the spatiotemporal evolution of the bubble–wall collision; (c i–iii) compare the results of a rising bubble shape with the experimental results of Bhaga & Weber (Reference Bhaga and Weber1981). The top subpanel shows the three-dimensional illustrations of the terminal bubble shapes, and the bottom subpanel illustrates a slice of the interface superposed on the experimental results of Bhaga & Weber (Reference Bhaga and Weber1981). The parameter values are (i) 
$(Ga, Bo)=(2.316, 29)$, (ii) (3.094, 29) and (iii) (4.935, 29).
Figure 3. Validation of the numerical technique: comparison of the temporal evolution of the particle settling velocity obtained from our predictions and published experimental data; (a–c) show results for 
$Re=11.6$ and 
$31.9$ (ten Cate et al. Reference ten Cate, Nieuwstad, Derksen and Van den Akker2002), and 
$Re=41$ (Mordant & Pinton Reference Mordant and Pinton2000), respectively.
3. Results and discussion
The purpose of this study is to gain insight into the interaction between a freely rising deformable bubble and a freely settling particle. We employed a hybrid level-set/front-tracking method to capture the interaction dynamics between the bubble and the particle. Unless otherwise stated, the particle and bubble radii are kept equal, 
$\varOmega =1$, and the particle and the bubble are released in-line, with a distance of 
$6R$ from their relative centres. The results in this section will investigate the effect of the bubble Bond and Galilei numbers on the interaction between a bubble and a particle. This work will also consider the effect of varying the solid-to-fluid density contrast 
$\zeta _p$, varying the bubble-to-particle size ratio 
$\varOmega$, and the initial bubble position on the interaction dynamics.
3.1. In-line bubble–particle interactions
The first part of the discussion will consider the interaction between a rising bubble and a settling particle with 
$Ar = 10 (\zeta _p=1)$ and 
$Ga = 10$. The bubble Bond number will be varied between 
$0.5$ and 
$5$. At these conditions, the bubble trajectory is expected to follow a rectilinear path as found by Zhang, Ni & Magnaudet (Reference Zhang, Ni and Magnaudet2021). Since the bubble and the particle are initially released in line, they will interact through head-on collisions.
Figures 4(a) and 4(b) show the spatiotemporal evolution of a rising bubble interacting with a settling particle with 
$Bo=0.5$ and 
$Bo=5$, respectively. As the bubble rises in a quiescent Newtonian fluid, several different bubble-shape regimes were found by Tripathi et al. (Reference Tripathi, Sahu and Govindarajan2015). Under the conditions considered in this study, the rising bubble has an axisymmetric shape, as seen in figure 4(a,b).
Figure 4. Spatiotemporal evolution of the rising bubble as it approaches the particle is shown in (a,b). The snapshots of the bubble contours are shown for 
$t=1, 2, 3, 4$, and the particle is only shown at 
$t=4$. In (a), a scale indicating the separation distance is provided between the bubble and the particle, indicating the bubble and the particle's initial centre positions at 
$z_b^0$ and 
$z_p^0$, respectively. The non-dimensional parameters are (a) 
$Bo=0.5, Ga=10, Ar=10$, (b) 
$Bo=5, Ga=10, Ar=10$; (d–f) show the change in energy of the system for three different 
$Bo$ numbers and maintaining 
$Ga=10$, and they are 
$Bo=0.5$, 
$Bo=2.0$ and 
$Bo=5.0$, respectively. Here, 
$E_m$ is the mechanical energy of the system, such that 
$\Delta E_m = \Delta E_p + E_k$. A new time scale is introduced here, 
$t^*$, which corresponds to the time for the reversal of kinetic energy due to the interaction. Panel (c) showcases the relative velocity of the bubble and the particle when varying 
$Bo$ and 
$Ga$ at 
$t^*$, and also showcases the difference in buoyancy force between the bubble and the particle when varying 
$\zeta _p$.
Prior to the interaction with the settling particle, the bubble is found to rise in a spherical shape when 
$Bo = 0.5$, and an ellipsoidal shape when 
$Bo = 5$. This is because when surface tension forces dominate at lower 
$Bo$, the deformation of the bubble will be weak (Qin, Ragab & Yue Reference Qin, Ragab and Yue2013). During the initial interaction process, the particle surface begins to deform the bubble as seen in figure 4(a,b), such that a dimple forms at the bubble apex, the part of the bubble closest to the particle; dimple formation at the bubble apex was observed for all the parameters studied in the present work, which is expected due to the increase of hydrodynamic pressure in the gap region between the bubble and the particle (Chan et al. Reference Chan, Klaseboer and Manica2011). The dimple can be characterised by a negative mean curvature at the centre, such that 
$(\kappa _1+\kappa _2)/2 <0$, indicating a concave deformation. A liquid film also begins to develop between the bubble and the particle surface, which drains due to capillary and gravitational forces, causing it to thin over time. The liquid film drainage described in this work during the interaction between the bubble and the particle is observed until they are at most one grid cell apart. In dimensional terms, one grid cell is of the order of 
${\sim }O (10\ \mathrm {\mu }\textrm {m})$, such that the drainage process is purely hydrodynamic and disjoining forces from the interaction can be neglected (Lee, Rudman & Slim Reference Lee, Rudman and Slim2023). Furthermore, the contact line dynamics between the bubble and the particle are considered for different particle wettability parameters.
Upon increasing 
$Bo$ to 
$Bo=2$ or 
$Bo=5$, the gravitational forces play a more significant role compared with the surface tension forces. As the bubble and the particle approach each other, two distinct observations can be made. Firstly, the onset of dimpling occurs at a position where the bubble and particle are farther away compared with when 
$Bo=0.5$. This can be seen in figure 4(a,b) and by comparing the bubble aspect ratio 
$\chi$ and the normalised distance 
$S$ between the bubble and the particle as seen in figure 5(a,b), respectively; here, 
$\chi$ is defined as the ratio of the bubble's height to its width, and is a measure of its deformation during its rise and interaction with the particle. The dimple formed at the apex is also more pronounced compared with the case when 
$Bo = 0.5$, and the deformation starts at the cap of the bubble, which leads to a change in the overall bubble shape. The increase in Bond number leads to flattening of the side of the bubble farthest away from the particle, as seen in the interaction between droplets and solid walls at low Reynolds numbers (Ascoli et al. Reference Ascoli, Dandy and Leal1990). It proves instructive to carry out an energy analysis of the system when considering the dynamics of the bubble–particle interactions. The four contributions to the total system energy are the potential energy 
$E_p$, surface energy 
$E_s$, kinetic energy 
$E_k$ and viscous dissipation 
$E_{\eta }$. The changes in the potential and surface energies, respectively, read
\begin{gather} \Delta E_p=\rho Vg(z_c-z_{0}), \end{gather}
\begin{gather}\Delta E_s=\gamma (A-A_0), \end{gather}while the kinetic energy and viscous dissipation are, respectively, expressed by
\begin{gather} E_k=\tfrac{1}{2}\rho \int_V (u^2+v^2+w^2) \, {\rm d} V, \end{gather}
\begin{gather}E_{\eta}=\int_0^t \int_V \xi \, {\rm d} V \, {\rm d} t, \end{gather}
where 
$A$ is the surface area, 
$V$ is the volume, 
$z_0$ is the bottom of the domain and 
$\xi$ is the viscous dissipation function. The temporal evolution of the energy budget during the initial interaction process when 
$Bo=0.5$, 
$Bo=2$ and 
$Bo=5$ are shown in figure 4(d–f), respectively. Here, 
$\Delta E_m$ is the mechanical energy of the system, such that 
$\Delta E_m=\Delta E_p^{total}+E_k^{total}$, and the results are normalised to the total energy of the system 
$E_t=E_p+E_k+E_s+E_\eta$. As the bubble rises and the particle settles, 
$E_k$ is found to increase up to a maximum, until the flow field of the particle and the bubble interact. At this point, 
$E_k$ starts to decrease as the gap distance between the bubble and the particle thins. We introduce a new time scale 
$t^*$ which specifies the time of kinetic energy reversal in figure 4(d–f). For 
$Bo=0.5$, 
$t^* \approx 3.0$ and for 
$Bo=2$ or 
$Bo=5$, 
$t^*\approx 3.5$. In figure 4(c), the relative interaction velocity at the peak kinetic energy is presented for the case of 
$Ga=10$ and 
$Ga=20$ when varying the 
$Bo$ number. When 
$Ga=10$, it can be seen that the relative velocity does not vary significantly when increasing 
$Bo$, however, the relationship between the relative velocity and 
$Bo$ remains the same for both 
$Ga$. When increasing 
$Bo$, the bubble rising velocity decreases, which allows the particle to attain a larger velocity as it settles. However, this increase in particle velocity does not compensate for the difference in bubble rise speed, which causes a decrease in the relative velocity prior to the interaction. A similar trend is observed for 
$Ga=20$, with a larger increase in relative velocity when decreasing 
$Bo$. Furthermore, the difference in buoyancy force between the bubble and the particle is also presented when varying 
$\zeta _p$, where the buoyancy force was taken as 
$F_b^i=Vg(\rho _i-\rho _l)$, where 
$i=(b, p)$ for the bubble or the particle, respectively. When 
$\zeta _p=1$, the buoyancy force of the particle is comparable to that of the bubble, however, decreasing 
$\zeta _p$ leads to a larger difference between the bubble and particle buoyancy force. A larger difference in buoyancy force between the bubble and the particle allows for the bubble to ‘push’ the particle upwards during the initial interaction, rather than remain in the same vertical position during film drainage. For low 
$Bo$, the particle kinetic energy is seen to decrease to zero at later times, which shows that the film drainage process causes the bubble and the particle to remain at approximately the same vertical position as the liquid film thins (see figure 4d). The formation of the liquid film between the bubble and the particle causes an outward radial velocity of the background fluid, which leads to the dissipation of kinetic energy, as seen at 
$t=3.5$ in figure 6 (Zawala & Dabros Reference Zawala and Dabros2013).
Figure 5. Temporal evolution of the aspect ratio, 
$\chi$, and the distance between the bottoms of the bubble and particle, 
$S$, shown in (a) and (b), respectively, for 
$Bo=0.5$, 
$Bo=2.0$ and 
$Bo=5$, with 
$Ga=10$, and 
$\zeta _p=1$. Panel (c) showcases the minimum aspect ratio when varying 
$Bo$ and maintaining 
$Ga=10$, while also plotting the change in interfacial surface area at that time step 
$(t_{\chi _{min}})$. The film thickness profile is presented for these cases at 
$t_{\chi _{min}}$ in (d), and the temporal evolution of the dimple radius 
$R_D$ is presented in (e).
Figure 6. Spatiotemporal evolution of the bubble–particle interaction dynamics shown for the same parameters used to generate figure 4; the rows are associated with (a) 
$Bo=0.5$, (b) 
$Bo=2.0$ and (c) 
$Bo=5$. The contours on the left-hand part of each panel show the logarithm of the viscous dissipation function 
$\log _{10}\xi$ in the fluid, and the contours on the right-hand part show the velocity magnitude in the fluid 
$V$. The bubble is shown with glyphs of the velocity to illustrate the relative motion of the bubble at different time steps.
As film drainage is a gradual process driven by pressure differences in the film region, this allows time for the surface tension forces to help the bubble recover a semispherical shape with a small dimple at the bubble apex. The dimple radius decreases during the drainage of the liquid film, which can be seen at 
$t \approx 5$ when 
$Bo = 0.5$ (see figure 6) and by comparing the liquid film thickness profiles and the bubble interface shape in figure 9(a,b). Yiantsios & Davis (Reference Yiantsios and Davis1990) have considered the interaction between a droplet and a solid surface and their numerical simulations showed that the terminal value of the dimple radius scales with 
$Bo$ as 
$R_{D}\sim (\frac {2}{3}Bo)^{{1}/{2}}R$, which is in excellent agreement with the dimple radius found in this work when 
$Bo=0.5 (R_{D}\approx 0.5)$. This result was first obtained by Derjaguin & Kussakov (Reference Derjaguin and Kussakov1939) by assuming the pressure in the dispersed phase is approximately equal to the pressure in the dimple, and balancing the pressure with the buoyancy force such that 
$({2\gamma }/{R}){\rm \pi} R_{D}^2 \sim \frac {4}{3}{\rm \pi} (\rho _l - \rho _b)g R^3$. In figure 5(e), the dimple radius is presented during the initial interaction period for different 
$Bo$ numbers when 
$Ga=10$ and 
$\zeta _p=1$, where the scaling of the dimple radius presents a good approximation of the results.
Upon increasing the 
$Bo$ number to 
$Bo=5$, a significantly different interaction process is observed at later times. As the particle deforms the bubble interface, the dimple that forms interacts with the southernmost part of the bubble, such that the bubble will rupture. Prior to its rupture, the bubble is seen to deform severely with a film thickness larger than that associated with the 
$Bo=0.5$ case (see figures 5d and 6 at 
$t=5$). The bubble topology resembles the splashing configuration found by Yoon & Shin (Reference Yoon and Shin2021) when considering the interaction between a droplet and a particle, with a liquid film between the bubble and the particle in the present work. The minimum aspect ratio attained during the interaction is presented when varying 
$Bo$ in figure 5(c). The time of maximum bubble deformation (minimum 
$\chi$) is taken as 
$t_{\chi _{min}}$, and the change in interfacial surface area 
$\Delta A/A_0$ is presented. At 
$t_{\chi _{min}}$, the film thickness profile between the bubble and the particle is plotted in figure 5(d). The profiles of the film thickness show that at larger 
$Bo$, the minimum gap distance is found at the bubble centre, and when decreasing 
$Bo$, the minimum gap distance shifts towards the bubble dimple. The aspect ratio 
$\chi$ reaches a minimum of 
$\sim 0.2$, and the normalised distance between the particle and the bubble, 
$S$, continues to decrease. The distortion of the bubble surface caused by the particle drives the increase in 
$\chi$ at later times, up to the point of bubble rupture. Figure 7(b) shows the passage of the particle post-bubble-rupture. The particle's kinetic (potential) energy continues to increase (decrease) as seen in figure 7(d), which shows that the particle continues to settle post-bubble-rupture. This bubble rupture mechanism was also found by the two-dimensional, axisymmetric simulations of the interaction between a rising bubble and a solid surface (Qin et al. Reference Qin, Ragab and Yue2013), where bubble rupture was observed at large 
$Bo$ numbers. As the bubble and the particle approach each other when 
$Bo=2$, the bubble is found to deform similarly to the 
$Bo=5$ case (see figure 5a), however, the bubble does not rupture. Therefore, the liquid film in the gap region between the bubble and the particle starts to drain. As the liquid film drains, the bubble's buoyancy drives it to shift away from the particle, as seen in figure 7(a), interrupting the film draining process. As the bubble moves away from the particle, we note that no contact was found between the bubble and the solid surface, as seen from the presented top view. The bubble then continues to rise, and the particle continues to settle, as seen by considering the change in mechanical energy in the energy budget plot presented in figure 7(c).
Figure 7. The spatiotemporal evolution of the bubble–particle interaction dynamics is shown in (a) for 
$Ga=10, Bo=2$ and 
$\zeta _p=1$, where the contours depict the velocity magnitude. Panel (b) shows the spatiotemporal evolution of the interaction dynamics when 
$Ga=10, Bo=5$ and 
$\zeta _p=1$. The time steps presented are relative to 
$t^*$, such that 
$t-t^*$ varies from 
$2.5$ to 
$7.0$. When 
$Bo=2$ and 
$t-t^*=7$, a top view of the interaction between the bubble and the particle is presented to showcase that no contact is found between them. In (c,d), the normalised energy budget 
$\tilde {E}/\tilde {E}_0$ is presented for 
$Bo=2$ and 
$Bo=5$, respectively, where 
$\tilde {E}_0=E_{p_0}^{particle}+E_{s_0}$.
In figure 8(a,b), the temporal evolution of 
$h(t)$, defined as the shortest bubble–particle distance at every time step, is presented for different 
$Bo$ numbers. As the bubble rises at early times, it assumes a spherical or ellipsoidal shape such that 
$h(t)$ is measured along centreline of the bubble and particle, which we define as 
$h_{0}(t)$. Figure 8(a) provides an illustration of 
$h_{0}(t)$ and 
$h(\boldsymbol {x},t)$, such that 
$h(t)= \min h(\boldsymbol {x},t)$. As the dimple forms on the bubble surface, 
$h(t)$ begins to deviate from 
$h_0(t)$. When comparing the minimum gap distance for different 
$Bo$, it can be seen that for smaller 
$Bo$, the film drainage process occurs fastest as the bubble deformations are less severe. This leads to the liquid in the film to drain over a smaller surface area (Manga & Stone Reference Manga and Stone1993). Furthermore, increasing 
$Bo$ leads to an increase in surface area in the gap region, which delays the film drainage process. This is found to occur up to 
$Bo=2$, because at 
$Bo=5$, the bubble ruptures. Comparing the difference between the gap thickness at the centre and the minimum gap thickness when increasing 
$Bo$, two different conclusions can be made. First, at low 
$Bo$, the development of the dimple occurs at earlier times when compared with the larger 
$Bo$ number cases. However, the dimple is much more pronounced with increasing 
$Bo$, at later times. Nevertheless, when the bubble ruptures at 
$Bo=5$, the minimum gap distance is always found at the centreline as the particle continues to settle. Previous studies (Jones & Wilson Reference Jones and Wilson1978; Yiantsios & Davis Reference Yiantsios and Davis1990; Chan et al. Reference Chan, Klaseboer and Manica2011; Quan Reference Quan2012; Denner Reference Denner2018) have shown that the film thickness during the interaction between a deformable interface and a solid surface scales as 
$t^{-\alpha }$ where 
$0\leq \alpha \leq 1$. During film drainage, we can assume that the outward radial velocity of the fluid is much greater than the vertical velocity of the bubble or the particle. Therefore, a new characteristic velocity scale 
$U$ can be obtained by balancing the buoyancy and viscous drag 
$(\Delta \rho R^3 g \sim \mu _b R U)$, such that 
$U=g\rho _l (2R)^2/6\mu _l\phi$ since 
$\rho _l \gg \rho _b$, and 
$\phi =(2+3/\lambda )/(1+1/\lambda )$ is a correction factor for the viscosity ratio 
$\lambda$ (Denner Reference Denner2018). We introduce a new draining time scale, 
$t_d =2R \sqrt {\mu _l/\gamma U}$, which can be obtained by assuming the film drainage is dominated by the balance of surface tension and viscous stresses (Klaseboer et al. Reference Klaseboer, Chevaillier, Gourdon and Masbernat2000; Chan et al. Reference Chan, Klaseboer and Manica2011; Denner Reference Denner2018). At smaller Bond numbers, the liquid film thickness is found to scale as 
$h\sim (t-t^* /t_d)^{-4/5}$, whereas when increasing the Bond number, the scale varies to 
$h\sim (t-t^* /t_d)^{-1/2}$. This result can be seen in figure 8(b). We note that when 
$Bo=2$, the film drainage result presented is prior to the bubble sliding away from the particle. Figure 9(c–f) presents the interaction between a bubble and a particle when 
$\varOmega =2$ when 
$Ga=10, Bo=0.5$ and 
$\zeta _p=1$. As the bubble and the particle interact, the approach velocity of the particle when 
$\varOmega =2$ is smaller, such that the kinetic energy from the larger bubble's motion is transferred to the particle, allowing it to reverse its direction and float. Here, 
$t^*$ was found to be around 
$t\approx 3.5$. As this kinetic energy from the bubble's motion is transferred to the particle, the bubble's second approach to the particle is at a substantially smaller velocity, such that the total kinetic energy stabilises around 
$t=5$ to 
$t=5.9$. The radius of the dimple on the bubble apex starts to decrease, as seen from the liquid film thickness profiles in figure 9(e), and the liquid film ruptures at 
$t-t^*\approx 2.5$. It can be seen that as opposed to the case where 
$\varOmega =1$, decreasing the particle size leads to the bubble continuing to rise during the film-draining process, which can be seen from the bubble apex position in figure 9(b, f) for the two cases. When the film ruptures, the larger bubble continues to push the particle upwards, until reaching the top boundary.
Figure 8. Panels (a) and (b) show the temporal evolution of the minimum gap thickness 
$h(t)$ for different 
$Bo$ when 
$Ga=10$ and 
$\zeta _p=1$. Two different scaling are also provided, such that 
$h\sim (t-t^* /t_d)^{-\alpha }$ where 
$\alpha$ here is either 
$-4/5$ or 
$-1/2$.
Figure 9. This figure showcases the film thickness and the bubble shape as it approaches the particle when 
$Ga=10, \zeta _p=1$ and 
$Bo=0.5$ in (a,b). In (c,d), the temporal evolution of the kinetic energy and particle velocity is presented for two cases, and they are when 
$\varOmega =1$, which is the base case, and when 
$\varOmega =2$. In (e, f), the film thickness profile and the bubble shape at different time steps during the interaction for 
$\varOmega =2$ are presented.
We further consider the effect of increasing 
$Ga$ and varying 
$\zeta _p$ on the interaction dynamics in figure 10. Seven different cases are presented, which correspond to 
$(Bo,Ga, \zeta _p)=(0.5,20,1$), (
$1,20,1$), (
$2,20,1$), (
$5,20,1$), (
$0.5,20,0$), (
$0.5,20,0.25$) and (
$0.5,20,0.75$). In figure 10(a), the temporal evolution of the particle velocity is presented. Previously, at a Bond number of 0.5, a Galilei number of 10 and a particle-to-fluid density ratio of 1, the bubble and the particle remained at relatively the same vertical position throughout the film draining process. This indicates that the bubble and particle experienced minimal vertical displacement relative to each other due to the difference in buoyancy force between the bubble and the particle. When 
$Ga$ is increased to 20 while maintaining 
$\zeta _p = 1$, the bubble approaches the particle with a greater velocity compared with when 
$Ga = 10$. However, the interaction dynamics remain consistent with the case when 
$Ga = 10$, where the magnitude of the particle velocity decreases, and the particle remains at relatively the same vertical position during film drainage (see figure 10c). In this investigation, we consider the effects of decreasing the particle-to-fluid density ratio (
$\zeta _p$), which introduces new dynamics into the interaction process. Specifically, a lower 
$\zeta _p$ makes the particle less dense relative to the fluid, increasing its susceptibility to the forces exerted by the bubble. This increased sensitivity results in more pronounced motion reversals and changes in velocity, demonstrating a critical dependence on the density ratio for the overall behaviour of the bubble–particle interaction. When 
$\zeta _p$ is decreased, the increased velocity of the bubble causes it to push the particle upwards as the gap distance between them decreases. The extent of flotation of the particle increases as 
$\zeta _p$ decreases. This can be verified by examining figure 10(a–c), which show the particle velocity, the total kinetic energy and the change in potential energy of the particle, respectively. The reversal of the particle's motion is similar to the behaviour observed when 
$\varOmega = 2$. As the particle reverses its motion due to its interaction with the bubble, the bubble's velocity decreases. This deceleration results from the transfer of energy from the bubble's motion to the particle, leading to a redistribution of momentum. The interaction process when 
$\zeta _p=[0,1]$, 
$Ga=20$ and 
$Bo=0.5$ can be seen in figure 10.
Figure 10. In (a–c), the particle settling velocity, total kinetic energy and the change in particle potential energy are presented, respectively. These results correspond to when 
$Ga=20, Bo=0.5$ and varying 
$\zeta _p$. Temporal evolution of the minimum gap thickness 
$h(t)$ is shown in (d) for 
$\zeta _p=1$, 
$Ga=20$ and 
$Bo=[0.5\unicode{x2013}5]$; an enlarged view of the plot in (d) is shown in ( f), for 
$Bo=0.5$ and 
$Bo=1$ as at these 
$Bo$ numbers, the bubble did not rupture. A scaling for the film drainage is provided such that 
$h\sim ( (t-t^*)/t_d )^{-4/5}$. Panel (e) showcases the film thickness profile at 
$t_{\chi _{min}}$ when varying 
$Bo$ and maintaining 
$Ga=20$ and 
$\zeta _p=1$.
When considering the temporal evolution of the gap thickness between the bubble and the particle for 
$Ga=20$ and 
$\zeta _p=1$, a similar general trend is found compared with that shown in figure 8. Increasing 
$Bo$ leads to an increase in the surface area between the bubble and the particle, delaying the process of film drainage. However, when comparing figure 10(d) with 8(a), the bubble ruptured for both 
$Bo=2$ and 
$Bo=5$ when increasing the impact relative velocity by varying 
$Ga$. Figure 10(e) shows the film thickness when varying 
$Bo$ and maintaining 
$Ga=20$ and 
$\zeta _p=1$ at 
$t_{\chi _{min}}$. As seen in figure 5(d), the film thickness profile indicates an increase in distance between the bubble and the particle when increasing 
$Bo$, and the minimum thickness shifts from the bubble dimple to the centreline. Following a similar scaling method to figure 8, the minimum film thickness when 
$Bo=0.5$ and 
$Bo=1$ scales as 
$h\sim (t-t^* /t_d)^{-4/5}$.
For head-on collisions at relatively small 
$Bo$ numbers, a contact line forms between the bubble and the particle after liquid film drainage. In figure 11, the evolution of the contact line dynamics and the effect of particle wettability on the overall dynamics is studied. Three different equilibrium contact angles are considered, and they are 
$\theta _e=(55, 105, 125)$ for 
$\varOmega =\zeta _p=1, Ga=10$ and 
$Bo=0.5$. First, figure 11(a) showcases the particle vertical velocity for the three cases. When 
$\theta _e=55$, the particle velocity is seen to increase significantly post-contact-line formation, which indicates that the bubble has attached to the particle and is pushing it upwards. When increasing 
$\theta _e$ to 
$105^\circ$ or 
$125^\circ$, the particle velocity is not as significantly affected, indicating the bubble is attempting to minimise its contact area with the particle. In figure 11(c), the radius of the contact line 
$R_{TPC}$ is shown when varying 
$\theta _e$ and 
$\varOmega$. When 
$\theta _e=55^{\circ }$, the contact line initially expands, and during the expansion, the bubble vertically deforms which leads to a peak in the aspect ratio 
$\chi$ at 
$t-t^*\approx 4.8$ (figure 11b). The contact line radius then maintains the same size as it floats with the particle until reaching the top boundary. When increasing 
$\theta _e$, the contact line radius initially expands, however, the bubble ‘slides’ across the particle and attempts to minimise the contact area. This leads to a decrease in the contact line radius as the bubble moves to the particle's apex, and the bubble then detaches due to buoyancy and continues to rise. The bubble detaches at 
$t-t^* \approx 11$ when 
$\theta _e=105^\circ$ and 
$\theta _e=125^\circ$, with a small part of the bubble still attached to the particle surface. When decreasing the particle size (
$\varOmega =2$) and maintaining 
$\theta _e=105^\circ$, the radius of the contact line expands during contact. However, due to the particle's upward velocity induced by the bubble during the initial interaction, the bubble and the particle continue to rise and the contact line remains relatively the same size. It is crucial to note that the initial radius of the contact line indicates the position at which the contact line forms, such that the dimple radius has a direct influence on the dynamics of the contact line.
Figure 11. This figure showcases the dynamics of the interaction after the formation of the contact line, when 
$Ga=10$, 
$Bo=0.5$ and 
$\zeta _p=1$. In panel (a), the particle's vertical velocity is presented for different particle wettability parameters, which correspond to a base case with 
$\theta _e=105^\circ$, and two additional cases with 
$\theta _e=55$, and 
$\theta _e=125$. Panel (b) showcases the aspect ratio 
$\chi$ for the three different cases, and the pentagram indicates reaching the top boundary. Lastly, (c) showcases the temporal evolution of the radius of the contact line 
$R_{TPC}$ when varying 
$\theta _e$ and 
$\varOmega$. Note that the 
$x$-axis is scaled with the contact line time 
$t_{TPC}$. The subpanels in (d) correspond to 
$t-t^*$ from 
$8$ to 
$11$ for the base case and 
$t-t^*$ from 
$5$ to 
$7$ when 
$\theta _e=55^{\circ }$.
3.2. Off-centre bubble–particle interactions
Here we consider bubble–particle interactions that do not correspond to head-on collisions. These are initiated by introducing a horizontal deviation 
$\delta$ in the initial bubble position. In this section, we study the dynamics associated with 
$\delta =[0.25, 0.5, 1]$ and keep 
$(Ga=10, \zeta _p=1)$ fixed. In figure 12, the spatiotemporal evolutions of the interaction between the bubble and the particle are presented for 
$Bo=0.5$ (figure 12a) and 
$Bo=5$ (figure 12b), and 
$\delta =1$. The background field shows the logarithm of the viscous dissipation function 
$\xi$ in the fluid. When initiating an offset between the bubble and the particle, it can be seen in figure 12 that for 
$Bo=0.5$, the bubble initially rises in a similar fashion as observed with the in-line collisions, but the dimple at the apex of the bubble forms on the side of the bubble rather than at the centreline. This causes the bubble to shift its velocity towards the other side (as seen in the glyphs) as it deforms during the interaction process. A liquid film forms between the bubble and the particle at the side of the bubble, at which point we have the largest dissipation of energy. The bubble tends to exhibit a ‘sliding’ motion, which is similar to the drafting–kissing–tumbling regime found by Zhang et al. (Reference Zhang, Ni and Magnaudet2021) between two bubbles rising in line. The bubble is then seen to continue its rise as it recovers its original spherical shape, which can be seen at 
$t=5.5$ and 
$t=6$. This can also be observed by the change in surface area plot in figure 13(b), where the change in surface area decreases after the interaction with the particle.
Figure 12. Spatiotemporal evolution of off-centre bubble–particle interaction dynamics with 
$Ar=10$, 
$Ga=10$, and an initial bubble location offset of 
$\delta =1$. The results shown in (a) and (b) were generated with 
$Bo=0.5$ and 
$Bo=5$, respectively. The contours are of the logarithm of the viscous dissipation function 
$\log _{10}\xi$ in the fluid. The bubble is shown with glyphs of the velocity to illustrate the relative motion of the bubble at different time steps. The snapshots are shown between 
$t=2.5$ and 
$t=6$ in time increments of 
$\delta t=0.5$.
Figure 13. Temporal evolution of the total kinetic energy, and the change in bubble surface area are shown in (a) and (b), respectively, for off-centre bubble–particle interactions with 
$(Bo,Ga,Ar,\delta ) = (0.5,10,10,1$) and (
$5,10,10,1$). Panels (c) and (d) show the temporal evolution of the settling velocity, and the 
$x$-component of the particle velocity, respectively. In (e, f), the side view of the particle trajectory with respect to the 
$x$-axis and the 
$y$-axis are presented, respectively.
Upon further increasing 
$Bo$ to 
$Bo=5$, a similar initial interaction process is observed such that the bubble rises with an ellipsoidal shape, and a dimple forms at the side of the bubble rather than at the centreline. However, due to the large deformations associated with this 
$Bo$ value, it can be seen that although a drafting–kissing–tumbling scenario is found, the side of the bubble associated with the impact forms a ‘tail’ (
$t=4.5$–5.5). As the tail forms, we have a sharp increase in the bubble interfacial area and energy dissipation, while the centre of mass of the bubble continues to rise as the total kinetic energy increases (see figure 13a). The tail of the bubble then retracts to its centre of mass (
$t=6$ in figure 12) and continues to rise, and the particle continues to settle.
From figure 13(d) it can be seen that upon increasing 
$Bo$, the deviation in the particle settling trajectory is not as affected as for low 
$Bo$. This is associated with the significant deformation of the bubble, which tends not to affect the 
$x$-component of the particle velocity as much as the sliding motion of the bubble observed for 
$Bo=0.5$. Furthermore, when fixing 
$Bo=0.5$ and varying 
$\delta$ between 0.25 and 0.5, figure 13(e, f) shows that the particle trajectory is affected more by decreasing the value of 
$\delta$. For 
$\delta =0.25$ and 
$\delta =0.5$, the particle's settling velocity (see figure 13c) decreases during the initial interaction process. Although the decrease in settling velocity is more pronounced as 
$\delta$ decreases, i.e. approaching the case of in-line bubble–particle interactions, the variation in the 
$x$-component of the velocity is more pronounced when 
$\delta =0.5$ compared with 
$\delta =0.25$ (see figure 13f). This leads to the particle trajectory being rather similar when 
$\delta =0.25$ and 
$0.5$ (see figure 13e).
4. Conclusion
This paper considered the interaction between a freely rising deformable bubble and a freely settling particle in a quiescent, Newtonian background fluid. The considered non-dimensional parameters were the bubble Bond and Galilei numbers, the particle Archimedes number, governed by the solid–fluid density contrast, and the size ratio between the bubble and the particle. The results highlighted the spatiotemporal evolution of the bubble shape during the interaction process, the bubble aspect ratio as a measure of its deformation, the liquid film thickness and an energy analysis of the system. When considering the bubble shape, it was found that increasing the Bond number led to more severe bubble deformations, leading up to bubble rupture and its penetration by the particle. When the 
$Bo$ number is relatively low, the interaction includes the formation of a liquid film between the bubble and the particle, which is associated with the dissipation of kinetic energy. When decreasing 
$\zeta _p$, the buoyancy force of the bubble exceeds that of the particle, and the bubble can lift the particle upwards if the approach velocity of the bubble is considerably larger than that of the particle. However, if the particle and the bubble approach velocities are similar in magnitude, the initial interaction process is associated with the dissipation of kinetic energy which causes the bubble and the particle to remain at their interaction position as the liquid film drains. The study also considered the temporal evolution of the minimum gap thickness 
$(h(t))$ between the bubble and the particle and found that increasing 
$Bo$ led to an increase in surface area in the gap region, further delaying the film-draining process. The effect of particle wettability on the contact line dynamics was also considered, such that the bubble remains attached to the particle surface at lower 
$\theta _e$ and slides away at larger 
$\theta _e$.
A complex situation arises when introducing an offset between the bubble and the particle initial positions. In this study, an offset is only considered for the bubble initial position in the horizontal direction at constant Archimedes and Galilei numbers. For low Bond numbers, a sliding motion of the bubble was observed as it approaches the particle, causing the bubble to push the particle away from its original settling trajectory. However, with increasing Bond number, the bubble deformation is seen to be much more severe such that the bubble centre of mass continues to rise with the formation of a tail in the particle's proximity. The tail then retracts, the bubble continues to rise and the particle trajectory is not as significantly affected.
This work has considered one viscosity and density ratio between the bubble and the background fluid 
$\lambda =\mu _l/\mu _b=100$ and 
$\beta =\rho _l/\rho _b=1000$. Possible avenues for future work may include the effects of varying the viscosity and density ratios, considering a larger difference between the particle and bubble sizes, introducing asymmetric particle shapes, and the effect of density stratification in the background fluid. Furthermore, the effects of surface mobility of the bubble on the interaction process, by introducing surfactants to the system for example, presents an interesting avenue for future work due to its relevance for industrial processes. All of these configurations may lead to different interaction dynamics between the bubble and the particle. Lastly, this work has considered the interaction between a single rising bubble and a single settling particle, however, the collective effect of the interaction between a rising bubble, or a group of bubbles, in the presence of a particle suspension remains an area for further exploration.
Acknowledgements
The numerical simulations were performed with code BLUE (Shin et al. Reference Shin, Chergui and Juric2017) and the visualisations were generated using ParaView.
Funding
This work was supported by the Engineering and Physical Sciences Research Council, UK, through the MEMPHIS (EP/K003976/1) and PREMIERE (EP/T000414/1) programme grant. O.K.M. acknowledges funding from PETRONAS and the Royal Academy of Engineering for a Research Chair in Multiphase Fluid Dynamics. A.A. acknowledges the Kuwait Foundation for the Advancement of Sciences (KFAS) for their financial support. A.A. and L.K. acknowledge HPC facilities provided by the Imperial College London Research Computing Service. A.A. also acknowledges fruitful discussions with Debashis Panda (Imperial, Chemical Engineering). D.J. and J.C. acknowledge support through HPC/AI computing time at the Institut du Developpement et des Ressources en Informatique Scientifique (IDRIS) of the Centre National de la Recherche Scientifique (CNRS), coordinated by GENCI (Grand Equipement National de Calcul Intensif) grant 2024 A0162B06721.
Declaration of interests
The authors report no conflict of interest.
Appendix A
In this section, a resolution test is presented for the case when 
$(Bo,Ga,Ar) = (0.5,10,10)$ in figure 14. Two different mesh resolutions are considered, and they correspond to 
$\Delta x/R =1/24$ and 
$\Delta x/R =1/32$. When considering the distance between the particle south pole and the bubble south pole, and the total kinetic energy of the system, an excellent agreement is found between the two mesh resolutions. The most difficult part to resolve is the film thickness. The lower resolution root mean square error compared with the higher resolution is approximately 
$0.5\,\%$, which is acceptable considering the higher resolution is the coarsest mesh utilised in this study. In figure 14(c), a domain size independence test is performed when 
$Bo=0.5$ and 
$Ga=10$, by increasing the domain size to 
$24R \times 24R \times 24R$. The results show that when increasing the domain size, there is no difference in the attained result, further increasing our confidence in the numerical set-up. In figure 14(d), a mesh size independence considering three different grid cell sizes is presented for the case of 
$Bo=5.0$ and 
$Ga=10$, where a resolution of up to 
$\Delta x/R =1/64$ is presented.
Figure 14. Panels (a,b) show the total kinetic energy of the system and the distance between the particle south pole and the bubble south pole (
$S$), at two different resolutions for 
$Ga=10, Bo=0.5$ and 
$\zeta _p=1$, respectively. Panel (c) shows the temporal evolution of the total kinetic energy when increasing the domain size for the same case. It can be seen that increasing the domain size leads to no difference in the results. Panel (d) shows the total kinetic energy for different mesh sizes for 
$Ga=10, Bo=5.0$ and 
$\zeta =1$, up to a resolution of 
$\Delta x/R = 1/64$. The results show that the domain size and mesh resolution utilised in this work is sufficient to capture the dynamics of the interaction between the bubble and the particle.
Appendix B
This section presents additional validation case studies for the numerical method utilised in this work. At first, figure 15 presents illustrations for the contact line methodology that is utilised in our in-house solver code BLUE. Furthermore, the adhesion of a bubble on a solid wall is presented and qualitatively compared with the experiments of Basarova et al. (Reference Basarova, Soušková and Zawala2018). The bubble is initialised close to the flat wall (
$0.1R$ from surface to surface), and as the contact line forms, the deformation of the bubble is accurately captured when compared with the results of the experiment. The bubble size is 
$0.705$ mm, and the viscosity and density ratios correspond to an air/water system. The equilibrium contact angle 
$\theta _e$ was set at 
$\theta _e \approx 100^{\circ }$, which corresponds to the surface wettability characteristics utilised in the experiments.
Figure 15. Panel (a) provides an illustration of the contact line numerical methodology between a bubble and a solid wall or solid particle. Panel (b) showcases the temporal evolution of the south pole of the bubble as it adheres to a solid wall with negligible inertia. The contact angle for this case was 
$\theta _A=102, \theta _R=99$. A qualitative comparison with the experimental results of Basarova et al. (Reference Basarova, Soušková and Zawala2018) is provided at different time steps.
Furthermore, additional validation cases are presented in figure 16 for a solid particle settling through an initially unperturbed interface with different fluid and particle properties. The particle was initially placed a distance of 
$6R$ away from the unperturbed interface (from the centre of the particle) and is allowed to settle due to gravity. The results attained from the numerical method are in close agreement with the results of the experiments of Pierson & Magnaudet (Reference Pierson and Magnaudet2018a,Reference Pierson and Magnaudetb), where a few differences were found in the case of a top fluid of V500 silicone oil with water as the bottom fluid and a steel sphere. The differences are mainly in the development of axisymmetric corollas near the particle surface as the particle penetrates the bottom fluid. The reason behind the difference may be attributed to the different initial position of the particle in the experiments as this parameter was not exactly specified, or due to the utilised resolution of the case study, which was 
$\Delta x/R=1/32$. Nonetheless, the quantitative results of the entrained volume of the first fluid into the second fluid 
$(V_e)$ is in close agreement with the experimental results, and the particle settling velocity is accurately captured. We note that in this case study, no contact line forms between the solid and the background fluids, such that an oil film coats the solid surface at all times, which is in agreement with the results found in the experiments of Pierson & Magnaudet (Reference Pierson and Magnaudet2018a).
Figure 16. This figure showcases the experimental results of Pierson & Magnaudet (Reference Pierson and Magnaudet2018a,Reference Pierson and Magnaudetb) for a spherical particle crossing an immiscible interface under different configurations. The plots showcase the performance of this numerical method compared with the experimental results in estimating the entrained fluid volume 
$V_e$ and the particle velocity. Qualitative results from the simulations are also provided.
Open access
