# Numerical Investigation of Two In-Line Two-Dimensional Bubbles Rising in a Two-Dimensional Quiescent Ambient Liquid by a Conservative Phase-Field Lattice Boltzmann Method

Authors: Bin Liao; Zaihe Yang; Shanqun Chen
Journal: Discrete Dynamics in Nature and Society
 (2022)

---

## Abstract

The coalescence and interaction of two in-line two-dimensional bubbles rising in viscous ambient liquids were studied using two-dimensional simulations. A mass-conserving lattice Boltzmann model that was combined by the phase-field lattice Boltzmann equation (LBE) and the pressure-evolution lattice Boltzmann method (LBM) with a multiple-relaxation-time (MRT) collision operator was used to solve the Navier–Stokes equations for multiphase flow. Starting from the circle shape, the interfaces of bubbles during the rising processes are captured by the calculation of the conservative phase-field LBE. The influence of liquid viscosity and interface tension on the coalescence and interaction of bubbles was studied by varying the Archimedes numbers (Ar) from 1 to 300 and several Bond numbers (Bo) from 5 to 200. It was found that the bubble–bubble interaction is enhanced with the decrease of the liquid viscosity, and the outcome of the coalescence is promoted by the two vortex rings (liquid circulations) around the two bubbles. A comprehensive map of the coalescence regime was obtained. Four distinct coalescence regimes (namely, Central I, Central II, Edge, and Hug) were identified, and three critical Ar that can distinguish the above four regimes were defined. Moreover, the effects of Ar and Bo on the coalescence time and the relative velocity between the two bubbles were also investigated. The current work enhances the understanding of the coalescence and interaction of bubbles.


---

## Body

## 1. Introduction

The phenomenon of bubbles rising in a quiescent ambient liquid is common both in nature and industrial processes, such as bubble column reactors, and heat exchangers in many chemical and petrochemical applications [1, 2]. An in-depth understanding of bubbles rising in quiescent ambient liquid not only helps to grasp or predict the industrial processes but also has a relatively high scientific significance in the field of multiphase flow.The research on a single bubble rising in ambient liquids originated in the work of Leonardo Da Vinci as early as the sixteenth century [3]. As the phenomenon is increasingly important in industry, the mechanism of a single bubble rising in ambient liquids has been extensively studied. Bhaga and Weber [4] systematically studied the shapes, wakes, and velocities of a single rising bubble under Morton numbers (Mo) greater than 4 × 10−3 and several Reynolds numbers (Re). It was found that the bubbles trailed and closed, and laminar toroidal wakes appeared under Re < 110, but the wakes became open and unsteady under Re > 110. Mário [5] reported that a set of bubbles raise velocity experiments in a liquid column using water or glycerol and found that the bubble’s terminal velocity is strongly dependent on the dynamic viscosity effect. Lalanne et al. [6] focused on shape-oscillations of a gas bubble rising in another liquid. A linear theory was proposed to distinguish the oscillation modes of lower and higher frequencies. Yan et al. [7] investigated the drag coefficient of a single bubble rising in deionized water, and a new correlation combined with nondimensional numbers was proposed to calculate the fluctuation of the drag coefficient. Because it was challenging to generate initial spherical bubbles in liquids, and there are some uncontrollable factors in experiments, many of the essential works in this topic were carried out by numerical simulations. The numerical methods used in the simulation of a single bubble rising in ambient liquids are the immersed-boundary method of Peskin [8], the front-tracking (FT) method of Tryggvason and collaborators [9, 10], the volume-of-fluid (VOF) method of Scardovelli and Zaleski [11] and Tripathi et al. [12], the level-set (LS) method of Sussman et al. [13], the lattice Boltzmann method (LBM) of Chen and Doolen [14], and the moving particle semi-implicit (MPS) method of Koshizuka and Oka [15]. According to Sahu’s review on rising bubble dynamics in viscosity-stratified fluids [16], we consider that the understanding of a single bubble rising in ambient liquids is relatively rich based on the existing experimental investigation and numerical simulations.By contrast, the understanding of two bubbles rising in a quiescent ambient liquid is highly limited since studies regarding this phenomenon are relatively rare; fewer works on this topic have been reported. Chen et al. [17] studied the coalescence of two identical bubbles rising in a stationary liquid pool. They found that interface tension and horizontal expansion of the bubble caused the bridge between the bubbles to reduce with the bubble size during the beginning of the coalescence. Islam et al. [18] focused on the effects of the initial horizontal bubble interval, oblique alignment, and rheological properties of non-Newtonian fluids on two side-by-side bubbles rising in a xanthan gum solution. Feng et al. [19] investigated the detailed coalescence and conjunction of two in-line bubbles rising in high-viscosity liquids. They discovered two coalescence modes, that is, conjunct coalescence and coalescence, without a conjunction. In addition, they proposed a hydrodynamic stability theory to predict the conjunct time. Tripathi et al. [20] systemically studied the three-dimensional motion of two bubbles side by side and revealed the destabilizing nature of the interaction between the wakes of the two bubbles. Zhang et al. [21] systemically studied the effects of the Galilei number (Ga) and Eötvös number (Eo) on two in-line bubbles rising in viscous ambient liquids. They identified four distinct coalescence regimes and defined three critical Ga, which distinguish the four regimes. It is necessary to note that most of the existing studies focused on the two bubbles rising under the axial-symmetric computer domain. Balla et al. [22] focused on the dynamics of a pair of initially three-dimensional spherical drops rising side by side in a surrounding liquid, and revealed the difference between the liquid–liquid system and gas–liquid system. However, the works on two bubbles rising in ambient liquids under two-dimensional conditions are still lacking, and the relevant mechanism of the coalescence and interaction of two bubbles is not fully understood.In this paper, a computational study on two in-line, two-dimensional bubbles rising in a two-dimensional quiescent ambient liquid is carried out. A numerical technique combining the conservative phase-field lattice Boltzmann equation (LBE) [23] and the pressure-evolution LBM with a multiple-relaxation-time (MRT) collision operator [24, 25] is used to address the complex interface topologies and the flow field during the two bubbles rising in ambient liquids. The effects of liquid viscosity and interface tension on bubble coalescence are systematically investigated by varying Ar from 1 to 300 and Bo from 5 to 200. Compared to the literature, the ranges of Ar and Bo are significantly extended. In particular, the bubble shape evolution, flow field, coalescence time, and relative velocity are explored for bubbles rising in ambient liquids.


## 2. Numerical Methods

### 2.1. Conservative Phase-Field Lattice Boltzmann Equation (LBE) [23]

The governing equation used for tracking the interface of bubbles rising in a two-dimensional quiescent ambient liquid can be written as follows:(1)∂ϕ∂t+∇⋅ϕu=∇⋅M∇ϕ−1−4ϕ−ϕm2ξn,where ϕ is the phase-field variable that indicates each phase such that it is zero in one phase (light fluid ϕl) and one in the other (heavy fluid ϕh), ϕm=ϕl+ϕh/2 indicates the interface location, t stands for time, M is the mobility, u stands for the macroscopic velocity vector, ξ is the interface thickness, and n is the unit vector normal to the interface:(2)n=∇ϕ∇ϕ.Under the equilibrium condition, for an interface located atr0=0, the phase-field distribution at r assuming a hyperbolic tangent profile is required: (3)ϕr=ϕm+ϕh−ϕl2tanh2r−r0ξ.Equation (1) can be solved by an inherently conservative LBE as follows:(4)∂hi∂t+Ci⋅∇hi=−hi−hieqλϕ,where hi is the phase-field distribution function, λϕ is the phase-field relaxation rate, and Ci is the microscopic velocity sets. For the 2D nine-velocity lattice (D2Q9), as i=0, Ci=0,0; as i=1−4, Ci=ccosθi,sinθi, where θi=i−1/2π; and as i=5−8, Ci=2ccosθi,sinθi, where θi=i−5/2π+π/4, c=δx/δt is the lattice speed, while δx and δt are, respectively, the lattice length and time scale.The phase-field distribution function can be written as follows:(5)hieq=ϕΓi+φwiCi⋅n⟶,in which(6)Γi=wi1+Ci⋅ucs+Ci⋅u22cs2−u22cs2,φ=Mcs21−4ϕ−ϕm2ξ,where wi is the weight coefficient set for the D2Q9 lattice, as i=0, wi=4/9; as i=1−4, wi=1/9; and as i=5−8, wi=1/36. cs=δx/3δt is the lattice speed of sound.Referring to the method proposed in [26], the phase-field LBE equation (4) can be solved in two steps:(1)
Collision:(7a)hi∗=hi−hi−hieqτϕ′+1/2.(2)
Advection:(7b)∂hi∗∂t+Ci⋅∇hi∗=0,whereτϕ′=λϕ/δt is the dimensionless relaxation time and the relationship between mobility and relaxation time is M=τϕ′cs2δt.To retain the second-order accuracy of the LBM in time and space, a Lax–Wendroff finite difference scheme [27] is used to solve the advection equation (7b):(8)hix,t+δt=hi∗x,t−αhi∗x,t−hi∗x−CiΔxc,t−12α1−αhi∗x+CiΔxc,t−2hi∗x,t+hi∗x−CiΔxc,t,where α is the Courant–Friedrichs–Lewy (CFL) number. Here, the Lax–Wendroff finite difference scheme is solved by an adaptive grid refinement method [26], and the value of α is equal to 0.25.The next step is to update the phase field. We take the zeroth moment of the distribution function to calculate it, and the equation is the following:(9)ϕ=∑i=08hi.Correspondingly, the densityρ is updated by(10)ρ=ρl+ϕ−ϕlϕh−ϕlρh−ρl,where ρh and ρl are the densities of the heavy and light fluids, respectively.The gradient of the phase-field, namely the vector normal to the interface in equation (2) can be calculated by(11)∇ϕ=32cΔx∑i=18Ciwiϕx+CiΔxc−ϕx−CiΔxc.


