12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

INLINE COALESCENCE OF TWO BUBBLES ASCENDING IN A
SHEAR-THINNING VISCOELASTIC LIQUID

Wenjun Yuan
School of Chemical Engineering and Technology
Xi’an Jiaotong University
Xi’an 710049, China
wenjun.yuan@xjtu.edu.cn

ABSTRACT

mum coalescence height of bubbles (Jiang et al., 2011; Sunet al., 2017). Especially, Lin & Lin (2003) studied two inline
oblate-cusped bubbles rising in a viscoelastic fluid using particle image analyzer, and noted that two approaching bubbles
result in a circulating flow back to the trailing bubble. Franket al. (2012) indicated that a viscoelasticity-dependent minimum height is required to induce the bubble coalescence. It
was observed that the coalescence for a large leading bubble
and a small trailing bubble coalescence is drastically promoted
in comparison with than that of two-equal bubble coalescence(Lin & Lin, 2009). According to the recent theoretical workof Zilonova et al. (2019), the elasticity can resist the motionof bubbles in space, and the relaxation time of the viscoelastic
medium significantly affects the bubbles’ translational movement. According to our recent work of Yuan et al. (2021a), thefluid elasticity can resist the motion of bubbles in space, andthe relaxation time of the viscoelastic medium significantly affects the bubbles’ translational movement. Unfortunately, due
to experimental limitations, the role of elasticity has not yet
been clearly clarified in the hydrodynamic interactions of abubble pair ascending in viscoelastic liquids.

Direct numerical simulations (DNS) are performed to investigate the inline rise of a pair of three-dimensional (3D)
air bubbles in a viscoelastic liquid using the volume-of-fluid(VOF) approach with an adaptive mesh refinement technique.
The exponential Phan-Thien Tanner (PTT) model is employed
as the non-linear viscoelastic constitutive equation for the liquid. The numerical model has been validated by comparisonswith previously published results, including the terminal velocity jump discontinuity of an isolated bubble rising in a viscoelastic fluid. Focusing on the inline rising bubble pair insuch a viscoelastic medium, we found that for pairs including aleading bubble and a trailing bubble, drafting-kissing scenario
occurs before the bubble-bubble coalescence. The long-range
repulsion and the short-range attraction due to fluid elasticity are critical to the bubble pair interactions. Interestingly,
the squeezing flow near the growing bubble neck seems to delay the coalescence process. The capillary wave propagating
down to the coalesced bubble tip together with the extensional
flow behind the stretched bubble, determines the generation ofsatellite microdroplets along the tail of the coalesced bubble.

Methodology

Introduction

We consider the isothermal, incompressible viscoelastic
two-phase flow. The basic equations for the mass conservationand momentum balance are as follows,

The dynamics of buoyancy-driven bubbles in liquids arenaturally and widely occurring in chemical industries, suchas as bubble column reactors (Besagni et al., 2018), heat exchangers (Singh et al., 2019) and petroleum furnaces (Joshi
et al., 2008), to name a few. Over the last decades, the investigations of bubbly flows (Liao et al., 2015) and rising bubble behaviors (Tripathi et al., 2015; Li & Gupta, 2019) havebeen extensively carried out for Newtonian fluids. However,many solutions in practical applications with the likes of emulsions, slurries and polymer melts, are non-Newtonian fluids(Chhabra, 2010), i.e., their viscosities are dependent on theshear strain applied to the liquid and they exhibit viscoelasticity. Although there are a considerable amount of literature available on the bubbles rising in non-Newtonian liquids(Zenit & Feng, 2018), the hydrodynamic interactions and coalescence of an inline, or tandem bubble pair in viscoelastic
liquids are seldom studied.
Using an arbitrary Lagrangian–Eulerian method, VélezCordero et al. (2011) examined the hydrodynamic interaction
between two bubbles rising for different initial alignments inshear-thinning power-law fluids, and reported that the wakeformed behind a leading bubble attracts the trailing bubble.
The shear-thinning effect was found to decrease the mini-∇ · u = 0,∂uρ+ u · ∇u = −∇p + ∇ · τ + f s − f g .
∂t

