Bubble rising and interaction in ternary fluid flow: a phase field study

Abstract

Bubble–droplet interaction is essential in the gas-flotation technique employed in wastewater treatment. However, due to the limitations of experimental methods, the details of the fluid flow involved have not been fully understood. Therefore, a phase field model for a three-phase flow was developed to study the rise of a single bubble and bubble–droplet interactions. The fluid–fluid interfaces are tracked by the Cahn–Hilliard equation, which is coupled with the Navier–Stokes equations with an equivalent volumetric force substituted for interfacial tensions. The model was discretized using an explicit finite difference method on a half staggered grid, and the pressure velocity coupling was tackled using the projection method. The in-house code was written in Fortran and run with the help of OpenMP, a shared memory parallelism. The model was validated against experiments with gratifying agreement achieved. Bubble–droplet interaction was simulated in two distinct situations: the first features a gas bubble crossing the interface between two other phases, and the second features a gas bubble chasing from behind an oil droplet in a surrounding fluid of the third phase.

---

1. Introduction

Bubble rising, a type of multiphase flow phenomenon, occurs in many industrial applications, such as boiling,^1–3 oil processing,^4–6 and exploitation of combustible ices.^7–9 A better understanding of bubble rising dynamics is thus crucial in developing optimal processing parameters for these processes. Previous studies, both numerical and experimental, have been mostly on the behavior of a single bubble rising in a heavier medium, and on the rising velocity, aspect ratio, terminal shape, and the like.^10–12 In different flow conditions, denoted by a couple of dimensionless parameters, like the Eötvös number (defined later on), bubbles display distinct final shapes and rising velocities. Besides, some empirical correlations for the rising velocity have been proposed.^13–15 The interactions of two or more bubbles of the same fluid have been investigated by Tripathi et al.,¹⁶ who found that the destabilizing nature of the wake leads to oscillatory trajectories.

Among numerical techniques to investigate bubble rising dynamics of a binary fluid system is the phase field method.^17–19 Distinct from traditional sharp interface methods, it embraces a thin but finite interface region instead of a sharp one of vanishingly small thickness.²⁰ In phase field modeling, the tracking of fluid–fluid interfaces via the Cahn–Hilliard equation is implicit, and its discretization and coding are easier to handle in comparison with other methods such as the volume of fluid method and the front tracking method. Moreover, the Cahn–Hilliard model boasts of sound conservation of mass, which is however a matter of concern in the level set method. Last but not least, with the nowadays GPU technology, the phase field method could be readily extended to three dimensions without the adaptive mesh refinement technique.²¹

The phase field model could be extended to more than three phases as well. Kim²² proposed a phase field model for ternary fluids, with the Cahn–Hilliard equation to track the interface between each pair of components. The interfacial tension is expressed as a volumetric term in the Navier–Stokes equations. This model conforms to the requirements of consistency and degeneracy.^23,24 Consistency means that when a phase is absent at the beginning it will be so all the time. Degeneracy states that when a phase is absent a ternary model degenerates into a binary one. Kim²² did not consider contact line dynamics for ternary fluid flow. Nevertheless, the issue has been studied by a number of scholars.^25–28 Other phase field models for ternary fluid flow can be found in a review paper.²⁹

Phase field modeling of the interaction between two immiscible drops could be found in a few studies.^30–32 They however mainly focused on the engulfing mechanism in a confined shear flow, and treated all the fluids as having equal densities and viscosities. For vertical head-on collisions of two immiscible droplets, Zhang et al.³³ investigated the effect of the ratio of interfacial tensions on the film thickness, the maximal deformation, and the contact time. It is noticed that the two droplets that were vertically aligned crashed into each other in the process.

To the authors' best knowledge, only a couple of studies concerning the gas-flotation technique in wastewater treatment were conducted numerically, as by Ming et al.³⁴ via a smooth particle hydrodynamics model and by Kalantarpour et al.³⁵ via a phase field Lattice Boltzmann method. Kalantarpour et al.³⁵ also simulated bubble–droplet interaction in water, but their model cannot access all the range of parameters, which is however not a matter of concern in the model by Kim.²²

Having surveyed the literature and found that little effort is paid to understand the details of fluid flow in the gas-flotation technique, this paper is therefore dedicated to that end, with a view to gaining better understanding into the mechanism behind the gas-flotation technique in wastewater treatment. The paper is organized as follows. First, mathematical statement and numerical schemes are given, followed by the mesh convergence study to determine a reasonable spatial step. Second, effect of the phase field mobility was examined. Third, single bubble rising across the interface between two other immiscible fluids was simulated. Fourth, bubble–droplet interaction in a heavier medium was studied. The phase field mobility was found to drastically influence bubble rising velocity and was adjusted based on the experimental work by Liu et al.³⁶

2. Mathematical statement

2.1 Mass and momentum equations

