Numerical Simulation of 3D Bubbles Rising in Viscous Liquids using

Numerical Simulation of 3D Bubbles Rising in Viscous Liquids

using a Front Tracking Method

Jinsong Hua

, Jan F. Stene

and Ping Lin

Institute of High Performance Computing, 1 Science Park Road, #01-01 The Capricorn,

Singapore 117528

Department of Mathematics, National University of Singapore, 2 Science Drive 2,

Singapore 117543

Submitted for Publication to

Journal of Computational Physics

May 2007

Manuscript correspondence should be directed to Dr. Jinsong Hua.

Postal address: Institute of High Performance Computing

1 Science Park Road, #01-01 The Capricorn

Singapore 117528

Email address: [email protected].edu.sg

Telephone number: +65-6419 1537

Fax number: +65-6778 0522

Numerical Simulation of 3D Bubbles Rising in Viscous Liquids

using a Front Tracking Method

Jinsong Hua

a ∗

, Jan F. Stene

and Ping Lin

Institute of High Performance Computing, 1 Science Park Road, #01-01 The Capricorn,

Singapore 117528

Department of Mathematics, National University of Singapore, 2 Science Drive 2,

Singapore 117543

Abstract

The rise of bubbles in viscous liquids is not only a very common process in many

industrial applications, but also an important fundamental problem in fluid physics. An

improved numerical algorithm based on the front tracking method, originally proposed by

Tryggvason and his coworkers, has been validated against experiments over a wide range

of intermediate Reynolds and Bond numbers using an axisymmetric model (Hua and

Lou, J. Comput. Phys. 222:769-795, 2007). In the current paper, this numerical algorithm

is further extended to simulate 3D bubbles rising in viscous liquids with high Reynolds

and Bond numbers and large density and viscosity ratios at the physical order of the

typical multi-fluid system of air bubbles in water. To facilitate the simulation, mesh

adaptation is implemented for both the front mesh and the background mesh, and the

governing Navier-Stokes equations for incompressible, Newtonian flow are solved in a

moving reference frame attached to the rising bubble. Specifically, the flow equations are

solved using a finite volume scheme based on the Semi-Implicit Method for Pressure-

Linked Equations (SIMPLE) algorithm, and it appears to be robust even in the range of

high Reynolds numbers and high density/viscosity ratios. 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. PARAMESH allows

background mesh adaptation as well as the solution of the governing equations in parallel

on a supercomputer. The interpolations between the front mesh and the background mesh

are done with Peskin’s distribution function. The current model has also been applied to

∗

Corresponding author, E-mail: [email protected]

simulate a number of cases of 3D gas bubbles rising in viscous liquids, e.g. air bubbles

rising in water. The simulation results are compared with experimental observations both

in aspect of terminal bubble shapes and terminal bubble velocities. In addition, we

applied this model to simulate the interaction between two bubbles rising in a liquid. The

simulation results provide us with more physical insights into the complex bubble rising

behavior in viscous liquids.

Keywords: Computational fluid dynamics; Incompressible flow; Multiphase flow; Bubble

rising; SIMPLE algorithm; Front tracking method; Adaptive mesh refinement; Moving

reference frame

1. Introduction

Multiphase flows are numerous in both everyday life and engineering practice [34].

Typical examples in nature include raindrops in air and gas bubbles in water, whereas

chemical reactions, combustion and petroleum refining are examples of multiphase flows

in industry. One very basic example of such flow is the rise of a single gas bubble in an

otherwise quiescent viscous liquid. The understanding of the flow dynamics of this

system is of great importance in engineering applications and to the fundamental

understanding of multiphase flow physics. Rising bubbles have long been studied

theoretically [8, 25], experimentally [1] as well as computationally through numerical

modeling [36]. While all these efforts have provided us with valuable insights into the

dynamics of bubbles rising in viscous liquids, there are still many questions that remain

unanswered due to the involvement of complex physics. The behavior of a bubble rising

in a viscous liquid is not only affected by the physical properties such as density and

viscosity of both phases [6], but also by the surface tension on the interface between the

two phases and by the bubble shape evolution [27, 2]. The difficulties in describing and

modeling the complex behavior of a rising bubble are to a large extent due to the strong

nonlinear coupling of factors such as buoyancy, surface tension, bubble/liquid

momentum inertia, viscosity, bubble shape evolution and rise history of the bubble. In

addition, the physics of the behavior of bubbles is of a three-dimensional nature. Hence,

most of the past theoretical works were done with a lot of assumptions, and the results are

only valid for certain flow regimes [25, 43]. The experimental works were limited by the

available technologies to monitor, probe and sense the moving bubbles without

interfering with their physics [1, 38, 44].

With the rapid advance of computing power and the development of robust numerical

methods, first principle based numerical simulations promise great potential in extending

our knowledge of the fundamental system of a single bubble rising in a viscous liquid.

However, there are still great challenges and difficulties in simulating such a system

accurately. This may be attributed to the following facts: (i) the sharp interface between

the gas bubble and the surrounding liquid should be tracked accurately without

introducing excessive numerical smearing; (ii) the surface tension gives rise to a singular

source term in the governing equations, leading to a sharp pressure jump across the

interface; (iii) the discontinuity of the density and viscosity across the fluid interface may

lead to numerical instability, especially when the jumps in these properties are high. For

example, the density ratio of liquid to gas could be as high as 1000; (iv) the geometric

complexity caused by bubble deformation and possible topological change is the main

difficulty in handling the geometry of interface; a large bubble may break up into several

small ones, and a bubble may also merge with other bubbles; (v) the complex physics on

the interface, e.g. the effects of surfactants, film boiling and phase change (heat and mass

transfer) and chemical reactions. Fortunately, various methods for multiphase flow have

been developed to address these difficulties, and each method typically has its own

characteristic strengths and weaknesses. Comprehensive reviews of numerical methods

for multiphase/interfacial flow simulation have been given by Scardovelli and Zaleski

[37] and Annaland et al. [41]. Most of the current numerical techniques applied in the

simulation of multiphase/interfacial flows have been developed with focus on the

following two aspects: (i) capturing/tracking the sharp interface, e.g. interface capturing,

grid fitting, front tracking or hybrid methods; and (ii) stabilizing the flow solver to handle

discontinuous fluid properties and highly singular interfacial source terms, e.g. the

projection-correction method [41] and the SIMPLE algorithm [6, 17].

The volume of fluid [14, 4], level-set [28, 42, 29] and phase-field [19] approaches fall

into the first category of front capturing methods. In these methods the interface is

captured using various volume functions defined on the grid used to solve the “one-fluid”

formulation of the governing equations for multiphase flow. Since interface capturing