(1)



(2)

where u, ρ, t and p are the velocity vector, fluid density, timeand pressure, respectively. τ is the extra-stress tensor. The
source terms f g = gρ j and f s = σ δs κn have been added on
the right hand side of the momentum equation, denoting thegravitational force acting in the negative y direction and thesurface tension force at the gas-liquid interface, respectively.
σ is the interfacial tension coefficient, and j indicates the unit
vector along the y axis direction. δs is the delta distributionfunction, and κ = ∇ · n is the interface curvature, wherein n is
the unit normal at the gas-liquid interface pointing towards thesurrounding liquid phase.
To simulate viscoelastic fluid flows, it is a common approach to split the total extra-stress tensor in asolvent contribution (τ s ) and a polymeric contribution τ p , τ = τ s + τ p .

1

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
where g(Θ) is a model-specific tensorial function dependingon Θ. The tensor Ω is anti-symmetric and generates a purerotation, while B is symmetric and generates a pure extension,
which have been given by Fattal & Kupferman (2004). After
solving Eq. 10, Θ is diagonalized in the form,

In order to have a closed set of equations, we use the constitutive equation presented by Phan-Thien & Tanner (1977) for
the simulation of polymer solutions and melts. Therefore, thesolvent and the polymeric stress tensors can be written as,τ s = 2ηs D(3)
Θ = RΛΘ RT

f tr τ pτ p + λ τfp = 2ηp D(4)

and the conformation tensor is recovered by the inverse relation of Eq. 9,

where ηs and ηp are the solvent viscosity and the polymericviscosity, respectively. λ is the relaxation time. The deformation gradient tensor D is given by,D=i1h∇u + (∇u)T
2

(11)

 
A = exp(Θ) = R exp ΛΘ RT(12)

so that the polymeric extra-stress tensor can be computed from
A (Eq. 7) and used in the momentum equation.

(5)

The Gordon-Schowalter derivative τfp is defined as:

Validation

∂τp
τfp =
+ u · ∇τ p − τ p · ∇u − ∇uT · τ p + ζ τ p · D + D · τ p∂t(6)where ζ ∈ [0,2] is a constant model parameter that controls theamount of movement between the fluid polymeric chains. The
polymeric extra-stress tensor is related with the conformation
tensor A by,τp =ηp(A − I)λ (1 − ζ )Fig. 1 shows the terminal mass centre rise velocity (Ub )of an isolated bubble as a function of its volume (Vb ). Note
that, both the coordinates are logarithmic. The obtained terminal rise velocities at different bubble volumes are quantitatively consistent with the experimental results (Pilz & Brenn,2007) and simulation data (Niethammer et al., 2019). Moreover, Ub slightly increases with Vb before the volume reaches
a critical value Vb,crit . There is an abrupt increase of Ub at thecritical bubble volume, corresponding to 40 mm3 <Vb <50
mm3 in our simulations. The value of Vb,crit and the magnitude of the rise velocity jump are both in good agreementswith the experimental results in Pilz & Brenn (2007). On theother hand, the velocity vector fields for subcritical (33 mm3 )and supercritical (51 mm3 ) volume bubbles are shown in theinset diagrams. The negative wake behind the rising bubble
with a cusped tail is visible for the supercritical state, which isalso in accordance with the observations in the literature (Hassager, 1979). Hence, our 3D simulations have quantitatively
predicted all the salient features for an isolated bubble risingin a viscoelastic fluid correctly (Yuan et al., 2021b), paving theroad for the present investigations on the hydrodynamic interaction and coalescence of two bubbles.
We define the above numerical simulations as Set A,
which indicates that there exists a critical bubble volume at(7)


For the function f tr τ p , we use its exponential form
given by Phan-Thien (1978),f tr τ p







ελ
= exp
tr τ pηp(8)