The governing equations for triple phase flows of incompressible Newtonian fluids turn up as follows, with eqn (1) being the continuity equation and eqn (2) the Navier–Stokes equation. Du/Dt in eqn (2) is the material derivative. The third term on the right hand side of eqn (2) represents interfacial tension, which will be defined later on.

∇·u = 0 (1)

In the equations above, u is velocity, t is time and T = −pI + σ is the total stress tensor. p is the mechanical pressure and σ = μ[∇u + (∇u)^T] is the Newtonian stress tensor, with μ being viscosity. g stands for the local gravitational acceleration. It is to be noted that all of the field parameters, like ρ = ρ₁c₁ + ρ₂c₂ + ρ₃c₃, are to be defined as continuous functions of the order parameters. Besides, c₁ + c₂ + c₃ = 1, where the subscripts represent different phases. c₁ = 1 represents water, c₂ = 1 air and c₃ = 1 Oliver oil.

2.2 Phase field equation for an evolving fluid–fluid interface

In this paper, the motion of a fluid–fluid interface is described by the Cahn–Hilliard equation, which tracks the interface implicitly with the phase indicator c_i.where M is the diffusion coefficient and φ_i is the chemical potential, with being the bulk free energy density. ξ stands for a measure of interfacial thickness. β(c₁,c₂,c₃) in eqn (4) is a Lagrange-multiplier to ensure . Since c₁ + c₂ + c₃ = 1, only two of them needs updating. Additionally, the diffusion strength M and the characteristic interface thickness ξ are assumed constant and equal for each fluid–fluid interface.

2.3 Interfacial tension

The interfacial tension between any two phases is converted to a volumetric force in phase field modeling. In this paper, a volumetric counterpart of capillary force according to Kim²² is adopted. It takes on the following form.where α = 6√2 and γ_i is the phase specific interfacial tension coefficient, satisfying the relationship σ_ij = γ_i + γ_j. σ_ij refers to the interfacial tension between phase i and phase j. Therefore, for triple phase systems, γ_i is uniquely determined.

2.4 Boundary/initial conditions and numerical procedures

Boundary conditions are needed to close a set of partial differential equations. The schematic of the problem is given in Fig. 1, with all the borders being treated as walls. Initially, the air bubble sits about 12Δx above the bottom in a quiescent fluid. The dashed line in Fig. 1(b) symbolizes an interface. The effect of computational domain has been ruled out. Besides, the diameters of the bubble and the droplet, unless otherwise stated, is set to 2.7 mm throughout the paper.

Fig. 1 Schematic of the problem. (a) For Section 3.1 to 3.3, (b) for Section 3.4 and (c) for Section 3.5.

Thanks to the symmetric motion of bubble rising, only half computational domain enters calculation. Zero Neumann boundary conditions are applied on all the walls and the axis of symmetry for the chemical potential, the order parameters, and pressure. As for velocity, the normal component of velocity vanishes and the tangential component is mirrored on the axis of symmetry. The no slip condition is imposed on all the walls. An explicit finite difference method on a half-staggered grid is adopted, as shown in Fig. 2.

Fig. 2 A half staggered grid. It is noted that only pressure is stored at the cell center, while all the others are stored at the cell vertices.

As for the discretization schemes, traditional central difference schemes are employed to discretize diffusion terms, and upwind schemes to approximate convection terms. The discretization procedure of eqn (5) could be found in Kim.²² During a time marching step, computation starts from the evolution of c_i, and then proceeds to the intermediate velocity u*, followed by the updating of the pressure pⁿ⁺¹ and by the renewal of the velocity uⁿ⁺¹, thus completing one marching. It will, however, not stop until the time duration set is reached. The solution step is summarized as follows in Fig. 3, where the superscript n flags the previous time step and n + 1 the current. Fⁿ contains all the other terms in the momentum equation. The time step Δt is on the order of magnitude of 10⁻⁶ s.

Fig. 3 Flow chart of the algorithm.

In addition, parallel programming based on OpenMP is utilized to accelerate computation. The pressure Poisson equation in step 3 in Fig. 3 is tackled using a Red/Black SOR algorithm, which is a parallel version of the traditional SOR algorithm. Fig. 4 describes a domain partition pattern of the Red/Black SOR algorithm. Notice that this partition pattern is used only for the updating of pressure.

Fig. 4 Checkboard partition for the Red/Black SOR algorithm of the pressure Poisson equation in step 3 in Fig. 3. Right is the pseudo Fortran code.

As indicated in Fig. 4, their update is divided into two steps. When the first group (red, for instance) is renewed using the values only at the black points, the second group (black) is updated using the newest values only at the red points. In this way, data dependency could be eliminated. Data dependency happens when a processor is reading the value and another is modifying the value of the same point in the meantime. The processor reading the value may or may not get the newest value of that point depending on the capabilities of the two processors. This would however not happen in the chessboard ordering. Besides, in each step, the successive over relaxation (SOR) method could be employed.

3. Results and discussion

