2. Statistical model

The collision frequency $K$ of particles with a bubble is determined by $K = \Gamma n_p$ , where $n_p$ is the particle number density, and the proportionality coefficient $\Gamma$ ( $\rm m^3\,s^{-1}$ ) is commonly referred to as the collision kernel and is determined by the flow (Pumir & Wilkinson Reference Pumir and Wilkinson2016). In quiescent flow, the collisions between a spherical bubble of radius $r_b$ rising at constant velocity $U_b$ with small inertial solid particles are deterministic. The particle grazing trajectory $\Psi _c$ determines a ‘collision tube’, where all the particles contained inside collide with the bubble, as illustrated in figure 1(a). In this case, the collision frequency is determined by the collision number flux, i.e. $K=Q_cn_p$ . The collision number flux $Q_c$ can be determined by the cross-sectional area of radius $r_c$ limited by the grazing trajectories upstream far from the bubble: $Q_c=\pi r_c^2U_b$ . Therefore, the collision kernel in quiescent flow can be expressed as $\Gamma _q=Q_c=\pi r_b^2U_bE_c$ , where the collision efficiency $E_c=r_c^2/r_b^2$ measures the ratio of collided particle number to the total encountering particle number. Notice here that $E_c$ depends on the bubble Reynolds number $Re_b=2r_bU_b/\nu$ , the ratio of the particle radius $r_p$ and $r_b$ , and the Stokes number $St_p = \tau _p/\tau _f$ , with $\tau _p$ the particle response time, and $\tau _f = 2r_b/U_b$ the time scale of the bubble–particle interaction (Dai, Fornasiero & Ralston Reference Dai, Fornasiero and Ralston2000; Sarrot et al. Reference Sarrot, Guiraud and Legendre2005; Huang et al. Reference Huang, Legendre and Guiraud2012).

Figure 1. (a) Sketch of the grazing trajectory (red dashed line) in quiescent flow. The shaded region indicates the collision tube where all particles collide on the bubble. (b) Sketch of the bubble–particle collision model under temporary bubble slip velocity $U'_{b}$ . (c) Mean flow streamlines around the bubble for the case of imposed velocity bubble (solid) and quiescent flow (dashed lines) at $\overline {Re_b}=120$ . (d) Trajectories of colliding particles ( $r_p/r_b=0.05,\ St_p=0.04$ ) for the case of imposed velocity bubble compared to the corresponding grazing trajectories (red lines) in quiescent flow at the same $Re_b =120$ . (e) Snapshot of the bubble–particle collision process in turbulence for the imposed-velocity bubble with flow from left to right in the bubble frame of reference. Incoming particles that end up colliding with the bubble are marked in red.

To provide an understanding of the relevant physical mechanisms of bubble–particle collision in turbulence, we propose a statistical model. We start from the assumption that the flow field in the vicinity of the bubble is approximately stationary and uniform during the bubble–particle interaction. Conceptually, this implies that the temporal scale $\tau _f$ is (much) shorter than the correlation time scale of flow fluctuations ( $O(\tau _\eta )$ ), but it remains to be tested from the data under what conditions exactly this assumption is valid in practice. Additionally, we consider the flow correlation length scale to be comparable to or larger than the bubble size. Within this ‘frozen turbulence’ assumption, the instantaneous bubble–particle collision process in turbulence can be related back to that observed in quiescent flow. Therefore, the equivalent steady flow problem is characterised by the magnitude of the instantaneous bubble slip velocity $U'_{b}$ with corresponding values $Re'_{b} = 2r_bU'_{b}/\nu$ and $St'_{p}=\tau _p/(2r_b/U'_{b})$ , as shown in figure 1(b). The entire bubble–particle collision process in a turbulent flow can then be viewed as a superposition of collision events in quiescent flow with varying parameters $Re'_{b}$ and $St'_{p}$ . In the framework, the collision number $N_c|_{\Delta \tau }$ during a short period $\Delta \tau$ ( $O(\tau _\eta )$ ) in turbulent flow, with mean particle number density $n_p$ , is given by $N_c|_{\Delta \tau }=\pi r_b^2E_c(Re'_{b},St'_{p})U'_{b}n'_p\,\Delta \tau$ , where $n'_{p}$ is the instantaneous encountering particle number density, which depends on $Re'_{b}$ and $St'_{p}$ . Importantly, $E_c(Re'_{b},\ St'_{p})$ denotes the collision efficiency in quiescent flow for varying flow parameters, which is deterministic and can be parametrised. Then the total expected collision number $N_c$ during a sufficiently long time period $T=\sum \Delta \tau$ is given by

(2.1) \begin{equation} N_c = \pi r_b^2T\int E_c(Re'_{b},\ St'_{p})\,U'_{b}f(U'_{b})n'_{p}\,{\rm d}U'_{b}. \end{equation}