with ε representing a parameter related to the elongational behaviour of the model.
Following our previous work Yuan et al. (2021b, 2020),
we adopt the finite volume adaptive solver originally developed by Pimenta & Alves (2018). The volume-of-fluid (VOF)
approach was used to represent the interface between the twophases Comminal et al. (2018). In particular, the simulations
were stabilized by reformulating the viscoelastic constitutive
equations with the log-conformation representation Fattal &Kupferman (2004) (LCR) and with both-sides diffusion Xue
et al. (1995) (BSD) in the momentum equation.
In the log-conformation tensor methodology, a new tensor Θ is defined as the natural logarithm of the conformation
tensor,Θ = ln(A) = R ln(Λ)RT(9)

In Eq. 9, the conformation tensor was diagonalized (A =RΛRT ) because it is positive definite, where R is a matrix containing in its columns the eigenvectors of A and Λ is a matrix
whose diagonal elements are the respective eigenvalues resulting from the decomposition of A. Eq. (4) written in terms ofΘ becomes,∂Θ1g(Θ) =+ u · ∇Θ − (ΩΘ − ΩΘ) − 2Bλ∂tFigure 1: Terminal rise velocity as a function of the bub-ble volume for an isolated bubble rising in a viscoelastic
fluid. Taken from (Yuan et al., 2021a).

(10)

2

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
ing the bubbles ascending, the leading one and trailing one getcloser and closer to each other due to the buoyancy-driven upward flow. However, in contrast with the observation in (Yuanet al., 2021a), the trailing bubble approaches the leading one
and then coalesces with it. This attraction seems to be causedby the strong buoyancy force of the trailing bubble, whichis acted on to generate the rapid initial acceleration (Yuan &Prosperetti, 1994).
Fig. 4 compares the detailed distributions of the polymeric stress components τxx and τyy in a vertical cutting plane
during the bubble-bubble coalescence process at different time
instants. For the short-range attraction in Figs. 4 before coalescence, τxx always has a positive large value between thetwo bubbles, while τyy is still modest. It is interesting to observe that as the bubbles are getting closer, the leading bubble
bottom takes a remarkable bulge at t= 0.096 s rather than the
flattened shape at the bubble lower hemisphere in Fig. 3 at t=
0.090 s. Moreover, it is obvious that at the lower raw of Fig.
4, the capillary wave propagates down from the growing coalesced neck to the bottom of the coalesced bubble. The comparison of these snapshots shows that the horizontal normal
stress component τxx in the neighbourhood of the coalescedneck gradually fades away, and both the normal stress components τxx and τyy come to emerge near the coalesced bubble
downstream pole. Recalling the results for an isolated bubble
rising in viscoelastic fluids (Yuan et al., 2021b), this indicatesthat polymeric stresses are squeezed in the small region near
the bubble trailing end, which will further force the ejection of
the satellite microdroplets. This is also the mechanism leading
to the formation of a cusp shape with a tip at the supercriticalbubble downstream pole, whereas the coalesced bubble has alarger volume, its tail is more likely to breakup with the propagating capillary wave.

Figure 2: The behaviors of two bubbles rising inline withdifferent volume configurations for initial separated distance h0 = 6 mm. Note that, circles and crosses indicate
the occurrence of coalescence and not, respectively.
Table 1: Rheological parameters for the polymer

solution used in the simulations.
β

λ [s]

γ[N/m]

ηs [Ns/m2 ]η p [Ns/m2 ]

Set A

0.02

0.207

0.076

0.031

1.511

Set B

0.02

0.207

0.038

0.002

0.098