uses the same grid as the flow solver, it is relatively easy to implement. However, the

accuracy of this approach is limited by the numerical diffusion from the solution of the

convection equation of the volume function. Various schemes have been developed to

advect, reconstruct / reinitialize the volume function to improve the accuracy in

calculating the interface position. One example is the high-order shock-capturing scheme

used to treat the convective terms in the governing equations [18]. Although the explicit

reconstruction of the interface is circumvented, the implementation of such high-order

schemes is quite sophisticated, and they do not work well for the sharp discontinuities

encountered in multiphase/interfacial flows. In addition, a relatively fine grid is needed in

the vicinity of the interface to obtain good resolution.

The second category of approaches tries to track the moving interface by fitting the

background grid points to the interface. The fitting is achieved through re-meshing

techniques such as deforming, moving, and adapting the background grid points. This

method is also well-known as “boundary-fitting approach”, and the “boundary” here

refers to the interface between the fluids. The grid-fitting approach is capable of

capturing the interface position accurately. Early development on this approach was done

by Ryskin and Leal [36]. Curvilinear grids were used to follow the motion of a rising

bubble in liquid. This method is suitable for relatively simple geometries undergoing

small deformations, and applications to complex, fully three-dimensional problems with

unsteady deforming phase boundaries are very rare. This is mainly due to difficulties in

maintaining the proper volume mesh quality and in handling complex interface geometry

such as topological change. In spite of these difficulties, recent work by Hu et al. [16]

showed some very impressive results on 3D simulations of moving spherical particles in

liquid.

The third category is the front tracking method. This approach solves the flow field

on a fixed grid and tracks the interface position in a Lagrangian manner by a set of

interface markers. These interface markers can be free particles without connection, or

they can be logically connected elements, possibly containing accurate geometric

information about the interface such as area, volume, curvature, deformation, etc. A front

tracking technique was proposed by pioneer researchers Glimm and his coworkers [11,

12, 13]. They represent the front interface using a set of moving markers and solve the

flow field on a separate background grid. The background grid is modified only near the

front to make background grid points coincide with the front markers of the interface. In

this case, some irregular grids are reconstructed and special finite difference stencils are

created for the flow solver, increasing the complexity of the method and making it more

difficult to implement. Independently, another front tracking technique was developed by

Peskin and collaborators [31, 10]. In their method, the interface is represented by a

connected set of particles which carry forces, either imposed externally or adjusted to

achieve a specific velocity at the interface. A fixed background grid is kept unchanged

even near the front interface, and the interface forces are distributed onto the background

to solve the “one-fluid” formulation of the fluid flow.

A number of combinations and improvements of these basic approaches have been

proposed to enhance the capabilities in dealing with the sharp, moving interface, where

complex physical phenomena and processes could occur. One of the most promising

approaches is arguably the front tracking method proposed by Tryggvason and his

collaborators [46, 45]. Actually, this method may be viewed as a hybrid of the front

capturing and the front-tracking techniques: a fixed background grid is used to solve the

fluid flow, while a separate interface mesh is used to track the interface position

explicitly. The tracked interface carries any jumps in the fluid properties, such as density

and viscosity, and any interfacial forces, such as surface tension. Fluid properties are then

distributed onto the fixed background grid according to the position of the interface. The

surface tension can be calculated according to the geometry of the interface and is also

distributed onto the background grid in the vicinity of the interface.

Besides the numerical techniques employed to capture/track the moving interface, it

is also very important to develop a stable numerical method to solve the governing

equations of the flow field. Some investigators have considered simplified models such

as Stokes flow [33], where inertia is completely ignored, and inviscid potential flow [15],

where viscous effects are ignored in. In both cases, the motion of deformable boundaries

can be simulated with boundary integral techniques. However, when considering the

transient Navier-Stokes equations for incompressible, Newtonian fluid flow, the so-called

“one-fluid” formulation for multiphase flow has proved most successful [4, 42, 46].

Popular modern methods that use the “one-fluid” formulation include the projection-

correction method [45, 41] and the SIMPLE algorithm [6, 17]. Various

multiphase/interfacial flow problems have been successfully simulated by the front

tracking method [45] with a projection-correction flow solver. It appears that previously

reported results have been limited to flows with low to intermediate Reynolds numbers

(<100) and small density ratios (<100) [5]. It is thus natural to re-examine the approach

and to make it more robust and applicable to wider flow regimes. Some revised versions

of the project-correction method have been proposed to improve its capability in handling

situations with large density and viscosity ratios [41]. Recently, Hua and Lou [17] tested

a SIMPLE-based algorithm to solve the incompressible Navier-Stokes equations. The

simulation results indicated that the newly proposed method could robustly solve the

Navier-Stokes equations with large density ratios up to 1000 and large viscosity ratios up

to 500.

Hua and Lou [17] presented extensive simulations and model validation on a single

bubble rising in a quiescent liquid. The comprehensive simulations show good results in

wide flow regimes with high density and viscosity ratios, and the algorithm is as such

promising in the direct numerical simulation of multiphase flow. Unfortunately, the

previous validation studies were limited to the 2D axisymmetric model where fluid flows

and bubble shapes are axisymmetric. Hence, it would be interesting to investigate the

robustness of the proposed numerical approach for multiphase flow in flow regimes of

higher Reynolds and Bond numbers where the bubble may not be axisymmetric anymore.

Therefore, a fully three-dimensional modeling approach is proposed in this paper. In

addition, other features such as mesh adaptation, moving reference frame and parallel

programming are introduced to enhance the model capability in simulating the rise of a

3D bubble in a viscous liquid.

The numerical algorithm proposed by Tryggvason and co-workers [46, 45], and

extended further by Hua and Lou [17], is adopted in this paper. The handling of the

moving interface in this method may be characterized as a hybrid of interface tracking

and capturing. The governing Navier-Stokes equations are solved on a fixed Cartesian

grid with an adaptive block structure, while the interface is represented by a set of

explicitly tracked front markers. These markers form an adaptive triangular surface mesh

that is advected with a velocity interpolated from the surrounding fluid. An illustration of

such a mesh system is shown in Figure 1. A single set of the governing equations are

solved in the entire computational domain by treating the two fluids as one single fluid

with variable fluid properties across the interface – often referred to as the “one-field” or

“one-fluid” approach. The interface is assumed to have a given finite thickness (normally

about two to four times the background grid size) so that jumps in the fluid properties

across the surface can be reconstructed smoothly by solving a Poisson equation.

A parallel adaptive mesh refinement (AMR) tool, PARAMESH [24], is integrated

with the modified SIMPLE flow solver, and the governing equations are solved in a non-