In this section, the mesh sensitivity study was conducted, followed by the validation of the model and by the investigation of the effect of the phase field mobility on bubble rising. Then, two distinct cases in triple phase flow were examined. The thermophysical quantities used are listed in Table 1, with those for Oliver oil found in Sahasrabudhe et al.³⁷ The word and/or numeral in the parenthesis indicates that the surface energy is about the interface between the substance of interest and that in the parenthesis. Air is denoted as component 2, water component 1 and Oliver oil component 3.

Table 1 Thermophysical quantities used in this paper

	Oliver oil (3)	Water (1)	Air (2)
Density (kg m^−3)	900	1000	1.18
Viscosity (mPa s)	75	1	0.0185
Surface energy (mN m^−1)	32 (air)	72 (air)
	16.82 (water)

---

3.1 Mesh sensitivity study

For phase field models, mesh convergence study means two things: how many nodes to resolve the thin interface and how thin the interface thickness should be compared with certain macroscopic length scales, for instance, bubble diameter. The first has been solved with the characteristic interface thickness ξ set equal to the spatial step Δx.³⁸ The second is to be decided with the following mesh sensitivity study, where an air bubble rising in water was simulated on grids of varied resolutions. A dimensionless number Cn = ξ/D, defined as the ratio of the characteristic interface thickness ξ to the bubble diameter D, is also employed to indicate the fineness of the grid.

For convenience of discussion, a couple of dimensionless number are defined.³⁹ They are listed as follows.

Eqn (9) defines the Reynolds (Re) number, measuring the relative importance between inertial and viscous forces, and eqn (10) defines the Eötvös (Eo) number, measuring the relative importance between gravitational and surface tension forces. The last is the Galilei (Ga) number. As the terminal velocity V_T is not a prior, the Re number is less used.

In the equations above, ρ_l denotes liquid density, V_T signals the terminal velocity of the bubble, μ_l stands for liquid viscosity, and σ surface tension coefficient. In this section, Re ∼ 108, Ga ∼ 155, and Eo ∼ 1 if the terminal velocity V_T is assumed to be on the order of magnitude of 0.04 m s⁻¹. The numerical outcome is shown in Fig. 5, where for instance Cn = Δx/(54Δx) = 1/54 means there are 54 cells across the bubble diameter.

Fig. 5 (a) Documents the history of the nose position and (b) depicts the bubble shape at 140 ms.

Fig. 5(a) plots the time evolution of the nose position of the bubble. Even with the coarsest grid, the captured position agrees quite well with the others. The relative error of the predicted top positions between the intermediate (Cn = 1/54) and finest grids (Cn = 1/108) at 140 ms is only 1.4%. The deviation in the finest grid around 40 ms is probably due to the stronger acceleration and to the choosing of the phase field mobility. As demonstrated in Section 3.2, the larger the phase field mobility, the higher the rising velocity. However, there is no consistent law to select the phase field mobility when performing the mesh refinement study for phase field modeling. Therefore, this is possibly caused by the fitting parameters.

Fig. 5(b) displays the bubble shape on various grids at 140 ms, where an evident loss of mass occurs if the coarsest grid is employed. As the mesh is refined, the problem is alleviated. Though the Cahn–Hilliard model can well conserve the total mass of a binary phase flow, the mass of one component may not be conserved, which has been also observed.^40,41 The shrinkage of drops can be reduced with the Cn number set below a critical value, typically on the order of magnitude of O(10⁻²) as suggested by Yue et al.⁴⁰

Given the result of the mesh sensitivity study here, the grid with Δx = 5 × 10⁻⁵ m or Cn = 1/54 is employed, unless otherwise stated, throughout the paper. Another issue is the choice of the phase field mobility, which has been adjusted according to the mesh size M ∼ 0.2Δx to achieve the sharp interface limit. Fig. 6 shows a sequence of bubble shapes and velocity field for the gird of Δx = 5 × 10⁻⁵ m.

Fig. 6 Bubble shape and velocity field at different instants on the grid of Δx = 5 × 10⁻⁵ m.

The leftmost column in Fig. 6 traces the history of bubble shapes, and the rest columns give velocity field distribution in the computational domain. Though a liquid jet below the bubble is induced as shown at 20 ms, it is not strong enough to pierce into the bubble. Moreover, the bubble, although slightly elliptical as shown at t = 20 ms and 140 ms, is almost spherical during the rising process.

The velocity of the top of the bubble is given in Fig. 7(b), which shows that it is nearly constant except at the very beginning when it experiences an acceleration. Chen et al.¹⁰ also noticed this.¹⁰ It is worth noting that the bubble is a prolate at t = 20 ms in Fig. 6, a shape that helps rise up, as was pointed out by Yang et al.¹⁸ who showed that a uniform vertical electric field elongates a bubble in the direction of rising, thereby speeding up the bubble. As the bubble passes the initial period of acceleration, the bubble shape appears to have been evolved as well. Inspection of Fig. 7(a) shows that the rising bubble, when reaching a steady state, becomes an oblate ellipsoid. This is consistent with the phase diagram proposed by Bhaga and Weber.⁴²