which the bubble terminal velocity sharply increases in viscoelastic fluids. For such a liquid, different bubble volume
configurations have been investigated for the leading and trailing bubbles with Vl and Vt = 14, 33, 51 and 65 mm3 for initial
separated distance h0 = 6 mm. The results are given in Fig.
2, where the circles denote the cases in which the two inline
bubbles can successfully coalesce during the rise process. We
can see that the requirement of “a successful coalescence case”
is that the trailing bubble is supercritical. However, it is noteworthy that the terminal rise velocity invariably increases withthe bubble volume (Yuan et al., 2021b). When the different
size bubbles exchange the initial positions, it is believed thatthe local flow patterns in the neighbourhood of rising bubbles
will make adjustment. For the different initial separated distance, the upward acceleration caused by the buoyancy may
also affect the rise behaviors of the two co-axial bubbles.
In addition, to contrast the volume effects on the bubble dynamics in wide flow regimes, we also perform threedimensional numerical simulations for different fluid physicalproperties. As the rise of an air bubble in viscoelastic fluids isa multi-physical phenomenon, we fix the polymer relaxation
parameter λ and the viscosity ratio β in the simulations. The
exponential PTT model for viscoelastic polymer solutions continues to be adapted and the extensibility parameter ε and theslip parameter ξ are set to constant 0.05 and 0, respectively.
Another different set is obtained by changing the interfacial
tension γ and the viscoelastic liquid viscosity ηB , as shown in
Table 1. In the following, we present the inline rise of two
bubbles for this Set B.

In Fig. 5, the vertical velocity fields as well as the distribution of the trace of the natural logarithm of the conformation tensor Θ are presented in a vertical cutting plane throughthe bubble center. To reduce image duplication, the snapshots
are visualized for the leading bubble and the trailing bubble
during the coalescence process, corresponding to the bubble
topology history in Fig. 4. From Figs. 5 , it is anticipatedthat the vertical fluid velocities around the coalesced bubble
lower and upper hemispheres are enhanced by the coalescence,
hence the coalesced bubble rises rapidly. For the instant before the coalescence in Fig. 5, the magnitudes of negative andpositive vertical velocity contours are increased at the leading bubble bottom and at the trailing bubble top, respectively,
providing evidence that the two bubbles attract each other andthen aggregate. This is at variance with the inline coalescence
between two bubbles in Newtonian flows (Chen et al., 2011;
Zhang et al., 2018), that the bubble-bubble interaction is only
determined by the influence of the leading bubble wake on the
trailing bubble. Regardless of the downward propagation of
the capillary wave in (Yuan et al., 2021a), the liquid circulation inside the growing neck is strong, reflecting the radicallychanging of the velocity at the top and the bottom of the coalesced bubble.
As shown in the lower raw of Figs. 5, the large value
regions of tr(Θ) exist in the narrow regions between the two
bubbles before coalescence. This indicates that the polymer
molecules tend to deform greatly as the two bubbles approach.
The squeezed film of the viscoelastic liquid is significantly
thicker than that in (Yuan et al., 2021a), so that the short-range
attraction is strong due to elasticity, related to the apparent
bulges at the bottom of the leading bubble and at the top of thetrailing bubble at t= 0.096 s. Owing to the fact that the large

ResultsFigs. 3 present the complex shapes of the inline rise of a
trailing bubble and a leading bubble (Vt = Vl = 51 mm3 ). Dur-3

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

Figure 3: Complex shapes variation with time for the inline rise of a leading bubble and a trailing bubble in a viscoelastic
fluid for Set B.

tr(Θ) occupies the region near the coalesced bubble neck, thepolymer deformation seems to hinder the film drainage during the coalescence process (Dekee et al., 1986). Conceivably,
the elasticity in the liquid delays the coalescence, as a resultof the squeezing flow in the film. The coalescence of the inline bubbles in this work is significantly different from the socalled drafting-kissing-tumbling scenario that occurs for bubble pairs in a Newtonian fluid and in a shear-thinning inelastic
fluid (Chen et al., 2011; Vélez-Cordero et al., 2011).