inertial moving reference frame attached to the rising bubble. The AMR feature allows a

relatively high-resolution mesh in the vicinity of the bubble surface. The non-inertial

moving reference frame technique translates the computational domain with the rising

bubble, allowing the computational domain to be relatively small and always centered

around the bubble. The latter feature is particularly useful for studying the path instability

of a rising bubble or the interaction of multiple bubbles, which may need very long

simulation periods. For example, it is observed in experiments that the paths of

millimeter-sized, rising air bubbles in water normally stabilize after a rise distance of 50-

100 times the initial bubble diameter. If a stationary frame was to be applied to simulate

this situation, the computational domain would be huge compared to the domain of

interest. Even though an AMR feature is adopted, the total number of grid points can still

be a big burden, slowing down the simulation.

The problem of a single bubble rising in a viscous liquid has been widely used as a

typical validation case for the development of new numerical methods for multiphase

flow [6, 40, 41]. Due to numerous experimental and numerical studies in the past, the

physical understanding of the bubble rise behavior in liquid has been well-established in

some flow regimes, e.g. regimes with lower Reynolds and Bond numbers [7, 1].

However, due to the complexity in multiphase flow physics and the difficulties in both

experiments and simulations, the behavior of a rising bubble with high Reynolds number

is not understood well [26, 21]. In this paper, we first validate our model through

comparing computational results of a single bubble rising with experimental results [1] in

aspects of both bubble shapes and terminal velocities. In addition, we apply the 3D model

to simulate air bubbles rising in water, and we compare the terminal velocities of the

bubbles predicted in our simulations with experimental results within a large range of

diameters: from 0.5 mm to 30 mm. There have been few numerical studies on this except

some recent ones [9, 21]. However, since air bubbles rising in water is such a common

process both in our daily life and in industrial applications, a better understanding of this

process is of great importance.

The rest of this paper is organized as follows. In Section 2 we present the governing

equations as well as the numerical method we apply to solve these equations. Numerical

results presented in Section 3 include a detailed sensitivity analysis of the computational

set-up as well as validation through a comparison of our numerical predictions with

available experimental data. Finally, we recapitulate our main findings in Section 4.

2. MATHEMATICAL MODEL AND NUMERICAL METHOD

2.1. Governing Equations

The problem of gas bubbles rising in liquids studied in this paper can be described as

an isothermal, multi-fluid system with two incompressible and immiscible Newtonian

fluids. We will use one single set of governing equations for the entire flow domain

where we treat the different fluids as one single fluid with material properties varying

across the interface. With this “one-fluid” approach, there is no need to deal with the

jump conditions across the interface when the governing equations are solved. However,

we will have to calculate the fluid property distributions and include the surface tension

as a singular source term in the solution domain through the use of a delta function before

the equations can be solved.

The mass conservation for the whole domain under the incompressibility condition

may be expressed in form of volume flux conservation,

⋅

∇

u . (1)

The momentum conservation (Navier-Stokes equations) takes the form,