### 2.2. Pressure-Evolution LBM with an MRT Collision Operator

To solve the Navier–Stokes equations for multiphase flow, the phase-field LBE, and the pressure-evolution LBM with an MRT collision operator are combined to propose a mass-conserving lattice Boltzmann model [24, 25].The general form of the discrete Boltzmann equation for nonideal fluids can be written as follows:(12)∂fi∂t+Ci⋅∇fi=−Λfi−fieq+Ci−u⋅Fcs2Γi,where fi is the particle distribution function, fieq=ρΓi is the equilibrium distribution function, and Λ is the generalized collision operator.The general forcing termF can be written as follows:(13)F=−∇p−ρcs2+Fb+Fs,where p is the pressure, Fb is the body force, Fs=μϕ∇ϕ is the interface tension force, μϕ=4ψϕ−ϕlϕ−ϕhϕ−ϕm−κ∇2ϕ is the chemical potential, while ψ and κ are related to the interface tension σ and interface width ξ by ψ=12σ/ξ and κ=3σξ/2, respectively.A new particle distribution function,gi=fics2+wip−ρcs2 is induced to replace the previous function, fi, and substituted into equation (12). Considering the incompressible limit, the Mach number (Ma) of the flow is much smaller Ma≪1. The term Γi−wi∇p∼ΟMa3 is negligible, leading to(14)∂gi∂t+Ci⋅∇gi=−Λgi−gieq+Ci−u⋅Γi−wi∇ρcs2+Γiμϕ∇ϕ+Fb,where gieq=fieqcs2+wip−ρcs2, ∇ρ=∂ρ/∂ϕ∇ϕ=ρh−ρl∇ϕ.Correspondingly, equation (14) can be further transferred to(15)∂gi∂t+Ci⋅∇gi=−Λgi−gieq+Γi−wiρh−ρlcs2+ΓiμϕCi−u⋅∇ϕ.+ΓiCi−u⋅Fb.Based on the method proposed in the previous section, equation (15) can be solved by splitting, and the resulting algorithm consists of collision and advection steps. The form of the collision step can be written as follows: (16)gi∗x,t=g¯ix,t−M−1S^Mg¯i−g¯ieq+δtΓi−wiρh−ρlcs2+ΓiμϕCi−u⋅∇Mϕ+δtΓiCi−u⋅Fb,where M is the orthogonal transformation matrix, M−1 is M’s inverse, S^ is the diagonal matrix, S^=diag1,1,1,1,1,1,1,sυ,sυ, where sυ=1/τ′+0.5, the kinematic viscosity υ=τ′cs2δt, and τ′ is the relaxation time explained by Lee and Lin [28].To maintain an explicit scheme,g¯i can be written as follows:(17)g¯i=gi−12S^g¯i−g¯ieq−δt2Γi−wiρh−ρlcs2+ΓiμϕCi−u⋅∇Cϕ−δt2ΓiCi−u⋅Fb,in which the modified equilibrium distribution function, g¯ieq, can be written as(18)g¯ieq=gieq−δt2Γi−wiρh−ρlcs2+ΓiμϕCi−u⋅∇Cϕ−δt2ΓiCi−u⋅Fb.The superscriptsC and M in equation (16)–(18) represent central and mixed finite differences, respectively, which can be written as follows:(19)1cCi⋅∇Cϕ=ϕx+CiΔx/c,t−ϕx−CiΔx/c,t2Δx,1cCi⋅∇Mϕ=−ϕx+2CiΔx/c,t+5ϕx+2CiΔx/c,t−3ϕx,t−ϕx−CiΔx/c,t4Δx.After the collision step, the advection step is carried out by the Lax–Wendroff scheme as follows:(20)g¯i∗x,t+δt=gi∗x,t−αgi∗x,t−gi∗x−CiΔx/c,t−12α1−αgi∗x+CiΔx/c,t−2gi∗x,t+gi∗x−CiΔx/c,t.Finally, the flow variables are calculated by(21)u=1ρcs2∑ig¯iCi+δt2μϕ∇Cϕ+Fb,(22)p=∑ig¯i+δt2ρh−ρlcs2u⋅∇Cϕ.


### 2.3. Numerical Validation

The benchmark study of a single bubble rising in ambient liquids is used to validate the reliability and accuracy of the proposed mass-conserving lattice Boltzmann model. A comparison between the present work and the MooNMD method [29] is made, and the results are shown in Figure 1. Two numerical simulations are conducted at different nondimensional parameters, which are Ar=35, Bo=10, and Ar=35, Bo=125, respectively. Here, Ar and Bo represent Archimedes number and Bond number, respectively. As shown in Figure 1(a), at the bubble diameter D=0.5 and the characteristic time t∗=3 (t∗=tgy/D), the interface topology of bubbles in present work is almost the same as the MooNMD results. Then, we compared the rising velocities, the definition of which is that the velocity of the bubble center of mass when the bubble rises, between the present work and the MooNMD method for t∗≤3, as shown in Figure 1(b). The rising velocities of bubbles in the present work are also in satisfactory agreement with the MooNMD results.Figure 1
Numerical results of a single bubble rising in ambient liquid for different nondimensional parameters atD=0.5, ρl/ρb=1000, μl/μb=100. (a) Interface topologies at t∗=3; (b) rising velocities for t∗≤3.
(a)(b)


## 2.1. Conservative Phase-Field Lattice Boltzmann Equation (LBE) [23]

The governing equation used for tracking the interface of bubbles rising in a two-dimensional quiescent ambient liquid can be written as follows:(1)∂ϕ∂t+∇⋅ϕu=∇⋅M∇ϕ−1−4ϕ−ϕm2ξn,where ϕ is the phase-field variable that indicates each phase such that it is zero in one phase (light fluid ϕl) and one in the other (heavy fluid ϕh), ϕm=ϕl+ϕh/2 indicates the interface location, t stands for time, M is the mobility, u stands for the macroscopic velocity vector, ξ is the interface thickness, and n is the unit vector normal to the interface:(2)n=∇ϕ∇ϕ.Under the equilibrium condition, for an interface located atr0=0, the phase-field distribution at r assuming a hyperbolic tangent profile is required: (3)ϕr=ϕm+ϕh−ϕl2tanh2r−r0ξ.Equation (1) can be solved by an inherently conservative LBE as follows:(4)∂hi∂t+Ci⋅∇hi=−hi−hieqλϕ,where hi is the phase-field distribution function, λϕ is the phase-field relaxation rate, and Ci is the microscopic velocity sets. For the 2D nine-velocity lattice (D2Q9), as i=0, Ci=0,0; as i=1−4, Ci=ccosθi,sinθi, where θi=i−1/2π; and as i=5−8, Ci=2ccosθi,sinθi, where θi=i−5/2π+π/4, c=δx/δt is the lattice speed, while δx and δt are, respectively, the lattice length and time scale.The phase-field distribution function can be written as follows:(5)hieq=ϕΓi+φwiCi⋅n⟶,in which(6)Γi=wi1+Ci⋅ucs+Ci⋅u22cs2−u22cs2,φ=Mcs21−4ϕ−ϕm2ξ,where wi is the weight coefficient set for the D2Q9 lattice, as i=0, wi=4/9; as i=1−4, wi=1/9; and as i=5−8, wi=1/36. cs=δx/3δt is the lattice speed of sound.Referring to the method proposed in [26], the phase-field LBE equation (4) can be solved in two steps:(1)
Collision:(7a)hi∗=hi−hi−hieqτϕ′+1/2.(2)
Advection:(7b)∂hi∗∂t+Ci⋅∇hi∗=0,whereτϕ′=λϕ/δt is the dimensionless relaxation time and the relationship between mobility and relaxation time is M=τϕ′cs2δt.To retain the second-order accuracy of the LBM in time and space, a Lax–Wendroff finite difference scheme [27] is used to solve the advection equation (7b):(8)hix,t+δt=hi∗x,t−αhi∗x,t−hi∗x−CiΔxc,t−12α1−αhi∗x+CiΔxc,t−2hi∗x,t+hi∗x−CiΔxc,t,where α is the Courant–Friedrichs–Lewy (CFL) number. Here, the Lax–Wendroff finite difference scheme is solved by an adaptive grid refinement method [26], and the value of α is equal to 0.25.The next step is to update the phase field. We take the zeroth moment of the distribution function to calculate it, and the equation is the following:(9)ϕ=∑i=08hi.Correspondingly, the densityρ is updated by(10)ρ=ρl+ϕ−ϕlϕh−ϕlρh−ρl,where ρh and ρl are the densities of the heavy and light fluids, respectively.The gradient of the phase-field, namely the vector normal to the interface in equation (2) can be calculated by(11)∇ϕ=32cΔx∑i=18Ciwiϕx+CiΔxc−ϕx−CiΔxc.


