Bubble–particle collisions in turbulence: insights from point-particle simulations

Abstract

Bubble–particle collisions in turbulence are central to a variety of processes such as froth flotation. Despite their importance, details of the collision process have not received much attention yet. This is compounded by the sometimes counter-intuitive behaviour of bubbles and particles in turbulence, as exemplified by the fact that they segregate in space. Although bubble–particle relative behaviour is fundamentally different from that of identical particles, the existing theoretical models are nearly all extensions of theories for particle–particle collisions in turbulence. The adequacy of these theories has yet to be assessed as appropriate data remain scarce to date. In this investigation, we study the geometric collision rate by means of direct numerical simulations of bubble–particle collisions in homogeneous isotropic turbulence using the point-particle approach over a range of the relevant parameters, including the Stokes and Reynolds numbers. We analyse the spatial distribution of bubble and particles, and quantify to what extent their segregation reduces the collision rate. This effect is countered by increased approach velocities for bubble–particle compared to monodisperse pairs, which we relate to the difference in how bubbles and particles respond to fluid accelerations. We found that in the investigated parameter range, these collision statistics are not altered significantly by the inclusion of a lift force or different drag parametrisations, or when assuming infinite particle density. Furthermore, we critically examine existing models and discuss inconsistencies therein that contribute to the discrepancy.

1. Introduction

Collisions between bubbles and particles in a turbulent flow are of significant technological relevance. In particular, such collisions are essential to the flotation process, which is a widely used separation technique especially in the mining industry (Nguyen & Schulze 2004). In this process, after grinding, small ore fragments are fed into a big water-filled cell that is agitated via a rotor and into which air bubbles are injected. Collisions between the ore particles and the bubbles then form the base for the decisive test: valuable mineral particles attach to the bubbles due to their hydrophobic surface and consequently rise to the top where they can be skimmed off as a froth, whereas the hydrophilic waste rock particles remain in suspension and are eventually discharged as tailings. This technology is already applied at staggering scales (Nguyen & Schulze (2004) estimated that a total of 2 billion tons of ore are treated annually). Given especially the relevance in the mining of copper, and in view of the strong push for electrification in response the the climate crisis, these numbers are likely to continue to rise in the future (Rogich & Matos 2008; World Bank Group 2017). The interest to understand the collision process better is driven by the demand for more reliable process modelling (Kostoglou, Karapantsios & Evgenidis 2020a) but also by the need for performance improvement. The latter is especially a concern for small particles with diameters smaller than $20\,\mathrm {\mu }{\rm m}$, where recovery is poor owing to their low collision rates (Nguyen, George & Jameson 2006; Miettinen, Ralston & Fornasiero 2010).

For modelling purposes, the collision process is generally separated into two components (Pumir & Wilkinson 2016): the ‘geometric collision rate’, which considers the collisions neglecting any interaction between the collision partners, and the ‘collision efficiency’, which quantifies how many of these collisions actually happen when taking the local modification of the flow field into account. The focus here is on the (ensemble-averaged) geometric collision rate between two species, ‘1’ and ‘2’, which when expressed per unit volume can be written as

(1.1)\begin{equation} Z_{12} = \varGamma_{12} n_1 n_2, \end{equation}

where $\varGamma _{12}$ is the collision kernel, and $n_1$ and $n_2$ denote the respective number densities of the two species. The collision kernel measures the rate at which the separation vector between particle centres crosses the collision distance. For spherical particles with a collision distance $r_c = r_1+ r_2$ (where $r_1,r_2$ denote the particle radii), $\varGamma _{12}$ can be expressed as (Sundaram & Collins 1997)

(1.2)\begin{equation} \varGamma_{12} = 4{\rm \pi} r_c^2\,g(r_c)\,S_-(r_c), \end{equation}

where besides the surface area $4{\rm \pi} r_c^2$ of the collision sphere, the other factors are the radial distribution function (RDF) at collision distance $g(r_c)$, which describes variations in the local particle concentration, and the effective radial approach velocity at contact,

(1.3)\begin{equation} S_-(r_c) ={-}\int_{-\infty}^{0} \Delta v_r\mathrm{p.d.f.}(\Delta v_r|r_c)\,\mathrm{d}(\Delta v_r). \end{equation}

Here, $\Delta v_r$ is the radial component of the relative velocity, which is positive when the pair separates, and $\mathrm {p.d.f.}(\Delta v_r|r_c)$ is the probability density function of $\Delta v_r$ conditioned on a pair with separation $r_c$.

The bidisperse collision kernel $\varGamma _{12}$ depends on a multitude of parameters characterising properties of the suspended particles and of the carrier flow. The discussion here is restricted to homogeneous isotropic turbulence for which the most relevant dependencies may be summarised in non-dimensional form as

(1.4)\begin{equation} \frac{\varGamma_{12}}{r_c^3/\tau_\eta} = f\left(St_1,St_2,\frac{\rho_1}{\rho_f},\frac{\rho_2}{\rho_f},Re_\lambda, \frac{r_1}{\eta},\frac{r_2}{\eta},Fr,\ldots \right). \end{equation}

Here, the Stokes number

(1.5)\begin{equation} St_i = \frac{\tau_i}{\tau_\eta} = \frac{r_i^2(2\rho_i/\rho_f + 1)}{9\nu\tau_\eta}, \end{equation}

($i = 1,2$) characterises how well particles follow the flow by relating the particle response time $\tau _i = {r_i^2(2\rho _i/\rho _f + 1)}{/(9\nu )}$ to the Kolmogorov time scale of the turbulence $\tau _\eta = (\nu /\varepsilon )^{1/2}$, with $\nu$ and $\varepsilon$ denoting the kinematic viscosity and the average rate of turbulent dissipation, respectively. Further, the ratios of particle ($\rho _i$) and fluid ($\rho _f$) densities are relevant as they characterise to what extent the particle motion is influenced by fluid (‘added mass’) inertia and buoyancy. The particle radii $r_i$ determine the collision radius $r_c$, and their size relative to the Kolmogorov length scale $\eta = (\nu ^3/\varepsilon )^{1/4}$ determines the range of turbulent scales relevant for their motion. In addition to turbulent driving, particle motion may also be affected by gravitational effects, and the relative importance of these two factors is captured by the Froude number $Fr = a_\eta /g$, where $a_\eta = \eta /\tau _\eta ^2$ and $g$ are the turbulence and gravitational accelerations, respectively. Finally, the intensity of the turbulence is measured by the Taylor Reynolds number $Re_\lambda = \sqrt {15/(\nu \varepsilon )}\,u'^2$, where $u'$ is the single-component root-mean-square (r.m.s.) fluid velocity. Obviously, the entire parameter space spanned by (1.4) is vast and cannot be studied comprehensively here. We therefore limit the present investigation to cases with $St_1 = St_2 = St$, and to the zero-gravity regime, i.e. $Fr \to \infty$. The benefit of these choices lies in the fact that they keep the problem simple enough to disentangle the relevant mechanisms. Similarly, these configurations are the most amenable to modelling approaches and therefore allow for their evaluation at the most basic level.

Our investigation is based on direct numerical simulations of bubbles and particles in statistically stationary homogeneous isotropic turbulence using the point-particle approach. Details of this approach will be described in § 3 after first reviewing relevant modelling approaches for the collision kernel in § 2. The results are shown in § 4, followed in § 5 by practical considerations in light of our results, and conclusions.

2. Theoretical background and existing models

2.1. The tracer limit $St \rightarrow 0$: shear mechanism

In the tracer limit of $St \to 0$, the suspended species follow the flow faithfully and distribute uniformly. This means that collisions occur only if particles of finite size are moved relative to each other due to shearing motions in the flow. Considering the dominant shear contribution of the smallest (Kolmogorov) scales of turbulence, and assuming local isotropy as well as Gaussian distribution of the flow velocity gradient, Saffman & Turner (1956) derived the classical result

