A numerical extension of the “smooth profile method” is presently suggested to simulate the attachment of a colloidal particle to the surface of an immersed bubble. In this approach, the two fluid-particle boundaries and the fluidic boundary are replaced with diffuse interfaces. The method is tested under various capillary numbers. Upon attachment to a stable bubble, it is found that the method is capable of reproducing the three microprocesses associated with the particle attachment. The change in the trajectory as the particle approaches the fluidic interface, the collision process, and the sliding down the bubble surface are all captured. Potential application of the present method shows great promise in the field of froth flotation, where the capture of hydrophobic particles by rising bubbles is of primary importance.

## I. INTRODUCTION

The transport of colloidal particles at the fluid-fluid interface of a binary fluid, henceforth simply referred to as the fluidic interface, is of significant importance in various processes, one of which is the flotation process.^{1} The valuable hydrophobic particles attach to the fluidic interface of rising bubbles while the undesired hydrophilic particles settle down to eventually be discharged. The simulation of the entire particle attachment process is complex. Current attachment models are still at an early stage of development. The fine attaching particles have so far been modelled as point particles,^{2,3} thereby neglecting the deformation of the fluidic interface and the coupling between the solid and the binary fluid. The use of a boundary-resolved simulation would overcome many of the shortcomings associated with point-particle models. Only a handful of studies dealing with the direct numerical simulation of particle motion in binary systems are found in the literature. To the best of the authors’ knowledge, none of these studies have yet been applied to simulate the dynamics of the particle attachment. This can probably be attributed to the great deal of difficulty in developing three-phase flow models for which the capillary effects are of chief importance. Aland *et al.*^{4} suggested a two-dimensional sharp interface model to simulate the fall of a colloidal particle across a horizontal fluidic interface. By changing the magnitude of the surface tension, the particle could either cross the fluidic interface or be trapped at the interface. Choi and Anderson^{5} developed a two-dimensional dynamic mesh procedure, termed the “extended finite element method,” to simulate the particle motion across a stratified binary fluid initially at rest. The method was later improved through the use of a mesh refinement across the diffuse fluidic interface.^{6} Fujita *et al.*^{7} combined an immersed boundary method and a level-set method to simulate the rearrangement of multiple particles during the drying of a colloidal film on a substrate. The lattice Boltzmann method was recently combined with a binary-fluid model to simulate the penetration of a spherical particle in a fluidic interface.^{8} Further reading on the recent development of lattice Boltzmann method for the transport of colloidal particles at a fluidic interface can be found in the Refs. 9–11. Here a novel combination of the smooth profile method^{12} with a diffuse binary-fluid model^{13} is suggested to directly simulate the attachment of a solid particle to an immersed bubble.

## II. THREE-PHASE MODEL

### A. Extended smooth profile method

The merit of the smooth profile method lies in the fact that the use of a dynamic mesh around the particle is not needed. In this method, the sharp particle boundary is replaced with a diffuse interface of length *ξ _{p}*, across which a particle concentration field

*ϕ*smoothly transitions from unity (particle region) to zero (binary fluid region). The following truncated hyperbolic function is used to define the solid field concentration:

_{p} The distance **r** = **x** − **X**_{p} is calculated from the particle centre of mass **X**_{p}. The terms *r _{p}* and

*r*, respectively, correspond to the particle radius and to the cutoff radius. The total fluid velocity field

_{c}**u**is decomposed into two components as

The first component (1 − *ϕ _{p}*)

**u**

_{f}corresponds to the velocity field of the binary fluid. The second component is the particle velocity field defined as

where **V**_{p} and **Ω**_{p} are, respectively, the translational and the rotational velocity of the particle. Note that the capital letters are used to distinguish the Lagrangian variables from their Eulerian counterparts. The evolution of the total field velocity, which satisfies the incompressibility condition ∇ ⋅ **u** = 0 over the entire computational domain, is governed by the following modified momentum equation:

where *σ*_{v} = *η*[∇**u** + (∇**u**)^{T}] is the viscous stress tensor, *ρ* the density of the binary fluid, *t* the time, *p* the pressure, and *η* the fluid viscosity. The term *ρϕ _{p}*

**f**

_{p}is introduced to enforce the particle rigidity. The exact formulation of the particle rigidity term, shown in Equation (18), is derived by assuming a momentum conservation between the particle and the fluid. Our newly defined capillary contribution