## 2.2. Pressure-Evolution LBM with an MRT Collision Operator

To solve the Navier–Stokes equations for multiphase flow, the phase-field LBE, and the pressure-evolution LBM with an MRT collision operator are combined to propose a mass-conserving lattice Boltzmann model [24, 25].The general form of the discrete Boltzmann equation for nonideal fluids can be written as follows:(12)∂fi∂t+Ci⋅∇fi=−Λfi−fieq+Ci−u⋅Fcs2Γi,where fi is the particle distribution function, fieq=ρΓi is the equilibrium distribution function, and Λ is the generalized collision operator.The general forcing termF can be written as follows:(13)F=−∇p−ρcs2+Fb+Fs,where p is the pressure, Fb is the body force, Fs=μϕ∇ϕ is the interface tension force, μϕ=4ψϕ−ϕlϕ−ϕhϕ−ϕm−κ∇2ϕ is the chemical potential, while ψ and κ are related to the interface tension σ and interface width ξ by ψ=12σ/ξ and κ=3σξ/2, respectively.A new particle distribution function,gi=fics2+wip−ρcs2 is induced to replace the previous function, fi, and substituted into equation (12). Considering the incompressible limit, the Mach number (Ma) of the flow is much smaller Ma≪1. The term Γi−wi∇p∼ΟMa3 is negligible, leading to(14)∂gi∂t+Ci⋅∇gi=−Λgi−gieq+Ci−u⋅Γi−wi∇ρcs2+Γiμϕ∇ϕ+Fb,where gieq=fieqcs2+wip−ρcs2, ∇ρ=∂ρ/∂ϕ∇ϕ=ρh−ρl∇ϕ.Correspondingly, equation (14) can be further transferred to(15)∂gi∂t+Ci⋅∇gi=−Λgi−gieq+Γi−wiρh−ρlcs2+ΓiμϕCi−u⋅∇ϕ.+ΓiCi−u⋅Fb.Based on the method proposed in the previous section, equation (15) can be solved by splitting, and the resulting algorithm consists of collision and advection steps. The form of the collision step can be written as follows: (16)gi∗x,t=g¯ix,t−M−1S^Mg¯i−g¯ieq+δtΓi−wiρh−ρlcs2+ΓiμϕCi−u⋅∇Mϕ+δtΓiCi−u⋅Fb,where M is the orthogonal transformation matrix, M−1 is M’s inverse, S^ is the diagonal matrix, S^=diag1,1,1,1,1,1,1,sυ,sυ, where sυ=1/τ′+0.5, the kinematic viscosity υ=τ′cs2δt, and τ′ is the relaxation time explained by Lee and Lin [28].To maintain an explicit scheme,g¯i can be written as follows:(17)g¯i=gi−12S^g¯i−g¯ieq−δt2Γi−wiρh−ρlcs2+ΓiμϕCi−u⋅∇Cϕ−δt2ΓiCi−u⋅Fb,in which the modified equilibrium distribution function, g¯ieq, can be written as(18)g¯ieq=gieq−δt2Γi−wiρh−ρlcs2+ΓiμϕCi−u⋅∇Cϕ−δt2ΓiCi−u⋅Fb.The superscriptsC and M in equation (16)–(18) represent central and mixed finite differences, respectively, which can be written as follows:(19)1cCi⋅∇Cϕ=ϕx+CiΔx/c,t−ϕx−CiΔx/c,t2Δx,1cCi⋅∇Mϕ=−ϕx+2CiΔx/c,t+5ϕx+2CiΔx/c,t−3ϕx,t−ϕx−CiΔx/c,t4Δx.After the collision step, the advection step is carried out by the Lax–Wendroff scheme as follows:(20)g¯i∗x,t+δt=gi∗x,t−αgi∗x,t−gi∗x−CiΔx/c,t−12α1−αgi∗x+CiΔx/c,t−2gi∗x,t+gi∗x−CiΔx/c,t.Finally, the flow variables are calculated by(21)u=1ρcs2∑ig¯iCi+δt2μϕ∇Cϕ+Fb,(22)p=∑ig¯i+δt2ρh−ρlcs2u⋅∇Cϕ.


## 2.3. Numerical Validation

The benchmark study of a single bubble rising in ambient liquids is used to validate the reliability and accuracy of the proposed mass-conserving lattice Boltzmann model. A comparison between the present work and the MooNMD method [29] is made, and the results are shown in Figure 1. Two numerical simulations are conducted at different nondimensional parameters, which are Ar=35, Bo=10, and Ar=35, Bo=125, respectively. Here, Ar and Bo represent Archimedes number and Bond number, respectively. As shown in Figure 1(a), at the bubble diameter D=0.5 and the characteristic time t∗=3 (t∗=tgy/D), the interface topology of bubbles in present work is almost the same as the MooNMD results. Then, we compared the rising velocities, the definition of which is that the velocity of the bubble center of mass when the bubble rises, between the present work and the MooNMD method for t∗≤3, as shown in Figure 1(b). The rising velocities of bubbles in the present work are also in satisfactory agreement with the MooNMD results.Figure 1
Numerical results of a single bubble rising in ambient liquid for different nondimensional parameters atD=0.5, ρl/ρb=1000, μl/μb=100. (a) Interface topologies at t∗=3; (b) rising velocities for t∗≤3.
(a)(b)


## 3. Problem Statement

Figure2 shows the schematic of the two in-line, two-dimensional bubbles rising in a quiescent ambient liquid problem. Simulations were conducted in a two-dimensional computational domain. The bounce back conditions [30] are imposed both on the top and bottom boundaries, and the periodic boundary conditions are both imposed on the left and right boundaries. Two initially identical circular bubbles are released at the centerline of the domain. The constant initial separation of S = 2.36R is considered, where R is the bubble radius. The initial velocity and pressure field of the domain are set to zero. In addition, an adaptive grid refinement method as mentioned above was used, and the grid setting are specified as the maximum and minimum cell sizes were Δmax=L/16 and Δmin=L/256, where L is the horizontal scale of the computational domain, respectively. The maximum refinement level was set to 8.Figure 2
Schematic of the two in-line, two-dimensional bubbles rising problem.The bubble rising can be strongly affected by four nondimensional parameters [3], that is, Archimedes number (Ar), (Ar=ρl−ρbgy1/2D3/2/μl, where ρl and ρb are the densities of ambient liquid and bubble, respectively. gy is the magnitude of the gravitational acceleration in the y-direction, D is the bubble diameter, μl is the liquid viscosity), Bond number (Bo), (Bo=ρl−ρbgyD2/σ, where σ is the interface tension between the bubble and ambient liquid), density ratio ρr (ρr=ρl/ρb) and viscosity ratio μr (μr=μl/μb). The Ar and Bo represent the ratio of external force to internal viscous force and gravitational force to interface tension force, respectively. The density ratio and viscosity ratio mainly represent the inertia effect and shear effect on the bubbles, respectively. Considering the common phenomenon of air bubbles rising in water, the density and viscosity ratios of ambient liquid to bubbles are fixed to 1000 and 100 in this paper, respectively. It is necessary to note that both ρr and μr remain constant except when stated otherwise.


## 4. Results and Discussions

### 4.1. Evolution of Shape and Flow Field Analysis

First, we conducted the simulations under a fixed Bo number (Bo = 50) to investigate the influence of the ambient liquid viscosity (Ar) on the rising behaviors of two in-line two-dimensional bubbles. Considering that the complexities of the bubbles rising in ambient liquids, the liquid viscosity condition (Ar) is separated into three regions: low Ar, medium Ar, and high Ar.