Based on this, the collision kernel can be derived as

(2.2) \begin{equation} \Gamma = \frac {K}{n_p}={\frac {\pi r_b^2}{n_p}\int E_c(Re'_{b},\ St'_{p})\,U'_{b}\,f(U'_{b})\,n'_{p}\,\textrm {d} U'_{b}}. \end{equation}

Note that collision efficiency has been studied extensively, and the general dependence of $E_c$ on the bubble Reynolds number and particle Stokes number has been established in the literature (Schulze Reference Schulze1989; Sarrot, Guiraud & Legendre Reference Sarrot, Guiraud and Legendre2005; Huang et al. Reference Huang, Legendre and Guiraud2012). The collision kernel is thus determined by the probability density function (PDF) $f(U'_{b})$ of the bubble slip velocity $U'_{b}$ , the modelling of which remains an open question. Furthermore, we note that the collision kernel will also be affected by the preferential sampling for particles with large $St_p$ , as this alters the incoming particle number density.

3. Numerical methods and simulations

We perform interface-resolved direct numerical simulations (DNS) to test this model. The turbulent flow is governed by the incompressible Navier–Stokes equations, which read

(3.1) \begin{eqnarray} \frac {\partial \boldsymbol{u}}{\partial t} + \boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla } \boldsymbol{u} &=& - \frac {1}{\rho _f}\,\boldsymbol{\nabla } p + \nu \, \nabla ^2 \boldsymbol{u} + \boldsymbol{f} + \boldsymbol{f}_b, \end{eqnarray}

(3.2) \begin{eqnarray} \boldsymbol{\nabla } \boldsymbol{\cdot} \boldsymbol{u} &=& 0, \end{eqnarray}