**f**

_{c}is defined in Subsec. II B. Only the salient feature of the smooth profile method was presented. Further reading may be found in the Ref. 14.

### B. Binary fluid model

The two concentration fields *ϕ*_{ℓ} and *ϕ _{g}*, which similarly smoothly transition from unity to zero, are used to distinguish between the two fluids. The two phases of the binary fluid are here introduced as “liquid” and “gas” for readability purposes. The concentration of the two fluids is driven by the minimisation of the free energy,

^{15,16}

where $V$ is the region of space occupied by the isothermal system, *k*_{B} the Boltzmann constant, *T*_{0} the temperature, *v*_{0} a constant unit volume, and *f* the free energy density. Without loss of generality it is postulated that the Ginzburg-Landau-type equation^{13} of the free energy density is made up of four terms,

The first term on the right hand side of Eq. (6) is the bulk energy density. Its formulation, based on the mean field theory,^{17} is expressed as

where *χ* is a parameter describing the affinity between the liquid and the gas. The second, the third, and the fourth term on the right-hand side of Equation (7) are the “capillary terms” and account for the three interfacial energies. The tunable interfacial length scales *ξ*_{ℓg}, *ξ*_{ℓp}, and *ξ _{gp}* are, respectively, introduced to control the interfacial energies. Because of the phase summation constraint

*ϕ*

_{ℓ}+

*ϕ*+

_{g}*ϕ*= 1, the order parameter

_{p}*ψ*(

**x**,

*t*) =

*ϕ*

_{ℓ}−

*ϕ*is now introduced. The bulk energy density

_{g}*f*

_{0}(

*ψ*,

*ϕ*) is shown in Figure 1. The logarithmic function is presently expanded to a fourth-order polynomial expression at the critical point and will emulate the fluid separation with a performance similar to its non-expanded counterpart.

_{p}^{18}The order parameter is updated in time according to the following modified “Cahn-Hilliard” equation:

^{19}

where *M* is the mobility, **I** the unit tensor, **n**_{p} = − ∇*ϕ _{p}*/|∇

*ϕ*| the local unit vector normal to the particle surface, and $\mu =\delta F/\delta \psi $ the chemical potential. The tensor operator (

_{p}**I**−

**n**

_{p}⊗

**n**

_{p}) prevents a mass flux normal to the diffuse boundary of the particle. It mimics the boundary condition

**n**

_{p}⋅ ∇

*μ*|

_{Sp}= 0 normally enforced at the sharp interface

*S*of a particle.

_{p}^{8}Far from the fluid-particle diffuse interface, where the magnitude of the gradient falls below an arbitrarily small value |∇

*ϕ*| <

_{p}*ϵ*, the outer product vector

**n**

_{p}⊗

**n**

_{p}equates the zero tensor. Finally, the formulation of the capillary term in Eq. (4) is inspired by a ternary fluid flow simulation.

^{20}It is here suggested as

where the second chemical potential is defined as $ \mu p =\delta F/\delta \varphi p $. The exact formulation of the two chemical potentials can be found in the supplementary material.

### C. Particle transport and numerical procedure

Let **u**^{n} and *ψ ^{n}* be the fields at the initial time

*t*. The field

^{n}*ψ*

^{n+1}is first determined by solving the Cahn-Hilliard Equation (8). In the next stage, the position of the particle centre of mass

**X**

_{p}and the particle orientation

**Θ**

_{p}are updated as

A fractional step approach is employed to solve momentum Equation (4). An intermediate velocity is first calculated as

The intermediate pressure term *p*^{∗} is calculated by solving a Poisson equation ∇ ⋅ **u**^{∗} = 0. Retaining only the hydrodynamic drag force **F**_{h}, the capillary force **F**_{c}, and an external force **F**_{e} representing the particle weight, the translational and rotational particle velocities are updated in time as^{21}

where *m _{p}* and

**I**

_{p}are the mass and the inertia tensor of the particle. The hydrodynamic force and torque exerted by the fluid on the particle are derived by assuming momentum conservation between the particle and the fluid. The time integral of the hydrodynamic force therefore equals the momentum change over the particle domain.

^{22}Since the capillary contribution is accounted for in the calculation of the intermediate velocity, the summation of the hydrodynamic force and of the capillary force is given by