#### 4.1.1. At Low Ar

The rising behaviors (mainly, the coalescence process) of bubbles in quiescent ambient liquid are shown in Figure3 at Bo = 50 and Ar = 1. We know that, at Bo = 50 and Ar = 1, the viscous force and interface tension both play important roles in the rising behaviors of bubbles. Therefore, the inertial force is relatively negligible, and the drag force induced by the relatively high-viscosity causes the bubbles to rise slowly. The interface tension constrains the bubble deformation. As shown in the first frame of Figure 3, the generated jet is rather weak, and only slight deformation is observed. Moreover, from the streamline distribution, it is found that the vortex cores are relatively far from the two bubbles, which causes the liquid circulations to be weak as well. Consequently, the upper interface of the following bubble alters slightly, and the lower interface of the following bubble is slightly elongated. As time increases, because of the shielding effect of the wake flow behind the leading bubble [31, 32], the upper interface of the following bubble is flattened by the dimple formed in the leading bubble. During the coalescence process, the following bubble remains elongated, and while the lower interface of the following bubble becomes slimmer; a smooth tail even appears at t∗=7.2.Figure 3
Coalescence process between two bubbles and corresponding flow field at Bo = 50, Ar = 1,ρr = 1000, μr = 100 and S/R = 2.36.On the other hand, the deformation of the leading bubble shows a reverse evolution during the rising. From left to right of Figure3, the hat shape of the leading bubble becomes more oblate. During the coalescence process, the leading bubble shows a slight deceleration before the coalescence. Consequently, the following bubble continuously approaches the leading bubble. Finally, the coalescence between two bubbles occurs at tc = 7.541, which represents the characteristic coalescence time. As shown in the third frame of Figure 3, the coalescence starts from a single contact point of the two bubbles, which is located almost at the symmetry axis of the two bubbles.Figure4 shows the rising behaviors of two bubbles at Bo = 50 and Ar = 10. As Ar increases to 10, the liquid viscosity decreases significantly. Compared to the case of Ar = 1, the influence of inertia increases a little based on the reduction of drag. Correspondingly, the generated jet becomes stronger, and the deformation of the following bubble increases slightly. Because the vortex cores are relatively near the two bubbles, the liquid circulations become stronger as well. Moreover, because the wake speed becomes faster, the shielding effect of the wake flow behind the leading bubble is strengthened. As shown in the second frame of Figure 4, the deformation of the following bubble becomes more significant, while the lower interface of the following bubble flattens rapidly. On the other hand, the leading bubble shows a similar deformation to the case of Ar = 1. As time increases, the following bubble approaches and contacts the leading bubble. The coalescence occurs at tc = 1.973, which is much earlier than the case of Ar = 1. After bubble coalescence occurs, the liquid circulations accelerate the lower face of the merged bubble, and a larger hat shaped bubble is formed at t∗=2.47.Figure 4
Coalescence process between two bubbles and corresponding flow field at Bo = 50, Ar = 10,ρr = 1000, μr = 100 and S/R = 2.36.It is worth mentioning that the detailed coalescence process of Ar = 10 is different from the case of Ar = 1. Although the coalescence region is located at the central part of both bubbles in the two cases, the contact area is different. In the case of Ar = 1, because of the weak liquid circulations, the impact between the two bubbles is very slight, and the contact area is small enough to be considered a single point. The liquid circulations become stronger in the case of Ar = 10. Under such conditions, there is a slightly stronger impact between the two bubbles during the coalescence process. Consequently, the contact area is enlarged, and the coalescence no longer starts from a contact point but from several points at the same time. It is necessary to note that the contact points are still located at the central part of the two bubbles.


#### 4.1.2. At Medium Ar

As Ar increases to 80, the rising behaviors of the two bubbles are shown in Figure5. According to References [3, 33], at Bo = 50 and Ar = 10, the inertial force and interface tension are of the same order. Namely, they both play significant roles in bubble dynamics. Because liquid viscosity decreases, the drag force has less influence on the rising behaviors of bubbles, and the generated jets are rather strong. Therefore, both the interfaces of the two bubbles are easier to deform. As shown in the first frame of Figure 5, an obvious deformation of the following bubble is observed. The lower interface of the following bubble is pushed upward, and an upward dimple is formed. Because the distance between the vortex cores to the lower interface of the following bubble is significantly smaller than the upper interface, the lower interface receives shearing forces induced by the strengthened liquid circulations. Consequently, the lower interface ascends much more rapidly than the upper interface. Then, the following bubble decelerates suddenly due to the obvious oblate deformation. On the other hand, the leading bubble shows a similar deformation as the smaller Ar cases. Namely, the two bubbles both experience an oblate deformation situation, but the deformation degree of the leading bubble is always greater than that of the following bubble. As shown in the second frame of Figure 5, the topology of the two bubbles seems to be a “hat” being put on a “head,” where the leading bubble is the “hat,” and the following bubble is the “head.” It is necessary to note that, compared to the above two smaller Ar cases, the degree of the “putting on” between the two bubbles is strengthened by the decreased liquid viscosity (increased Ar).Figure 5
Coalescence process between two bubbles and corresponding flow field at Bo = 50, Ar = 80,ρr = 1000, μr = 100 and S/R = 2.36.As time increases, under the continuous action of the two vortex rings (liquid circulations), the oblate deformation of the two bubbles becomes more significant. That is, the dimple of the lower interface of the following bubble gradually expands and squeezes most of the bubble to both sides. Then, with the horizontal development of two bubbles, the probability of potential collision on the central part of the two bubbles is decreased. Moreover, as shown in the third frame of Figure5, as the vortex cores are gradually closed to the edge of the bubbles, the shearing effect increases the probability of bubble coalescence in peripherals. Consequently, the coalescence between the two bubbles starts from the peripheral part instead of the central part. Interestingly, the coalescence time is tc = 2.097, which is slightly larger than the case of Ar = 10. As shown in the fourth frame of Figure 5, after the contact area first appears in the peripheral part of the two bubbles, the coalescence subsequently expands to the central part.


#### 4.1.3. At High Ar

Figure6 gives the results of the rising behaviors of the bubbles at Bo = 50 and Ar = 300. Under such a condition, with the liquid viscosity reduced significantly, the bubble dynamics are almost unaffected by the drag force. Instead, the bubbles mostly face inertia and gravity. Because the scale of the interface tension is much smaller than the inertia, the interface tension plays a relatively unimportant role in this case. Therefore, the liquid circulations and the generated jets become quite strong. As shown in the first frame of Figure 6, during the approach of the two bubbles, because of the strengthened liquid jet behind the leading bubble, the following bubble significantly deforms and elongates in the vertical direction. Meanwhile, the lower interface of the following bubble is driven by the liquid jet, and an upward dimple is formed. At this time, the shape of the lower interface resembles a skirt. Compared to the first frame of Figure 5, compared to the previous case (Ar = 80), the skirt topology is more significant. Subsequently, the following bubble rising behaviors are pushed by the continuous action of the two vortex rings. As shown in the second frame of Figure 6, most of the gas inside the following bubble moves upward rapidly. Correspondingly, the gas in the tail of the following bubble decreases, and the “neck” of the skirt becomes thinner. Ultimately, because the pressure induced by the ambient liquid flow is rather strong, the interface tension is unable to maintain the existing equilibrium. As shown in the third frame of Figure 6 (t∗ = 1.78), two satellite bubbles are shed from the tail of the following bubble. The phenomenon can be captured for a single bubble under the condition of a rather large Eo [3, 33].Figure 6
Coalescence process between two bubbles and corresponding flow field at Bo = 50, Ar = 300,ρr = 1000, μr = 100 and S/R = 2.36.On the other hand, for the leading bubble, with the continuous action of the two strengthened vortex rings, its axial symmetry thickness decreases continuously until central breakup occurs. As shown in the second frame of Figure6, the leading bubble is divided into two parts. It is necessary to note that the thinning and central breakup of the leading bubble is promoted by the inertia effect induced by the rather fast approach of the following bubble. After the breakup of the leading bubble, the following bubble continues to rise through the hollow space of the leading bubble. Meanwhile, the two parts of the leading bubble spread out on both sides. Then, as shown in the third frame of Figure 6, similar to the previous case (Ar = 80), coalescence occurs in the peripheral part of the two bubbles. However, we consider that the coalescence processes are different between the two cases (Ar = 80 and Ar = 300). The following two aspects can explain the difference. First, the leading bubble is divided into two parts in the case of Ar = 300, which does not occur in the case of Ar = 80. The rising behaviors in the two cases are different. Second, in the case of Ar = 300, coalescence occurs after the breakup of the leading bubble. Although the coalescence area is located at the peripheral part of two bubbles, the coalescence process looks more like the leading bubble’s two arms, giving the following bubble a “hug.” To better interpret this type of coalescence, Figure 7 gives detailed results of the coalescence process at Bo = 200 and Ar = 200. As shown in the third frame of Figure 7, coalescence occurs when the inner interface of the leading bubble is in contact with the lower interface of the following bubble, namely, the contact area is located on the backside of the following bubble, rather than its peripheral part. Consequently, we consider that the coalescence process at Bo = 50 and Ar = 300 is more appropriately described by using the word “hug.”Figure 7
Coalescence process between two bubbles at Bo = 200, Ar = 200,ρr = 1000, μr = 100 and S/R = 2.36.