Fig. 7 (a) Deviation of bubble shape from spherical to elliptical and (b) history of the bubble rising velocity.

3.2 Validation of the model

Having performed mesh sensitivity studies, attention is paid to the comparison of computed results with experimental data to validate the phase field model. In the case above, surface tension dominates as is evident by the Eötvös number. In this section, inertia dominated cases were considered, where the experiment conducted by Sharaf et al.⁴³ was used for comparison. The drop radius is 19.27 mm and the numerical outcome is provided in Fig. 8. In this section, Ga ∼ 2960, and Eo ∼ 50, meaning strong inertia but weak surface tension effect.

Fig. 8 Comparison of bubble shape with experimental results by Sharaf et al.⁴³

Fig. 8 displays a bubble resembling a skirted cap at t = 6. Overall speaking, the comparison between numerical and experimental data is reasonable. As the bubble starts moving upwards due to buoyancy around t = 3, a vortex flow as in Fig. 6 will develop around the periphery, resulting in a large dynamic pressure on the bottom of the bubble. Since the dynamic pressure inside the bubble is much less than outside the bubble, the pressure difference dimples the bottom, as demonstrated in Fig. 8. The gas inside would be accelerated, rushing to the upper surface, where the pressure gradient is initially lower than on the bottom. As the bottom continues deforming, the gas inside pushes further the upper surface, rendering it to move.

Fig. 9 gives the speed distribution at particular instants. At t = 3 when the bottom has been deformed due to the liquid jet from the difference between buoyancy and drag resistance, it is clear that the vortex core comes near the skirt, facilitating its growing. The maximum speed shows up in the wake at t = 3. The uprising jet then brings the air in motion via its viscosity. As the jet pushed by pressure gradient continues forging ahead, the distance between the top and bottom however would cease increasing at the center line, because of a higher pressure developed in front of the jet, or more precisely beneath the top surface of the bubble. Thus a negative pressure gradient forms, reducing the velocity of the jet, as shown at t = 4. Meantime, the wake velocity also dwindles due to viscosity. Since the vortex core sits near the skirt, the skirt is elongated, which in turn thins the vortex pattern and leads to an expanded and flattened bottom, as displayed at t = 5 and t = 6.

Fig. 9 Speed distribution at different instants for the validation case.

3.3 Effect of the phase field mobility on an air bubble rising in water

In a phase field model, an important parameter is the phase field mobility. In general, it should not be too large lest it dampens convective motion, nor can it be too small since it can increase unwanted deformation. Not easily amenable to experimental observation, its choice is mostly empirical, often on the order of magnitude of χΔx² with χ being a tuning parameter. In the current model, M ∼ χΔx is chosen due to the incorporation of ξ⁻¹ into M and to the scale of ξ ∼ Δx. Different phase field mobilities may lead to distinct fluid behaviors. Therefore, it is of practical importance to check its effect on bubble rising.

A set of computed results is given in Fig. 10, where numerical configurations are the same as in Section 3.1, except for the phase field mobility, which is allowed to vary. The history of the nose position under a variety of phase field mobilities is provided in Fig. 10(a), and the bubble shape at 100 ms is shown in Fig. 10(b).

Fig. 10 (a) Effect of the phase field mobility on an air bubble rising in water. (b) Depicts bubble shapes under various phase field mobilities at 100 ms.

The effect of the phase field mobility on the rising speed is evident in Fig. 10(a), if not drastic. Basically, the rising velocity is roughly proportional to the phase field mobility. The little plateau formed by the orange diamond in Fig. 10(a) shows that the bubble has been in contact with the top awhile. Another matter of concern can be identified in Fig. 10(b), where a loss of mass is aggravated as the phase field mobility is increased. The problem can be alleviated through a reasonable choice of the phase field mobility and via a proper volume ratio between the bubble and the whole computational domain.

Though the mesh convergence study helps choose a proper spatial step, the terminal rising velocity has not been validated against existing experimental or theoretical data. Liu et al.³⁶ conducted experiments on air bubbles rising in water, and graphed the rising velocity against the dimensionless number, Eo = gD²(ρ_l − ρ_g)/σ. The definition here is slightly different from eqn (10), which neglects gas density. The rising velocity, derived from the nose position in Fig. 10(a) when M = 5 × 10⁻⁵ m³ s kg⁻¹, approaches that predicted by Liu et al.,³⁶ which is on the order of magnitude of O(10⁻¹). Consequently, the phase field mobility is fixed to 5 × 10⁻⁵ m³ s kg⁻¹ throughout the following sections.

3.4 An air bubble crossing the oil–water interface driven by buoyancy force