(2.1)\begin{equation} \varGamma_{12}^{(ST)} = \sqrt{\frac{8{\rm \pi}}{15}}\,\frac{r_c^3}{\tau_\eta}, \end{equation}

predicting the rate of shear-driven collisions in a turbulent flow. Note that here we use the spherical formulation for the collision kernel, which was shown to be the appropriate form by Wang, Ayala & Xue (2005). Employing the concept of a collision cylinder, instead, results in a slightly different value of the prefactor ($\sqrt {8{\rm \pi} }/3\approx 1.67$ instead of $\sqrt {8{\rm \pi} /15} \approx 1.29$).

2.2. Intermediate $St$: preferential sampling and velocity decorrelation

For non-zero $St$, the suspended species no longer completely follow the flow. Such inertial effects influence the collision rate via two different pathways. First, even if the drift is small, its accumulated effect leads to preferential concentration, inducing clustering in the particle field. Additionally, the increasing decorrelation between the local particle and fluid velocities affects the collision velocities.

Preferential concentration of inertial particles is widely observed experimentally (Aliseda et al. 2002; Monchaux, Bourgoin & Cartellier 2010; Obligado et al. 2014; Petersen, Baker & Coletti 2019; Li et al. 2021) and numerically (Bec et al. 2007; Goto & Vassilicos 2006; Calzavarini et al. 2008b; Voßkuhle et al. 2014; Ireland, Bragg & Collins 2016). It is rather straightforward to extend (2.1) to account for this effect by simply multiplying it with the RDF $g_{12}(r_c)$ (Voßkuhle et al. 2014), yielding

(2.2)\begin{equation} \varGamma_{12}^{(STc)} = \sqrt{\frac{8{\rm \pi}}{15}}\,\frac{r_c^3}{\tau_\eta}\,g_{12}(r_c). \end{equation}

Various approaches have been proposed to explain the phenomenon of preferential concentration (Maxey 1987; Bec et al. 2005, 2007; Chen, Goto & Vassilicos 2006; Goto & Vassilicos 2006; Coleman & Vassilicos 2009; Zaichik & Alipchenkov 2009; Fouxon 2012). Among these, the most intuitive one is the ‘centrifuge picture’ (Maxey 1987), according to which heavy particles are ejected out of eddies due to their inertia and hence accumulate in regions of low vorticity and high strain. For collisions between heavy particles (e.g. cloud droplets), it is therefore found that $g_{12}(r_c)\geq 1$ generally (Zhou, Wexler & Wang 2001), i.e. clustering enhances the collision rate in some cases even by multiple orders of magnitude (Voßkuhle et al. 2014; Ireland et al. 2016; Pumir & Wilkinson 2016).

For the collision velocity, inertial effects can be split into a local mechanism and a non-local mechanism. The local mechanism results from the fact that particles react differently to the same fluid forcing provided that they have different properties. Hence this effect contributes an additional relative velocity for bidispersed collisions only, and plays no role in monodispersed cases. For this to occur, the particle trajectories should not deviate significantly from the pathlines (i.e. $St$ not too large). An extension of the Saffman–Turner approach accounting for this effect has been reported by Yuu (1984).

In contrast, the non-local mechanism refers to the situation in which particles arrive at the same location with different particle velocities. Illustratively, one can think of these particles as being ‘slung out’ of neighbouring eddies, and the effect is therefore also known as the ‘sling effect’ (Falkovich, Fouxon & Stepanov 2002). Unlike the local mechanism, the sling effect is also active for monodisperse suspensions, and so far has been studied almost exclusively in this context (e.g. Falkovich et al. 2002; Wilkinson, Mehlig & Bezuglyy 2006; Falkovich & Pumir 2007; Ijzermans, Meneguz & Reeks 2010; Bewley, Saw & Bodenschatz 2013; Voßkuhle et al. 2014). From these studies, it has become clear (see e.g. Pumir & Wilkinson 2016, for an overview) that the sling effect can significantly enhance monodisperse collision rates by increasing $S_-$. A widely used parametrisation is $S_- \sim u_\eta \,F(St,Re_\lambda )$, where $u_\eta = \eta /\tau _\eta$ is the Kolmogorov velocity, such that the sling-induced collision rate is given by

(2.3)\begin{equation} \varGamma_{11}^{(slg)} = 4{\rm \pi} r_c^2 u_\eta\,F(St,Re_\lambda). \end{equation}

It has then been proposed (Wilkinson et al. 2006; Voßkuhle et al. 2014) to obtain the overall collision rate from the sum

(2.4)\begin{equation} \varGamma_{11}^{(tot)} = \varGamma_{11}^{(STc)} + \varGamma_{11}^{(slg)}. \end{equation}

The underlying idea for this decomposition is that particles that are clustered close to each other collide with low velocities, whereas those with high relative velocities can be assumed to be uniformly distributed. Note that both (2.3) and (2.4) are given for the monodisperse case only, since related results for the bidisperse case have not been reported yet.

2.3. Large $St$ limit: kinetic gas behaviour

At very large (but finite) $St$, the velocities of the suspended particles arriving at the same point are increasingly uncorrelated. Assuming entirely random and isotropic particle velocities in the spirit of the kinetic gas theory, Abrahamson (1975) derived the collision kernel

(2.5)\begin{equation} \varGamma_{12}^{(A)} = \sqrt{8{\rm \pi}}\,r_c^2\,\sqrt{v_1'^2 + v_2'^2}, \end{equation}

where the mean-square single-component particle velocity $v_i'^2$ is related to flow properties via

(2.6)\begin{equation} v_i'^2 = A_iu'^2 = \frac{\dfrac{T_{fL}}{\tau_i}+\gamma_{i}^2}{\dfrac{T_{fL}}{\tau_i}+1}\,u'^2, \end{equation}

with $T_{fL}$ denoting the fluid Lagrangian integral time scale, and $\gamma _i=3\rho _f/(2\rho _i + \rho _f)$. It has been pointed out (Voßkuhle et al. 2014) that (2.5) is not strictly valid for turbulence as it fails to account for the multiscale structure of the flow. Alternative derivations (Völk et al. 1980; Mehlig, Uski & Wilkinson 2007) based on the Kolmogorov (1941) phenomenology arrived at $F(St,\infty ) \sim K\sqrt {St}$, with $K$ being a universal dimensionless constant, at the limit of intense turbulence in the context of (2.3).

2.4. Modelling approaches for bubble–particle collisions in turbulence

Next, we outline briefly the different approaches to modelling bubble–particle collisions reported in particular in the mining literature so far. For additional details, we refer the reader to the reviews on the topic (e.g. Nguyen et al. 2016; Hassanzadeh et al. 2018; Kostoglou, Karapantsios & Oikonomidou 2020b). Note that here and in the following we use subscripts $i = b,p$ to denote bubbles and heavy particles, respectively.

The most commonly adopted approach is to assume the high-$St$ limit and to base the collision models on (2.5). The differences from the theory of Abrahamson (1975) are related to the expressions used to model the r.m.s. velocities. Instead of (2.6), Schubert (1999) and later Bloom & Heindel (2002) used the relation given by Liepe & Möckel (1976), which reads

(2.7)\begin{equation} w_i'^{(LM)} = 0.57\,\frac{\varepsilon^{4/9}r_i^{7/9}}{\nu^{1/3}}\left(\frac{\rho_f-\rho_i}{\rho_f}\right)^{2/3}. \end{equation}