### 4.2. Ar-Bo Map for the Coalescence Regime and Critical Ar

In the previous section, we established that distinct rising behaviors and coalescence processes occurred under varying Ar and fixed Bo values. To fully understand the relationship between the coalescence regime and the condition parameters, a full picture of the effect of Ar-Bo on the coalescence regime should be plotted.Above all, we classify the coalescence regime into three categories: Central, Edge, and Hug regimes based on the bubble shape topologies. The Central regime is the most common regime in our simulations. In this regime, the coalescence starts from single or several contact points at the central part of the two bubbles (e.g., the third frames of Figures3 and 4). The Edge regime refers to the coalescence that starts from the peripheral parts of the two bubbles (e.g., the third frames of Figure 5). In the Hug regime, the leading bubble breaks into two parts and gives the following bubble a “hug” before the coalescence (e.g., the third frames of Figures 6 and 7). Moreover, the Central coalescence regime was subdivided into two regimes, that is, Central I and Central II. In the Central I regime, the coalescence starts from a single point, which is located at the symmetry axis of the two bubbles. While in the Central II regime, the coalescence starts from several contact points, which are located in the central part near the symmetry axis.Correspondingly, an Ar-Bo map of coalescence regimes is plotted, as shown in Figure8. Obviously, the Central I regime generally occurs at low Ar or low Bo, and the Hug regime always occurs under high-condition parameters (Bo ≥ 50 and Ar ≥ 150). The Central II regime occurs under most of the condition parameters. Interestingly, when Bo is rather small (Bo ≤ 10), even if Ar is as large as 500, the Hug regime is still not observed. Accordingly, we consider that the Edge regime may not translate to the Hug regime even if the liquid viscosity is small enough to be ignored. Moreover, the Edge regime only exists under certain condition parameters (Bo ≤ 75 and Ar ≥ 60) in our simulations. This means that the interface tension and the liquid viscosity are both essential for the Edge regime. However, when the interface tension is rather small (Bo ≥ 100), the Edge regime might not always occur.Figure 8
Ar-Bo map of coalescence regimes plotted as Bo versus Ar.To distinguish the four regimes, three critical Ar values (namely, Ar1, Ar2, and Ar3) are defined as the boundaries among them. As shown in Figure 8, several relationships between the critical Ar and Bo are obtained. First, with the increase of Bo, Ar1 decreases continuously and Ar2 decreases first and then increases, which is different from the results under the axisymmetric computational domain in Reference [21]. The changes in Ar1 or Ar2 suggest that the reduction in interface tension promotes the gas in bubbles to transfer faster in the vertical direction, which makes the contact area between the two bubbles larger. Second, because the leading bubble is more likely to break up into two parts under such conditions, Ar3 decreases with increasing Bo. We consider that the Hug regime always occurs under relatively low-interface tension.


### 4.3. Effect of Ar on Coalescence Time at Different Bo Values

Interface tension is an important factor in the coalescence time of two rising bubbles [34–36]. However, in the two-dimensional domain, the influence of liquid viscosity on the coalescence time under a single interface tension condition has not yet been systematically studied. To obtain a comprehensive understanding, we investigated the changes in coalescence time with varying Ar at several Bo values (Bo = 5, 10, 25, 50, 100, 200) in this section, where Ar corresponds to the liquid viscosity and several Bo values corresponds to different degrees of interface tension.Figure9 gives the results of coalescence time tc versus Ar at different Bo values. In general, with the increase in Ar, tc decreases basically under all Bo values, suggesting that irrespective of Bo, Ar shows a significant effect on tc. With the increase in Bo, tc shows a certain extent of decline, suggesting that the decreasing interface tension promotes coalescence. When the liquid viscosity is rather high (Ar ≤ 3), bubbles coalescence occurs much later. As Ar increases (Ar ≤ 20), according to the comparative interpretation of Figures 3 and 4, the coalescence is considerably promoted by the strengthened inertia effect. However, when Ar ≥ 30, tc shows little change and basically appears to stabilize. According to the interpretation of Figure 5, this implies that the influence of liquid viscosity on tc can be negligible under such conditions. It is necessary to note that, under Bo ≥ 100 and Ar ≥ 125, tc suddenly increases slightly, which is caused by the transition of the Edge to the Hug coalescence regime. According to the interpretation of Figures 3, 5, and 6, compared to the Edge regime, tc in the Hug regime is always considerably larger.Figure 9
Effect of Ar on coalescence time at different Bo.


### 4.4. Effect of Ar and Bo on Relative Rising Velocity

The influence of interface tension (Bo) and liquid viscosity (Ar) on bubble relative rising velocity Ur are illustrated in Figures10 and 11, respectively. In Figure 10, the effect of Ar on Ur is investigated by using three representative Ar values at Bo = 50. In Figure 11, the effect of Bo on Ur is investigated by using three representative Bo values at Ar = 25. In these figures, the ends of the curves all correspond to the beginning of bubble coalescence. The relative rising velocity is specifically analyzed from two views, namely, versus time history (Figures 10(a) and 11(a)) and dependence on instantaneous distance (S/R) between the two bubbles (Figures 10(b) and 11(b)). In general, the time evolution of Ur depends significantly on Ar, Bo, and the instantaneous distance, S/R, between the two bubbles.Figure 10
The effect of Ar on the instantaneous relative velocity at Bo = 50. (a) Relative velocity versus time; (b) relative velocity versus bubble center distance.
(a)(b)Figure 11
The effect of Bo on the instantaneous relative velocity at Ar = 25. (a) Relative velocity versus time; (b) relative velocity versus bubble center distance.
(a)(b)The relationship among the relative velocity Ur, Ar, and instantaneous distance,S/R, can be included as follows: first, as time increases and S/R decreases, Ur first increases and then decreases, which is a result of gas transfers inside the leading and following bubbles. Second, the slopes of the Ur curves are enlarged by increasing Ar, which indicates that low-liquid viscosity promotes gas transfer, increasing the relative velocity. It is necessary to note that, the relationships among the relative velocity Ur, Bo, and instantaneous distance, S/R, in Figure 11 are similar to the results in Figure 10.


## 4.1. Evolution of Shape and Flow Field Analysis

First, we conducted the simulations under a fixed Bo number (Bo = 50) to investigate the influence of the ambient liquid viscosity (Ar) on the rising behaviors of two in-line two-dimensional bubbles. Considering that the complexities of the bubbles rising in ambient liquids, the liquid viscosity condition (Ar) is separated into three regions: low Ar, medium Ar, and high Ar.


### 4.1.1. At Low Ar