The foregoing sections deal only with binary phase flows, for which phase field models can be very useful. This section examines an air bubble rising in two different phases, with the lighter Oliver oil floating on the heavier water. Since the bubble may cross the interface between the other two phases, the effect of interfacial tensions may bring about some interesting results. Computational configuration is the same as in Fig. 1(b) and Fig. 11 depicts a sequence of snapshots showing an air bubble rising across the oil–water interface.

Fig. 11 A sequence of snapshots showing an air bubble rising across the oil–water interface.

Fig. 11 is obtained by extracting the contours of c₁ = c₂ = c₃ = 0.5 and by mirroring the right half with respect to the axis of symmetry. The air bubble at t = 30 ms is depicted intentionally in both green and red, since the bubble at this instant is completely submerged in water so that the contours of c₁ = c₂ = 0.5 coincide with each other. As the air bubble rises close to the oil–water interface, the flow in front of the bubble will deform the interface, driving it upwards. At t = 30 ms, the bubble appears to be spherical in shape, having an aspect ratio of around unity. As the bubble continues rising upwards, making its way through the oil–water interface around t = 40 ms, there is a net upward Laplace pressure driving it onwards, since the lower part of the bubble in contact with water has a larger interfacial tension than does the upper part in contact with oil. The bubble thus accelerates more in the rear.

This is worth elaborating upon. Integrating the surface tension force along the lower arc AB counterclockwise at t = 40 ms in Fig. 11, one has

In eqn (12) n is the unit normal perpendicular to but pointing away from the arc AB, and t is the unit tangent normal. It is to be noted that t_a and t_b are collinear with σ₁₂. Since surface tension σ₁₂ makes an angle θ₁ with the vertically upward z axis, the resultant lifting force of eqn (12) comes out as 2 × σ₁₂ × cos θ₁. In a similar manner, the effective drag force from surface tension σ₂₃ would be 2 × σ₂₃ × cos θ₂. Given σ₁₂ > σ₂₃ and a comparable θ₁ and θ₂, one deduces that the combined force, pointing vertically upwards, would pull the bubble up. In addition, the bottom experiences stages of shape evolution, from semicircular in the water to flattened across the interface, and then to semicircular again if it is completely wetted by oil as shown at t = 140 ms in Fig. 11.

As the bubble gradually departs from the oil–water interface, the extra lifting force caused by the interfacial tension dwindles, since the interfacial tension gradually becomes uniform and the integral in eqn (12), with the kernel σ₂₃dt, could be integrated along the whole periphery, meaning that it would vanish eventually. As a result, the bubble is subjected only to gravity, pressure, and viscous force. Completely submerged in oil, the bubble takes on again a spherical shape in accordance with the phase diagram proposed by Bhaga and Weber⁴² all the time from t = 70 ms to 140 ms.

Li et al.⁴⁴ proposed four regimes for a bubble ascending in two media, with the lighter fluid being upon the heavier fluid. The results depend on the relative magnitude of interfacial tensions. The bubble could cross the interface between the two media if σ₁₂ > σ₂₃ + σ₁₃ and σ₁₃ < σ₂₃ + σ₁₂, which is satisfied in Fig. 11.

Fig. 12 depicts the distribution of speed for particular instants. At t = 30 ms, the bubble is making its way through the interface, the process of which resembles an impact on a soft substrate. Two lateral flows develop around the bubble nose, as shown therein due to mass conservation, with the speed being the largest. A pair of vortexes are also seen at t = 30 ms. As time progresses to 40 ms, the lateral flow moves downward, nearing the equator of the bubble. Meantime, the maximum speed within the bubble is located around the center, since the air velocity must be larger than the rising velocity of the bubble, so that the stagnation pressure beneath the top surface is higher than on it, generating a positive pressure gradient to drive the bubble to move. As the bubble is about to leave the interface as displayed at t = 50 ms, a cusp in the speed contour within the bubble is found, which is caused the extra acceleration induced by surface tension, as explained above.

Fig. 12 Flow fields around the air bubble as it is crossing the oil–water interface. The lower scale bar is only for t = 70 ms.

Fig. 13 monitors the positions of the nose and the rear. These two lines are generally parallel to each other. Nevertheless, from t = 40 ms on, the distance between the nose and rear starts to decrease, reaching a minimum around t = 50 ms. Afterwards, it begins to increase until t = 70 ms and later remains essentially constant.

Fig. 13 History of the positions of the nose and rear as a bubble crosses the oil–water interface.

3.5 An air bubble rising behind an oil droplet in water

A pair of bubbles or a combination of a bubble and a droplet, both of different fluids, interact with each other when they encounter during the process of rising up. There may be a number of possibilities when several bubbles rise and interact. If one from behind is moving faster than the other, it may push up the slower-moving one without direct contact, or may penetrate into the bigger one and then emerge out from the opposite side. Here the rising and interaction of an air bubble and an oil droplet in water is considered. Numerical configurations could be found in Fig. 1(c). Computed results of bubble–droplet interactions are plotted in Fig. 14.