This result is based on an analogy to gravitational settling with the fluid acceleration in the inertial range replacing gravity. The resulting apparent weight is balanced by Allen's drag, which scales with the particle slip velocity $\boldsymbol {w}_i = \boldsymbol {v}_i - \boldsymbol {u}$ as $|\boldsymbol {w}_i|^{3/2}$ (with boldface denoting vectors). Here, $\boldsymbol {v}_i$ is the bubble/particle velocity, and $\boldsymbol {u}$ is the flow velocity. Later work by Nguyen & Schulze (2004) used (2.7) with a different constant, 0.83, for bubbles. For particles, they replaced the inertial subrange acceleration with that in the dissipation range, and Allen's drag with Stokes drag, in order to account for the small size of typical particles. This resulted in

(2.8)\begin{equation} w_p'^{(NS)} = \frac{2r_p^3\varepsilon}{135\nu^2}\left(\frac{\rho_p-\rho_f}{\rho_f}\right). \end{equation}

Importantly, both (2.7) and (2.8) are expressions for the relative velocity between bubbles/particles and the surrounding fluid. Their use with Abrahamson's theory is therefore inconsistent since (2.5) contains velocities in a fixed frame of reference. This was already pointed out in Kostoglou et al. (2020a). We further note that the quasi-static assumption underlying (2.7) and (2.8) is equivalent to a low-$St$ approximation and therefore not valid for the high-$St$ limit in which (2.5) applies. In fact, (2.8) is consistent with the rigorously derived small-$St$ limit (Fouxon 2012)

(2.9)\begin{equation} \boldsymbol{w}_i=\beta_i\tau_i \boldsymbol{a}_f, \end{equation}

with

(2.10)\begin{equation} \beta_i = \frac{2(\rho_f - \rho_i)}{\rho_f + 2\rho_i}, \end{equation}

if the fluid acceleration is approximated by the dissipative scaling $|\boldsymbol {a}_f| \approx \varepsilon r_i/(15 \nu )$. Note, however, that it appears more appropriate to use $|\boldsymbol {a}_f| \approx a_\eta$ for small particles ($r_i/\eta \lessapprox 1$) because $|\boldsymbol {a}_f| \rightarrow 0$ otherwise. In the general case, the fact that the fluid acceleration experienced by the particle may vary considerably over the particle response time precludes a simple relation between $\boldsymbol {w}_i$ and $\boldsymbol {a}_f$.

In a different approach, Ngo-Cong, Nguyen & Tran-Cong (2018) extended the model by Yuu (1984) to the bubble–particle case. The resulting expression takes the form