The rising behaviors (mainly, the coalescence process) of bubbles in quiescent ambient liquid are shown in Figure3 at Bo = 50 and Ar = 1. We know that, at Bo = 50 and Ar = 1, the viscous force and interface tension both play important roles in the rising behaviors of bubbles. Therefore, the inertial force is relatively negligible, and the drag force induced by the relatively high-viscosity causes the bubbles to rise slowly. The interface tension constrains the bubble deformation. As shown in the first frame of Figure 3, the generated jet is rather weak, and only slight deformation is observed. Moreover, from the streamline distribution, it is found that the vortex cores are relatively far from the two bubbles, which causes the liquid circulations to be weak as well. Consequently, the upper interface of the following bubble alters slightly, and the lower interface of the following bubble is slightly elongated. As time increases, because of the shielding effect of the wake flow behind the leading bubble [31, 32], the upper interface of the following bubble is flattened by the dimple formed in the leading bubble. During the coalescence process, the following bubble remains elongated, and while the lower interface of the following bubble becomes slimmer; a smooth tail even appears at t∗=7.2.Figure 3
Coalescence process between two bubbles and corresponding flow field at Bo = 50, Ar = 1,ρr = 1000, μr = 100 and S/R = 2.36.On the other hand, the deformation of the leading bubble shows a reverse evolution during the rising. From left to right of Figure3, the hat shape of the leading bubble becomes more oblate. During the coalescence process, the leading bubble shows a slight deceleration before the coalescence. Consequently, the following bubble continuously approaches the leading bubble. Finally, the coalescence between two bubbles occurs at tc = 7.541, which represents the characteristic coalescence time. As shown in the third frame of Figure 3, the coalescence starts from a single contact point of the two bubbles, which is located almost at the symmetry axis of the two bubbles.Figure4 shows the rising behaviors of two bubbles at Bo = 50 and Ar = 10. As Ar increases to 10, the liquid viscosity decreases significantly. Compared to the case of Ar = 1, the influence of inertia increases a little based on the reduction of drag. Correspondingly, the generated jet becomes stronger, and the deformation of the following bubble increases slightly. Because the vortex cores are relatively near the two bubbles, the liquid circulations become stronger as well. Moreover, because the wake speed becomes faster, the shielding effect of the wake flow behind the leading bubble is strengthened. As shown in the second frame of Figure 4, the deformation of the following bubble becomes more significant, while the lower interface of the following bubble flattens rapidly. On the other hand, the leading bubble shows a similar deformation to the case of Ar = 1. As time increases, the following bubble approaches and contacts the leading bubble. The coalescence occurs at tc = 1.973, which is much earlier than the case of Ar = 1. After bubble coalescence occurs, the liquid circulations accelerate the lower face of the merged bubble, and a larger hat shaped bubble is formed at t∗=2.47.Figure 4
Coalescence process between two bubbles and corresponding flow field at Bo = 50, Ar = 10,ρr = 1000, μr = 100 and S/R = 2.36.It is worth mentioning that the detailed coalescence process of Ar = 10 is different from the case of Ar = 1. Although the coalescence region is located at the central part of both bubbles in the two cases, the contact area is different. In the case of Ar = 1, because of the weak liquid circulations, the impact between the two bubbles is very slight, and the contact area is small enough to be considered a single point. The liquid circulations become stronger in the case of Ar = 10. Under such conditions, there is a slightly stronger impact between the two bubbles during the coalescence process. Consequently, the contact area is enlarged, and the coalescence no longer starts from a contact point but from several points at the same time. It is necessary to note that the contact points are still located at the central part of the two bubbles.


### 4.1.2. At Medium Ar

As Ar increases to 80, the rising behaviors of the two bubbles are shown in Figure5. According to References [3, 33], at Bo = 50 and Ar = 10, the inertial force and interface tension are of the same order. Namely, they both play significant roles in bubble dynamics. Because liquid viscosity decreases, the drag force has less influence on the rising behaviors of bubbles, and the generated jets are rather strong. Therefore, both the interfaces of the two bubbles are easier to deform. As shown in the first frame of Figure 5, an obvious deformation of the following bubble is observed. The lower interface of the following bubble is pushed upward, and an upward dimple is formed. Because the distance between the vortex cores to the lower interface of the following bubble is significantly smaller than the upper interface, the lower interface receives shearing forces induced by the strengthened liquid circulations. Consequently, the lower interface ascends much more rapidly than the upper interface. Then, the following bubble decelerates suddenly due to the obvious oblate deformation. On the other hand, the leading bubble shows a similar deformation as the smaller Ar cases. Namely, the two bubbles both experience an oblate deformation situation, but the deformation degree of the leading bubble is always greater than that of the following bubble. As shown in the second frame of Figure 5, the topology of the two bubbles seems to be a “hat” being put on a “head,” where the leading bubble is the “hat,” and the following bubble is the “head.” It is necessary to note that, compared to the above two smaller Ar cases, the degree of the “putting on” between the two bubbles is strengthened by the decreased liquid viscosity (increased Ar).Figure 5
Coalescence process between two bubbles and corresponding flow field at Bo = 50, Ar = 80,ρr = 1000, μr = 100 and S/R = 2.36.As time increases, under the continuous action of the two vortex rings (liquid circulations), the oblate deformation of the two bubbles becomes more significant. That is, the dimple of the lower interface of the following bubble gradually expands and squeezes most of the bubble to both sides. Then, with the horizontal development of two bubbles, the probability of potential collision on the central part of the two bubbles is decreased. Moreover, as shown in the third frame of Figure5, as the vortex cores are gradually closed to the edge of the bubbles, the shearing effect increases the probability of bubble coalescence in peripherals. Consequently, the coalescence between the two bubbles starts from the peripheral part instead of the central part. Interestingly, the coalescence time is tc = 2.097, which is slightly larger than the case of Ar = 10. As shown in the fourth frame of Figure 5, after the contact area first appears in the peripheral part of the two bubbles, the coalescence subsequently expands to the central part.


### 4.1.3. At High Ar

Figure6 gives the results of the rising behaviors of the bubbles at Bo = 50 and Ar = 300. Under such a condition, with the liquid viscosity reduced significantly, the bubble dynamics are almost unaffected by the drag force. Instead, the bubbles mostly face inertia and gravity. Because the scale of the interface tension is much smaller than the inertia, the interface tension plays a relatively unimportant role in this case. Therefore, the liquid circulations and the generated jets become quite strong. As shown in the first frame of Figure 6, during the approach of the two bubbles, because of the strengthened liquid jet behind the leading bubble, the following bubble significantly deforms and elongates in the vertical direction. Meanwhile, the lower interface of the following bubble is driven by the liquid jet, and an upward dimple is formed. At this time, the shape of the lower interface resembles a skirt. Compared to the first frame of Figure 5, compared to the previous case (Ar = 80), the skirt topology is more significant. Subsequently, the following bubble rising behaviors are pushed by the continuous action of the two vortex rings. As shown in the second frame of Figure 6, most of the gas inside the following bubble moves upward rapidly. Correspondingly, the gas in the tail of the following bubble decreases, and the “neck” of the skirt becomes thinner. Ultimately, because the pressure induced by the ambient liquid flow is rather strong, the interface tension is unable to maintain the existing equilibrium. As shown in the third frame of Figure 6 (t∗ = 1.78), two satellite bubbles are shed from the tail of the following bubble. The phenomenon can be captured for a single bubble under the condition of a rather large Eo [3, 33].Figure 6
Coalescence process between two bubbles and corresponding flow field at Bo = 50, Ar = 300,ρr = 1000, μr = 100 and S/R = 2.36.On the other hand, for the leading bubble, with the continuous action of the two strengthened vortex rings, its axial symmetry thickness decreases continuously until central breakup occurs. As shown in the second frame of Figure6, the leading bubble is divided into two parts. It is necessary to note that the thinning and central breakup of the leading bubble is promoted by the inertia effect induced by the rather fast approach of the following bubble. After the breakup of the leading bubble, the following bubble continues to rise through the hollow space of the leading bubble. Meanwhile, the two parts of the leading bubble spread out on both sides. Then, as shown in the third frame of Figure 6, similar to the previous case (Ar = 80), coalescence occurs in the peripheral part of the two bubbles. However, we consider that the coalescence processes are different between the two cases (Ar = 80 and Ar = 300). The following two aspects can explain the difference. First, the leading bubble is divided into two parts in the case of Ar = 300, which does not occur in the case of Ar = 80. The rising behaviors in the two cases are different. Second, in the case of Ar = 300, coalescence occurs after the breakup of the leading bubble. Although the coalescence area is located at the peripheral part of two bubbles, the coalescence process looks more like the leading bubble’s two arms, giving the following bubble a “hug.” To better interpret this type of coalescence, Figure 7 gives detailed results of the coalescence process at Bo = 200 and Ar = 200. As shown in the third frame of Figure 7, coalescence occurs when the inner interface of the leading bubble is in contact with the lower interface of the following bubble, namely, the contact area is located on the backside of the following bubble, rather than its peripheral part. Consequently, we consider that the coalescence process at Bo = 50 and Ar = 300 is more appropriately described by using the word “hug.”Figure 7
Coalescence process between two bubbles at Bo = 200, Ar = 200,ρr = 1000, μr = 100 and S/R = 2.36.


## 4.1.1. At Low Ar