where $\boldsymbol{u}$ is the fluid velocity, $p$ denotes the pressure, and the parameters are the kinematic viscosity $\nu$ and the reference liquid density $\rho _f$ . The vector $\boldsymbol{f}$ denotes an external random large-scale volume force, which is statistically homogeneous and isotropic, with constant-in-time global energy input (Perlekar et al. Reference Perlekar, Biferale, Sbragaglia, Srivastava and Toschi2012). This force is used to generate and maintain the turbulent flow. The bubble-free turbulent flow intensity is characterised by the Reynolds number based on the Taylor microscale, $Re_\lambda =\sqrt {15u'/(\nu \epsilon )}$ , where $u'$ denotes the root mean square velocity of the turbulence, and $\epsilon$ is the mean energy dissipation rate. The vector $\boldsymbol{f}_b$ accounts for the bubble–fluid two-way coupling.

Table 1. Parameter of the numerical simulations and relevant turbulence scales: $Re_\lambda$ is the Taylor–Reynolds number; $\eta =(\nu ^3/\epsilon )^{1/4}$ is the Kolmogorov dissipation length scale in grid space units $\Delta x$ ; $\tau _\eta$ is the Kolmogorov time scale in time-step units $\Delta t$ ; $L=u'^3/\epsilon$ is the integral scale; $T_L=L/u'$ is the large-eddy turnover time; $\lambda = u'\sqrt {15 \nu / \epsilon }$ is the Taylor micro-scale; $\overline {Re_b}=2r_b\overline {U_b}/\nu$ is the bubble Reynolds number based on the mean rising velocity $\overline {U_b}$ ; $T_i=u'/\overline {U_b}$ is the turbulent intensity.

In the practice of flotation, commonly a large amount of surfactants is present in the liquid. Therefore, it is reasonable to assume that the bubbles with a typical Reynolds number below 200 are fully contaminated and approximately spherical, resulting in a boundary condition that is nearly no-slip (Nguyen & Schulze Reference Nguyen and Schulze2003; Huang et al. Reference Huang, Legendre and Guiraud2012). The Weber number $We=2r_b\rho _f U_b^2/\chi$ based on the surface tension $\chi$ of water and the bubble rise velocity is $O(0.1)$ for our simulations, and even lower based on turbulent fluctuations. Therefore, bubble deformations remain negligible even if surface tension is lowered by the presence of surfactants. Under these conditions, the bubble behaves similarly to buoyant spheres. In this case, the translation and rotation of the bubble are governed by the Newton–Euler equations

(3.3) \begin{eqnarray} m_b\, \frac {{\rm d} \boldsymbol{U}_b}{\textrm {d}t} &=& \oint _{S_b} \boldsymbol {\sigma } \boldsymbol{\cdot} \boldsymbol{n}\, \textrm {d}S + V_b(\rho _b-\rho _f)\boldsymbol{g}, \end{eqnarray}

(3.4) \begin{eqnarray} \frac {\textrm {d} \boldsymbol {{I}} \boldsymbol{\Omega }}{\textrm {d}t} &=& \oint _{S_b} (\boldsymbol{x}-\boldsymbol{x}_b) \times (\boldsymbol {\sigma } \boldsymbol{\cdot} \boldsymbol{n})\, \textrm {d}S, \end{eqnarray}

where $\boldsymbol{U}_b(t)=\textrm {d}\boldsymbol{x}_b/\textrm {d}t$ and $\mathbf {\Omega }(t)$ , respectively, are the bubble velocity and angular velocity vectors at position $\boldsymbol{x}_b(t)$ . The bubble mass is given by $m_b=\rho _b V_b$ (where $\rho _b$ is the bubble density, and $V_b$ is the volume), and $\boldsymbol {{I}}$ is the moment of inertia tensor. Here, $\boldsymbol {\sigma } = -p \boldsymbol {I}+\rho \nu (\boldsymbol{\nabla} \boldsymbol{u} + \boldsymbol{\nabla} \boldsymbol{u}^{\rm T})$ is the fluid stress tensor, $\boldsymbol{x}-\boldsymbol{x}_b$ is the position vector relative to the bubble centre, and $\boldsymbol{n}$ is the outward-pointing normal to the bubble surface $S_b$ . This term $\boldsymbol {\sigma }$ is solved by the immersed boundary method (IBM) to account for the two-way coupling between the flow and the bubble. In the last decades, the IBM has been used widely for studying the multi-phase flow (Peskin Reference Peskin2002; Mittal & Iaccarino Reference Mittal and Iaccarino2005; Uhlmann Reference Uhlmann2005). In the IBM, the Euler grids of the fluid are fixed, and the surface of the bubble is represented by Lagrangian nodes that move with the bubble motion. To avoid the formation of a mean flow due to the buoyancy driving exerted by the bubble, we compensate the average force applied from the bubble to the liquid to attain a statistically stationary state (Höfler & Schwarzer Reference Höfler and Schwarzer2000; Chouippe & Uhlmann Reference Chouippe and Uhlmann2015). To achieve this, the spatial average of the IBM volume force term needs to be subtracted. More explicitly, we compute the average at each simulation time step:

(3.5) \begin{eqnarray} \langle \boldsymbol{f}^{ibm}\rangle _\Omega (t) = \frac {1}{\|\Omega \|}\int _\Omega \boldsymbol{f}^{ibm}(\boldsymbol{x},t)\, \rm {d}\boldsymbol{x}, \end{eqnarray}

where $\Omega$ denotes the entire computational domain, $\langle \cdots \rangle _\Omega$ indicates the spatial average over the entire $\Omega$ region, and $\|\Omega \|$ is the volume of the $\Omega$ region. Then the bubble-related contribution to the volume force, $\boldsymbol{f}_b$ , is obtained from

(3.6) \begin{eqnarray} \boldsymbol{f}_b(\boldsymbol{x},t)=\boldsymbol{f}^{ibm}(\boldsymbol{x},t)-\langle \boldsymbol{f}^{ibm}\rangle _\Omega (t). \end{eqnarray}

We adopt an implicit method (Tschisgale, Kempe & Fröhlich Reference Tschisgale, Kempe and Fröhlich2017) to solve the bubble dynamics as the conventional explicit method is numerically unstable when the bubble is light (Uhlmann Reference Uhlmann2005). To avoid the build-up of momentum in the triply periodic domain over time, the force applied by the bubble on the liquid is compensated to attain a statistically stationary state (Chouippe & Uhlmann Reference Chouippe and Uhlmann2015).

We use a code based on the lattice Boltzmann method to solve the Navier–Stokes equations (Calzavarini Reference Calzavarini2019; Jiang et al. Reference Jiang, Wang, Liu, Sun and Calzavarini2022). Two sets of simulations are carried out: (1) a simplified ideal configuration where a constant bubble velocity is imposed, and (2) simulations with a freely rising bubble, where the bubble–fluid density ratio is $\rho _b/\rho _f=10^{-3}$ . Matching practically relevant conditions (Wang et al. Reference Wang, Hoque, Evans and Mitra2022), we consider a moderately turbulent flow ( $Re_\lambda =32$ and 64) and a bubble Reynolds number $Re_b\sim O(100)$ . Detailed parameters of the simulations are summarised in table 1. There are two computational conditions that should be satisfied in the simulations. First, the wake of the fully contaminated bubble presents a steady axisymmetric vortex at such a high $Re_b$ in a quiescent flow (Johnson & Patel Reference Johnson and Patel1999; Sarrot et al. Reference Sarrot, Guiraud and Legendre2005). As the flow domain is periodic, the incoming flow ahead of the bubble might be disturbed by the remnants of this bubble wake. To avoid this issue, the flow domain in the bubble rising direction should be sufficiently large. The whole flow domain adopted in this work is rectangular, with uniform grid sizes $N_x \times N_y\times N_z=3072\times 512\times 512$ , where the bubble rises in the $x$ -direction. We found that using this domain length was sufficient to avoid wake effect. Second, the number of lattices within the boundary layer of the bubble should be sufficient. The approximate boundary layer thickness is estimated as $\delta =1.13/Re_b^{1/2}$ (Johnson & Patel Reference Johnson and Patel1999), where $\delta$ denotes the boundary layer thickness normalised by $2r_b$ . Additionally, the minimum particle size ( $r_p/r_b=0.025$ ) should be larger than the lattice unit to ensure the accuracy of interpolation when the particle is close to the bubble. To this end, the bubble diameter is resolved by 80 grids to make sure that the boundary layer is well resolved, as well as the minimum distance of the point-like particle to the bubble being larger than one grid unit. The ratio of turbulent dissipation time scale $\tau _\eta$ to $\tau _f$ in our simulations is close to 1.2. We note that the typical auto-correlation time scale of the turbulent flow is longer than $\tau _\eta$ , which implies the validity of the assumption of frozen turbulence in our model. In the flotation process, the mineral particles are significantly smaller than the turbulent dissipation scale $\eta$ . We represent these as one-way coupled point particles, considering the non-Stokesian drag force and the added mass force. We do not account for potential alterations to the drag force when a particle is close to the bubble surface. The particle response time $\tau _p=r_p^2/(3\beta \nu )$ includes the density coefficient $\beta =3\rho _f/(\rho _f+2\rho _p)$ , where $\rho _p$ denote the density of the particle. The particle dynamics is driven by the instantaneous turbulent flow. The collision frequency $K=N_c/T$ , and thus the collision kernel $\Gamma = K/n_p$ , is measured by counting the number of particles ( $N_c$ ) that collide with the bubble over a long time period $T$ . A collision is detected when the distance between the particle and bubble is smaller than $r_p+r_b$ . Additionally, we conduct simulations of bubble–particle collisions in a quiescent flow at $Re_b$ ranging from 80 to 210 to obtain the dependence of $E_c$ on $Re_b$ and $St_p$ . In these simulations, a constant-velocity inflow is applied at the inlet, while a homogeneous Neumann boundary condition is imposed at the outlet. The bubble is fixed at the centre of the domain, and the spatial resolution is kept the same as for the turbulent flow case and the turbulence forcing term $\boldsymbol{f}$ is switched off.