gxxnuuuu

)()()]([

)(

dsp

ρρδσκμρ

−+−+∇+∇⋅∇+−∇=⋅∇+

∂

∫

, (2)

where u is the fluid velocity,

is the fluid density,

is the density of the liquid phase,

is the pressure,

is the fluid viscosity,

is the surface tension coefficient,

is the

interface curvature,

is the unit normal vector to the interface,

is the gravitational

acceleration, and )

x(x −

is a delta function that is defined as the product of three one-

dimensional delta functions:

)(z)()( yx)(

x , ),zy,(x

x . The subscript refers to a

point on the interface

Γ . It is worth pointing out that the material properties density

and viscosity

will be discontinuous across the interface, and there will generally be a

jump in the pressure

across the interface as well. 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 dimensionless characteristic

variables as follows,

x =

u =

, t

g =

where is the diameter of a sphere with the same volume as the bubble and

D g=g

Thus we may re-express the Navier-Stokes equations as

gxxnuuuu

)1()(

)]([

1)(

−+−+∇+∇⋅∇+−∇=⋅∇+

∂

∫

ρδκμρ

BoRe

, (3)

in which the superscript * has been omitted for convenience. Note that the non-

dimensional Reynolds and Bond numbers used here are thus defined as

2/32/1

= and

= .

By studying the non-dimensional formulation, it can be noticed that the flow is entirely

characterized by the following four dimensionless parameters: The density and viscosity

ratios of the fluids, the Reynolds number and the Bond number. It is also noted that the

definition of the non-dimensional Reynolds number in the experimental works is

different from the one defined here. Normally, experimental works prefer the following

non-dimensional numbers: Eotvos number (

, also known as Bond number); Morton

number (

) and Reynolds number ( ), defined as Re

= ,

σρ

= , and

∞

=Re ,

where is the terminal rise velocity of the bubble measured in the experiments.

∞

2.2. Treatment of the Discontinuities across the Interface

When the governing Navier-Stokes equations are solved numerically on a fixed grid,

the values of the density and viscosity on these grid points are required. It is a reasonable

assumption that each fluid is incompressible, and fluid properties such as density and

viscosity are constant in each fluid phase. Hence the density and viscosity are physically

discontinuous across the interface between the two immiscible fluids, and this abrupt

jump at grid points adjacent to the interface has traditionally caused great problems in

many numerical methods. In Tryggavson et al. [46, 45] a fixed background mesh is

adopted to solve the governing flow equations, and a separated mesh is applied to track

the position of the interface as well as the discontinuities across the front. An illustration

of such a mesh system is shown in Figure 1. The discontinuities across the front are

distributed from the front mesh to the background mesh, and continuous distributions of

the fluid properties on the fixed background mesh can be reconstructed. The singular

source term on the front is distributed to the background grid similarly, and the governing

equations can thus be solved on the fixed background grid using any preferred numerical

approach.

Let us first assume zero interface thickness. Consider the associated reconstruction of

the field distribution of material properties at time

t in the whole domain through

a certain indicator function . Let the indicator function be zero in the liquid phase

and one in the gas phase. We may then write,

),( tb x

),(

),()(),(

tIbbbtb

lbl

⋅

−

, (4)

where is either fluid density or viscosity, and the subscripts

and

refer to liquid

and gas phase, respectively. Further let

),( tb x

be the domain of the gas phase and let

the interface between the two phases. The indicator function may then be expressed as

∫

′′

−=

)(

)(),(

dtI xxxx

. (5)

Taking the gradient of the indicator function and applying Stokes’ theorem, we get

∫∫

ΓΩ

′

−=

′′

−∇=∇

)()(

)()(),(

ddtI sxxnxxxx

δδ

, (6)

where

n is the outer unit normal vector of the interface. Taking the divergence yields a

Poisson equation for the indicator function:

∫

′

−⋅∇=∇

)(

)(),(

dtI sxxnx

. (7)

We solve this equation and then calculate the distribution of material properties from

equation (4).

Unverdi and Tryggvason [46] addressed the sharp jump in fluid properties across the

interface in their front tracking algorithm. They introduced an artificial thickness of the

interface inside which the material properties vary continuously from one fluid to the

other. According to this idea, a distribution function is introduced to approximate

the delta function

)(xD

)(x

with the assumption of an artificial thickness of the interface. We

here adopt the traditional Peskin distribution function [31],

⎪

⎩

⎪

⎨

⎧

<−−+

=−

−

otherwise. 0

2 if )

cos(1)4(

)(

xxxxΠ

. (8)

Substituting the Peskin distribution function into equation (7), the indicator function can

be reconstructed by solving the Poisson equation. The resulting indicator function will

then be zero in the pure liquid phase, vary continuously from zero to one in the artificial

thickness region, and one in the gas phase.

Besides the discontinuity in fluid properties across the bubble interface, the surface

tension, , a singular source term on the bubble interface, brings

another great challenge for numerical methods in multiphase flow. In the current study

the net force caused by surface tension on the surface elements is calculated, thus

circumventing the high-order derivatives involved in curvature calculations. Figure 2

shows the surface tension force exerted on a central surface element (E0) by its

neighboring elements (E1, E2 and E3). The surface tension force acting on an edge

shared between the central element and a neighboring element can be calculated by

∫

−= dsF

)( xxn

δσκ

3 and 2 1,i )(

iii

ntF

, (9)

where is the vector of edge

i and its unit outer normal. Hence, for a central

element E0, the net force caused by surface tension can be expressed as,

0,i

∑∑

×==

)(

ntFF

. (10)

According to equation 10, the net surface tension force on all surface elements can be

calculated and then distributed to the background mesh for solving the momentum

equations:

)()(

∑

−=

XD xFxF

σσ

(11)

where represents the mass centre of the k-th element used to triangulate the bubble

surface.

2.3. Tracking the Moving Interface

With the techniques introduced in Sections 2.1 and 2.2, the governing equations can

be solved on a fixed background grid to obtain the flow field. An adaptive, unstructured

triangular mesh (front markers) is used to represent the interface between the two fluid

phases. Hence, the velocity of the moving front markers can be obtained by interpolation

from the flow field on the background mesh, and then the front mesh points can be

advected in a Lagrangian manner. Thus the front moves with the same velocity as the

surrounding fluid, and the so-called no-slip condition of the interface is satisfied. In this

paper, the interpolation is carried out using the same distribution function as the one used

for the transfer of fluid properties to the background grid:

∑

−=

xxxuxu )(),(),(

fff

Dtt , (12)

Δ+=

uxx . (13)

As the front marker points are advected, the mesh size and quality may consequently

change. The resolution of the front mesh has a strong effect on the information exchange

with the fixed background grid, which may eventually affect the accuracy of the

simulation results. Therefore, it is of key importance that the front mesh has a more or

less constant quality and uniform size throughout the duration of the simulation. To

ensure this, the front quality is examined at each time step and adapted when necessary.

In this paper, the resolution of the triangular mesh for the 3D surface of the bubble is

maintained more or less uniform through adaptation as the interface evolves.

2.4. Mesh Adaptation

As two sets of mesh are applied in the current front tracking method, the resolution of

the front mesh and the background mesh near the front plays an important role in

resolving the interfacial physics of the multiphase flow. From physical principles it is

known that fluid particles on the bubble interface will move downwards towards the

bottom of the bubble as the bubble rises. Similarly, the mesh points also move

downwards on the bubble surface in the front tracking method as the bubble rises. As a

result, the mesh on the upper part of the bubble becomes coarser. On the other hand, the

mesh at the lower part of the bubble becomes increasingly dense. Figure 3(a) shows the

variation of front mesh quality as the bubble rises in liquid without front mesh adaptation.

It is obvious that the accuracy will be affected when the mesh on the top is too coarse,

and that the dense fine mesh at the bottom of the bubble will consume excessive

computing power without much benefit in accuracy. Thus, the front mesh adaptation as

shown in Figure 3(b) is essential to ensure the accuracy and efficiency of the simulation.

In this aspect, three basic operations are adopted to adapt the front mesh, namely edge

swap, edge split and edge deletion. For long edges, the edge swap operation as shown in

Figure 4(a) is a simple and easy operation to improve the mesh quality. In the edge split

operation as shown in Figure 4(b), a new point is generated by surface fitting of the

existing neighboring mesh points. This new point is then inserted into the two associated

meshes, and new and finer triangles are generated to replace the old ones. Thus edge split

is important to refine the front mesh. For the deletion of short edges as shown in Figure

4(c), the triangles associated with the short edge will be deleted, and the resulted gap will

be sealed through merging the old nodes of the short edge. Hence, edge deletion is

important to coarsen the front mesh. In addition, consistent checking of the mesh

connectivity is also important to ensure the accuracy in calculating the surface tension.

Accordingly, the resolution of the background mesh also plays an important role in

capturing the flow behavior – particularly so in the vicinity of the interface. If the

background mesh resolution is too low, then the detailed flow dynamics will not be

captured reasonably well, resulting in unreliable and inaccurate simulations. Therefore, it

is desirable to have relatively high-resolution grids, particularly near the interface, while

coarse grids may be used away from the interface. This is achieved in our model by the

use of the block-based adaptive mesh refinement (AMR) tool PARAMESH [24]. In this

study, the refining and coarsening of the grid blocks is based on whether there exists a

bubble front within the blocks. An example of the block-wise Cartesian mesh refinement

generated by PARAMESH can be seen in Figure 5. It can be seen that fine background

grids are located in the vicinity of the bubble, while coarser background grids are applied

in regions further away from the bubble front. This feature makes it more efficient to

solve the governing equations and to capture the flow physics near the interface

accurately. An excellent feature of PARAMESH is that the Cartesian grid at all levels of

blocks has the same structure. Hence, once the flow solver is developed for one grid

block, it can be easily applied to all levels of blocks. In addition, the different blocks can

be distributed to different CPUs in an MPI parallel environment, which speed up the

problem solving cycle.

2.5. Flow Solver

The numerical method for the governing flow equations is one of the key components

in the simulation. Traditionally, an explicit projection-correction method based on

second-order central differences on a regular, staggered Cartesian grid has been used

along with the front tracking approach [46]. However, it seems that this approach is

unable to handle the large density ratios typical of systems in industrial applications as

well as in nature. In search of a more robust solver, Hua and Lou [17] implemented a

modified version of the classical SIMPLE method [30] for axisymmetric multiphase

flow. Further details about the modified SIMPLE algorithm can be found in the past

works [17, 6]. Computational results indicate that this approach is more robust than the

projection correction method – especially for multiphase flows with large density ratios

such as the air-water system. This improvement is most probably due to the fact that the

SIMPLE algorithm avoids solving the problematic pressure equation directly. Instead, the

pressure and the velocity are corrected iteratively based on the governing equations.

2.6. Moving Reference Frame

In many applications of multiphase flow it is often desirable to study the long-term

behavior and evolution of the moving interface between the fluids. In such applications

the front may move a considerable distance, and the study of the rise path of an air bubble

in water is a typical example. The computational domain must then be correspondingly

large to accommodate such extensive movement. However, in a three-dimensional model

with a high-resolution grid, a large computational domain is computationally very

expensive. Computational cost may therefore limit the domain size and thus also long-

time simulations.

To remedy this problem, we have therefore incorporated a moving reference frame

into our numerical algorithm. The idea is to move the reference frame together with the

front such that the front (e.g. a bubble) remains more or less fixed in the computational

domain. The size of the computational domain may then be chosen independently of the

duration of the simulation, and this will in turn reduce the computational cost

significantly. As a result, we may carry out long-time simulations of moving interface

problems which could not be done in a stationary reference frame.

Figure 6 illustrates a moving reference frame. There, the frame

stands for a

stationary reference frame and the frame

′

for a moving reference frame. The

positions of a monitoring point in the frames

and

′

are represented as and

, respectively, which are correlated with the position of the moving reference frame

( ) according to equation (14). The velocity of the monitoring point

′

is in the

frame

)t,(xu

and in the moving frame

),'( txu

′

, and the velocity of the moving frame is

. The following correlation can be obtained: )(

u t

pmp

xxx

′

(14)

),()(),(

ttt

xuuxu

′

. (15)

Allowing translational, but not rotational, movement of the frame in the present

study, the following is the updated governing flow equations in the moving reference

frame:

gxxn

)1()(

)]([

1)(

−+

′

−

′

∇

′

∇

′

⋅∇

′

+∇

′

−=+

′′

⋅∇

′

∂

′

∂

∫

ρδκ

μρρ

. (16)

′

⋅

∇

′

u (17)

The moving front will generally be accelerating and so will the moving reference

frame. Thus the frame of reference in which we solve the governing equations is no

longer an inertial frame, and we must therefore modify the momentum equations to take

into account the acceleration of the frame. In addition, according to equation 15, when

the governing equations are solved in a moving reference frame, the velocity condition on

the boundary ( ) should be modified as

x )(),(),'(' ttt

mBB

uxuxu

−

Notice that the additional term on the left-hand side

dtd

u , which denotes the

acceleration of the moving reference frame, is added in Equation (16). We aim to

choose so that the rising bubble remains as fixed as possible in the moving frame, i.e.

ideally the acceleration of the frame is equal to the acceleration of the bubble. The bubble

acceleration is of course unknown, so we need to approximate this acceleration at each

time step. We shall adopt the prediction as presented by Rusche [35], namely

−=

, (18)

where

−

=Δ

−

xxxx

λλ

. (19)

Here is the position of the centre of mass of the bubble relative to the moving frame

at time step

, , and

)(

−

−=Δ

nnn

ttt

and

are appropriate under-relaxation factors.

It was found that

= 1.0

gave good results.

2.7. 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 , is calculated through interpolation of

the fluid velocity field

according to Equation (12).

(2)

The front is advected to its new position

1+n

x by using the normal interface

velocity

u found in step (1) – see Equation (13). The front elements are then

subject to examination for adaptation and topological change. Meanwhile, volume

conservation is enforced.

(3)

The indicator function

)(

11 ++ n

is computed based on the interface position

1+n

x . This is done by solving the Poisson problem in Equation (7) with the

discrete delta distribution from Equation (8). Subsequently, the distribution of the

I x

density

1+n

, the viscosity

1+n

and the surface tension

1+n

F is updated on the

flow solver grid points.

(4)

We find the velocity field

1+n

u and the pressure

1+n

p 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

2+n

t .

3. RESULTS AND DISCUSSION

In this Section, we report various numerical results for different purpose, e.g. model

sensitivity analysis, validation case, and model capabilities exploration. A summary of all

the simulations cases can be found in Table I, in which all simulation parameters are

listed.

3.1. Sensitivity Analysis - Size of the Computational Domain

Extensive experiments have been performed in the past to study the rise and

deformation of single bubbles in quiescent liquids. Often these experiments have been

done in large containers with a size of at least 20 bubble diameters in each spatial

direction to avoid wall containment effects [1]. In this paper, we intend to validate

simulation results against such experiments. To achieve this, the computational domain

should also be rather large to avoid any significant effects caused by the wall

confinement. On the other hand, if the domain is chosen too large, excessive computing

time is needed to complete the simulation. To analyze the influence of the domain size, a

number of numerical tests were run with domain sizes ranging from two to twelve bubble

diameters in each of the spatial dimensions. The grid resolution was kept constant and

accommodated approximately twenty cells inside the bubble in each direction. All

computations were carried out in a moving reference frame.

The terminal rise velocity and the terminal bubble shape were used to assess the wall

confinement effects for various domain sizes. The aim of this sensitivity analysis is to

find the smallest possible computational domain in which wall containment effects have

negligible impact on the bubble terminal velocity and shape. Figure 7 presents the

simulation results of the terminal bubble shape and the rise velocity using different

domain sizes under the conditions of

0.243Bo*

, 24.15*R

e , 1000/ =

and

100/ =

(Case A3). There is a notable change in the terminal bubble shape and rise

velocity as the domain size is increased from two to six bubble diameters. When the

domain size is increased beyond six bubble diameters, no significant change in the

simulation results is observed. Moreover, it is noted that the wall confinement has a

strong effect on the terminal velocity for small domain sizes from two to six bubble

diameters. However, the change in terminal velocity is only around 1% when increasing

the computational domain size from eight to ten bubble diameters in each spatial

dimension. Based on these observations we conclude that a domain size with side length

of eight bubble diameters should be sufficient in our simulations. Actually, experimental

results by Krishna et al. [23] also indicate that wall effects on a rising bubble is negligible

when the diameter of the liquid container is larger than eight bubble diameters.

3.2. Sensitivity Analysis - Grid Resolution

PARAMESH divides the computational domain into a number of blocks in each

spatial direction. Each block consists of a certain number of grid cells in each direction,

and the governing equations are discretized and solved numerically on these grid cells.

For a computational domain of fixed size, there are therefore two ways to change the grid

resolution using PARAMESH: either by changing the number of blocks, which is

determined by the maximum number of refinement levels, or by changing the number of

cells in each block. In our grid sensitivity analysis we kept the maximum refinement level

fixed and changed the number of cells in each block. All simulations were done in a

moving, cubic computational domain with side length equal to eight bubble diameters

using equal grid spacing in each of the spatial directions.

The results from the grid sensitivity analysis can be found in Figure 8 under the

conditions of ,

0.243Bo* = 24.15*R

e , 1000/

and 100/

(Case A3). It is

noted that the terminal bubble velocity is highly sensitive to the mesh resolution up until

16 cells per bubble in each space direction. However, increasing the number of cells from

16 to 20 yields less than 1% change in terminal velocity. The terminal bubble shape also

has a very strong dependence on the grid resolution - especially when the grid resolution

is relatively low. On the other hand, we can even see a slight change in shape when we

increase the number of cells per bubble diameter from 16 to 20. However, increasing the

number of cells per bubble diameter further to 32 does not lead to any visually detectable

change in bubble shape. Based on the above we find it sufficient to use a resolution of 20

cells per bubble diameter in our simulations.

3.3. Comparison of Results Obtained in a Stationary and a Moving Reference Frame

Since a moving reference is introduced in this paper to perform long-time simulations

of single bubbles rising as well as two-bubble interactions, it is important to evaluate the

accuracy and impact of using a moving reference as opposed to a stationary frame. For

this purpose, simulation results for Case A4 obtained in a stationary reference frame are

compared with those obtained in a moving reference frame.

The bubble shapes predicted in both stationary and moving frames at different time

steps are shown in Figure 9. There are no visually observable differences between the

results from the two frames. A comparison of the velocity profiles of bubbles rising in a

stationary and a moving reference frame is shown in Figure 10, and they are in

reasonable agreement. Note that the velocity of the moving frame has been added to the

rise velocity obtained in the moving reference frame to enable a direct comparison with

the rise velocity in the stationary frame. In addition to comparing bulk behavior of rising

bubbles, we would also like to compare the detailed flow patterns around the bubbles.

A comparison of the streamlines predicted in a stationary and a moving frame can be

found in Figure 11. Since one frame is at rest and one is moving, comparing streamlines

would only be sensible if we modify the velocity in one of the frames to account for the

different velocities of the frames. Here we have chosen to subtract the velocity of the

moving reference frame from the velocity field computed in the stationary frame before

comparison. It can be seen from Figure 11 that an excellent agreement of the streamline

patterns in the stationary and the moving frame is obtained. Furthermore, pressure

distributions predicted in a stationary and a moving frame are compared in Figure 12, and

again the results agree well.

Based on the various tests and comparisons of different flow characteristics carried

out and described above, we can conclude that the use of a moving reference frame yields

numerical results equivalent to those obtained in a stationary reference frame.

3.4. Model Validation with Experiments

In this section we will compare experimental results available in the literature [1]

with predictions obtained by our numerical method. In Figure 13 we compare observed

and predicted terminal bubble shapes for a range of Reynolds and Bond numbers, and the

results agree very well. In Table II we compare the associated terminal rise velocities,

and again there is reasonable agreement between experiments and our numerical

predictions. However, note that the relative deviation in Case A1 is a little bit high. This

may be due to the low rise velocity where the relative error will be high even though

there is no change in the simulation accuracy.

Air bubbles rising in water are common in many industrial processes. Examples in

chemical engineering include bubble columns, loop reactors, agitated stirred reactors,

flotation, or fermentation reactors. For the design of efficient two-phase reactors, detailed

knowledge of bubble sizes and shapes, slip velocities, internal circulations, swarm

behaviors, bubble induced turbulence and mixing, and bubble size distributions

(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 studies have been performed to

study air bubbles rising in water [7, 44]. Their measurements of the terminal rise velocity

of air bubbles in water are presented in Figure 14 as a function of the bubble size. It is

found that the measurements of the 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 [7], but more recently both

Wu and Gharib [47] and Tomiyama et al. [44] attributed this variation to the manner in

which the initial bubbles were generated. The issue continues to be a matter of discussion

- refer to Yang and Prosperetti [48].

Due to difficulties in measuring the physical properties on the bubble, a fundamental

understanding of the system of a single bubble rising in high Reynolds number regimes is

not well-established. With the recent rapid increase in computing power, numerical

simulations of two-phase flows based on continuum mechanics models with moving free

interfaces have become feasible and proved extremely useful for a better understanding

of fundamental processes and phenomena. However, numerical modeling of the multi-

fluid system of air bubbles rising in water is still quite challenging due to the large

density ratio of water to air, the low liquid viscosity of water, high Reynolds numbers,

and large bubble deformations. Koebe et al. [21] started early trials of 3D direct

numerical simulation of air bubbles rising in water at high Reynolds number using the

volume of fluid (VOF) method. They studied bubbles with diameters from 0.5 mm to 15

mm, and their numerical predictions on the terminal rise velocity of the bubbles agree

reasonably with experimental data. However, they introduced some initial white noise in

the simulations, which may introduce non-physical perturbations to the simulation

system. The recent work by Dijkhuizen et al. [9] reported their trial on simulation of

single air bubbles rising in initially quiescent pure water using both a 3D front tracking

method and a 2D VOF method for bubble diameters ranging from 1 mm to 8 mm. The

calculated terminal rise velocites by the 3D front tracking method are quite close to the

experimental observations by Tomiyama et al. [44], but they over-predicted the velocity

for bubble diameters larger than 3 mm.

In this paper, we use the front tracking method with features of mesh adaptation and

moving reference frame, allowing a finer mesh in the region of the bubble surface.

Consequently, better accuracy is obtained in the current simulations. We simulate a single

air bubble rising in initially quiescent pure water with the bubble diameters ranging from

0.5 mm to 30 mm. The numerically predicted rise velocities of the bubbles agree well

with the upper bound of the experimental measurements by Tomiyama et al. [44] within

the whole range of different bubble sizes. When the bubble diameter is in the range from

2.0 mm to 10 mm, oscillation of the bubble rise velocity and the bubble shape is also

predicted in the simulations. The terminal bubble rise velocity is calculated through

averaging the instantaneous rise velocity over a period of time. Since we assume the

initial bubble shape to be spherical and the surface tension coefficient to be constant, the

bifurcation of the bubble rise velocity is not revealed in the current simulation. However,

this is an interesting topic to be explored in the future.

3.5 Numerical Studies on the Interaction between Two Rising Bubbles in Viscous Liquid

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 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 would like to extend this model to explore the complex interaction between

two bubbles rising in a liquid. Figures 15 and 17 illustrate the simulation of the

interaction of two initially spherical bubbles rising in a quiescent liquid due to buoyancy.

In the first simulation, one smaller bubble is initially located 2.5D above a bigger

bubble in vertical direction, and 0.5D axis-off from the bigger bubble in the horizontal

direction of Y. Here, D represents the effective diameter of the bigger bubble. The

diameter of the smaller bubble is half that of the bigger bubble. The flow conditions for

the bigger bubble are as follows:

6.134*R

0.115Bo*

1181/

and (Case

B1). Figure 15 shows the temporal bubble shape evolution of two rising bubbles. As the

bigger bubble has a higher rise velocity, it will catch up with the smaller bubble (Figure

15(

5000/ =

μμ

0.4=

)). When they are close enough, the trailing bigger bubble is significantly

affected by the low-pressure zone in the wake of the leading smaller bubble. The trailing

bubble undergoes large deformations and moves towards the bottom wake zone of the

leading bubble (Figure 15(

0.6=

)). Finally, the trailing big bubble merges with the

leading smaller bubble, and a toroidal bubble ring is formed (Figure 15(

0.10=

)).

Similar bubble shape evolution patterns have also been predicted by other numerical

predictions [40]. In addition, Figure 16 shows the temporal variation of the position of the

bubbles in both vertical and horizontal directions. It can be seen from Figure 16(a) that

the trailing bigger bubble has a higher rise speed than the smaller 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 returns to a normal situation of a single bubble. The lateral

movement of the trailing bubble caused by the leading bubble can be seen in Figure

16(b). Even though initially the leading bubble 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.

In the second simulation, the smaller bubble is initially located 2.5D above the bigger

bubble in vertical direction, and 1.0D center-off from the big bubble in the horizontal

direction of Y. Here, D represents the effective diameter of the bigger bubble. The

diameter of the smaller bubble is half that of the bigger bubble. The flow conditions for

the bigger bubble are as follows:

24.15*R

0.243Bo*

, and

1181/ =

ρρ

5000/ =

(Case B2). Figure 17 shows the temporal evolution of the two rising

bubbles. As a result, in this case the bigger bubble finally overtakes the smaller one, since

the big bubble is further horizontally off-set the small one. It can be seen from Figures 17

and 18(b) that the leading bubble first starts moving laterally away from the trailing

bubble, whereas the trailing bubble then starts moving towards the leading bubble before

the overtaking occurs (Figure 17(

0.4

, 0.6

, 0.8

)) . After the bigger bubble has

overtaken the smaller one, the smaller bubble is significantly affected by the wake of the

bigger bubble. In fact, the smaller bubble is attracted to the wake of the bigger bubble,

resulting in a highly deformed and elongated bubble shape as shown in Figure 17

(

0.12=

0.14=

0.16=

). It is also noticed from Figure 18(a) that the smaller bubble

is accelerated and rises fast in the wake of the bigger bubble, and finally it catches up and

merges with the bigger bubble (Figure 17(

0.18

, 0.20

, 0.22

)). In this case, the

smaller trailing bubble has little effect on the rising speed of the bigger leading bubble.

4. CONCLUSION

The numerical algorithm used in this paper is a further extension of the algorithm

given in [17] for simulating 3D gas bubbles rising in viscous liquids at high Reynolds and

Bond numbers for systems with large density and viscosity ratios such as air/water. To

achieve this, mesh adaptation is implemented for both the front mesh and the background

mesh, and a moving reference frame attached to the rising bubble is used to solve the

governing incompressible Navier-Stokes equations. The solution method is a finite

volume scheme based on the Semi-Implicit Method for Pressure-Linked Equations

(SIMPLE), and the solver appears to be robust even in the range of high Reynolds

numbers and high density/viscosity ratios. The 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. It includes features such as background mesh adaptation and

parallel implementation of solvers for the governing equations on supercomputers. The

interpolations between the front mesh and the background mesh are done with Peskin’s

traditional approximation of the delta function. The current model has been applied to

simulate a number of examples of 3D gas bubbles rising in viscous liquids, e.g. air

bubbles rising in water. The simulation results are compared with experimental

observations in aspect of both the terminal bubble shape and the terminal bubble velocity.

In addition, we use this model to simulate the interaction between two bubbles rising in

liquid. The simulation results provide us with some physical insights into the complex

behavior of bubbles rising in viscous liquids.

REFERENCES

D. Bhaga and M.E. Weber, Bubbles in viscous liquids: shapes, wakes and

velocities,

J. Fluid Mech. 105, 61 (1981).

T. Bonometti and J. Magnaudet, Transition from spherical cap to toroidal bubbles,

Phys. Fluids. 18, 052102 (2006).

L.A. Bozzi, J.Q. Feng, T.C. Scott, and A.J. Pearlstein, Steady axisymmetric

motion of deformable drops falling or rising through a homoviscous fluid in a

tube at intermediate Reynolds number,

J. Fluid Mech. 336, 1 (1997).

J.U. Brackbill, D.B. Kothe, and C. Zenmach, A continuum method for modeling

surface tension,

J. Comput. Phys. 100, 335 (1992).

B. Bunner, and G. Tryggvason, Dynamics of homogenous bubbly flows. Part I:

Rise velocity and microstructure of the bubbles,

J. Fluid Mech. 466, 17 (2002).

L. Chen, S.V. Garimella, J.A. Reizes, and E. Leonardi, The development of a

bubble rising in a viscous liquid,

J. Fluid Mech. 387, 61 (1999).

R. Clift, J.R. Grace, and M.E. Weber, Bubbles, Drops, and Particles, (Academic

Press, San Diego, 1978).

R.M. Davies, and F.I. Taylor, The mechanism of large bubbles rising through

extended liquids and through liquids in tubes,

Proc. R. Soc. Lond. A 200, 375

(1950).

W. Dijkhuizen, E.I.V. van den Hengel, N.G. Deen, M. van Sint 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 method

models,

Chem. Eng. Sci. 60, 6169 (2005).

10.

A.L. Fogelson, and C.S. Peskin, A fast numerical method for solving the three-

dimensional Stokes equations in the presence of suspended particles,

J. Comput.

Phys.

79, 50 (1988).

11.

J. Glimm, O. McBryan, R. Menikoff, and D.H. Sharp, Front tracking applied to

Rayleigh–Taylor instability,

SIAM J. Sci. Statist. Comput. 7, 230 (1986).

12.

J. Glimm, J.W. Grove, B. Lindquist, O. McBryan, and G. Tryggvason, The

bifurcation of tracked scalar waves,

SIAM J. Sci. Stat. Comput. 9, 61 (1988).

13.

J. Glimm, J.W. Grove, X.L. Li, and D.C. Tan, Robust computational algorithms

for dynamic interface tracking in three dimensions,

SIAM J. Sci. Comput. 21,

2240 (2000).

14.

C.W. Hirt and B.D. Nichols, The Volume of Fluid (VOF) method for the

Dynamics of Free Boundaries,

J. Comput. Phys. 39, 201 (1981).

15.

T.Y. Hou, J.S. Lowengrub, and M.J. Shelley, Boundary integral methods for

multicomponent fluids and multiphase materials,

J. Comput. Phys. 169, 302

(2001).

16.

H.H. Hu, N.A. Patankar, and M.Y. Zhu, Direct numerical simulation of fluid-solid

system using the arbitrary Lagrangian-Eularian technique,

J. Comput. Phys. 169,

427 (2001).

17.

J. Hua and J. Lou, Numerical simulation of bubble rising in viscous liquid. J.

Comput. Phys.

222, 759 (2007).

18.

M. Ida, An improved unified solver for compressible and incompressible fluids

involving free surfaces. Part I. Convection.

Computer Physics Communications

132, 44 (2000).

19.

D. Jacqmin, Calculation of two-phase Navier-Stokes flows using phase-field

modeling,

J. Comput. Phys. 155, 96 (1999).

20.

S.I. Kang, and G.L. Leal, Numerical solution of axisymmetric, unsteady free-

boundary problem at finite Reynolds number. I. Finite difference scheme and its

application to the deformation of a bubble in a uniaxial straining flow,

Phys.

Fluids

30, 1929 (1987).

21.

M. Koebe, D. Bothe, J. Pruess, and H.-J. Warnecke, 3D direct numerical

simulation of air bubbles in water at high Reynolds number, In Proceedings of the

2002 ASME Fluids Engineering Division Summer Meeting (FEDSM2002-

31143), Montreal, Quebec, Canada, July 14-18, 2002.

22.

A. Koynov, J.G. Khinast, and G. Tryggvason, Mass transfer and chemical

reactions in bubble swarms with dynamics interface,

AIChE J. 51, 2786 (2005).

23.

R. Krishna, M.I. Urseanu, J.M. van Baten, and J. Ellenberger, Wall effects on the

rise of single gas bubbles in liquids,

Int. Comm. Heat Mass Transfer 26, 781

(1999).

24.

P. MacNeice, K.M. Olson, C. Mobarry, R. deFainchtein, and C. Packer,

PARAMESH: A parallel adaptive mesh refinement community toolkit.

Comput.

Phys. Commun.

126, 330 (2000).

25.

D.W. Moore, The rise of a gas bubble in viscous liquid, J. Fluid Mech. 6, 113

(1959).

26.

G. Mougin and J. Magnaudet, Wake-induced forcesand torqueson a zigzagging /

spiraling bubble,

J. Fluid Mech. 567,185 (2006).

27.

M. Ohta, T. Imura, Y. Yoshida, and M. Sussman, A computational study of the

effect of initial bubble conditions on the motion of a gas bubble rising in viscous

liquids,

Int. J. Multiphase Flow 31, 223 (2005).

28.

S.J. Osher and J.A. Sethian, Fronts propagating with curvature dependent speed.

Algorithms based on Hamilton-Jacobi formulations,

J. Comput. Phys. 79, 12

(1988).

29.

S. Osher and R.P. Fedkiw, Level set methods: An overview and some recent

results,

J. Comput. Phys. 169, 463 (2001).

30.

S.V. Patankar, Numerical Heat Transfer and Fluid Flow, (Hemisphere, New

York, 1980).

31.

C.S. Peskin, Numerical analysis of blood flow in the heart, J. Comput. Phys. 25,

220 (1977).

32.

C.S. Peskin and B.F. Printz, Improved volume conservation in the computation of

flows with immersed boundaries,

J. Comput. Phys. 105, 33 (1993).

33.

C. Pozrikidis, Interfacial dynamics for stokes flow, J. Comput. Phys. 169, 250

(2001).

34.

A. Prosperetti, Bubbles, Phys. Fluids 16, 1852 (2004).

35.

H. Rusche, Computational Fluid Dynamics of Dispersed Two-Phase Flows at

High Phase Fractions, Ph.D. Thesis, Imperial College of Science, Technology &

Medicine, 2002.

36.

G. Ryskin and L.G. Leal, Numerical simulation of free-boundary problems in

fluid mechanics. Part 2. Buoyancy-driven motion of a gas bubble through a

quiescent liquid,

J. Fluid Mech. 148, 19 (1984).

37.

R. Scardovelli and S. Zaleski, Direct numerical simulation of free-surface and

interfacial flow.

Annu.Rev. Fluid. Mech. 31, 567 (1999).

38.

W.L. Shew and J.-F. Pinton, Viscoelastic effects on the dynamics of a rising

bubble,

J. Stat. Mech. P01009, (2006).

39.

S. Shin, and D. Juric, Modelling three-dimensional multiphase flow using a level

contour reconstruction method for front track with connectivity.

J. Comput. Phys.

180, 427 (2002).

40.

M. van Sint 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, 2999 (2005).

41.

M. van Sint Annaland, W. Dijkhuizen, N.G. Deen, and J.A.M. Kuipers,

Numerical simulation of behaviour of gas bubbles using a 3-D front-tracking

method,

AIChE J. 52, 99 (2006).

42.

M. Sussman, P. Smereka, and S. Osher, A level set approach for computing

solutions to incompressible two-phase flows,

J. Comput. Phys. 114, 146 (1994).

43.

T.D. Taylor, and A. Acrivos, On the deformation and drag of a falling viscous

drop at low Reynolds number,

J. Fluid Mech. 18, 466(1964).

44.

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

28, 1497 (2002).

45.

G. Tryggvason, B.B. Bunner, A. Esmaeeli, D. Juric, N. Al-Rawahi, W. Tauber, J.

Han, S. Nas, and Y.J. Jan, A front-tracking method for the computations of

multiphase flow.

J. Comput. Phys. 169, 708 (2001).

46.

S.O. Unverdi and G. Tryggvason, A front-tracking method for viscous,

incompressible, multi-fluid flows.

J. Comput. Phys. 100, 25 (1992).

47.

M.M. Wu and M. Gharib, Experimental Studies on the Shape and Path of Small

Air Bubbles Rising in Clean Water,

Phys. Fluids 14, 49 (2002).

48.

B. Yang, A. Prosperetti, and S. Takagi, The Transient Rise of a Bubble Subject to

Shape or Volume Changes,

Phys. Fluids 15, 2640 (2003).