The rising behaviors (mainly, the coalescence process) of bubbles in quiescent ambient liquid are shown in Figure3 at Bo = 50 and Ar = 1. We know that, at Bo = 50 and Ar = 1, the viscous force and interface tension both play important roles in the rising behaviors of bubbles. Therefore, the inertial force is relatively negligible, and the drag force induced by the relatively high-viscosity causes the bubbles to rise slowly. The interface tension constrains the bubble deformation. As shown in the first frame of Figure 3, the generated jet is rather weak, and only slight deformation is observed. Moreover, from the streamline distribution, it is found that the vortex cores are relatively far from the two bubbles, which causes the liquid circulations to be weak as well. Consequently, the upper interface of the following bubble alters slightly, and the lower interface of the following bubble is slightly elongated. As time increases, because of the shielding effect of the wake flow behind the leading bubble [31, 32], the upper interface of the following bubble is flattened by the dimple formed in the leading bubble. During the coalescence process, the following bubble remains elongated, and while the lower interface of the following bubble becomes slimmer; a smooth tail even appears at t∗=7.2.Figure 3
Coalescence process between two bubbles and corresponding flow field at Bo = 50, Ar = 1,ρr = 1000, μr = 100 and S/R = 2.36.On the other hand, the deformation of the leading bubble shows a reverse evolution during the rising. From left to right of Figure3, the hat shape of the leading bubble becomes more oblate. During the coalescence process, the leading bubble shows a slight deceleration before the coalescence. Consequently, the following bubble continuously approaches the leading bubble. Finally, the coalescence between two bubbles occurs at tc = 7.541, which represents the characteristic coalescence time. As shown in the third frame of Figure 3, the coalescence starts from a single contact point of the two bubbles, which is located almost at the symmetry axis of the two bubbles.Figure4 shows the rising behaviors of two bubbles at Bo = 50 and Ar = 10. As Ar increases to 10, the liquid viscosity decreases significantly. Compared to the case of Ar = 1, the influence of inertia increases a little based on the reduction of drag. Correspondingly, the generated jet becomes stronger, and the deformation of the following bubble increases slightly. Because the vortex cores are relatively near the two bubbles, the liquid circulations become stronger as well. Moreover, because the wake speed becomes faster, the shielding effect of the wake flow behind the leading bubble is strengthened. As shown in the second frame of Figure 4, the deformation of the following bubble becomes more significant, while the lower interface of the following bubble flattens rapidly. On the other hand, the leading bubble shows a similar deformation to the case of Ar = 1. As time increases, the following bubble approaches and contacts the leading bubble. The coalescence occurs at tc = 1.973, which is much earlier than the case of Ar = 1. After bubble coalescence occurs, the liquid circulations accelerate the lower face of the merged bubble, and a larger hat shaped bubble is formed at t∗=2.47.Figure 4
Coalescence process between two bubbles and corresponding flow field at Bo = 50, Ar = 10,ρr = 1000, μr = 100 and S/R = 2.36.It is worth mentioning that the detailed coalescence process of Ar = 10 is different from the case of Ar = 1. Although the coalescence region is located at the central part of both bubbles in the two cases, the contact area is different. In the case of Ar = 1, because of the weak liquid circulations, the impact between the two bubbles is very slight, and the contact area is small enough to be considered a single point. The liquid circulations become stronger in the case of Ar = 10. Under such conditions, there is a slightly stronger impact between the two bubbles during the coalescence process. Consequently, the contact area is enlarged, and the coalescence no longer starts from a contact point but from several points at the same time. It is necessary to note that the contact points are still located at the central part of the two bubbles.


## 4.1.2. At Medium Ar

As Ar increases to 80, the rising behaviors of the two bubbles are shown in Figure5. According to References [3, 33], at Bo = 50 and Ar = 10, the inertial force and interface tension are of the same order. Namely, they both play significant roles in bubble dynamics. Because liquid viscosity decreases, the drag force has less influence on the rising behaviors of bubbles, and the generated jets are rather strong. Therefore, both the interfaces of the two bubbles are easier to deform. As shown in the first frame of Figure 5, an obvious deformation of the following bubble is observed. The lower interface of the following bubble is pushed upward, and an upward dimple is formed. Because the distance between the vortex cores to the lower interface of the following bubble is significantly smaller than the upper interface, the lower interface receives shearing forces induced by the strengthened liquid circulations. Consequently, the lower interface ascends much more rapidly than the upper interface. Then, the following bubble decelerates suddenly due to the obvious oblate deformation. On the other hand, the leading bubble shows a similar deformation as the smaller Ar cases. Namely, the two bubbles both experience an oblate deformation situation, but the deformation degree of the leading bubble is always greater than that of the following bubble. As shown in the second frame of Figure 5, the topology of the two bubbles seems to be a “hat” being put on a “head,” where the leading bubble is the “hat,” and the following bubble is the “head.” It is necessary to note that, compared to the above two smaller Ar cases, the degree of the “putting on” between the two bubbles is strengthened by the decreased liquid viscosity (increased Ar).Figure 5
Coalescence process between two bubbles and corresponding flow field at Bo = 50, Ar = 80,ρr = 1000, μr = 100 and S/R = 2.36.As time increases, under the continuous action of the two vortex rings (liquid circulations), the oblate deformation of the two bubbles becomes more significant. That is, the dimple of the lower interface of the following bubble gradually expands and squeezes most of the bubble to both sides. Then, with the horizontal development of two bubbles, the probability of potential collision on the central part of the two bubbles is decreased. Moreover, as shown in the third frame of Figure5, as the vortex cores are gradually closed to the edge of the bubbles, the shearing effect increases the probability of bubble coalescence in peripherals. Consequently, the coalescence between the two bubbles starts from the peripheral part instead of the central part. Interestingly, the coalescence time is tc = 2.097, which is slightly larger than the case of Ar = 10. As shown in the fourth frame of Figure 5, after the contact area first appears in the peripheral part of the two bubbles, the coalescence subsequently expands to the central part.


## 4.1.3. At High Ar

Figure6 gives the results of the rising behaviors of the bubbles at Bo = 50 and Ar = 300. Under such a condition, with the liquid viscosity reduced significantly, the bubble dynamics are almost unaffected by the drag force. Instead, the bubbles mostly face inertia and gravity. Because the scale of the interface tension is much smaller than the inertia, the interface tension plays a relatively unimportant role in this case. Therefore, the liquid circulations and the generated jets become quite strong. As shown in the first frame of Figure 6, during the approach of the two bubbles, because of the strengthened liquid jet behind the leading bubble, the following bubble significantly deforms and elongates in the vertical direction. Meanwhile, the lower interface of the following bubble is driven by the liquid jet, and an upward dimple is formed. At this time, the shape of the lower interface resembles a skirt. Compared to the first frame of Figure 5, compared to the previous case (Ar = 80), the skirt topology is more significant. Subsequently, the following bubble rising behaviors are pushed by the continuous action of the two vortex rings. As shown in the second frame of Figure 6, most of the gas inside the following bubble moves upward rapidly. Correspondingly, the gas in the tail of the following bubble decreases, and the “neck” of the skirt becomes thinner. Ultimately, because the pressure induced by the ambient liquid flow is rather strong, the interface tension is unable to maintain the existing equilibrium. As shown in the third frame of Figure 6 (t∗ = 1.78), two satellite bubbles are shed from the tail of the following bubble. The phenomenon can be captured for a single bubble under the condition of a rather large Eo [3, 33].Figure 6
Coalescence process between two bubbles and corresponding flow field at Bo = 50, Ar = 300,ρr = 1000, μr = 100 and S/R = 2.36.On the other hand, for the leading bubble, with the continuous action of the two strengthened vortex rings, its axial symmetry thickness decreases continuously until central breakup occurs. As shown in the second frame of Figure6, the leading bubble is divided into two parts. It is necessary to note that the thinning and central breakup of the leading bubble is promoted by the inertia effect induced by the rather fast approach of the following bubble. After the breakup of the leading bubble, the following bubble continues to rise through the hollow space of the leading bubble. Meanwhile, the two parts of the leading bubble spread out on both sides. Then, as shown in the third frame of Figure 6, similar to the previous case (Ar = 80), coalescence occurs in the peripheral part of the two bubbles. However, we consider that the coalescence processes are different between the two cases (Ar = 80 and Ar = 300). The following two aspects can explain the difference. First, the leading bubble is divided into two parts in the case of Ar = 300, which does not occur in the case of Ar = 80. The rising behaviors in the two cases are different. Second, in the case of Ar = 300, coalescence occurs after the breakup of the leading bubble. Although the coalescence area is located at the peripheral part of two bubbles, the coalescence process looks more like the leading bubble’s two arms, giving the following bubble a “hug.” To better interpret this type of coalescence, Figure 7 gives detailed results of the coalescence process at Bo = 200 and Ar = 200. As shown in the third frame of Figure 7, coalescence occurs when the inner interface of the leading bubble is in contact with the lower interface of the following bubble, namely, the contact area is located on the backside of the following bubble, rather than its peripheral part. Consequently, we consider that the coalescence process at Bo = 50 and Ar = 300 is more appropriately described by using the word “hug.”Figure 7
Coalescence process between two bubbles at Bo = 200, Ar = 200,ρr = 1000, μr = 100 and S/R = 2.36.


## 4.2. Ar-Bo Map for the Coalescence Regime and Critical Ar