ulation of the planar extrudate swell of pseudoplastic andviscoelastic fluids with the streamfunction and the VOF
methods. J. Non-Newton Fluid Mech. 252, 1–18.
Dekee, D, Carreau, PJ & Mordarski, J 1986 Bubble velocity and coalescence in viscoelastic liquids. Chem. Eng. Sci.
41 (9), 2273–2283.
Fattal, Raanan & Kupferman, Raz 2004 Constitutive laws forthe matrix-logarithm of the conformation tensor. J. NonNewton Fluid Mech. 123 (2-3), 281–285.
Frank, Xavier, Charpentier, Jean-Claude, Ma, Youguang, Midoux, Noël & Li, Huai Z 2012 A multiscale approach formodeling bubbles rising in non-Newtonian fluids. Ind. Eng.
Chem. Res. 51 (4), 2084–2093.
Hassager, O. 1979 Negative wake behind bubbles in nonNewtonian liquids. Nature 279 (5712), 402–403.
Jiang, Shaokun, Ma, Youguang, Fan, Wenyuan, Yang, Ke &Li, Huaizhi 2011 Chaotic behavior of in-line bubbles risingwith coalescences in non-Newtonian fluids: A multiscale
analysis. Korean J. Chem. Eng. 28 (1), 56–63.
Joshi, Jyeshtharaj B, Pandit, Aniruddha B, Kataria, Kamal L,
Kulkarni, Rohit P, Sawarkar, Ashish N, Tandon, Deepak,
Ram, Yad & Kumar, Man Mohan 2008 Petroleum residue
upgradation via visbreaking: A review. Ind. Eng. Chem.
Res. 47 (23), 8960–8988.
Li, Weihua & Gupta, Nivedita R 2019 Buoyancy-driven motion of bubbles in the presence of soluble surfactants in a
Newtonian fluid. Ind. Eng. Chem. Res. 58 (18), 7640–7649.
Liao, Yixiang, Rzehak, Roland, Lucas, Dirk & Krepper, Eckhard 2015 Baseline closure model for dispersed bubbly
flow: Bubble coalescence and breakup. Chem. Eng. Sci.
122, 336–349.
Lin, Tsao-Jen & Lin, Gen-Ming 2003 The mechanisms of bubble coalescence in a non-Newtonian fluid. Can. J. Chem.
81 (3-4), 476–482.
Lin, Tsao-Jen & Lin, Gen-Ming 2009 Mechanisms of in-line
coalescence of two-unequal bubbles in a non-Newtonian
fluid. Chem. Eng. J. 155 (3), 750–756.
Niethammer, Matthias, Brenn, Günter, Marschall, Holger &Bothe, Dieter 2019 An extended volume of fluid method
and its application to single bubbles rising in a viscoelastic
liquid. J. Comput. Phys. 387, 326–355.
Phan-Thien, Nhan 1978 A nonlinear network viscoelastic
model. J. Rheol. 22 (3), 259–283.
Phan-Thien, Nhan & Tanner, Roger I 1977 A new constitutive
equation derived from network theory. J. Non-Newton Fluid

Conclusions
To improve the understanding of bubble dynamics in
non-Newtonian media, we have carried out direct numericalsimulations (DNS) to explore the ascent of a pair of threedimensional bubbles in a viscoelastic liquid. The adaptive
mesh refinement approach is used to refine the mesh near theregion of interest, and the exponential PTT model is employed
to describe the shear-thinning viscoelasticity. On the basis of
the good agreement between the present simulation results andexperimental data by Pilz & Brenn (2007), we focus on the inline coalescence of a tailing bubble and a leading bubble. By
probing the polymeric stresses and deformation, our results
demonstrated that the short-range bubble-bubble attraction isstrong due to the presence of elasticity, while the squeezing
flow hinders the film drainage and delays the coalescence. In
addition, the coalescence of large bubbles is more prevailingin highly elastic flows. The capillary wave propagating fromthe growing neck to the apex of the coalesced bubble and theextensional flow along the coalesced bubble tail, determine thegeneration of satellite microdroplets.

REFERENCES
Besagni, Giorgio, Inzoli, Fabio & Ziegenhein, Thomas 2018
Two-phase bubble columns: A comprehensive review.
ChemEngineering 2 (2), 13.
Chen, RH, Tian, WX, Su, GH, Qiu, SZ, Ishiwatari, Yuki &Oka, Yoshiaki 2011 Numerical investigation on coalescenceof bubble pairs rising in a stagnant liquid. Chem. Eng. Sci.
66 (21), 5055–5063.
Chhabra, Rajendra P 2010 Non-Newtonian fluids: an introduction. In Rheology of complex fluids, pp. 3–34. Springer.
Comminal, Raphaël, Pimenta, Francisco, Hattel, Jesper H,
Alves, Manuel A & Spangenberg, Jon 2018 Numerical sim-