In the second stage of the fractional step, the particle velocity field is enforced onto the total fluid velocity field as

The pressure *p _{p}*, originating from the rigidity constraint, is obtained from the incompressibility condition ∇ ⋅

**u**

^{n+1}= 0. The time integral of the force density field needed to enforce the particle rigidity is calculated as

The above equations are discretised in space using a second-order central differencing scheme and in time using a forward Euler method. The authors are aware of more advanced schemes.^{23} We believe however that the above discretisation schemes are justified by the motivation of this work, which primarily involves the development of a new boundary-resolved method to simulate the dynamics of the attachment process.

## III. MODEL VALIDATION

It is assumed that the two fluids have equal viscosity and that the three phases have equal density. This assumption greatly simplifies the numerical modelling and, at low capillary number, should not cause a significant loss of precision in the simulation of the particle trajectories. The equations were implemented in their non-dimensional form using the Reynolds number Re, the Peclet number Pe, and the capillary number Ca given by

The reference length is defined as *L*_{0} = *ξ*_{ℓg}, the reference time as $ t 0 = L 0 2 / D 0 $, and the reference velocity as *U*_{0} = *γ*/*η*. The two terms *D*_{0} = *Mk*_{B}*T*_{0}/*v*_{0} and *γ* = *L*_{0}*k*_{B}*T*_{0}/*v*_{0} correspond to the diffusion coefficient and to the interfacial energy between the two fluids. The parameters used in the following four simulation sets S1, S2, S3, and S4 are summarised in Table I.

. | χ
. | ξ/Δ
. _{p} | ξ_{ℓg}/Δ
. | r/Δ
. _{p} | r/_{c}ξ
. _{p} | $ C \u2113 2 $ . | $ C g 2 $ . | Pe . | Re . | Ca . | |F_{e}|/γL_{0}
. | N
. _{grid} | c × 10^{3}
. |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|

S1 | … | 3 | … | 5 | 0.5 | … | … | … | 0.01 | … | … | 64^{2}, 64^{3} | 19.1, 2.0 |

S2 | 2.5 | 3, 6, 3 | 1, 2, 1 | 10, 20, 10 | 0.5, 1, 0.5 | 1 | 1..2 | 1 | 0.01 | 1 | 0 | 64^{2}, 128^{2}, 64^{3} | 79.7, 79.7, 15.9 |

S3 | 2.5 | 2 | 1 | 5 | 1 | 1 | 1 | 1 | 0.01 | 0.1, 1, 10 | 1 | 128^{2} | 4.8 |

S4 | 2.5 | 2 | 1 | 5 | 1 | 1 | 2 | 1 | 0.01 | 0.1 | 0.5, 1, 1.5 | 256^{2} | 1.2 |

. | χ
. | ξ/Δ
. _{p} | ξ_{ℓg}/Δ
. | r/Δ
. _{p} | r/_{c}ξ
. _{p} | $ C \u2113 2 $ . | $ C g 2 $ . | Pe . | Re . | Ca . | |F_{e}|/γL_{0}
. | N
. _{grid} | c × 10^{3}
. |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|

S1 | … | 3 | … | 5 | 0.5 | … | … | … | 0.01 | … | … | 64^{2}, 64^{3} | 19.1, 2.0 |

S2 | 2.5 | 3, 6, 3 | 1, 2, 1 | 10, 20, 10 | 0.5, 1, 0.5 | 1 | 1..2 | 1 | 0.01 | 1 | 0 | 64^{2}, 128^{2}, 64^{3} | 79.7, 79.7, 15.9 |

S3 | 2.5 | 2 | 1 | 5 | 1 | 1 | 1 | 1 | 0.01 | 0.1, 1, 10 | 1 | 128^{2} | 4.8 |

S4 | 2.5 | 2 | 1 | 5 | 1 | 1 | 2 | 1 | 0.01 | 0.1 | 0.5, 1, 1.5 | 256^{2} | 1.2 |

### A. Hydrodynamic drag (S1)

To validate the in-house implementation of the smooth profile method, the Stokes flow of a single liquid phase through a periodic array of spheres is first simulated. A sphere, subject to a constant external force, is initially placed at the centre of the system. After a short period of time, the sphere reaches its terminal velocity $ V p \u221e $. The sedimentation velocity is calculated as

The magnitude of the external force is varied and the corresponding drag coefficient calculated. Figure 2 shows an excellent agreement with the theory for two- and three-dimensional tests. The theoretical formulation of the drag coefficient presently used for comparison purposes is given in Appendix. The effect of the interfacial thickness and of the particle resolution on the drag has been extensively discussed.^{22,24} With *r _{p}* = 5 and

*ξ*/Δ = 3, the hydrodynamic drag was found to deviate by about 10%, a value in close agreement with previously reported data.