In the previous section, we established that distinct rising behaviors and coalescence processes occurred under varying Ar and fixed Bo values. To fully understand the relationship between the coalescence regime and the condition parameters, a full picture of the effect of Ar-Bo on the coalescence regime should be plotted.Above all, we classify the coalescence regime into three categories: Central, Edge, and Hug regimes based on the bubble shape topologies. The Central regime is the most common regime in our simulations. In this regime, the coalescence starts from single or several contact points at the central part of the two bubbles (e.g., the third frames of Figures3 and 4). The Edge regime refers to the coalescence that starts from the peripheral parts of the two bubbles (e.g., the third frames of Figure 5). In the Hug regime, the leading bubble breaks into two parts and gives the following bubble a “hug” before the coalescence (e.g., the third frames of Figures 6 and 7). Moreover, the Central coalescence regime was subdivided into two regimes, that is, Central I and Central II. In the Central I regime, the coalescence starts from a single point, which is located at the symmetry axis of the two bubbles. While in the Central II regime, the coalescence starts from several contact points, which are located in the central part near the symmetry axis.Correspondingly, an Ar-Bo map of coalescence regimes is plotted, as shown in Figure8. Obviously, the Central I regime generally occurs at low Ar or low Bo, and the Hug regime always occurs under high-condition parameters (Bo ≥ 50 and Ar ≥ 150). The Central II regime occurs under most of the condition parameters. Interestingly, when Bo is rather small (Bo ≤ 10), even if Ar is as large as 500, the Hug regime is still not observed. Accordingly, we consider that the Edge regime may not translate to the Hug regime even if the liquid viscosity is small enough to be ignored. Moreover, the Edge regime only exists under certain condition parameters (Bo ≤ 75 and Ar ≥ 60) in our simulations. This means that the interface tension and the liquid viscosity are both essential for the Edge regime. However, when the interface tension is rather small (Bo ≥ 100), the Edge regime might not always occur.Figure 8
Ar-Bo map of coalescence regimes plotted as Bo versus Ar.To distinguish the four regimes, three critical Ar values (namely, Ar1, Ar2, and Ar3) are defined as the boundaries among them. As shown in Figure 8, several relationships between the critical Ar and Bo are obtained. First, with the increase of Bo, Ar1 decreases continuously and Ar2 decreases first and then increases, which is different from the results under the axisymmetric computational domain in Reference [21]. The changes in Ar1 or Ar2 suggest that the reduction in interface tension promotes the gas in bubbles to transfer faster in the vertical direction, which makes the contact area between the two bubbles larger. Second, because the leading bubble is more likely to break up into two parts under such conditions, Ar3 decreases with increasing Bo. We consider that the Hug regime always occurs under relatively low-interface tension.


## 4.3. Effect of Ar on Coalescence Time at Different Bo Values

Interface tension is an important factor in the coalescence time of two rising bubbles [34–36]. However, in the two-dimensional domain, the influence of liquid viscosity on the coalescence time under a single interface tension condition has not yet been systematically studied. To obtain a comprehensive understanding, we investigated the changes in coalescence time with varying Ar at several Bo values (Bo = 5, 10, 25, 50, 100, 200) in this section, where Ar corresponds to the liquid viscosity and several Bo values corresponds to different degrees of interface tension.Figure9 gives the results of coalescence time tc versus Ar at different Bo values. In general, with the increase in Ar, tc decreases basically under all Bo values, suggesting that irrespective of Bo, Ar shows a significant effect on tc. With the increase in Bo, tc shows a certain extent of decline, suggesting that the decreasing interface tension promotes coalescence. When the liquid viscosity is rather high (Ar ≤ 3), bubbles coalescence occurs much later. As Ar increases (Ar ≤ 20), according to the comparative interpretation of Figures 3 and 4, the coalescence is considerably promoted by the strengthened inertia effect. However, when Ar ≥ 30, tc shows little change and basically appears to stabilize. According to the interpretation of Figure 5, this implies that the influence of liquid viscosity on tc can be negligible under such conditions. It is necessary to note that, under Bo ≥ 100 and Ar ≥ 125, tc suddenly increases slightly, which is caused by the transition of the Edge to the Hug coalescence regime. According to the interpretation of Figures 3, 5, and 6, compared to the Edge regime, tc in the Hug regime is always considerably larger.Figure 9
Effect of Ar on coalescence time at different Bo.


## 4.4. Effect of Ar and Bo on Relative Rising Velocity

The influence of interface tension (Bo) and liquid viscosity (Ar) on bubble relative rising velocity Ur are illustrated in Figures10 and 11, respectively. In Figure 10, the effect of Ar on Ur is investigated by using three representative Ar values at Bo = 50. In Figure 11, the effect of Bo on Ur is investigated by using three representative Bo values at Ar = 25. In these figures, the ends of the curves all correspond to the beginning of bubble coalescence. The relative rising velocity is specifically analyzed from two views, namely, versus time history (Figures 10(a) and 11(a)) and dependence on instantaneous distance (S/R) between the two bubbles (Figures 10(b) and 11(b)). In general, the time evolution of Ur depends significantly on Ar, Bo, and the instantaneous distance, S/R, between the two bubbles.Figure 10
The effect of Ar on the instantaneous relative velocity at Bo = 50. (a) Relative velocity versus time; (b) relative velocity versus bubble center distance.
(a)(b)Figure 11
The effect of Bo on the instantaneous relative velocity at Ar = 25. (a) Relative velocity versus time; (b) relative velocity versus bubble center distance.
(a)(b)The relationship among the relative velocity Ur, Ar, and instantaneous distance,S/R, can be included as follows: first, as time increases and S/R decreases, Ur first increases and then decreases, which is a result of gas transfers inside the leading and following bubbles. Second, the slopes of the Ur curves are enlarged by increasing Ar, which indicates that low-liquid viscosity promotes gas transfer, increasing the relative velocity. It is necessary to note that, the relationships among the relative velocity Ur, Bo, and instantaneous distance, S/R, in Figure 11 are similar to the results in Figure 10.


## 5. Concluding Remarks

The detailed coalescence and interaction of two in-line, two-dimensional bubbles rising in a two-dimensional quiescent ambient liquid have been studied by using a numerical technique, which combines the conservative phase-field lattice Boltzmann equation (LBE) and the pressure-evolution LBM with a multiple-relaxation-time (MRT) collision operator. Bubble shape evolution, flow field, coalescence time, and relative rising velocity have been analyzed in detail. In particular, we have studied the effect of liquid viscosity and interface tension on the coalescence and interaction of the two bubbles by varying Ar (1 ≤ Ar ≤ 300) for several Bo values (5 ≤ Bo ≤ 200).Four distinct coalescence regimes, that is, Central I, Central II, Edge, and Hug, were revealed by the numerical simulations under the conditions of Ar increasing while other parameters remain constant. The coalescence regimes are the results of a combination of vertical elongation of the following bubble and horizontal expansion of the leading bubble. In particular, with increasing Ar, the bubble shape showed a further complexity, and the possibility for bubble coalescence in the central part was obviously reduced.The relationships between coalescence dynamics, Ar, and Bo are illustrated. Three critical Ar values (namely, Ar1, Ar2, Ar3), which distinguish the boundaries between different regimes, are marked clearly. The changes in the three critical Ar all dependent on Bo. In particular, the Edge regime only exists under certain conditions, and the Hug regime is absent for low Bo and Ar. The coalescence time significantly depends on both Ar and Bo, and the coalescence is promoted by the decreasing interface tension.The motions of bubbles under different condition parameters were studied. With time and a decrease in the instantaneous distance between the two bubble centers, the relative velocity curves first increase and then decrease, and the slope of the Ur curves is enlarged by increasing Ar or increasing Bo.In general, the interaction between the two bubbles is determined by two vortex rings around the two bubbles, which is strengthened by the reduction of liquid viscosity. Coalescence is an outcome of the interaction between the two rising bubbles, and the effects of Ar and Bo on coalescence are revealed as follows. First, in the range of 1 ≤ Ar < Ar1, the increase in interface tension results in very late coalescence. Second, the interaction between the two bubbles is mainly dominated by the interface tension and shearing forces induced by the two vortex rings (liquid circulations) around the two bubbles. The bubbles collide and coalesce in the central part in most cases. Third, under Ar > Ar3, the decrease in Bo postpones the occurrence of bubble coalescence. Under low Bo conditions, the leading bubble shows drastic deformation and breaks into two parts. The coalescence of bubbles is delayed by the transfer of gas in the leading bubble.In the current work, the initial distance between two bubbles was restricted to a relatively constant small value. We consider that a larger initial distance may cause more drastic interactions. Thus, the effect of the initial distance on the coalescence and interaction between the two bubbles is the next step to be explored. Moreover, the initial setup of the two bubbles exhibits mirror symmetry with respect to the central axis, how the features change if this assumption is lifted? Obviously, the bubble rising dynamics in non-Newtonian liquids is a problem worthy of investigation. The above scientific issues have a high probability of becoming our next research projects.



---

*Source: 4090324-2022-06-25.xml*