Contemporary Mathematics
Numerical Simulation of Gas Bubbles Rising in Viscous
Liquids at High Reynolds Number
Jinsong Hua, Ping Lin, and Jan F Stene
Abstract. A robust algorithm for direct numerical simulation of fully threedimensional, incompressible two-phase flow is presented. The method is introduced in the context of gas bubbles rising in viscous liquids, e.g. air bubbles
rising in water. Key strengths of the simulation approach include the ability
to simulate flows in an extended, wide range of Reynolds and Bond numbers
and also the ability to handle large, realistic ratios of the density and viscosity of the fluids. The aptitude of the algorithm is due to the combination of
several powerful components: The interface between the phases is tracked explicitly by an unstructured, adaptive, triangular mesh (front tracking), while
the equations governing incompressible, Newtonian flows are solved using a
finite volume scheme based on the Semi-Implicit Method for Pressure-Linked
Equations (SIMPLE) algorithm. Further, the SIMPLE flow solver is integrated
with PARAMESH: a block-based, adaptive mesh refinement (AMR) tool for
multi-processing, allowing the solution of the governing equations in parallel
on a supercomputer. Finally, the use of a non-inertial, moving reference frame
attached to the rising bubble facilitates long-time simulations without requiring a larger computational domain. The methodology just described has been
applied to simulate 3D gas bubbles rising in viscous liquids. In particular, air
bubbles rising in water were studied, and the numerically predicted terminal
rise velocities agree well with experimental measurements for bubbles with initial diameters ranging from 0.5 mm to 30 mm. Simulation results also reveal
that bubbles rising in viscous liquids at high Reynolds numbers may rise on a
zigzag or/and a spiral path. At last, we tried to use this model to simulate the
interaction between two bubbles in liquid. The simulation results provide us
with more physical insights into the complex bubble rising behavior in viscous
liquids.
2000 Mathematics Subject Classification. Primary 76T10; Secondary 76D05, 76D45.
Key words and phrases. Two-phase flow, bubble dynamics, air-water, path instability, front
tracking method, PARAMESH, bubble-bubble interaction, moving reference frame.
This paper was also presented on the 6th International Conference on Multiphase Flow,
ICMF 2007, Leipzig, Germany, July 9 - 13, 2007.
Authors were supported in part by Singapore Academic Research Funds R-146-000-053-112
and R-146-000-064-112.
c
°0000
(copyright holder)
1
2
JINSONG HUA, PING LIN, AND JAN F STENE
1. Introduction
Multiphase flows occur in a wide range of situations such as many important
biological and industrial processes, and they have been extensively studied by researchers both theoretically and experimentally. One very fundamental example of
multiphase flows is that of an isolated gas bubble rising and deforming in a viscous
liquid. It is clear that a thorough knowledge of this basic system is of major importance to understand more complex flows such as bubbly flows. However, despite
its apparent simplicity, important aspects of the behavior of a single bubble rising
still remain uncertain. This incomplete knowledge is certainly not due to a lack
of research effort. Numerous experimentalists have studied single bubbles rising
in otherwise quiescent liquids ([Cl78],[WG02] and [To02]). One major obstacle
they face is the limitations in measuring techniques. For example, the measuring
techniques that are well-suited for one phase is disturbed by the presence of the
other phase, and it is difficult to get complete information about the flow field and
pressure distribution in the entire domain. It is also non-trivial to retrieve fully
three-dimensional representations of the evolving bubble shape.
In addition to the limitations related to measurements, there have often been
inconsistencies in the physical quantities that are indeed measurable. One important reason for these inconsistencies is probably the effect of different degrees of
purity of the liquids used. The resulting variations in surfactants have then in
certain flow regimes influenced the results (See [Cl78]). Another reason may be
attributed to the way the bubbles were formed since initial shape perturbations
have been reported to greatly affect the evolution of the flow dynamics of a rising
bubble (See [WG02] and[To02]).
In recent years, improved numerical algorithms and more computing power
have also seen direct numerical simulations contributing to a better physical understanding of multiphase flows. The great advantage of numerical simulations
is of course that they avoid some of the aforementioned practical difficulties encountered in experiments: important data such as the flow field and the pressure
is readily available throughout the computational domain; truly pure water without any surfactants can be assumed throughout the simulations; and there are no
uncertainties associated with the initial bubble formation. Thus direct numerical
simulations have great potential as a complement to experiments in the pursuit of
a better understanding of flow problems such as a single bubble rising.
Many numerical studies on bubbles rising in viscous liquids have been reported.
Unfortunately, they are limited to the cases with either simple bubble shapes, e.g.
axisymmetric bubbles and bubbles with prescribed fixed shapes, low density and
viscosity ratios, or lower Reynolds number regimes. Axisymmetric bubble shapes
are assumed in the papers [KB99], [RR00], [Da04] and [BM06]. This restricts
the applicability of the numerical method to a limited flow regime only. Mogin
and Magnaudet ([MM02] and [MM06]) are amongst those who have assumed a
non-deformable, fixed bubble shape to simplify the simulations of a bubble rising in
the high Reynolds number regime. Low values for the density and viscosity ratios,
which are smaller than those of real gas-liquid systems such as the air-water system,
are used in the numerical studies on multiphase flow by Unverdi and Tryggvason
[UT92], Chen et al. [Ch99] and Annaland et al. [An05].
Fully three-dimensional simulations of a single bubble rising with realistic parameters have been attempted by a number of researchers. Sussman et al. [Su99]
GAS BUBBLES RISING IN VISCOUS LIQUIDS
3
introduced an adaptive level set approach and presented some preliminary results
for an inviscid, three-dimensional air bubble rising in water. The same system
was also studied in [Ko02] by applying a three-dimensional volume of fluid (VOF)
method. Initial spherical bubbles with diameters ranging from 0.5 mm to 15 mm
were studied and reasonable results were obtained for their terminal velocities.
However, they reported that bubbles with an initial diameter larger than 15 mm
had a tendency to break up before terminal velocity was reached. Lrstad and Fuchs
[LF04] presented an improved 3D VOF method and applied it successfully to study
ellipsoidal air bubbles in water with equivalent diameters ranging from 1.82 mm
to 6 mm. Yet another VOF approach was proposed in [Oh05] when they studied
three-dimensional air bubbles rising in silicone oils with various Morton numbers.
They reported good agreement of predicted bubble shapes in most flow regimes, but
they were unable to reproduce reasonable shapes in the spherical cap regime with
high Reynolds numbers. Both a 2D VOF method as well as a 3D front tracking
approach were used in [Di05]. The front tracking algorithm reasonably predicted
the rise velocity and aspect ratio of a single air bubble rising in water for diameters
in the range of 1 mm to 7 mm. Further, a 3D VOF method and a 3D front tracking
method were used in [An05] and [An06], respectively. According to their results,
the latter approach appears more robust in the sense that it can simulate bubbles
rising in a far wider flow regime. However, as was the case with [Oh05], no results
were reported for the spherical cap regime characterized by high Reynolds numbers.
This paper is to present a direct numerical simulation algorithm that is capable of simulating the complex bubble rising behavior in viscous liquids when the
Reynolds number is high. The numerical algorithm that allows these predictions
is a state-of-the-art front tracking method. The robustness and versatility of the
method is due to the combination of several powerful components: the moving interface is handled with the front tracking method (cf. [UT92], [Tr01] and [HL07]),
while the governing Navier-Stokes equations are solved with a modified SIMPLE
scheme [Pa80] on an adaptive Cartesian grid using PARAMESH (See [Ma00]). To
allow long-term simulations, the system is solved in a moving reference frame that
moves together with the rising bubble.
In addition, we present numerical results for a single air bubble rising in water
for bubble diameters from 0.5 mm to 30 mm, which is a much wider range than
other simulations reported so far. We predict terminal velocity and bubble shape
in the spherical, wobbling and also the challenging spherical cap regime. Moreover,
our simulations also reveal the path instability when a single bubble rises at high
Reynolds numbers in the wobbling regime. Both zigzag and spiral paths have been
reproduced in numerical simulations under certain conditions. At last, the current
numerical algorithm is applied to investigate the interaction mechanism of two
different-sized bubbles rising in liquid. These simulations provided some physical
insights into the bubble-bubble interaction dynamics.
2. Numerical Method
Governing Equations
We can reasonably assume that the multi-fluid systems studied in this paper are
isothermal systems of two Newtonian, incompressible and immiscible fluids. Using
4
JINSONG HUA, PING LIN, AND JAN F STENE
the ’one-fluid’ approach (cf. [SZ99]), we get one single set of governing equations
for the entire flow domain, namely mass conservation,
(2.1)
∇·u=0
and momentum conservation
Z
∂(ρu)
(2.2)
+∇·ρuu = −∇p+∇·[µ(∇u+∇T u)]+ σf κf nf δ(x−xf )ds+(ρ−ρl )g.
∂t
Γ
It is worth pointing out that the material properties density ρ and viscosity µ will
be discontinuous across the interface. Similarly, there will generally be a jump in
the pressure p across the interface as well. Here x is the point at which the equation
is evaluated and xf is a point on the interface f . Note that the surface tension
term is a singular term that only comes into effect on the interface between the two
fluids.
We non-dimensionalize the equations by introducing the dimensionless characteristic variables
x∗ =
x
u
, u∗ = √
, τ∗ =
D
gD
r
g
ρ
p
µ
g
t, ρ∗ = , p∗ =
, µ∗ = , κ∗ = Dκ, g∗ =
,
D
ρl
ρl gD
µl
kgk
where D is the diameter of a sphere with the same volume as the bubble. Consequently, we may re-express the Navier-Stokes equations as
(2.3)
Z
∂(ρ∗ u∗ )
1
1
+∇·ρ∗ u∗ u∗ = −∇p∗ + ∇·[µ∗ (∇u∗ +∇T u∗ )]+
κ∗f nf δ(x∗ −x∗f )ds+(ρ∗ −1)g∗ .
∂t
Re
Bo Γ
Note that the dimensionless Reynolds and Bond numbers used here are thus defined
as
(2.4)
Re =
ρl g 1/2 D3/2
ρl gD2
and Bo =
.
µl
σ
By studying the non-dimensional formulation, one can observe that the flow may be
entirely characterized by the following four dimensionless parameters: The density
and viscosity ratios of the fluids, the Reynolds number and the Bond number.
Front Tracking Method for Multi-phase Flow
A relatively brief overview of the numerical method that has been used in
the current simulations is presented in this paper. A more detailed description of
the axisymmetric version of this method can be found in [HL07]. The algorithm
basically consists of the following main parts: front tracking of the interface which
is represented by an unstructured triangular mesh; a SIMPLE-based scheme to
solve the flow equations on a Cartesian grid; an adaptive mesh refinement tool,
PARAMESH, used to generate the flow solver grid; and all computations are carried
out in a moving reference frame.
Discontinuities across the Interface
We represent the interface with an unstructured triangular mesh, often referred
to as the front mesh. This mesh is a completely separate data structure from the
three-dimensional Cartesian grid on which the governing flow equations are solved
- often referred to as the background grid. Figure 1 shows a schematic diagram of
GAS BUBBLES RISING IN VISCOUS LIQUIDS
5
Figure 1. Schematic model of the front tracking method for
multi-phase flow simulations.
the front and background mesh used for the multiphase flow simulations. It can
be noted from Equation (2.2) that the density and viscosity distributions on the
background grid are needed to solve the governing flow equations. However, these
fluid properties are discontinuous across the interface. The density and viscosity is
therefore reconstructed on the background grid based on the position of the front
mesh. This is done by introducing an indicator function I(x, t) which has the value
of one when x is in the gas phase and zero in the liquid phase. Letting Γ(t) be the
interface between the phases, can be found by solving
Z
(2.5)
∇2 I(x, t) = ∇ ·
nδ(x − x0 )ds.
Γ(t)
To obtain a smooth indicator function, we replace the delta distribution δ(x) with
the following discrete delta distribution introduced in [Pe77]:
½
Q3
π
(4h)−3 i=1 (1 + cos( 2h
|xi − xi,f |)) if |xi − xi,f | < 2h
(2.6) D(x − xf ) =
,
0
otherwise
where h is the grid size of background mesh, and i represents the three directions
of the coordinate system.
It is also noted from Equation (2.2) that the value of the singular surface tension
term on the background grid is needed. In this study, the surface tension force is
calculated by considering the net force caused by surface tension on each triangular
surface element. The surface tension force is then distributed to the background
grid:
(2.7)
Fσ (x) =
X
D(x − xf )
f
3
X
σ(ti × ni,0 ),
i=1
where t and n are the edge vector and outer normal of the surface element, respectively.
Front Mesh Advection and Adaptation
The triangular front mesh consists of unstructured points often referred to as
front markers. When the flow field is solved on the background mesh, the velocity
on these front markers can be obtained through interpolation from the background
6
JINSONG HUA, PING LIN, AND JAN F STENE
(a)
(b)
Figure 2. Mesh quality variation of the triangular mesh representing the rising bubble surface: (a) without, and (b) with mesh
adaptation.
(a)
(b)
Figure 3. An illustration of the multi-blocks of the background
mesh for the flow field in the vicinity of the bubble surface: (a)
before, and (b) after refinement.
flow field. Subsequently, the interface at time tn is advected to its new position at
time tn+1 based on the velocity field at time tn .
Moving the front markers with the surrounding fluid flow will change the resolution and quality of the front mesh. This change will influence the exchange
of information with the background grid and thus eventually also the simulation
results. The front mesh is therefore examined and dynamically adapted to ensure
a uniform resolution throughout the duration of the simulation – see Figure 2 for
an illustration of the effect of this adaptation. This triangular mesh adaptation
consists of the three basic operations of swapping, splitting, and deleting edges.
Further, it is of critical importance that the background grid has an adequate
resolution. In particular, a fine grid is required near the interface due to the smallscale dynamics in this region, whereas the resolution requirements are more relaxed
away from the interface. It may thus be desirable to have a fine grid near the
interface and a coarser grid away from the bubble surface. To achieve this, the
block-based adaptive mesh refinement tool PARAMESH (See [Ma00]) is adopted
in this study to manage the background mesh. Figure 3 presents an illustration of
the different levels of refinement of the background grid.
GAS BUBBLES RISING IN VISCOUS LIQUIDS
7
Figure 4. Schematic diagram of the moving reference frame attached to the rising bubble.
Obviously there must be a correspondence between the resolution of the front
mesh and the background grid. The resolutions are of the same order, but typically
the front mesh is slightly finer than the background grid.
Flow Solver
The choice of numerical method to solve the governing flow equations can
be made independent of how we deal with the interface. Traditionally, explicit
projection-correction methods have frequently been used in conjunction with front
tracking of the interface (See [UT92]). However, this approach has encountered
problems for the large density ratios typical of many natural and industrial processes.
An alternative solver presented in [HL07] seems to be able to resolve this problem. Their method is basically a modified version of the classical SIMPLE method
(See [Pa80]). The SIMPLE algorithm avoids solving the pressure equation directly,
instead correcting the pressure and velocity iteratively based on the governing equations.
A Moving Reference Frame
Many applications in multiphase flow involve moving fronts, e.g. a gas bubble
rising in a liquid. It is often desirable to study the long-term behavior and the
temporal evolution of the front, in which case it may move a relatively long distance. Thus the computational domain must be correspondingly large. However,
in a three-dimensional setting such a large domain is likely to be computationally
infeasible, and this constraint may therefore prevent the desired long time simulation. To remedy this problem, we are solving the flow equations in a non-inertial
reference frame that moves with the front as shown in Figure 4. Here, the frame
XY is a stationary reference frame and the frame X 0 Y 0 a moving reference frame.
The velocity of the monitoring point P is u(x, t) in the frame XY and u0 (x0 , t) in
the moving frame X 0 Y 0 . The velocity of the moving frame is um (t). The relation
8
JINSONG HUA, PING LIN, AND JAN F STENE
between the velocities is thus
(2.8)
u(x, t) = um (t) + u0 (x0 , t).
Using this non-inertial frame, the front remains more or less fixed in the computational domain as time goes by. Importantly, this allows the duration of the
simulation to be chosen independently of the domain size. Allowing translational,
but not rotational movement of the frame, we then get the following updated flow
equations:
(2.9)
Z
∂(ρu0 )
1
1
+∇·ρu0 u0 +ρam = −∇p+
∇·[µ(∇u0 +∇T u0 )]+
κf nf δ(x−xf )ds+(ρ−1)g.
∂t
Re
Bo Γ
Notice the additional term on the left-hand side, where am denotes the acceleration
of the moving reference frame. This quantity must be predicted at each time step,
and in this study we make this prediction following the approach presented in
[Ru02]. The negative of the frame velocity, −um (t), is then used as the boundary
condition for the velocity field in the moving reference frame.
Solution Procedure
We may now summarize the main steps in advancing the solution from one
time step to the next as follows:
(1) The velocity of the front marker points, u∗f , is calculated through interpolation of the fluid velocity field u∗ .
(2) The front is advected to its new position xn+1
by using the normal inf
n
terface velocity uf found in (1). The front elements are then subject to
examination for adaptation and topological change, and volume conservation is enforced.
(3) The indicator function I n+1 (xn+1
) is computed based on the interface
f
position xn+1
.
Subsequently,
the
distribution of the density ρn+1 , the
f
viscosity µn+1 and the surface tension Fn+1
is updated on the flow solver
σ
grid points.
(4) We find the velocity field un+1 and the pressure pn+1 by solving the
mass continuity and momentum equations using a modified version of the
SIMPLE algorithm. Appropriate boundary conditions are applied.
(5) Repeat steps (1) to (4) to advance the solution to time tn+2 .
3. Results and Discussion
Numerical Model Validation
Despite the tremendous improvement in computing power, accurate numerical
simulations of three-dimensional two-phase flows are still computationally challenging. Hence, the ideal solution to this is to minimize the computational requirements
while still ensuring satisfactory simulation results. Achieving this in the numerical
simulations of a bubble rising in liquid requires answering questions like: How large
must the computational domain be to avoid wall confinement effects? What mesh
resolution is required to capture the relevant physics? For a detailed sensitivity
analysis addressing these questions, the reader is referred to the previous work in
GAS BUBBLES RISING IN VISCOUS LIQUIDS
(a)
9
(b)
Figure 5. Comprison case A: (a) experimentally observed bubble
with experimental conditions Bo = 116, M = 1.31 and Re = 20.4,
and (b) numerically predicted bubble with modeling conditions
and the predicted terminal velocity Bo = 116, Re∗ = 33.02 and
U∗ = 0.59.
(a)
(b)
Figure 6. Comprison case B: (a) experimentally observed bubble
with experimental conditions Bo = 339, M = 43.1 and Re = 18.3,
and (b) numerically predicted bubble with modeling conditions
and the predicted terminal velocity Bo = 339, Re∗ = 30.83 and
U∗ = 0.57.
[HL07]. In the current simulations, the computational domain is cubic with side
lengths equal to eight times the diameter of the rising bubble. Further, the resolution of the background grid in the vicinity of the bubble corresponds to 20 grid cells
per bubble diameter, while the resolution of the front mesh is around twice as fine.
The dimensionless time step is 0.005 throughout the duration of the simulations.
The modeling conditions and results for the two cases studied here are given
below. The conditions are similar in these two cases except when it comes to the
Bond number: Case B has a significantly higher Bond number than case A. This
implies that the surface tension is much lower in Case B than in Case A, and
the deformation of the bubble is indeed much more pronounced in B than in A.
Importantly, Figures 5 and 6 show good agreement between the experimentally
observed (See [BW81]) and numerically predicted bubble shapes for both cases.
Table 1 below compares the experimental Reynolds number Re with the computed Reynolds number ReC based on the numerically predicted terminal velocity.
There is an excellent agreement in both cases with deviations less than 5%.
Air Bubbles Rising in Water
Air bubbles rising in water may be widely observed in many industrial processes.
Examples in chemical engineering include bubble columns, loop reactors, agitated
stirred reactors, flotation, and fermentation reactors. For the design of efficient
two-phase reactors, detailed knowledge of bubble sizes and shapes, slip velocities,
10
JINSONG HUA, PING LIN, AND JAN F STENE
Table 1. Comparison of terminal rise velocities found in experiment and predicted by numerical simulations.
Test cases
A
B
Experiment
Re
Re∗
20.4
33.02
18.3
30.83
Simulation
U∗
ReC = U∗ · Re∗
0.59
19.5
0.57
17.6
Deviation
|Re − ReC |/Re · 100
4.4%
3.8%
Figure 7. Comparison of the numerically predicted and experimentally measured terminal rise.
internal circulation, swarm behavior, bubble induced turbulence and mixing, and
bubble size distribution (including coalescence and breakup) is of fundamental importance. In such industrial applications, bubbles often have non-spherical and
even dynamic shapes as well as asymmetric wake structures. Extensive experimental work has been carried out to study the behavior of an air bubble rising in water
(See [Cl78] and [To02]). Their measurements of the terminal velocity of air bubbles rising in water are presented in Figure 7 as a function of the bubble size. It
is found that the measurements of the bubble terminal velocity vary significantly
(or bifurcation) when the bubble size is greater than 0.5 mm and smaller than 10
mm. Traditionally this variation has been explained by the presence of surfactants
(See [Cl78]), but more recently both Wu and Gharib [WG02] and Tomiyama et
al. [To02] attributed this variation to the manner in which the initial bubbles
were generated. The matter continues to be a matter of discussion, as reported in
[YP03].
We simulate a single air bubble rising in initially quiescent pure water with the
bubble size ranging from 0.5 mm to 30 mm. The numerically predicted terminal
rise velocities are plotted in Figure 7 for comparison with experimental results.
We see a relatively good agreement between the numerical predictions and the
upper bound of the experimental measurements in [To02] within the whole range
of different bubble sizes. Do note that when the bubble sizes range from 2.0 mm to
10 mm, oscillation of the bubble rise velocities and the bubble shapes was observed
in the simulations. For such cases, the terminal bubble rise velocity was calculated
through averaging the rising velocity over a period of time.
GAS BUBBLES RISING IN VISCOUS LIQUIDS
(a)
11
(b)
Figure 8. The trajectory of the mass centre of a bubble rising on
a zigzag path, (a) 3D view and (b) XY plane view.
The bifurcation of the bubble rise velocity found experimentally is not observed
in the numerical predictions. This is as expected since the two probable causes of
this bifurcation are not present in our numerical simulations: i) There are no surfactants present in the water, and ii) We do not simulate the bubble generation process,
but rather assume the existence of a perfectly spherical bubble at the start of our
simulations, eliminating any initial shape disturbances and their effects. However,
including such effects numerically with the aim of reproducing the bifurcation in
the rise velocities would be a most interesting study.
Path Instability of a Rising Bubble
Many experiments (e.g. [WG02][ME00]) have demonstrated that millimetresized bubbles rising in low viscosity liquids do not generally follow a straight trajectory. In pure water, the transition from a straight rise path to zigzag occurs when
the equivalent diameter of the bubble exceeds 1.8 mm. In this regime, the bubbles
exhibit approximately oblate spheroidal shapes, and they rise in zigzag within a
vertical plan or they spiral around a vertical axis. In this paper we applied the
front tracking method to simulate the rise behavior of a single, initially stationary and spherical bubble in a quiescent viscous liquid. Both zigzag and spiral rise
patterns were revealed using the current simulation algorithm.
Figure 8 shows that the numerically predicted rise path of the bubble is a zigzag
when the simulation parameters are set to be Bo = 10.0, Re = 1000, ρl /ρb = 1000
and µl /µb = 100. A 3D view of the trajectory of the bubble mass center is shown in
Figure 8(a). It clearly indicates that the bubble starts moving towards one side as
it begins to rise. When the velocity of the bubble becomes high enough, the bubble
starts following a zigzag path. A projection of the trajectory into the xy-plane
is shown in Figure 8(b), revealing that the zigzag path lies almost entirely in one
single vertical plane.
The mechanism of the bubble rising behavior is closely related to the wake
structure created by the rising bubble. A careful study of the wake structure of
a rising bubble is illustrated in Figure 9, which shows the flow stream line path
around the bubble and a pressure contour in the bubble wake. When the bubble
rising speed is low, a bubble wake with symmetric, closed recirculation rings is
formed, and a low pressure zone is generated at the recirculation centre as shown
12
JINSONG HUA, PING LIN, AND JAN F STENE
(a)
(b)
(c)
(d)
(e)
(f)
Figure 9. Variations of the bubble wake structure for a bubble
rising on a zigzag path.
in Figure 9(a). As the rise velocity increases, the flow instability is amplified and
the bubble wake starts to detach from one side of the bubble bottom as shown in
Figure 9(b). Due to the asymmetric wake structure, the drag and lift forces acting
on the bubble will also become unbalanced, and the bubble is tilted as shown
in Figure 9(b). As the bubble speed increases further, the bubble wake becomes
more asymmetric and the bubble tilting becomes more pronounced. As a result,
the recirculation ring of the bubble wake is fully broken on one side as shown in
Figure 9(c), and the other end of the recirculation ring starts attaching itself to one
bottom side of the oblate bubble as shown in Figure 9(d). Consequently, two open
recirculation rings attached to one side of the bubble bottom are finally formed,
resulting in tilting and lateral movement of the bubble. The lateral movement
makes the open recirculation rings in the bubble wake switch from one side to the
other side of the bubble bottom as shown in Figures 9(e) and (f). The oscillation
of the two open wake recirculation rings from one side of the bubble bottom to the
other causes the bubble to rise on a zigzag path.
GAS BUBBLES RISING IN VISCOUS LIQUIDS
(a)
13
(b)
Figure 10. The trajectory of the mass centre of a bubble rising
on a spiral path, (a) 3D view and (b) XY plane view.
(a)
(b)
Figure 11. Variations of the bubble wake structure for a bubble
rising on a spiral path.
Figure 10 shows a spiral pattern for the bubble rise path as predicted by simulations when the parameters are set to be Bo = 32.0, Re = 1400, ρl /ρb = 1000 and
µl /µb = 100. A 3D view of the trajectory of the mass center of the bubble is shown
in Figure 10(a). It clearly indicates that the bubble starts rising in a spiral with
increasing radius until an almost constant radius is obtained. The top view of the
bubble rise trajectory shown in Figure 10(b) indeed reveals an almost constant spiral radius as the bubble rises. A careful study of the bubble wake structure is shown
in Figure 11. It is found that the wake structure for a spiraling bubble is totally
different from the one observed for a zigzagging bubble. For the spiraling bubble,
only one strong, open recirculation ring is attached to one side of bubble bottom,
the point of attachment being the lowest point of the bubble. This is maybe due to
the high deformability of the spiraling bubble with higher Bond number. The point
of attachment of the open recirculation ring moves along the bottom side edge of
the oblate bubble and causes the bubble to follow a spiral trajectory.
Interaction of Two Rising Bubbles
The problem of a single bubble rising in a viscous liquid is an ideal case for
numerical model validation. However, the final goal when developing a numerical model for multiphase flow is not just investigating the flow behavior of single
14
JINSONG HUA, PING LIN, AND JAN F STENE
(a) τ = 2
(b) τ = 4
(d) τ = 8
(c) τ = 6
(e) τ = 10
Figure 12. The interaction of two initially spherical bubbles rising due to buoyancy.
bubbles rising in viscous liquids, but also investigating multi-fluid systems with
multiple bubbles. With the confidence from validating the current model for a single bubble rising in a viscous liquid, we extend the model to explore the complex
interaction between two bubbles rising in a liquid. It is here worth mentioning
that in the present study, the bubble coalescence process is simplified and based
purely on geometric criteria rather than criteria related to the complex interface
physics. This means that coalescence occurs when nodes on the two bubble surfaces
get close enough, e.g. less than one fifth of the average triangle side length in the
surface mesh. The associated topological change and volume conservation is also
dealt with in a geometric manner.
Figure 12 shows the interaction of two initially spherical bubbles rising in a
quiescent liquid due to buoyancy. The smaller bubble is initially located 2.5D
above the big bubble in vertical direction, and 0.5D axis-off from the big bubble
in the horizontal direction of Y . Here, D represents the effective diameter of the
big bubble. The diameter of the small bubble is half that of the big bubble. The
flow conditions are: Bo = 115.0 and Re = 134.6 for the big bubble, Bo = 28.75 and
Re = 47.6 for the smaller bubble, while ρl /ρb = 1181 and µl /µb = 5000 for both
bubbles. Figure 12 shows the temporal bubble shape evolution of two ring bubbles.
As the big bubble has a higher rising velocity, it will catch up with the small bubble.
When they are close enough, the trailing big bubble is significantly affected by the
low-pressure zone in the wake of the leading small bubble. The trailing bubble
therefore undergoes large deformations and moves towards the bottom wake zone
GAS BUBBLES RISING IN VISCOUS LIQUIDS
(a)
15
(b)
Figure 13. Numerical predictions of (a) the vertical position and
(b) the lateral position of the rising bubbles before and after coalescence.
of the leading bubble as shown in Figure 12(c). Finally, the trailing big bubble
merges with the leading smaller bubble, and a toroidal bubble ring as shown in
Figure 12(e) is formed.
In addition, Figure 13 shows the temporal variation of the position of the
bubbles in both vertical and horizontal directions. It can be seen from Figure 13(a)
that the trailing big bubble has a higher rise speed than the small leading bubble.
The interesting finding is that when the two bubbles are close enough, then the
rise speed of both bubbles increases significantly. After the coalescence of the two
bubbles, the resulting merged bubble represents the familiar situation of a single
bubble rising in a liquid. The lateral movement of the trailing bubble caused by
the leading bubble can be seen in Figure 13(b). Even though the leading bubble
initially moves slightly away from the trailing bubble laterally, this distance is quite
small. However, the trailing bubble, despite its big size, is significantly affected by
the leading bubble and moves towards it.
4. Conclusions
An improved numerical algorithm for the front tracking method (See [HL07])
has been proposed and validated against experiments within a wide range of intermediate Reynolds and Bond numbers using an axisymmetric model. In this paper,
the numerical algorithm is further extended to simulate 3D bubbles rising in viscous
liquids at high Reynolds and Bond numbers with large density and viscosity ratios
typical of the air-water multi-fluid system. Certain mesh adaptation techniques are
implemented for both the front mesh and the background mesh, and the governing
Navier-Stokes equations for incompressible, Newtonian flows are solved in a moving reference frame attached to the rising bubble. The governing flow equations
are solved using a finite volume scheme based on the Semi-Implicit Method for
Pressure-Linked Equations (SIMPLE) algorithm. The 3D bubble surface is tracked
explicitly using an adaptive unstructured triangular mesh. The model is integrated
with the software package PARAMESH, a block-based adaptive mesh refinement
(AMR) tool developed for parallel computing, allowing background mesh adaptation and the solution of the flow equations in parallel on a supercomputer. With
16
JINSONG HUA, PING LIN, AND JAN F STENE
these new features, the modeling algorithm can offer high simulation accuracy and
track the rising bubble for a long duration.
The current model has also been applied to simulate a number of cases of
3D gas bubbles rising in viscous liquids, e.g. air bubbles rising in water. The
simulation predictions on bubble terminal velocity agree well with experimental
measurements in [To02] over a wide range of bubble sizes from 0.5 mm to 30 mm.
The simulation results also reveal that bubbles may rise on a zigzag or/and a spiral
path in a viscous liquid at high Reynolds numbers. The study of the bubble wake
structure qualitatively indicates the mechanisms that lead to a zigzag or spiral
bubble path. It is found that the zigzagging bubble has lower deformability and
two open recirculation rings attached to its bottom in the bubble wake. On the
other hand, the spiraling bubble has a higher deformability and only one open
recirculation ring attached to the side edge of its bottom. In addition, we tried to
use this model to simulate the interaction between two rising bubbles in a liquid.
In conclusion, the simulation results provide us with more physical insights into the
complex behavior of bubbles rising in viscous liquids.
Appendix: Nomenclature
Bo
D
F
g
h
I()
n
p
Re
t
t
u
Bond number (dimensionless)
Effective/Equivalent bubble diameter (m)
Surface tension force(N )
Gravitational acceleration (ms−2 )
Grid size of background mesh
Indicator function (dimensionless)
Unit outer normal on bubble surface
Pressure (N m−2 )
Reynolds number (dimensionless)
Time (s)
Edge vector of bubble surface element
Flow velocity (ms−1 )
Greek letters
∆t Time step size
δ() Delta function (dimensionless)
Γ() Fluid interface
κ Curvature of the bubble
µ Fluid viscosity (Pa s)
ρ Fluid density (kg m−3 )
σ Surface tension (N m−1 )
τ
Non-dimensional time
Subscripts
b Bubble
f
Gas-liquid interface
l
Liquid
∗ Non-dimensionalized parameter
GAS BUBBLES RISING IN VISCOUS LIQUIDS
17
References
[An05] M. van S. Annaland, N. G. Deen and J. A. M. Kuipers, Numerical simulation of gas
bubbles behaviour using a three-dimensional volume of fluid method, Chem. Eng. Sci. 60
(2005), 2999-3011.
[An06] M. van S. Annaland, W. Dijkhuizen, N. G. Deen, and J. A. M. Kuipers, Numerical simulation of behavior of gas bubbles using a 3-D front-tracking method, AIChE J. 52 (2006),
99-110.
[BW81] D. Bhaga and M. E. Weber, Bubbles in viscous liquids: shapes, wakes and velocities, J.
Fluid Mech. 105 (1981), 61-85.
[BM06] J. Bonometti and J. Magnaudet, Transition from spherical cap to toroidal bubbles, Phys.
Fluids 18 (2006).
[Ch99] L. Chen, S. V. Garimella, J. A. Reizes and E. Leonardi, The development of a bubble
rising in viscous liquid, J. Fluid. Mech. 387 (1999), 61-96.
[Cl78] R. C. Clift, J. R. Grace and M. E. Weber, Bubbles, Drops, and Particles, Academic, New
York, 1978.
[Da04] J. Dai, J. D. Sterling and A. Nadim, Numerical simulation of air bubbles in water using
an axissymmetric VOF method, Advances in Fluid Mechanics: Computational Methods in
Multiphase Flow II. WIT Press, 2004.
[Di05] W. Dijkhuizen, E. I. V. van den Hengel, N. G. Deen, M. van S. Annaland and J. A. M.
Kuipers, Numerical investigation of closures for interface forces acting on single air bubbles
in water using Volume of Fluid and Front Tracking models, Chem. Eng. Sci. 60 (2005),
6169-6175.
[HL07] J. Hua and J. Lou, Numerical simulation of bubble rising in viscous liquid, J. Comput.
Phys. 222 (2007), 769-795.
[Ko02] M. Koebe, D. Bothe, J. Pruess and H. J. Warnecke, 3D direct numerical simulation of
air bubbles in water at high Reynolds number, Proceedings of the ASME Fluids Engineering
Division Summer Meeting, 2002.
[KB99] R. Krishna and M. V. Baten, Simulating the motion of gas bubbles in a liquid, Nature
398 (1999), 208.
[LF04] D. Lrstad and L. Fuchs, High-order surface tension VOF-model for 3D bubble flows with
high density ratio, J. Comput. Phys. 200 (2004), 153-176.
[Ma00] P. MacNeice, K. M. Olson, C. Mobarry, R. deFainchtein and C. Packer, PARAMESH: A
parallel adaptive mesh refinement community toolkit, Comput. Phys. Commun. 126 (2000),
330-354.
[ME00] J. Magaudet and I. Eames, The motion of high-Reynolds-number bubbles in inhomogeneous flows, Annu. Rev. Fluid Mech. 32 (2000), 659.
[MM02] G. Mougin and J. Magnaudet, The generalized Kirchhoff equations and their application
to the interaction between rigid body and an arbitrary time-dependent viscous flow, Int. J.
Multiphase Flow 28 (2002), 1837-1851.
[MM06] G. Mougin and J. Magnaudet, Wake-induced forces and torques on a zigzagging/spiraling
bubble, J. Fluid. Mech. 567 (2006), 185-194.
[Oh05] M. Ohta, S. Haranaka, Y. Yoshida and M. Sussman, Three-dimensional numerical simulations of the effect of initial bubble conditions on the motion of a bubble rising in viscous
liquids, J. Chem. Eng. Jpn 38 (2005), 878-882.
[Pa80] S.V. Patankar, Numerical Heat Transfer and Fluid Flow, Hemisphere, 1980.
[Pe77] S. Peskin, Numerical analysis of blood flow in the heart, J. Comput. Phys. 25 (1977),
220-252.
[RR00] F. Raymond and J. M. Rosant, A numerical and experimental study of the terminal
velocity and shape of bubbles in viscous liquids, Chem. Eng. Sci. 55 (2000), 943-955.
[Ru02] H. Rusche, Computational Fluid Dynamics of Dispersed Two-phase Flows at High Phase
Fractions, Ph.D. Thesis, Imperial College of Science, Technology and Medicine, 2002.
[SZ99] R. Scardovelli and S. Zaleski, Direct numerical simulation of free-surface and interfacial
flow, Annu. Rev. Fluid Mech. 31 (1999), 567-603.
[Su99] M. Sussman, A. S. Almgren, J. B. Bell, P. Colella, L. H. Howell and M. L. Welcome, An
Adaptive Level Set Approach for Incompressible Two-Phase Flows, J. Comput. Phys. 148
(1999), 81-124.
18
JINSONG HUA, PING LIN, AND JAN F STENE
[To98] A. Tomiyama, Struggle with computational bubble dynamics, Third International Conference on Multiphase Flow, ICMF’98 Proceedings, Lyon, France, 1998.
[To02] A. Tomiyama, G. P. Celata, S. Hosokawa and S. Yoshida, Terminal velocity of single
bubbles in surface tension force dominant regime, Int. J. Multiphase Flow, Vol. 28, 2002, pp.
1497-1519.
[Tr01] G. Tryggvason, B. B. Bunner, A. Esmaeeli, D. Juric, N. Al-Rawahi, W. Tauber, J. Han,
S. Nas and Y. S. Jan, A front-tracking method for the computations of multiphase flow, J.
Comput. Phys., Vol. 169, 2001, pp. 708-759.
[UT92] S. O. Unverdi and G. Tryggvason, A front-tracking method fror viscous, incompressible,
multi-fluid flows, J. Comput. Phys., Vol. 100, 1992, pp. 25-37.
[YP03] b. Yang and A. Prosperetti, The transient rise of a bubble subject to shape or volume
changes, Phys. Fluids, Vol. 15, 2003, pp. 2640-2648.
[WG02] M. Wu and M. Gharib, Experimental studies on the shape and path of small air bubbles
rising in clean water, Phys. Fluids, Vol. 14, 2002, L49.
Institute of High Performance Computing, 1 Science Park Road, #01-01 The Capricorn Singapore Science Park II, Singapore 117528
E-mail address: [email protected]
Department of Mathematics, National University of Singapore 2, Science Drive 2,
Singapore 117543
Current address: Division of Mathematics, University of Dundee, Dundee, Scotland DD1
4HN, UK
E-mail address: [email protected]
Department of Mathematics, National University of Singapore 2, Science Drive 2,
Singapore 117543
E-mail address: [email protected]