_{p}^{12}

### B. Particle hydrophobicity (S2)

The particle hydrophobicity can be changed by raising or lowering the interfacial energies.^{25} A change in the particle hydrophobicity is here achieved by varying the Cahn numbers,

The first Cahn number *C*_{ℓ} is arbitrarily set to unity while the second Cahn number *C _{g}* is allowed to vary. When

*C*

_{ℓ}=

*C*= 1, the particle is neither hydrophobic nor hydrophilic. This is because the liquid-particle interfacial energy

_{g}*γ*

_{ℓp}and the gas-particle interfacial energy

*γ*are equal and so the contact angle, given by

_{gp}*θ*= (

*γ*

_{ℓp}−

*γ*)/

_{gp}*γ*, equates

*θ*= 90°. For

*C*≥ 1, the contact angle is first solved analytically using the hyperbolic function shown in Equation (1). The two cutoff values shown in Figure 3 are also used for the analytical determination of

_{g}*θ*. The two theoretical curves slightly differ by a few degrees. Two and three dimensional tests are then performed to estimate the geometrical contact angle. In line with the analytical solution, it is found that an increase in $ C g 2 $ causes a linear decrease in the contact angle. A finer grid resolution brings the linear fit closer to the analytical solution. A greater number of grid nodes could not be achieved because of the high computational cost. It can also be seen that the slope of the linear fit is highly sensitive to the cutoff radius

*r*so that an exact comparison to the analytical solution is difficult to achieve. With

_{c}*r*=

_{c}*ξ*/2, only the far end of the hyperbolic tail is cut off.

_{p}## IV. RESULTS AND DISCUSSIONS

### A. Capillary effect on bubble deformation (S3)

A set of simulations is now performed to qualitatively assess the effect of the capillary number on the particle attachment. All the following simulations are performed in a two dimensional space. Three dimensional tests would lead to a qualitatively similar conclusion. At the initial time *t* = 0, a circular bubble is placed at the centre of the computational domain. The position of the bubble centre of mass is maintained throughout the simulation by subtracting the contribution of the bubble velocity from the total momentum. This step is performed at each iteration right after the calculation of the intermediate velocity. The bubble velocity is calculated as

The condition *ψ*/|*ψ _{g}*| < − 0.8 corresponds to the bubble region excluding both the interfacial fluidic area and the particle area. The particle is released at an arbitrarily upstream position and is brought in motion by a fictitious gravitational force, arbitrarily set (Table I). The chosen low Reynolds number indicates here a Stokes flow around the moving particle. As seen in Figure 4, an increase in the capillary number causes a larger deformation of the fluidic interface. The resolved flow field around the particle is also shown. Setting the capillary number to Ca = 0.1 results in a realistic attachment to an immersed gas bubble. The bubble shape remains fairly stable and the three micro-processes,

^{1}which are the particle approach, the collision, and the downward sliding, are all captured. A capillary number lower than unity is in fact expected for the present simulation since the capillary force is the dominant force during the sliding process.

^{26,27}

. | Present sim. . | Reference 27 . | Reference 28 . |
---|---|---|---|

r/_{p}r _{b} | 0.125 | 0.05 ± 0.025 | 0.2 ± 0.1 |

θ | 60° | ∼65° | ∼65° |

Re | 10^{−2} | Re_{p} ∼ 10^{−2} | Re_{p} ∼ 10^{−2} |

Ca | 0.1 | ∼10^{−4} | ∼10^{−4} |

|F_{e}|/(γr) _{p} | ∼10^{−1} | ∼10^{−3} | ∼10^{−4} |

### B. Particle attachment (S4)