(2.11)\begin{equation} \varGamma_{bp}^{(NC)} = (r_b + r_p)^2\sqrt{\frac{8{\rm \pi}}{3}\left\{3(A_b + A_p - 2B)u'^2 + (A_br_b^2 + A_pr_p^2 + 2Br_pr_b)\,\frac{\varepsilon}{3\nu}\right\}}, \end{equation}

where the term proportional to $u'^2$ represents the ‘local’ inertial effect on the relative velocity that adds to the shear-driven collisions that are accounted for by the $\varepsilon$ term. Aside from bubble/particle properties, the coefficients $A_i$ and $B$ depend also on $T_{fL}$, similar to (2.6). The bubble–particle velocity correlation is determined by $B$. Ngo-Cong et al. (2018) additionally incorporated (2.7) and (2.8) into this model. However, doing so suffers from the same inconsistencies outlined above. As a consequence, this expression featured a negative radicand when evaluated for the parameters in this study, and is therefore not included further. It is worth mentioning that the term proportional to $u'^2$ in (2.11) does not approach (2.5) even at large $St$, since $A_i=B$ for identical particles. We have therefore extended the model along the lines of the theory of Kruis & Kusters (1997), which attempted to reconcile this issue by marrying the concept of Yuu (1984) for the small $St$ case and that of Williams & Crane (1983) (which does not consider shear-driven collisions) for the large $St$ case, to bubble–particle collisions (see Appendix A). The original formulation of Kruis & Kusters (1997) accounted for collisions of particles of equal density only. We further note the work of Fayed & Ragab (2013), who employed the model of Zaichik, Simonin & Alipchenkov (2010), which was developed originally for collisions between particles with arbitrary but equal density assuming a joint-normal fluid–particle velocity distribution, to the bubble–particle case. Recently, Kostoglou et al. (2020a) proposed another model that is specific to the case of very fine particles that essentially follow the flow, such that the relative velocity is dominated by the bubble slip velocity.

A common feature of almost all these models is that they are direct adaptations of concepts developed for collisions between heavy particles. They therefore fail to account for fundamental differences in how bubbles and heavy particles react to a turbulent flow. Most strikingly, for example, the response to fluid accelerations given in (2.9) is in opposite directions as $\beta$ switches sign from negative to positive from $\rho _i>\rho _f$ to $\rho _i<\rho _f$. Similarly, applying the centrifuge picture to light particles such as bubbles, one expects them to concentrate in regions of high vorticity as they travel to the centre of the eddies. This implies that heavy and light particles segregate in a turbulent flow, as has indeed been observed (Calzavarini et al. 2008a; Fayed & Ragab 2013). As a consequence, it is expected that $g_{bp}(r_c)< 1$ in these cases, such that preferential concentration is expected to lead to a decrease in $\varGamma _{bp}^{(STc)}$ for collisions between heavy and light particles. Such aspects cannot be quantified easily from laboratory-scale flotation set-ups (e.g. Darabi et al. 2019) since it is difficult to disentangle the many factors influencing the overall flotation rate. Numerically, Wan et al. (2020) observed no reduction of the RDF, but increased relative velocities for bubble–particle pairs based on a point-particle simulation while taking gravitational effects into account. For cases without gravity, Fayed & Ragab (2013) reported $g_{bp}(r_c) < 1$ and relative velocities matching those predicted by the theory of Zaichik et al. (2010) for simulations of the extreme case $\rho _b/\rho _f = 0$ and $\rho _p/\rho _f\to \infty$ across a range of $\tau _i$. It is the goal of the present study to add to this a systematic investigation of how bubble–particle relative behaviour affects their collision statistics in the low- to moderate-$St$ regime, and to assess how this affects modelling outcomes.

3. Methods

3.1. Fluid phase

To obtain the background turbulence, we solve the Navier–Stokes equation and the continuity equation for incompressible flow:

(3.1a)$$\begin{gather} \frac{\mathrm{D}\boldsymbol{u}}{\mathrm{D}t} ={-}\frac{1}{\rho_f}\,\boldsymbol{\nabla} P + \nu\, \nabla^2\boldsymbol{u} + \boldsymbol{f}_\varPsi, \end{gather}$$

(3.1b)$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} = 0, \end{gather}$$

where $\mathrm {D}/\mathrm {D}t$ is the material derivative following a fluid element, $t$ is the time, and $P$ is the pressure. The forcing $\boldsymbol {f}_\varPsi$, which is non-zero only for the wavenumbers $|\boldsymbol {\kappa }|/|\boldsymbol {\kappa }_0| < 2.3$ (i.e. the largest scales), with $|\boldsymbol {\kappa }_0|$ being the smallest wavenumber along each direction, is added to counter dissipation and maintain statistical stationarity. We employ the widely used Eswaran & Pope (1988) forcing scheme (Chouippe & Uhlmann 2015; Spandan et al. 2020). In brief, a complex vector is generated in Fourier space for the forced wavenumbers by multiple Uhlenbeck–Ornstein processes. This vector is then projected onto the plane normal to the wavevector, thereby ensuring that $\boldsymbol {f}_\varPsi$ is divergence-free. To simulate the fluid motion, a second-order finite-difference solver is implemented on a staggered grid (Verzicco & Orlandi 1996; van der Poel et al. 2015). All spatial derivatives, including the nonlinear terms, are discretised by a second-order central finite-difference method. The simulation domain is a cubic box with length $L_{box}=1$ and triply periodic boundary conditions to eliminate boundary effects. Time marching is performed using a fractional step third-order low-storage Runge–Kutta scheme and the implicit Crank–Nicolson scheme for all viscous terms at a maximum Courant–Friedrichs–Lewy (CFL) number 1.2. The CFL number is $\max (|u_1| + |u_2| + |u_3|)\,\Delta t/\Delta x$ over all cells, where the grid spacing $\Delta x$ is identical in all the three dimensions. Here, $u_1,u_2,u_3$ are the $x$-, $y$-, $z$-components of $\boldsymbol {u}$, and $\Delta t$ is the time step. The simulation is parallelised via slab decomposition along the $z$-direction.

Turbulence with $Re_\lambda = 72$ and $175$ is generated from an initially quiescent fluid following the tuning method proposed by Chouippe & Uhlmann (2015). The simulations are allowed to run until the pseudo-dissipation $\bar {\varepsilon }$ and $Re_\lambda$ are statistically stationary over at least $30T_L$. These statistically stationary flow fields are then used as the starting fields for the point-particle simulations. The flow statistics are listed in table 1. For validation, the longitudinal and transverse energy spectra are plotted in figure 1. Excellent agreement with the literature is shown.

Table 1. Statistics of the homogeneous isotropic turbulence: the grid size ($\mathcal {N}$), pseudo-dissipation ($\bar {\varepsilon }$), Kolmogorov length scale ($\eta$), maximum wavenumber ($k_{max}$), Kolmogorov velocity ($u_\eta$) scale, r.m.s. velocity fluctuations ($u'$), large-scale isotropy ($u_x'/u_y'$), and large eddy turnover time ($T_L$) relative to the Kolmogorov time scale ($\tau _\eta$). Here, $N_b,N_p$ are the numbers of bubbles and particles, respectively.

Figure 1. The (a) longitudinal and (b) transverse energy spectra in single-phase statistically stationary homogeneous isotropic turbulence. Dashed lines show the data from Jiménez et al. (1993); triangles represent the $-5/3$ power law. The agreement with the literature is excellent, so the dashed lines can be obscured, especially in (b).

3.2. Suspended phases

Bubbles and particles in the system are modelled using the point-particle approximation, where forces act on point masses. The bubble and particle dynamics are governed by (Maxey & Riley 1983; Tchen 1947)

(3.2)\begin{equation} \frac{4}{3}\,{\rm \pi} r_i^3\rho_i\,\frac{{\rm d}\boldsymbol{v}_i}{{\rm d}t} = 6{\rm \pi}\mu r_i\,f_i(\boldsymbol{u}-\boldsymbol{v}_i) + \frac{4}{3}\,{\rm \pi} r_i^3\rho_f\,\frac{\mathrm{D}\boldsymbol{u}}{\mathrm{D}t} + \frac{2}{3}\,{\rm \pi} r_i^3\rho_f\left(\frac{\mathrm{D}\boldsymbol{u}}{\mathrm{D}t} - \frac{{\rm d}\boldsymbol{v}_i}{{\rm d}t}\right), \end{equation}

where $\mu =\nu \rho _f$ is the absolute viscosity, $r_i$ is determined from $St_i$, and $\boldsymbol {u}$ is evaluated at particle position in this equation. The three terms on the right-hand side of (3.2) are the drag force, the pressure gradient force and the added mass force, respectively. Note that history forces, lift and reverse coupling are not considered to render it easier to disentangle the relative behaviour of bubble and particles. Nevertheless, as the lift force can be expected to play a role for bubbles, its effect will be discussed in § 4.3. Furthermore, a one-way coupled system is a consistent choice to study the geometric collision rate where hydrodynamic interactions are neglected. It entails the assumption of the dilute limit in which turbulence modifications by the suspended species are negligible (Brandt & Coletti 2022). The correction factor $f_i = 1 + 0.169\,Re_i^{2/3}$ accounts for finite bubble or particle Reynolds number $Re_i = 2r_i\,|\boldsymbol {w}_i|/\nu$, and implies the assumption of rigid spheres that obey the no-slip boundary condition for both species (Nguyen & Schulze 2004). This is realistic since liquids in flotation cells typically contain significant amounts of surfactants, so that the bubble surfaces would likely be contaminated (Nguyen & Schulze 2004; Huang, Legendre & Guiraud 2012). Although other commonly used expressions for $f_i$ are available, the difference is minimal, as shown in Appendix C.

Equation (3.2) is solved for each bubble and particle using a finite-difference scheme. To determine the flow velocities and the velocity gradients at bubble and particle positions, these quantities are interpolated from the Eulerian grid of the fluid solver to the particle positions using tri-cubic Hermite spline interpolation with a stencil of four points per direction. This choice is made as Hermite splines are comparable in accuracy (van Hinsberg, Clercx & Toschi 2017) and computationally cheaper to implement than B-splines (Ostilla-Monico et al. 2015). Time marching of $\boldsymbol {v}_i$ is performed with the explicit forward Euler method, and that of the positions of the suspended phases is done using the second-order Adams–Bashforth scheme. For stability, the time step is restricted such that neither the fluid nor the particle CFL number ($\max (v_{i1},v_{i2},v_{i3})\,\Delta t/\Delta x$) exceeds the value 1.2. This limit is enforced in both the suspended- and fluid-phase solvers, which run with the same simulation time step. We have compared particle statistics to data from the literature in order to verify our code, and those results are included in Appendix B.

Collisions between particles and bubbles are treated as ‘ghost collisions’. Under this approach, a ‘collision’ occurs once the centres of the members of an approaching pair reach the collision distance and the colliding pair pass each other without interaction. The collision radii are determined through the virtual radii $r_i$ as computed from $St_i$. This scheme is often employed by simulations of particles in turbulence (Bec et al. 2005, 2007; Goto & Vassilicos 2008; Voßkuhle et al. 2014; Ireland et al. 2016) and has been shown to be consistent with the formulation of $\varGamma _{pp}^{(ST)}$ (Wang, Wexler & Zhou 1998). To suppress the effect of different length scales when comparing collision pairs, we take $r_c = r_b + r_p$ for every type of collision. This is unlike studies examining solely particle–particle collisions, which define $r_c = 2r_p$ (e.g. Voßkuhle et al. 2014; Ireland et al. 2016). Numerically, these collisions are detected with the ‘proactive’ detection scheme in Sundaram & Collins (1996).

The details of the simulation of the suspended phases are as follows. Particles ($\rho _p/ \rho _f = 5$) and bubbles ($\rho _b/ \rho _f = 1/1000$) with $St_i = 0.1\unicode{x2013}3$ were seeded randomly and homogeneously ($10\,000\unicode{x2013}140\,000$ each) into the turbulent flow after it had reached a statistically stationary state as described in § 3.1. The density ratios correspond to sulphide minerals colliding with air bubbles in water, while the simulated $St$ range falls within $0.1 \lesssim St \lesssim 100$, which corresponds to particle sizes yielding the highest mineral recovery rate in conventional flotation cells ($10\,\mathrm {\mu }{\rm m}\lesssim r_p\lesssim 100\,\mathrm {\mu }{\rm m}$ when $\nu = 10^{-6}\,{\rm m}^{2}\,{\rm s}^{-1}$ and $\varepsilon = 5\,{\rm m}^{2}\,{\rm s}^{-3}$, following Ngo-Cong et al. 2018). For each case, bubbles and particles have identical $St_i$ to focus on the effect of different densities and to keep the parameter space manageable. The upper limit $St_i \leq 3$ is mandated by the fact that the increasingly large (virtual) bubble radius for larger $St$ violates the point-particle approximation. For $St_b = 3$ and at $Re_\lambda = 175$, we have $r_b/\eta \approx 5$, which is already marginal (Homann & Bec 2010). We monitor the p.d.f.s of bubble and particle positions for statistical stationarity. Once this has been reached, collision statistics are collected over at least $7.7T_L$. Two kinds of statistics are evaluated: in a fixed reference frame considering all particles/bubbles (Eulerian statistics), and along individual trajectories of the pairs that collide (Lagrangian statistics). The former are computed at least every $\sim 0.06T_L$, while the latter are calculated every time step (with only 25 % of all bubbles for the bubble–bubble case as the number of colliding bubbles is high). All the statistics presented are time- and ensemble-averaged unless indicated otherwise.

4. Results

4.1. Eulerian statistics

4.1.1. Collision kernel

Figure 2(a) shows simulation results for the dimensionless bubble–particle (bp), bubble–bubble (bb) and particle–particle (pp) collision kernels $\varGamma$. In addition to determining $\varGamma$ based on counting the number of collisions per time step (solid triangles), the collision kernels for $Re_\lambda = 175$ are determined indirectly via the RDF and the effective radial collision velocity according to (1.2) (shown as dots). Both results match closely, verifying our analysis procedure. With the normalisation by $\tau _\eta /r_c^3$ suggested by the Saffman–Turner framework, the results are insensitive to the change in $Re_\lambda$ for the pp and bp cases, but not for bb collisions. From the data, it is further evident that the relative behaviour of bubbles and particles is distinct from that of identical particles. For collisions between identical species, $\varGamma \tau _\eta /r_{c}^3$ is maximum when $St \sim 1$, while $\varGamma _{bp}\tau _\eta /r_{c}^3$ exhibits a minimum for this value of $St$. This trend is not captured by any of the models for the bp case discussed in § 2.4, for which the predictions are included as lines in figure 2(a). Generally, the model predictions are also significantly higher than the actual collision rates obtained from the simulations, the exception being the Saffman & Turner (1956) model, which best captures the magnitude yet fails to predict the proper $St$-trend, as is illustrated more clearly in figure 2(b).

Figure 2. (a) The dimensionless bubble–particle (bp), bubble–bubble (bb) and particle–particle (pp) collision kernels at collision distance. The triangle symbols denote $\varGamma$ determined directly from the collision rate, while $\varGamma _{4{\rm \pi} }$ (shown only for the $Re_\lambda = 175$ cases) is calculated according to (1.2). Unless specified otherwise, the colour conventions for all figures follow this figure. (b) Zoomed-in version of (a) plotted on linear scale.

The fact that the Abrahamson (1975) and large-$St$ Kruis & Kusters (1997) predictions do not match the data is to be expected as the present $St$ range does not match the assumptions made in these frameworks. Naturally, this also transfers to all approaches based on $\varGamma _{bp}^{(A)}$, and the somewhat better agreement with our data for models employing (2.7) and (2.8) instead of (2.6) is rather an artefact of the inconsistencies discussed earlier. This is emphasised by figure 3(a), where the large difference between the modelled slip velocities and the mean-square bubble/particle velocities $v_i'^2$ is obvious. Notably from the same figure, $v_p'^2$ is well predicted by the models of Abrahamson (1975), Kruis & Kusters (1997) (small $St$) and Zaichik et al. (2010). On the other hand, model estimates for $v_b'^2$ are generally too high. It is insightful to note the stark overprediction by (2.6) as the underlying framework by Abrahamson (1975) is also employed in the models by Ngo-Cong et al. (2018). The resulting overprediction of $v_b'^2$ might therefore explain in part why the estimate of $\varGamma _{bp}^{(NC)}$ is too large even though their concept (which is based on Yuu (1984)) is in principle more suitable for the moderate values of $St$ here. Unlike Abrahamson (1975), Kruis & Kusters (1997) (small $St$) include dissipation range scaling in their modelling, which is seen to improve the prediction at $St\lessapprox 1$ but still overestimates $v_b'^2$ for larger $St$. Finally, there is another factor, which affects all models. Bubbles especially do not sample the flow randomly, such that the mean-square fluid velocity at bubble locations $\langle u'^2\rangle _b$ is almost 10 % lower than $u'^2$, as shown in figure 3(b). This implies that bubbles sample flow regions with weaker fluid velocity fluctuations such that using $u'^2$ as model input may overestimate $v'^2_b$. Figure 3(b) also shows that this effect is less pronounced for heavy particles. The underlying cause for these observations is preferential sampling of flow regions, which we will discuss in § 4.1.2.

Figure 3. (a) The mean-square velocity of bubbles and particles at various $St$. Models are shown for $Re_\lambda = 175$ only. (b) The mean-square fluid velocity conditioned at bubble/particle positions for $Re_\lambda = 175$.

4.1.2. Bubble/particle spatial distribution

In order to elucidate in particular the $St$ trends for $\varGamma _{bp}$, we first investigate the distribution of bubbles and particles in the flow. While all the models assume homogeneous distributions, this does not hold at intermediate $St$, as figure 4 confirms, where instantaneous snapshots of the bubble and particle fields at different $St$ and $Re_\lambda =175$ are shown. For $St = 1$ and 3, bubbles and particles are seen to cluster but do so in different regions of the flow. This behaviour has been observed previously in the literature (Calzavarini et al. 2008a; Fayed & Ragab 2013; Wan et al. 2020) and is additionally shown by the different mean-square fluid velocity at bubble/particle positions in figure 3(b). To investigate this preferential concentration, we plot the norm of the rotation $\langle R^2 \rangle _r$ and strain $\langle \mathcal {S}^2 \rangle _r$ rates of the flow at the particle/bubble positions in figure 5, with $\langle \cdot \rangle _r$ denoting an ensemble average over particles in pairs with separations smaller than $r$. For tracers ($St=0$), we obtain $\tau _\eta ^2\langle R^2 \rangle _r=\tau _\eta ^2\langle \mathcal {S}^2 \rangle _r =0.5$ with $r\to \infty$, consistent with the analytical result in statistically stationary homogeneous isotropic turbulence. For the other cases, $\langle R^2\rangle _r$ and $\langle \mathcal {S}^2\rangle _r$ are conditioned on pairs with $r\leq 2r_c$. Figure 5(a) shows that bubbles (particles) cluster in regions of high (low) rotation rate, which is consistent with the centrifuge picture. Due to the clustering, conditioning has little effect for monodisperse collisions, and the results are therefore very close to single particle statistics. This is different for bp pairs, where $R^2$ is consistently lower (higher) for bubbles (particles) close to a particle (bubble). These observations imply that bp collisions occur for a subset of bubbles/particles that is located outside of their respective ‘preferred’ location within the flow. We note that correlations between particle/bubble locations and the strain rate (figure 5b) are much weaker and do not display the same qualitative trends observed for $R^2$. This is consistent with observations by e.g. Ireland et al. (2016) and Wang et al. (2020), who also found only a weak correlation between the positions of heavy particles and regions of the flow with low strain rate.

Figure 4. Instantaneous snapshots of bubbles and particles in a slice with width $\times$ height $\times$ depth = $L_{box} \times L_{box} \times 20\eta$ in the non-transient state for $Re_\lambda = 175$ at (a) $St = 0.1$, (b) $St = 1$, and (c) $St = 3$. The size of the $St = 0.1$ bubbles and particles is tripled for visibility.

Figure 5. Average value of the norm of (a) the rotation rate $R^2$, and (b) the strain rate $\mathcal {S}^2$, of the flow at bubble/particle/tracer positions. The bubble and particle data are conditioned on pairs with separation $r \leq 2r_{c}$. Here, bp,b (bp,p) refer to bubbles (particles) in bubble–particle pairs.

A consequence of the preferential concentration in different flow regions is that bubbles and particles become segregated. This effect is quantified by the RDF $g(r)$. Essentially, $g(r)$ relates the actual number of pairs with separation $r$ to that expected for uniformly distributed particles. Therefore, $g(r) = 1$ when particles are uniformly distributed, $g(r)>1$ implies clustering, and $g(r)<1$ implies segregation (Saw 2008). Figure 6(a) shows results for the RDF at collision distance from our simulations. Bubble–particle pairs indeed segregate, and the segregation is strongest when $St = 1$, as reflected by the local minimum of $g(r_{c})$. This is compatible with the corresponding behaviour of bubbles and particles, which always cluster for the tested parameters, with maximum clustering when $St\sim 1$. Increasing $Re_\lambda$ increases segregation slightly, but the effect is weak.

Figure 6. (a) RDF at collision distance. (b) The bubble–particle RDF as a function of $r$ and the segregation length scale $r_{seg}$ at $Re_\lambda = 175$. The number above each line is the corresponding $St$, and the dotted segments denote the RDF below the collision distance. Inset shows $r_{seg}$ at various $St$.

Figure 6(b) displays the bubble–particle RDF as a function of $r$ at $Re_\lambda =175$, which helps us to understand the effect of different length scales. The shape of the $g(r)$ curves is distinctly different from the power-law behaviour reported for the monodisperse case (Ireland et al. 2016). For all $St$, there is an essentially flat region for $g(r)$ at small scales, followed by a transition region with the steepest gradient at intermediate scale ${\sim }10\eta$ before approaching 1 for large separations. This behaviour suggests that the relevant length scale for the segregation, i.e. a typical distance between bubble and particle clusters, is at intermediate scales. To quantify this more precisely, we define the separation corresponding to the point of inflexion in $g_{bp}(r)$ as the segregation length scale $r_{seg}$. The values of $r_{seg}$ are marked by crosses in figure 6(b), and plotted against $St$ in the inset. The figure shows that $r_{seg}/\eta$ increases approximately linearly with $St$ and increases slightly with $Re_\lambda$. These findings are consistent with those by Calzavarini et al. (2008a), who studied segregation based on a concept inspired by Kolmogorov's distance measure (Kolmogorov 1963). In particular, these authors also report segregation scales ${\sim }10\eta$ with an increasing trend for higher $Re_\lambda$.

The trends observed for $g(r_c)$ in figure 6(a) remarkably resemble those discussed for the $St$ dependence of $\varGamma$ earlier, in the context of figure 2. It is therefore suggestive to think that the effects of segregation may explain in particular the discrepancy between the data and $\varGamma ^{(ST)}$. We can check this by plotting $\varGamma ^{(STc)}$ as defined in (2.2), which is shown as hollow symbols in figure 7(a). From this plot, it can be seen that correcting $\varGamma ^{(ST)}$ with the RDF leads to an almost perfect match with the pp collision kernel. This differs from, but is only seemingly at odds with, the results in Voßkuhle et al. (2014) due to the larger collision distance considered here as a consequence of keeping $r_c= r_b + r_p$ constant for all collision types. This therefore implies that while pp collisions at moderate $St$ are governed by the sling mechanism at small separations, they remain shear-dominated at larger ones. In contrast, $\varGamma ^{(STc)}$ overcorrects the bp collision kernel and undercorrects that of bb collisions. In figure 7(b), we replot the same data in the form of the ratio $\varGamma ^{(STc)}/\varGamma$, which can be interpreted as the relative contribution of the shear mechanism (compensated for segregation/clustering) to the overall collision rate. For pp collisions, this ratio is very close to 1 throughout. The situation is different for the bb case, where $\varGamma ^{(STc)}/\varGamma$ approaches 1 only for the lowest $St$ considered, and the value drops significantly for the higher $St$. Consistently the lowest values for $\varGamma ^{(STc)}/\varGamma$ are observed for bp collisions where the value quickly drops to around 0.5, which implies that a significant part of the relative velocities cannot be explained by the shear mechanism in this case. It is important to stress here that keeping $r_c$ constant for all collision types for a given $St$ removes this as a factor, such that the observed trends across species can only be rooted in differences in the relative approach velocities, which will be studied next.

Figure 7. (a) The dimensionless collision kernel (solid symbols) and the Saffman–Turner prediction after accounting for the RDF $\varGamma ^{(STc)}$ (hollow symbols). (b) The ratio between $\varGamma ^{(STc)}$ and the collision kernel.

4.1.3. Relative velocities

We show results for the non-dimensional effective radial approach velocity $S_-/u_\eta$ at collision distance in figure 8(a). The non-dimensional approach velocity increases with $St$ and also slightly with $Re_\lambda$ for collisions involving bubbles, whereas it has only a very weak $Re_\lambda$ dependence for pp collisions, as also reported in Ireland et al. (2016). Consistent with figure 7, the approach velocities are highest for bp collisions and lowest in the pp case.

Figure 8. (a) The effective radial collision velocity. (b) The effective bubble–particle radial approach velocity against pair separation $r$ at $Re_\lambda = 175$. (c) The effective radial approach velocity of each type of collisions across $St = 0.1$ to $St = 3$ at $Re_\lambda = 175$. The dashed lines show the prediction by Zaichik et al. (2010).

To understand these trends, it is instructive to consider $S_-$ as a function of $r$, which is presented in figure 8(b) for bp collisions for different $St$, and in the plots of figure 8(c) for different collision types at constant $St$. Also shown in these figures is the effective approach velocity for tracer particles, for which $S_-$ is proportional to $r$ in the dissipation range, as expected. In particular, for low $St$, where $S_-(r)$ remains close to the tracer curve, this in part explains the $St$ dependence of $S_-(r_c)$ as $r_c$ increases with increasing $St$. However, the curves also ‘peel off’ the linear scaling at increasingly larger distances and to a larger extent as $St$ increases. For the pp case, this is a well-known manifestation of the ‘sling effect’, where due to path history effects, particles arrive at the same location with different velocities such that their relative velocities exceed those of the fluid. Note that while the bb case has a larger $S_-$ than the tracer case, the linear scaling remains mostly unchanged there, such that the difference presumably is rather due to preferential concentration effects and the fact that the slip velocity at the same $St$ is higher for this case (see (2.9)). The deviations from the tracer behaviour are strongest for the bp case throughout, and already for $St\geq 2$ $S_-$ are essentially independent of $r$ for $r\lessapprox 10\eta$. The decorrelation between bubble and particle velocities implicit in these observations is the reason why $S_-$ is largest for all $St$ in figure 8(a). The model of Zaichik et al. (2010) fails to reproduce this feature. Instead, $S_-$ from the model is highest for bb collisions, such that the trend falls in line with that of $v_i'^2$ (see figure 3). Plotting the model prediction against the separation distance in figure 8(c) confirms that this trend persists over all $r$. Approach velocities are also generally overpredicted by the model, which is in part also related to the differences observed for $v_b'^2$ in figure 3 and to the significant peel off the linear scaling especially for the bb case. Consistent with the results in figure 7, $S_-^{pp}$ at both $Re_\lambda$ is matched very well by the Saffman–Turner model $S_-^{(ST)} = 0.5r_c\sqrt {2\varepsilon /15\nu {\rm \pi}}$, where the prefactor 0.5 accounts for the fact that only separating pairs are considered here.

In modelling approaches (e.g. Abrahamson 1975; Yuu 1984; Kruis & Kusters 1997; Zaichik et al. 2010; Ngo-Cong et al. 2018), it is common to deduce $S_-$ based on a Gaussian distribution of the relative radial velocity $\Delta v_r$. In this case, the ratio $S_-/\sqrt {S_{2\parallel }}$ can be determined to be $1/\sqrt {2{\rm \pi} } \approx 0.4$, where the variance is defined as

(4.1)\begin{equation} S_{2\parallel}(r) = \int_{-\infty}^{\infty} (\Delta v_r)^2\mathrm{p.d.f.}(\Delta v_r|r)\,\mathrm{d}(\Delta v_r). \end{equation}

It was already reported in Ireland et al. (2016) that $S_-/\sqrt {S_{2\parallel }}$ can drop significantly below the Gaussian value for pp collisions at $St \sim 1$ and for small $r_c$. This is confirmed by our results in figure 9, where also the bb case is seen to follow similar trends. For bp collisions, the ratio remains much closer to the Gaussian value, almost indifferent to $St$ and close to the tracer result. Relative velocities are therefore better approximated by the Gaussian assumption for the bp case, while the overprediction of $S_-$ resulting from doing so is larger for the monodisperse collisions.

Figure 9. The ratio of the effective radial approach velocity to the standard deviation of the radial component of the relative velocity $\Delta v_r$ at $Re_\lambda = 175$.

For a more quantitative investigation of variation of $S_-$ across collision types, we plot the excess bubble–particle effective radial collision velocity $\Delta S_- = S_-^{bp} - (S_-^{bb}+S_-^{pp})/2$ in figure 10(a). At small $St$, such an additional approach velocity arises for the bp case from the fact that, due to the change in sign of $\beta _i$ in (2.9), bubbles and particles react differently when subjected to the same fluid acceleration. Indeed, figure 10(b) shows that $\Delta v_r^{bp}$ conditioned on pairs with $r\in [r_c-\eta /2,r_c+\eta /2]$ at small $St$ closely follows the resulting expression for the relative velocity $\Delta v_r/u_\eta = St\,(\beta _b - \beta _p)a_{fr}/a_\eta$, where $a_{fr}$ is the average radial fluid acceleration at bubble and particle positions. We therefore conceptualise the resulting excess bp collision velocity as the local ‘turnstile mechanism’, since bubbles and particles respectively move in opposite directions, as illustrated in figure 10(c). An estimate for this effect can be obtained using $a_f \sim a_\eta$ in (2.9) and assuming a random orientation between the particle separation vector and the fluid acceleration, resulting in an average angle $\phi _{rdm} = 1\,\mathrm {rad} \approx 57^\circ$. This is shown by the dashed line in figure 10(a). While not strictly a prediction for $\Delta S_-$, this estimate is largely consistent with the data for $St<1$, suggesting that the enhanced approach velocity in this range is indeed due to the local turnstile effect. This is further consistent with the fact that for $St\geq 1$, the correlation of $\Delta v_r^{bp}$ with the local fluid acceleration is seen to vanish almost completely in figure 10(b). A positive $\Delta S_-$ is also predicted by the model of Zaichik et al. (2010), but underestimates the value in the low-$St$ regime. Beyond $St = 1$, $\Delta S_-$ flattens off at high values. Along with the observations made in figure 8(c), this suggests that the ‘path history’ effect active at these $St$ values is enhanced for the bp case. We can rationalise this in terms of a non-local ‘turnstile mechanism’ (see figure 10c), where differences in the particle velocities arise as bubbles and particles interact differently with larger and more energetic eddies in the flow. As a consequence, their local velocities differ more than for monodisperse cases at equal $St$. Both the local and non-local inertial effects discussed here enhance the relative approach velocity and therefore explain why $\varGamma ^{(STc)}_{bp} < \varGamma _{bp}$ as observed before.

Figure 10. (a) The excess bubble–particle effective radial collision velocity $\Delta S_- = S_-^{bp} - (S_-^{bb}+S_-^{pp})/2$ from simulations over a range of $St$. The factor $1/2$ for the dashed line accounts for the fact that $S_-$ is obtained by averaging only over the negative portion of $\mathrm {p.d.f.}(\Delta v_r|r_c)$ as defined in (1.3). (b) The radial component of the relative velocity conditioned on pairs with $r\in [r_c-\eta /2,r_c+\eta /2]$ binned by the local radial fluid acceleration at $Re_\lambda =175$. (c) Schematic illustrations of the local and non-local turnstile mechanisms.

4.2. Lagrangian statistics

To better examine the collision process, we adopt the Lagrangian point of view, where colliding pairs are tracked individually for separations $r \leq 100\eta$ up to collision. The corresponding statistics averaging over all pairs with the same $r$ are denoted by the suffix $\mathcal {L}$.

In figure 11, we show how $R^2|_\mathcal {L}$ and $\mathcal {S}^2|_\mathcal {L}$ vary as functions of $r$ for the colliding pairs. Consistent with the results in figure 5, variations are significantly stronger for $R^2|_\mathcal {L}$ (figure 11a) compared to $\mathcal {S}^2|_\mathcal {L}$ (figure 11b). In contrast to the monodisperse case, bubbles and particles need to leave their preferred flow regions in order to collide. While at $r/\eta \sim 100$ the unconditioned preferential concentration (indicated by the symbols in figure 11) is mostly recovered, a pronounced drift generally sets in for $r/\eta \lessapprox 50$. This might suggest that interaction with eddies at this intermediate range plays a role in bringing bubbles and particles together from their segregated locations. At $St = 0.1$, bubbles and particles appear to occupy similar flow regions for $r/\eta \lessapprox 10$ before colliding, i.e. in particular, the $R^2|_\mathcal {L}$ curve flattens at these scales, whereas the same is not the case at $St = 1$ and $St = 3$, where the collision distance is larger. These observations are consistent with the considerations regarding the local and non-local turnstile mechanisms made above. It is noteworthy that $R^2|_\mathcal {L}$ at $St \leq 1$ varies mostly for bubbles as they approach particles, while at $St = 3$, $R^2|_\mathcal {L}$ is almost constant for the bubbles but varies more for the particles. The latter case also stands out in the $\mathcal {S}^2|_\mathcal {L}$ plot as the slight decrease in strain values towards collision is not observed there.

Figure 11. Norm of (a) the rotation rates and (b) the strain rates of the flow at bubble and particle positions against pair separation conditioned on colliding pairs at $Re_\lambda = 175$ for $r \geq r_{c}$. Only the $St = 1$ case is shown for bubble–bubble and particle–particle pairs. Also plotted at $r\rightarrow +\infty$ is the unconditioned $R^2$ and $\mathcal {S}^2$ at the respective $St$.

In figure 12, we present velocity statistics conditioned on the colliding pairs. In particular, we consider the approach velocity $S_-|_\mathcal {L}$ (figure 12a) as well as the r.m.s. relative speed $(\Delta v)'|_\mathcal {L}$ (figure 12b). The approach velocities differ slightly in magnitude from the unconditioned results in figure 8, but the trends remain consistent. In particular, $S^{bp}_-|_\mathcal {L}$ is again largest for separations close to the respective collision distances for all $St$ shown. This is remarkable in view of the fact that $(\Delta v)'|_\mathcal {L}$ is actually largest for the bb case. This implies that $\Delta \boldsymbol {v}$ and $\boldsymbol {r}$ tend to be more anti-aligned for the bp case to yield the higher relative approach velocity.

Figure 12. (a) The effective approach velocity, and (b) the r.m.s. relative speed, conditioned on colliding pairs at $Re_\lambda = 175$.

To measure this, the angle between the separation vector and the relative velocity vector $\Delta \phi$ (see the inset of figure 13a) is used. For a head-on approach, $\Delta \phi =0^\circ$, while $\Delta \phi = 90^\circ$ means that the other particle is circling around in the rest frame of the collision partner. Note that as colliding particles must be approaching each other, $0^\circ \leq \Delta \phi <90^\circ$. As demonstrated in figure 13(a), the bubble–particle pairs are indeed the most anti-aligned, i.e. $\Delta \phi |_\mathcal {L}$ is lowest, relative to the other types of collisions near $r_{c}$. Figure 13(b) shows the ratio of the cosines of $\Delta \phi |_\mathcal {L}(r_c)$ as a measure for how much differences in alignment lead to a relative enhancement of $S_-(r_c)$. The ratio is up to 10 for bp relative to bb collisions, and up to ${\sim }1.5$ compared to particle–particle collisions, and values for both cases are slightly higher at the higher $Re_\lambda$. When interpreting the smaller difference between pp and bp collisions here, it should be kept in mind that based on the discussion around figure 7, the collision mechanism differs between these two cases: whereas pp collisions appear predominantly shear-driven, where a good alignment may be expected, inertial effects, for which this is not necessarily the case, play a much more significant role for bp collisions.

Figure 13. (a) The angle between the separation and relative velocity vectors $\Delta \phi$ for $r \geq r_{c}$ at $Re_\lambda = 175$. Bubble–particle/bubble–bubble/particle–particle pairs are represented by green/blue/red lines as in figure 12. The inset shows the definition of $\Delta \phi$ in the rest frame of the particle on the right. (b) The extra angular contribution to the radial component of the bubble–particle collision velocity relative to the bubble–bubble and particle–particle cases.

4.3. Effect of lift, finite particle density and nonlinear drag

Our simulations are performed without the lift force using a finite particle density $\rho _p/\rho _f=5$ and a nonlinear drag law with $f_i = 1+0.169\,Re_i^{2/3}$. To test how sensitive our results are to these parameters, we ran addition simulations for $Re_\lambda =175$ changing one parameter at a time: adding the lift force $\boldsymbol {F_L} = -2/3\rho _f{\rm \pi} r_i^3(\boldsymbol {v}_i-\boldsymbol {u})\times (\boldsymbol {\nabla }\times \boldsymbol {u})$ to the right-hand side of (3.2), setting $\rho _p/\rho _f=\infty$, or using Stokes drag ($f_i = 1$) for both phases. For all of these simulations, $r_c$ is chosen to be identical to the case when $\rho _p/\rho _f=5$. We compare the results for $\varGamma$ in figure 14, and those for $g(r_{c})$ and $S_-(r_{c})$ in figure 15. The plots show that the collision statistics are not sensitive to the lift force. Also more generally, the effect of the lift force is limited, with the most noticeable change occurring for the preferential sampling where we observed a slight decrease in $\tau _\eta \langle R^2\rangle _r$ at bubble clusters, consistent with Mazzitelli, Lohse & Toschi (2003). It is also evident that the particle density affects these quantities only marginally. Similarly, the results remain largely unchanged for the $f_i = 1$ case, with the exception of a slight increase in $\varGamma _{bp}$ due to a higher $S_-$ at the larger $St$. We furthermore note that other expressions for $f_i$ have been proposed in the literature (Schiller & Naumann 1933; Clift, Grace & Weber 2005; Wan et al. 2020). These are compared in Appendix C, which shows that they do not differ significantly when $Re_i \lesssim 100$. Hence we conclude that our results would not be sensitive to whichever one of these parametrisations is chosen. Based on the findings in this section, it therefore appears appropriate to neglect lift, and to use the simplifications of linear drag and $\rho _p/\rho _f \to \infty$ for simulations and modelling approaches in the present parameter range.

Figure 14. The dimensionless collision kernel for different particle densities simulated, with the lift force included and for the case when the drag correction factor is $f_i = 1$.

Figure 15. (a) The RDF at collision distance, and (b) the effective radial collision velocity for different particle densities, with the lift force included and for the case when the drag correction factor is $f_i = 1$. The symbols follow figure 14.

5. Discussion and conclusion

We studied bubble–particle collisions in turbulence using a point-particle approach. Besides a critical appraisal of current models for the problem, our results highlight that the difference in the density ratios between bubbles and particles critically affects the collision statistics. An overview of the physical picture that emerges is presented in figure 16. For $St \to 0$, the shear mechanism is applicable, meaning that the collision velocity is well-described by the local fluid velocity gradient. Once $St$ reaches finite values, bubbles and particles concentrate preferentially in different regions of the flow, which leads to their segregation. This effect reduces the bubble–particle collision kernel $\varGamma _{bp}$ and is maximal for $St \sim 1$. The correction remains ${O}(1)$, which is comparable to the enhancement of the collision kernel due to clustering for particle–particle cases, but much less than that for the bubble–bubble cases, respectively. For particle–particle collisions (evaluated at $r_c = r_b + r_p$), the relative approach velocities are found to be entirely consistent with the shear mechanism for the full range of $St$ explored here, such that the collision kernel in these cases is very well approximated by extending the shear-induced collision kernel to account for non-uniform particle concentration $\varGamma ^{(STc)}$. The same does not hold for bubble–particle and bubble–bubble collisions, where inertial effects additionally play a role and enhance the effective relative approach velocities beyond the shear scaling. This enhancement is strongest for the bubble–particle case, and we conceptualise this effect in terms of a ‘turnstile’ mechanism, which is related to the fact that fluid accelerations lead to opposing drift velocities relative to the fluid for bubbles and particles. This happens locally, i.e. on the same fluid element, if $St$ is sufficiently small for bubbles/particles to follow pathlines ($St\leq 0.5$ based on the results in figure 10), but also appears to play a role at higher $St$, where the path history becomes more relevant. For $St \gtrapprox 2$, the ratio $\varGamma ^{(STc)}/\varGamma$ and the values of the effective radial approach velocity at contact $S_-(r_c)$ are comparable for bubble–particle and bubble–bubble collisions, which indicates that the ‘non-local’ turnstile is less prevalent there, in line with the fact that local particle/bubble velocities become increasingly random at larger $St$.

Figure 16. A sketch of the bubble–particle collision mechanisms with their approximate $St$ dependence.

Our method does not allow us to assess the high $St$ (kinetic-gas-like) regime since bubble sizes violate the point-particle assumption. Applying the Abrahamson (1975) model pertinent to this regime to the present data vastly overpredicts the collision kernel, and the apparent improvement by adaptations of this framework is merely an artefact of using inconsistent velocity expressions. Our data for bubble–particle collisions are better approximated by the models of Saffman & Turner (1956) and Zaichik et al. (2010) when taking an additional correction for the inhomogeneous distribution into account. The predictions based on the theory of Zaichik et al. (2010) generally overestimate $S_-$ and fail to reproduce the finding that $S_-^{bp}$ is largest in the present parameter range when comparing across species. These observations are not affected when using a linear drag relation, i.e. $f_i = 1$ in (3.2), or in the limit $\rho _p/\rho _f \to \infty$ as shown in § 4.3. We therefore could not reproduce the good agreement with the Zaichik et al. (2010) theory for $S_-^{bp}$ reported in Fayed & Ragab (2013). Potentially, this is due to different choices for $St_b$, which we set to equal to $St_p$, whereas this parameter is kept constant in Fayed & Ragab (2013).

In our simulations, we employ the non-dimensional control parameters $St$ and $Re_\lambda$. In practice, the fluid and particle properties are controlled separately so using $r_b$, $r_p$ and $\varepsilon$ is more common. However, by definition, changing $\varepsilon$ affects both $St$ and $Re_\lambda$. To assess the overall effect, consider the region between $St = 0.1$ and $St = 0.5$, where the non-dimensional collision kernel $\varGamma _{bp}\tau _\eta /r_c^3$ decreases most sharply and drops by a factor of approximately 1.3. As the particle parameters are constant, $\tau _i = \textrm {const.}$, this is associated with a 5-fold decrease in $\tau _\eta$. According to the shear scaling, which appears to be applicable to our data, this increase in turbulence intensity implies also a 5-fold increase in $\varGamma _{bp}$. Hence there remains a net benefit from increasing the turbulence intensity even in this range, although the effect is reduced somewhat due to segregation.

Our results show that considering finite particle density, nonlinear drag and the lift force have limited effect on the bubble–particle collision rate. However, other relevant factors are yet to be explored. Most importantly, this concerns the buoyancy force and finite size effects in particular for the bubble motion and non-identical $St$ for bubbles and particles, as well as modelling the collision behaviour of bubbles and particles. More complex simulations as well as experiments are certainly needed to make precise predictions of the bubble–particle collision rate in realistic settings. However, as far as modelling approaches and a physical understanding are concerned, we believe that the most promising approach to disentangle the effect of additional factors (such as buoyancy) is to treat them as modifications to a simpler base case like the one studied here.