Fig. 14 (a) A sequence of snapshots showing bubble–droplet interaction in water (upper). (b) History of the positions of the nose and rear as a bubble rising behind an oil droplet (lower).

Fig. 14 presents one possible outcome of bubble–droplet interaction. Initially, the air bubble is below the oil droplet. The air bubble strives to catch up with the slower oil droplet by squeezing the water in its front. At t = 30 ms, the bubble and the oil droplet have been in contact with each other and the oil droplet is deformed by the rising air bubble. Since the water–air interfacial tension is larger than the oil–air one, a resultant extra force helps further deform the oil droplet or help the bubble squeeze into the oil droplet, as explained at t = 30 ms in Fig. 11. At t = 40 ms, the air bubble appears to be surrounded only partially, but at t = 64 ms is swallowed up completely by the oil droplet. Subsequently, the air bubble, while being contained in the oil droplet, continues to ascend and attempts to separate from it because of a higher rising velocity due to buoyancy.

The spreading coefficient, defined as S = σ₁₂ − (σ₂₃ + σ₁₃), determines whether the droplet can spread on the bubble or whether the bubble can be wetted by the droplet. The spreading would occur if S is positive, which is the case at t = 64 ms as shown in Fig. 14(a). Inspection of the results in Fig. 14(b) shows that the distance between the nose and rear of the air bubble is reduced when the air bubble is to be swallowed up by the oil droplet, because the lower part of the air bubble is subjected to a higher interfacial tension, hence a larger acceleration rate. Afterwards, the distance changes little since an equilibrium state has been established. A similar case was studied by Kalantarpour et al.³⁵ using a three-component phase field Lattice Boltzmann method. However, their model has some limitations: for instance, physical quantities like air viscosity are not assigned real values, but tailored for numerical convenience. Fig. 15 depicts the speed distribution for particular instants in Fig. 14. It is to be noted that the lower scale bar is only for t = 100 ms.

Fig. 15 Speed distribution for the bubble–droplet interaction at particular instants. The lower scale bar is only for t = 100 ms.

Fig. 15 indicates that when the bubble and the drop are to come in touch with each other, a higher pressure region would develop beneath the drop bottom to form a pressure gradient, driving the drop to rise. As they come in real touch, the bubble bottom would be deformed, resulting in an oblate. A large portion of air inside the bubble is also experiencing an thrust. As the bubble is being swallowed up by the drop, two pairs of vortexes develop as shown at t = 40 ms. As time proceeds to t = 100 ms, an equilibrium state established.

Computed results with the air–water and the air–oil interfacial tensions both set to 0.072 N m⁻¹ are plotted in Fig. 16. In this case, the air bubble experiences a uniform interfacial tension as it makes its way through the water–oil mixture. It is apparent from the figure that having dimpled the bottom of the oil droplet at t = 40 ms, the air bubble continues ascending, piercing into the oil droplet from below, and yet the droplet is unable to enclose the air bubble. Clearly, this is because the interfacial tension between air and oil is the same as between air and water herein. At t = 120 ms, the oil droplet evolves into a toroidal shape, signaling that the bubble is about to separate from the bubble. Besides, the spreading coefficient S herein is negative, indicating that full wetting of the bubble by the droplet is impossible, as shown in Fig. 16. It is worth noting that for the above two cases in Fig. 14 and 16, the rising of oil droplets is being accelerated by the air bubble, a process that can be of benefit to remove oil droplets in wastewater treatment.⁴⁵

Fig. 16 Bubble–droplet interaction when the air–water and air–oil interfacial tensions are the same.

Fig. 17 compares the nose positions of the air bubble under different conditions, where the air bubble is released at the same position. In general, the air bubble moves fastest when it is behind the oil droplet, and slowest when it is to cross water–oil interface, since the bubble becomes an oblate when crossing the interface, a geometric shape that hinders rising as verified by Yang et al.¹⁸

Fig. 17 Comparison of the nose position of the air bubble under different circumstances.

4. Concluding remarks

In this paper, a phase field model for ternary fluid flow was employed to examine bubble rising dynamics and bubble–droplet interaction in the gas-flotation technique. Three phase fields are updated each time step to track fluid–fluid interfaces. In single bubble rising dynamics, the phase field mobility was found to drastically influence the rising velocity. In the triple phase flow, the results indicate that by integrating the interfacial tension along the bubble perimeter, one would find an effective lifting force, pushing the bubble upwards as it tries to cross the oil–water interface. This resultant force could cause additional deformation on the bottom, inducing another clockwise vortex but retarding the upward motion. Besides, when chasing an oil droplet in water, the air bubble could merge with the oil droplet, floating together upwards at a faster pace. Though the 2D simulation could shed some light on the inner physics, the real situation is 3D. Therefore, further work could be done on fully 3D simulations in terms of interactions among a swarm of bubbles and droplets in a wide range of parameters, which may be made possible by the emerging GPU technology and may resemble more the factual situation in the gas-flotation technique.

