{'fraction_non_alphanumeric': 0.02900132275132275, 'fraction_numerical': 0.006779100529100529, 'mean_word_length': 4.7558051008755236, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Numerical modeling of bubble-particle interaction in a volume-of-fluid framework

A numerical method is presented to simulate gas-liquid-solid flows with bubble-particle interaction, including particle collision, sliding, and attachment. Gas-liquid flows are simulated in an Eulerian framework using a volume-of-fluid method. Particle motions are predicted in a Lagrangian framework. Algorithms that are used to detect collision and determine the sliding or attachment of the particle are developed. An effective bubble is introduced to model these bubble-particle interaction. The proposed numerical method is validated through experimental cases that entail the rising of a single bubble with particles. Collision and attachment probabilities obtained from the simulation are compared to model and experimental results based on bubble diameters, particle diameters, and contact angles. The particle trajectories near the bubble are presented to show differences with and without the proposed bubble-particle interaction model. The sliding and attachment of the colliding particle are observed using this model.


Introduction
Bubble-particle interaction is a common physical phenomenon that can be observed in various industrial applications, such as in the recovery of coal and minerals from ores, de-inking of paper for the recycling of waste paper, and inclusion removal in steel industries [1]. Specifically, an inclusion removal refers to a technique that is widely used to manufacture ultra-clean steel. An inclusion refers to a micro-sized particle created from the oxidation of molten steel, which degrades steel quality. To enhance steel quality, the removal technique utilizes bubbles to collect inclusions by attaching them on the bubble surface [2]. Thus, understanding the physics of inclusion attachment is directly related to collision efficiency. However, inclusion removal is a complex phenomenon that is challenging to be described by a numerical simulation. Therefore, a simple method that reflects the physical phenomena of the removal process must be used. Mechanisms of bubble-particle interaction are related to three different phenomena, namely, oscillating collision, liquid film creation and sliding. When the flow-driven particle is in contact with the bubble, it oscillates because of the force balance between particle inertia and surface tension, i.e., oscillating collision.
Then, a liquid film is created while the particle remains on the bubble surface, i.e., liquid film creation. During existence of the liquid film, the particle slides on the bubble surface, i.e., sliding. Furthermore, the particle is attached if the liquid film is ruptured. These attachment processes are characterized by three time intervals, namely, collision, sliding, and induction time intervals.
The collision time interval is characterized by particle oscillation on the bubble surface. This oscillating motion is described as the harmonic oscillation [3]. In the same sense, the sliding time is characterized by the interval in which the particle slides over the bubble surface. This phenomenon was modeled by assuming inertialess particle and potential flow fields [4,5] and adapting the experimental results [6,7]. Lastly, the induction time indicates the lasting duration of lasting the liquid film. According to Nguyen et al. [5], these time intervals are interpreted as a relative position between a particle and a bubble. They adapted a concept of a critical angle to describe the induction and sliding time, which will be explained in detail in Section 3. This study uses characteristic time and geometrical relations.
Recently, bubble-particle interaction is solved numerically. Bubbles and nearby flow fields are solved in the volume-of-fluid (VOF) framework, whereas particles are solved in a Lagrangian framework [8,9,10,11,12]. In most works, bubble-particle interactions, such as collision, sliding, and attachment, are not considered. However, such phenomena are critical factors in determining particle collecting efficiency or probability. Some work consider actual physical phenomena. They modeled hydrophobic forces that act on particles near a stationary bubble and concluded that considering the acting forces on particles and bubbles is critical in predicting bubble-particle interaction [13,14]. However, known flow fields are used, and the dynamics of a phase interface are ignored.
In the present study, we target to model complex bubble-particle interaction by considering particle attachment mechanisms. Flow fields and bubble interfaces are solved in an Eulerian framework using the VOF method. Particle motions are solved in a Lagrangian framework. To model bubble-particle interaction, we adapt the theoretical model developed by Nguyen et al. [5], which reflects the actual physical phenomena. In this model, the position of particles relative to the bubble surface is a key factor. However, in the VOF framework, a bubble interface is represented in an Eulerian framework. Thus, the exact collision point cannot be determined. To address this problem, an effective bubble is adopted to model particle motions on the bubble surface accurately, which will be demonstrated in Section 4.
The present paper begins with a description of the numerical methods for simulating gas-liquid-solid flows in Section 2. The physical background and numerical algorithms of bubble-particle interaction are presented in Sections 3 and 4, respectably. Section 5 provides the validation results of the present numerical method and compared with the experimental results by Hewitt et al. [15]. The conclusions are drawn in Section 6.

Governing equations
Gas-liquid flows are governed by the Navier-Stokes equations with surface tension. The one-fluid formulation is adopted using total density and dynamic viscosity instead of the solving governing equations of each phase. The incompressible variable-density Navier-Stokes equations with surface tension are given as follows: where , , , , , , and are the Cartesian coordinates, the corresponding velocity components, pressure, total density, total dynamic viscosity, gravitational acceleration, and surface tension, respectively. The total density and dynamic viscosity are defined as follows: where , , , , and are the densities, dynamic viscosities of the gas and liquid phases, and the volume fraction of the gas phase.
The numerical schemes for solving Eq. (2) are similar to those of Choi and Moin [16]. The Crank-Nicolson time advancement scheme is employed for Eq. (2).
A Newton-iterative method is applied to solve the nonlinear terms. A second-order central difference scheme is used to discretize spatial derivative terms. Then, a fractional step method and an approximate factorization technique [17] are applied to the discretized equations. Eq. (1) leads to the Poisson equation for pressure, which is solved by an algebraic multigrid method [18].
The solid particle motions are predicted in the Lagrangian framework. The Lagrangian equations that govern the particle motions are expressed as follows: where x p , u p , u, , and p are the particle position and velocity, fluid velocity, drag coefficient, and particle Stokes number. Schiller and Naumann [19] suggested that the drag coefficient to be modeled as follows: where p is the particle Reynolds number. The particle Stokes number p is defined as follows: where p and p are the particle density and the particle diameter, respectably.
Eqs (5) and (6) are solved using a third-order Runge-Kutta method. The fluid velocity u at the particle position x p is obtained by employing a fourth-order Lagrange interpolating polynomial.

Volume-of-fluid method
The phases are distinguished by volume fraction and it is advected by the VOF method. The advection equations for the volume fraction is written as follows: A piecewise linear interface calculation (PLIC) approach [20] is used to solve Eq. (9). In this approach, the phase interface of each cell is represented as a line or a plane equation. The VOF/PLIC method proceeds in two steps. First, the interface is reconstructed by the volume fraction . The reconstruction step requires the estimation of the normal vector and the plane constant of the interface.
The interface of each cell is represented as follows: where , , and are the local normal vector of the interface, the position vector, and the plane constant of the interface, respectably. In the present study, the normal vector is calculated by a mixed Youngs-centered (MYC) method, as suggested by Aulisa et al. [21]. Then, the plane constant can be obtained analytically using the geometrical relation between the normal vector and the volume fraction [22]. After the reconstruction step, the volume fraction is advected by Eq. (9) using a geometrical flux computation. The conservative advection scheme suggested by Weymouth and Yue [23] is adopted to guarantee that the error of mass conservation is within the range of machine accuracy. The total density and dynamic viscosity are calculated with Eqs. (3) and (4) using the updated volume fraction, respectably.
Furthermore, surface tension needs to be estimated at the phase interface.
This force is calculated by the continuum surface force method of Brackbill [24] : where and are the surface tension coefficient and the curvature of the phase interface , respectably. The curvature is estimated by the height function method combined with the paraboloid-fitted method developed by Popinet [25]. The sharpness of the surface tension is conserved by implementing a balanced force algorithm proposed by Francois et al. [26].

Mechanisms of sliding and attachment of the colliding particle
The numerical modeling of bubble-particle interaction requires understanding on how a particle moves on a bubble surface. In this study, we focus on the case where the particle size is much smaller than the bubble size that is the particle diameter p is less than 0.1 times of the minimum bubble diameter b , which is 0.75 mm. Moreover, a particle shape is assumed to be a rigid sphere. Under this condition, the particle momentum is negligible compared with the bubble momentum. In other words, the particle motions do not affect those of the bubble.
In addition, these size differences cause capillary force to become dominant at the bubble surface. Consequently, the particle cannot transform the bubble surface.
The particle motions on the bubble surface can be illustrated using three characteristic time intervals. These time intervals determines whether the particle is attached or slides away after collision (Fig. 1). First, the collision time c refers to the duration of oscillation before sliding. Second, the sliding time s corresponds to the period when the particle slides over the bubble surface. Lastly, the induction time i covers the period from thin liquid film formation to rupture. The thin liquid film plays an important role in determining sliding or attachment of the colliding particle. When the thin liquid film exists, the particle keeps sliding because the thin film interrupts the attraction force between the particle and the bubble surface.
In other words, the particle is attached on the condition that c or s > i . When c > i , the thin liquid film is ruptured during particle oscillation, and then the particle is attached. If s > i , then the thin liquid film is ruptured during sliding.
When s i , the thin liquid film exist during sliding. Finally, the particle slides away. The relationships among the time criteria are summarized as follows: i. If c > i , then the particle is attached during oscillation.
ii. If s > i , then attachment occurs during sliding.
iii. If s i , then the particle slides away.
Particle motions on the bubble surface can also be described according to relative position of the bubble and the particle. Relative position is represented by the angle measured from a stagnation point to the point of particle contact. angle criteria. Three characteristic angles are used to determine particle motions.
First, the bubble-particle angle o is the angle at the particle position from the stagnation point. Secondly, the critical angle cr is the angle beyond which no attachment occurs. Lastly, the collision angle m is the angle beyond which no collision occurs. From the definition, no collision occurs when o > m . If cr o m , then the particle collides but the particle slides away. The reason is that the collision point is beyond the critical angle cr . When o < cr , the particle is attached. The relationships among the angle criteria are summarized as follows: i. If o > m , then no collision occurs.
ii. If cr o m , then the particle slides away.
iii. If o < cr , then the particle is attached during oscillation or sliding.
In the present model, the time and angle criteria are used to determine the particle

Collision and attachment probability
According to particle motions on the bubble surface, the number of the colliding and attached particles are determined. From these numbers, a collision and attachment probability are defined. These probabilities can be used to estimate the The collision probability c is defined as the fraction of an actual number of colliding particles and an ideal number of colliding particles [5] as follows: where ci is the ideal number of colliding particles, and cr is the actual number of colliding particles. Fig. 3 schematically shows that how particles collide to the bubble surface. Ideal collision implies that all particles within the bubble sweep path collide with the bubble surface. In other words, if the location of the particles is within the bubble radius b , then the particles ideally collide with the bubble surface. However, the actual collision occurs when the particle location is within the collision radius OC . From this definition, the collision probability c can be calculated geometrically as follows: where 2 b is the area where the ideal collision occurs, and 2 OC is the area where the actual collision occurs.
The collision probability c also can be modeled using the material properties and flow conditions [27] as follows: where the dimensionless parameters , , , and are given by where b , p , and b = b b / are the bubble velocity, particle velocity, and bubble Reynolds number, respectably. To calculate Eq (14), the bubble velocity and particle velocity need to be obtained by using material properties. The bubble velocity b has an empirical relation in the case of a rising single bubble [28] as The particle velocity is given by Oeters [29] as follows: where p is the particle density. Eqs. (16) and (17) are substituted into Eq. (14) to obtain the collision probability. Then, the collision angle m can be calculated using Eqs. (13) and (14).
The attachment probability a is defined by the fraction of as actual number of attached and colliding particles [5] as follows: where a is the number of attached particles. The actual attachment occurs when the particle location is within the attachment radius OA (Fig. 3). a of 1 denotes that all colliding particles are attached on the bubble surface. In the same manner of defining the collision probability, we can define the attachment probability a as follows: where 2 OA is the area where the particles are attached. In this study, the attachment probability a is fitted from the experimental measurements conducted by Hewitt et al. [15]. An exponential form is adopted as follows: where model constants 1 and 2 are written by where cont is the contact angle. Then, the critical angle cr can be calculated using the collision angle m , Eqs. (19) and (20).

Numerical modeling collision, sliding, and attachment
The present study mainly aims to model the physics of the bubble-particle interaction explained in Section 3 numerically. Collision occurs when the particle is in contact with the bubble surface. Then, the particle is attached or slides away according to the initial collision position of the particle. Particle motions including collision, sliding, and attachment, are predicted using the time and angle criteria obtained from empirical models. The present numerical model consists of two procedures. The first one detects the collision of the particle near the bubble surface. The second procedure determines the sliding or attachment of the colliding particle. Fig. 4 schematically illustrates the total procedures of modeling of the bubble-particle interaction. Detailed algorithms and methodologies are presented in following subsections.

Relative position and angle between bubble and particle
As mentioned in Section 3.1, the angle criteria are used to determine the particle motions on the bubble surface. To apply the angle criteria, the bubble-particle angle o is required. This angle can be calculated from the position relationships of the bubble surface and the particle in a single cell (Fig 5). Each vector described by the bold letter is obtained using the VOF method and simple vector calculation.
The following descriptions show how to obtain the bubble-particle angle o .
n is the normal vector of the bubble surface. This vector is calculated from the geometric VOF method, as described in Section 2.2. p is the vector from the coordinate origin to the particle center. a is the vector from the coordinate origin to the center of the bubble surface, which is equal to the center of the  Figure 5: Schematic of position relationships between the bubble and particle. a is the vector from the coordinate origin to the bubble surface center. n is the normal vector of the bubble surface. b is the vector from the coordinate origin to the bubble center. p is the vector from the coordinate origin to the particle center. p − b is the vector from the bubble center to the particle center. plane equation that represents the bubble surface. This location can be found by calculating center of the plane equation. Lastly, b is the vector from the coordinate origin to the bubble center. Fig. 5 shows that the direction of the normal vector of the bubble surface n corresponds to a − b. The distance between the bubble surface and the bubble center is the bubble radius b . This length is calculated from the curvature of the bubble obtained from the VOF method, that is b = 2/ . By using this relation, the position vector of the bubble center b can be calculated where = b /|n|. Finally, we can obtain the particle position vector relative to (a) * Figure 6: Schematic of (a) the effective bubble and (b) the actual bubble. b-p is the distance between the bubble surface to the particle center. The actual bubble has the diameter b and the effective bubble has the diameter b + 2 × * . the bubble center p − b. The bubble-particle angle o denotes the angle between the bubble velocity b and the particle position vector relative to the bubble center On the basis of this definition, the bubble-particle angle o can be obtained using the inner product formula as follows:

Effective bubble
Determining particle motions on the bubble surface begins from detecting bubble-particle collision because the sliding and attachment are only examined for the colliding particle. The collision of bubble-particle implies that the particle is in a physical contact with the bubble surface. In other words, a collision condition refers to the bubble-particle distance is equal to the particle radius. However, a new definition of collision is required because the exact collision position cannot be determined for the following reason. The bubble surface represented on the computational domain is not exactly defined in the VOF framework. This interface is reconstructed from the volume fraction in an Eulerian framework. However, the particle size is smaller than the grid size. Thus, the exact distance between the particle center and the bubble surface is unable to be calculated.
An effective bubble is introduced due to this difficulty in finding the exact position of collision. The collision is detected on the effective bubble surface.
This phenomenon implies that the sliding and attachment are determined on the effective bubble surface. Fig. 6 shows the regions of two bubbles. The dark-shaded region represents the actual bubble with the diameter b . The light-shaded region denotes the effective bubble with the diameter b + 2 × * . The bubble-particle distance b-p is defined by the distance between the particle center and the bubble surface. The bubble-particle distance b-p of the -th particle can be calculated where p denotes the position vector of the -th particle. The size of the effective bubble is controlled by * .

Effective collision and critical angle
As mentioned in the previous section, the particle motions on the bubble surface need to be determined on the effective bubble. The collision point shifts from the actual bubble to the effective bubble. Accordingly, the collision and critical angles need to be tuned. Fig. 7 shows the geometrical definition of the effective collision angle * m and the effective critical angle * cr . In this research, an effective factor is multiplied to the original collision angle m and critical angle cr . The effective collision angle * m and critical angle * cr are defined as follows: * Appendix A presents the detailed form of effective factors m and cr and tuning process.

Detecting collision
The algorithm for detecting collision is presented in the VOF framework.
Algorithm 1 summarizes the algorithm for detecting collision. The interface cell corresponds to the cell with the volume fraction of 0 < < 1. First, the existence of the interface is examined for each particle after updating the particle position.
A 3 × 3 × 3 stencil is the candidate cell for examination. If the interface cell exists, then the normal vector n of the interface and the position of the bubble center are obtained using the geometrical VOF method. Next, the bubble-particle distance b-p is calculated using Eq. (24). The condition of b-p * implies that the particle is within the effective bubble. For this case, the bubble-particle angle o and the effective collision angle * m are calculated. The final collision is determined by comparing these two angles. If o * m , then the colliding particle is tagged.

Determination of sliding or attachment
After detecting collision, the last step is to determine sliding or attachment of the colliding particle using the time and angle criteria. First, the collision time c and the induction time i are computed. The collision time c modeled by Evans [3] is used. This model is nearly identical to the mean value of other collision time models [2]. The collision time c is modeled as follows: This time refers to the duration of oscillation after particle collision. The induction time i modeled by Nguyen et al. [5] is written as follows: where dimensionless parameters A and B are given by where and are given by Eqs. (15a) and (15b), respectably. The effective critical angle * cr is obtained by Eq. (25b). The procedures for determining sliding or attachment is summarized in Algorithm 2. As discussed in Section 3.1, i < c indicates that the colliding particle is attached to the bubble surface because the thin liquid film has been ruptured during oscillation. Otherwise, the bubble-particle angle o and the effective critical angle * cr are compared. When o < * cr , the particle is turned out attachment and eliminated in the computational domain.

Rising of a single bubble with particles
Numerical simulations of the rising of a single bubble with particles are conducted to verify the present gas-liquid-solid flow solver using the bubble-particle interaction model. Fig. 8 shows the computational configuration. The simulation results are compared with the experimental results obtained by Hewitt et al. [15].   Before validating the bubble-particle interaction model, the terminal velocity of a rising single bubble without a particle is compared with that in the empirical relation given by Eq. (16). Terminal velocity is estimated by the rising velocity of the bubble center. The simulation results are in good agreement with the empirical relation (Fig. 9). The following subsections show results of the rising of the single bubble with particles.

Collision probability and diameter
The collision probabilities obtained from the present simulation using Eq. (12) are compared with the model written in Eq. (14). The ideal number of colliding particles ci is 201 2 . The particles are squared because they are linearly distributed.
Similarly, the actual number of colliding particles cr is obtained by computing the square of the total sum of the tagged particles. Fig. 10 shows the collision probabilities as a function of particle and bubble diameters. The simulation results are in good agreement with the model results. This finding shows that the present model can predict the collision probability using the effective bubble.  The collision diameters OC are further compared with the results calculated using the model equation. Eqs. (12) and (13) show that the square of the collision diameter is proportional to the total number of colliding particles. This length can be obtained using Eq. (13), that is OC = b √ c . Fig. 11 shows that the collision diameters as a function of particle and bubble diameters. The collision diameter OC linearly increases with the particle diameter p . This trend implies that the total number of colliding particles is proportional to the square of the particle diameter p .

Attachment probability and diameter
The attachment probability and diameter are compared with the model and experimental results. The attachment probability in the present simulation is estimated by using Eq. (18). The total number of attached particle a is obtained by computing the square of the total sum of eliminated particles. Fig. 12 shows that The attachment diameter is calculated by OA Eqs. (13) and (19). The square of the attachment diameter can be interpreted as the total number of attached particles. Fig. 13 shows the attachment diameters as a function of particle diameters, bubble diameters, and contact angles. Overall, the simulation results are well matched with the model and experimental results.
When the bubble diameter is 0.75 mm, the attachment diameter OA increases with the particle diameter p . However, for the condition of the bubble diameter of 1.2 mm, the attachment diameter does not change much as the particle diameter p increases.
This tendency is due to the product of the collision and attachment probability as follows: where is called a flotation probability [5] or total attachment probability [2]. This probability corresponds to the fraction of the number of attached particles and the ideal number of colliding particles. This relation also shows that the attachment diameter OA is proportional to sin cr when the bubble diameter b is fixed. The model results of Fig. 13 is obtained from b √ . Moreover, the attachment diameter OA increases with the contact angle cont because the increment of the contact angle cont causes the attachment probability a to rise (Fig. 12).

Particle trajectory
The present model directly distinguishes the sliding or attachment of the colliding particle, thereby enabling us to observe the particle motions near the bubble surface. Figs. 14 and 15 show the particle trajectories near the bubble. The bubbles are represented using the iso-surfaces of the volume fraction of 0.5.
The particle trajectories are compared with and without the bubble-particle      Fig. 15(b) is larger than that in Fig. 15(c). The reason is that the attachment diameter of the former case is 90 μm, which is larger than the latter case of 60 μm (Figs. 13(c) and 13(d)).

Rising of multiple bubbles with particles
The present bubble-particle interaction model can be applied to multiple bubbles. The computational domain is equivalent to the single bubble case, as presented in Section 5.1. Periodic boundary conditions are applied to the side walls instead of using the slip wall conditions. Fig. 16 shows that the 3 × 3 stationary    that the present bubble-particle interaction model can predict particle motions on the bubble surface regardless of the number of bubbles.

Conclusions
Numerical algorithms for simulating gas-liquid-solid flows, including bubbleparticle interaction, have been proposed . Phase interfaces are described and advected by the geometrical conserving VOF method. The kinematics of the particles is solved in the Lagrangian framework. The bubble-particle interaction is modeled in the VOF framework to predict the particle motions on the bubble surface, including collision, sliding, and attachment. The new algorithms consist of two procedures, namely, detecting collision and determining sliding or attachment for the colliding particle. These algorithms require the angle between two vectors, the bubble velocity, and the vector of the particle position from the bubble center.
The calculation of this angle is proposed using the geometric VOF method. In addition, the use of the effective bubble is suggested because the exact point of collision cannot be obtained in the VOF framework.

Significance and novelty
• A novel numerical framework for simulating gas-liquid-solid flows is proposed with the bubble-particle interaction model.
• New algorithms for detecting collision and determining sliding or attachment in a VOF framework are suggested using the geometrical relationship between a bubble and a particle and by adopting the concept of an effective bubble.
• The present numerical model for bubble-particle interaction is validated and visualized using simulation cases of rising of a single bubble and rising of multiple bubbles with particles in terms of bubble diameters, particle diameters, and contact angles.

Highlights
• Simulating gas-liquid-solid flows by combining an Eulerian approach and a Lagrangian approach • Proposing algorithms for detecting collision and determining sliding or attachment • Validating the numerical method based on the collision and attachment probability • Visualizing bubbles through particles with and without using the bubbleparticle interaction model