The number of grid points is here increased. In this way, the particle is released from a higher altitude and the trajectories are analysed in greater detail. Figure 5 shows the effect of the external force on the particle trajectory. The present two-dimensional particle trajectories are compared with three-dimensional data taken from the literature. The data from the first reference were obtained from a particle-point model^{3} and proved accurate in the simulation of finely attaching particles.^{27} In the second comparative experiment,^{28} the liquid was brought in slow laminar motion around a pinned gas bubble. Although the microspheres later detached after a temporary sliding, the experimental trajectories are suitable for comparison during the particle approach. Table II compares the present simulation parameters with those of the two references. The experimental external force includes the particle weight and the buoyancy force exerted by the water. The present particle trajectories, illustrated in Figure 5, show a qualitatively good agreement with the reference data. It is seen that the slowly increasing trajectory deviation as the particle approaches the fluidic interface is captured by the present model. The trajectory deviation, defined as |*X _{p}*(

*t*) −

*X*(0)|, is the distance from the current particle position to the unperturbed vertical straight path, taken by the particle in the far-bubble region. The model also predicts an increasing deviation with increasing external force. This tendency is in-line with experimental data.

_{p}^{29}Figure 6(a) shows the evolution of the particle velocity magnitude as a function of the gap

*h*. The velocity is normalised with particle terminal velocity. The gap

*h*is defined as the shortest distance from the particle surface to the fluidic interface. An exact formulation of the gap is however difficult because of the deformation and of the diffusivity of the fluidic interface. It is therefore assumed that the gap drops to zero as soon as the collision process starts. It is found that, during the approach, the particle velocity magnitude decreases as it comes closer the interface. This tendency is in-line with the reference data.

^{27}The same evolution of the particle velocity is also shown in Figure 6(b) as a function of the polar angle and brings a more natural visualisation of the collision and sliding processes. The polar angle ranges from

*α*≈ 0° (upstream pole of the bubble) to

*α*= 180° (downstream pole of the bubble). The bubble equator is indicated by

*α*= 90°. It is seen that, during the collision process, the particle velocity sharply peaks over a short period of time. The particle velocity then reaches a second maxima close to the bubble equator, where the tangential component of the external force is largest.

### C. Limitations

The presently suggested method has proved successful in the simulation of a particle attachment to an immersed bubble under various capillary numbers. The preliminary tests showed that the three micro-processes (approach, collision, and sliding) observed during the attachment to a stable bubble could be numerically reproduced. Some limitations are however associated with our method. In the close vicinity of the fluidic interface, the particle trajectory suddenly changes even though the sharp particle boundary and the bubble isoline *ψ* = 0 do not yet meet. Upon this close encounter, the diffuse fluid-particle interfacial region and the diffuse fluidic interfacial region start overlapping. The capillary force therefore acts as a short-range attraction force that rapidly pulls the particle towards the fluidic interface. This effect could be reduced by diminishing the interfacial thickness. An exact comparison with experimental data is also difficult. The experimental capillary number of a typical air-water mixture is smaller than that of the present study. The use of such a small value is however not practicable numerically.

## V. CONCLUSIONS

A boundary-resolved model has been presented to simulate the attachment of a particle to the fluidic interface of an immersed bubble under various capillary numbers. The smooth profile method, originally developed for the simulation of particle transport, has successfully been combined with a new binary fluid model. Our preliminary tests showed that, upon attachment to a stable bubble, the three micro-processes associated with the particle attachment could be numerically reproduced. The model could predict the change in trajectory as the particle approaches the fluidic interface, the collision process, and the sliding down the bubble surface. Extension of the present method shows great promise in the field of froth flotation, where the attachment of mineral hydrophobic particles to rising air bubbles is of utmost importance. Future improvements will involve a lower capillary number and fluids with different densities and viscosities.

## SUPPLEMENTARY MATERIAL

See supplementary material for the mathematical formulation of two chemical potentials *μ* and *μ _{p}* shown in Equation (9).

## Acknowledgments

This work was supported by a Marie Curie International Outgoing Fellowship with the European Union Seventh Framework Program for Research and Technological Development (2007-2013) under the Grant Agreement No. 623518.

### APPENDIX: HYDRODYNAMIC DRAG

The Stokes flow through a periodic cubic array of spheres has been studied theoretically by Hasimoto.^{30} The corrected drag formulation as a function of the solid fraction, defined as $c=4\pi r p 3 / ( 3 L 3 ) $, is given by

The above formula is valid for *c* < 0.25.^{31} The extension to the flow through a periodic array of infinitely long circular cylinders also exists. The drag force per length unit is given by

where $c=\pi r p 2 / L 2 $.