Author contributions

Conceptualization, M. S.; methodology, M. S.; software, M. S.; validation, M. S.; formal analysis, M. S.; investigation, M. S.; resources, M. S.; data curation, M. S.; writing—original draft preparation, M. S.; writing—review and editing, B. Q. L.; visualization, M. S.; supervision, B. Q. L.; project administration, B. Q. L.; funding acquisition, M. S. All authors have read and agreed to the published version of the manuscript.

Conflicts of interest

None declared.

Acknowledgements

This research was funded by the Natural Science Foundation of the Jiangsu Higher Education Institutions of China (Grant No. 21KJB460034).

References

1. H. J. Ivey, Relationships between bubble frequency, departure diameter and rise velocity in nucleate boiling, Int. J. Heat Mass Transfer, 1967, 10(8), 1023–1040.
2. J. F. Klausner, R. Mei, D. M. Bernhard and L. Z. Zeng, Vapor bubble departure in forced convection boiling, Int. J. Heat Mass Transfer, 1993, 36(3), 651–662.
3. A. Khorram and S. Mortazavi, Direct numerical simulation of film boiling on a horizontal periodic surface in three dimensions using front tracking, Phys. Fluids, 2022, 34(5), 052117.
4. G. F. Bennett and R. W. Peters, The removal of oil from wastewater by air flotation: A review, Crit. Rev. Environ. Control, 1988, 18(3), 189–253.
5. Z. Bai, H. Wang and S. T. Tu, Oil–water separation using hydrocyclones enhanced by air bubbles, Chem. Eng. Res. Des., 2011, 89(1), 55–59.
6. M. Shirzadi, M. Ueda, K. Hada, T. Fukasawa, K. Fukui, Y. Mino, T. Tsuru and T. Ishigami, High-resolution numerical simulation of microfiltration of oil-in-water emulsion permeating through a realistic membrane microporous structure generated by focused ion beam scanning electron microscopy images, Langmuir, 2022, 38(6), 2094–2108.
7. E. Suess, G. Bohrmann, J. Greinert and E. Lausch, Flammable ice, Sci. Am., 1999, 281(5), 76–83.
8. X. L. Wang, G. E. Li and J. Du, Triaxial test study on marine sediments sample with simulated combustible ice, Appl. Mech. Mater., 2014, 580–583, 376–379.
9. C. Yan, Y. Li, Y. Cheng, J. Wei, W. Tian, S. Li and Z. Wang, Multifield coupling mechanism in formations around a wellbore during the exploitation of methane hydrate with CO₂ replacement, Energy, 2022, 245, 123283.
10. L. Chen, S. V. Garimella, J. A. Reizes and E. Leonardi, The development of a bubble rising in a viscous liquid, J. Fluid Mech., 1999, 387, 61–96.
11. F. Ren, B. W. Song and M. C. Sukop, Terminal shape and velocity of a rising bubble by phase-field-based incompressible Lattice Boltzmann model, Adv. Water Resour., 2016, 97, 100–109.
12. M. Vahabia, H. Hadavandmirzaeib, B. Kamkaric and H. Safarid, Interaction of a pair of in-line bubbles ascending in an Oldroyd-B liquid: A numerical study, Eur. J. Mech. B Fluids, 2021, 85, 413–429.
13. S. G. G. Stokes, Mathematical and Physical Papers, Cambridge University Press, London, 1980.
14. R. Davies and G. Taylor, The mechanics of large bubbles rising through extended liquids and through liquids in tubes, Proc R Soc Lond A Math Phys Sci., 1950, 375–390.
15. A. Tomiyama, G. Celata, S. Hosokawa and S. Yoshida, Terminal velocity of single bubbles in interfacial tension force dominant regime, Int. J. Multiphase Flow, 2002, 28, 1497–1519.
16. M. K. Tripathi, A. R. Premlata, K. C. Sahu and R. Govindarajan, Two initially spherical bubbles rising in quiescent liquid, Phys. Rev. Fluid, 2017, 2, 073601.
17. A. Fakhari and M. H. Rahimian, Phase-field modeling by the method of lattice Boltzmann equations, Phys. Rev. E, 2010, 81(3), 036707.
18. Q. Z. Yang, B. Q. Li, J. Y. Shao and Y. C. Ding, A phase field numerical study of 3D bubble rising in viscous fluids under an electric field, Int. J. Heat Mass Transfer, 2014, 78, 820–829.
19. X. Mao, V. Joshi and R. Jaiman, A variational interface-preserving and conservative phase-field method for the surface tension effect in two-phase flows, J. Comput. Phys., 2021, 433, 110166.
20. D. Jacqmin, Contact-line dynamics of a diffuse fluid interface, J. Fluid Mech., 2000, 402, 57–88.
21. S. Mirjalili, C. B. Ivey and A. Mani, Comparison between the diffuse interface and volume of fluid methods for simulating two-phase flows, Int. J. Multiphase Flow, 2019, 116, 221–238.
22. J. Kim, Phase field computations for ternary fluid flows, Comput. Methods Appl. Mech. Eng., 2007, 196(45–48), 4779–4788.
23. F. Boyer and C. Lapuerta, Study of a three component Cahn–Hilliard flow model, ESAIM: Math. Model. Numer. Anal., 2006, 40, 653–687.
24. F. Boyer and S. Minjeaud, Hierarchy of consistent n-component Cahn–Hilliard systems, Math. Models Methods Appl. Sci., 2014, 24, 2885–2928.
25. C. Y. Zhang, H. Ding, P. Gao and Y. L. Wu, Diffuse interface simulation of ternary fluids in contact with solid, J. Comput. Phys., 2016, 309, 37–51.
26. Q. Zhang and X. P. Wang, Phase field modeling and simulation of three-phase flow on solid surfaces, J. Comput. Phys., 2016, 319, 79–107.
27. S. Dong, Wall-bounded multiphase flows of N immiscible incompressible fluids: Consistency and contact-angle boundary condition, J. Comput. Phys., 2017, 338, 21–67.
28. H. Liang, J. Xu, J. Chen, Z. Chai and B. Shi, Lattice Boltzmann modeling of wall-bounded ternary fluid flows, Appl. Math. Model., 2019, 73, 487–513.
29. J. Kim, Phase-field models for multi-component fluid flows, Commun. Comput. Phys., 2012, 12(03), 613–661.
30. J. M. Park and P. D. Anderson, Diffuse-interface modeling of three-phase interactions, Appl. Phys. Lett., 2016, 108(19), 191604.
31. W. Liu and J. M. Park, Ternary modeling of the interaction between immiscible droplets in a confined shear flow, Phys. Rev. Fluid, 2022, 7(1), 013604.
32. W. Liu and J. M. Park, Numerical study on the engulfing behavior between immiscible droplets in a confined shear flow, Chem. Eng. Sci., 2023, 266, 118265.
33. J. T. Zhang, H. R. Liu and H. Ding, Head-on collision of two immiscible droplets of different components, Phys. Rev. Fluid, 2020, 32(8), 082106.
34. F. R. Ming, P. N. Sun and A. M. Zhang, Numerical investigation of rising bubbles bursting at a free surface through a multiphase SPH model, Meccanica, 2017, 52, 2665–2684.
35. R. Kalantarpour, A. Ebadi, S. M. Hosseinalipour and H. Liang, Three-component phase-field Lattice Boltzmann method with high density ratio and ability to simulate total spreading states, Comput Fluids., 2020, 204, 104480.
36. L. Liu, H. Yan, G. Zhao and J. Zhuang, Experimental studies on the terminal velocity of air bubbles in water and glycerol aqueous solution, Exp. Therm. Fluid Sci., 2016, 78, 254–265.
37. S. N. Sahasrabudhe, V. Rodriguez-Martinez, M. O'Meara and B. E. Farkas, Density, viscosity, and surface tension of five vegetable oils at elevated temperatures: Measurement and modeling, Int. J. Food Prop., 2017, 20(S2), S1965–S1981.
38. C. Zhou, P. Yue, J. J. Feng, C. F. Ollivier-Gooch and H. H. Hu, 3D phase-field simulations of interfacial dynamics in Newtonian and viscoelastic fluids, J. Comput. Phys., 2010, 229(2), 498–511.
39. J. Zhang, M. J. Ni and J. Magnaudet, Three-dimensional dynamics of a pair of deformable bubbles rising initially in line. Part 1. Moderately inertial regimes, J. Fluid Mech., 2021, 920, A16.
40. P. Yue, C. Zhou and J. J. Feng, Spontaneous shrinkage of drops and mass conservation in phase-field simulations, J. Comput. Phys., 2007, 223(1), 1–9.
41. C. Zhang and Z. Guo, Spontaneous shrinkage of droplet on wetting surface in phase-field model, Phys. Rev. E, 2019, 100(6–1), 061302.
42. D. Bhaga and M. E. Weber, Bubbles in viscous liquids: shapes, wakes and velocities, J. Fluid Mech., 1981, 105, 61–85.
43. D. M. Sharaf, A. R. Premlata, M. K. Tripathi, B. Karri and K. C. Sahu, Shapes and paths of an air bubble rising in quiescent liquids, Phys. Fluids, 2017, 29, 122104.
44. E. Li, S. Al-Otaibi, I. Vakarelski and S. Thoroddsen, Satellite formation during bubble transition through an interface between immiscible liquids, J. Fluid Mech., 2014, 744, R1.
45. M. Dudek and G. Oye, Microfluidic study on the attachment of crude oil droplets to gas bubbles, Energy Fuels, 2018, 32(10), 10513–10521.

---