4

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

Figure 4: Distributions of the polymeric stress tensor components τxx (left) and τyy (right) in a vertical cutting plane throughthe centre of the rising bubbles during the inline coalescence process.

Mech. 2 (4), 353–365.
Pilz, C. & Brenn, G. 2007 On the critical bubble volume at
the rise velocity jump discontinuity in viscoelastic liquids.
J. Non-Newton Fluid Mech. 145 (2-3), 124–138.
Pimenta,F&
Alves,
MA2018
Rheotool.
https://github.com/fppimenta/rheoTool.
Singh, Jogender, Montesinos-Castellanos, Alejandro &Nigam, Krishna DP 2019 Process intensification for compact and micro heat exchangers through innovative technologies: A review. Ind. Eng. Chem. Res. 58 (31), 13819–
13847.
Sun, Wanpeng, Zhu, Chunying, Fu, Taotao, Yang, Hui, Ma,
Youguang & Li, Huaizhi 2017 The minimum in-line coalescence height of bubbles in non-Newtonian fluid. Int. J.
Multiph. Flow 92, 161–170.
Tripathi, Manoj Kumar, Sahu, Kirti Chandra & Govindarajan,
Rama 2015 Dynamics of an initially spherical bubble risingin quiescent liquid. Nat. Commun. 6 (1), 6268.
Vélez-Cordero, J Rodrigo, Sámano, Diego, Yue, Pengtao,
Feng, James J & Zenit, Roberto 2011 Hydrodynamic interaction between a pair of bubbles ascending in shearthinning inelastic fluids. J. Non-Newton Fluid Mech. 166 (12), 118–132.
Xue, S-C, Phan-Thien, N & Tanner, RI 1995 Numerical studyof secondary flows of viscoelastic fluid in straight pipes
by an implicit finite volume method. J. Non-Newton Fluid
Mech. 59 (2-3), 191–213.
Yuan, H & Prosperetti, A 1994 On the in-line motion of two
spherical bubbles in a viscous fluid. J. Fluid Mech. 278,
325–349.
Yuan, Wenjun, Zhang, Mengqi, Khoo, Boo Cheong & PhanThien, Nhan 2020 Dynamics and deformation of a threedimensional bubble rising in viscoelastic fluids. J. NonNewton Fluid Mech. 285, 104408.
Yuan, Wenjun, Zhang, Mengqi, Khoo, Boo Cheong & PhanThien, Nhan 2021a Hydrodynamic interaction and coalescence of two inline bubbles rising in a viscoelastic liquid.
Phys. Fluids 33 (8), 083102.
Yuan, Wenjun, Zhang, Mengqi, Khoo, Boo Cheong & PhanThien, Nhan 2021b On peculiar behaviours at critical volumes of a three-dimensional bubble rising in viscoelastic
fluids. J. Non-Newton Fluid Mech. 293, 104568.
Zenit, R. & Feng, J.J. 2018 Hydrodynamic interactions among
bubbles, drops, and particles in non-Newtonian liquids.
Annu. Rev. Fluid Mech. 50, 505–534.
Zhang, Yang, Chen, Ke, You, Yunxiang & Ren, Wei 2018 Coalescence of two initially spherical bubbles: Dual effect ofliquid viscosity. Int. J. Heat Fluid Flow 72, 61–72.

5

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

Figure 5: Vertical fluid velocity contours and distributions of the trace of the natural logarithm of the conformation tensor
Θ in a vertical cutting plane through the rising bubbles during the inline coalescence.

coelastic drag. Phys. Rev. E 99 (2), 023109.

Zilonova, Ekaterina, Solovchuk, Maxim & Sheu, TWH 2019
Dynamics of bubble-bubble interactions experiencing vis-

6