Group Title: 7th International Conference on Multiphase Flow - ICMF 2010 Proceedings
Title: 4.4.2 - A Novel Lattice Boltzmann Method for Direct Simulation of Bubble Dynamics
Full Citation
Permanent Link:
 Material Information
Title: 4.4.2 - A Novel Lattice Boltzmann Method for Direct Simulation of Bubble Dynamics Computational Techniques for Multiphase Flows
Series Title: 7th International Conference on Multiphase Flow - ICMF 2010 Proceedings
Physical Description: Conference Papers
Creator: Yu, Z.
Fan, L.-S.
Publisher: International Conference on Multiphase Flow (ICMF)
Publication Date: June 4, 2010
Subject: lattice Boltzmann method
adaptive mesh refinement
Abstract: The traditional lattice Boltzmann method (LBM) for two-phase flow simulation is often hindered by numerical instability at low viscosities and insufficient resolution at the gas-liquid interface. As a result, they usually lack the capability of simulating bubbles with high Reynolds numbers or large deformation. To overcome such a problem, a novel LBM technique is developed with multi-relaxation-time (MRT) and adaptive mesh refinement (AMR) scheme. The new LBM approach is based on the improved interaction potential model, which utilizes a new force incorporation scheme that is able to maintain grid-independent fluid properties and can be easily integrated into the MRT scheme. The MRT scheme is found to greatly enhance the numerical stability at low viscosities, and consequently enables the study of high-Reynolds number bubble flows. On the other hand, as bubbles move in the computational domain, fine grid resolution is dynamically adapted to regions near the bubble surface in order to accurately capture the bubble shape. The novel algorithm that incorporates MRT and AMR is presented in detail in this study. Numerical simulations of the buoyant rise of bubbles are conducted under various conditions. The shape and velocity of the simulated bubbles are compared with experimental results and correlations in the literature. The improvement of the current method in accuracy and stability may significantly enhance the capability of LBM in the simulation of complex multiphase flows under realistic conditions.
General Note: The International Conference on Multiphase Flow (ICMF) first was held in Tsukuba, Japan in 1991 and the second ICMF took place in Kyoto, Japan in 1995. During this conference, it was decided to establish an International Governing Board which oversees the major aspects of the conference and makes decisions about future conference locations. Due to the great importance of the field, it was furthermore decided to hold the conference every three years successively in Asia including Australia, Europe including Africa, Russia and the Near East and America. Hence, ICMF 1998 was held in Lyon, France, ICMF 2001 in New Orleans, USA, ICMF 2004 in Yokohama, Japan, and ICMF 2007 in Leipzig, Germany. ICMF-2010 is devoted to all aspects of Multiphase Flow. Researchers from all over the world gathered in order to introduce their recent advances in the field and thereby promote the exchange of new ideas, results and techniques. The conference is a key event in Multiphase Flow and supports the advancement of science in this very important field. The major research topics relevant for the conference are as follows: Bio-Fluid Dynamics; Boiling; Bubbly Flows; Cavitation; Colloidal and Suspension Dynamics; Collision, Agglomeration and Breakup; Computational Techniques for Multiphase Flows; Droplet Flows; Environmental and Geophysical Flows; Experimental Methods for Multiphase Flows; Fluidized and Circulating Fluidized Beds; Fluid Structure Interactions; Granular Media; Industrial Applications; Instabilities; Interfacial Flows; Micro and Nano-Scale Multiphase Flows; Microgravity in Two-Phase Flow; Multiphase Flows with Heat and Mass Transfer; Non-Newtonian Multiphase Flows; Particle-Laden Flows; Particle, Bubble and Drop Dynamics; Reactive Multiphase Flows
 Record Information
Bibliographic ID: UF00102023
Volume ID: VID00107
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: 442-Yu-ICMF2010.pdf

Full Text

Paper No 7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30-June 4, 2010

A Novel Lattice Boltzmann Method for Direct Numerical Simulation of Bubble Dynamics

Zhao Yu and L-S Fan*

William G. Lowrie Department of Chemical and Biomolecular Engineering
The Ohio State University
140 West 19th Avenue, Columbus, Ohio 43210, USA and fan.1

Keywords: lattice Boltzmann method, bubble, adaptive mesh refinement


The traditional lattice Boltzmann method (LBM) for two-phase flow simulation is often hindered by numerical instability at
low viscosities and insufficient resolution at the gas-liquid interface. As a result, they usually lack the capability of simulating
bubbles with high Reynolds numbers or large deformation. To overcome such a problem, a novel LBM technique is developed
with multi-relaxation-time (MRT) and adaptive mesh refinement (AMR) scheme. The new LBM approach is based on the
improved interaction potential model, which utilizes a new force incorporation scheme that is able to maintain
grid-independent fluid properties and can be easily integrated into the MRT scheme. The MRT scheme is found to greatly
enhance the numerical stability at low viscosities, and consequently enables the study of high-Reynolds number bubble flows.
On the other hand, as bubbles move in the computational domain, fine grid resolution is dynamically adapted to regions near
the bubble surface in order to accurately capture the bubble shape. The novel algorithm that incorporates MRT and AMR is
presented in detail in this study. Numerical simulations of the buoyant rise of bubbles are conducted under various conditions.
The shape and velocity of the simulated bubbles are compared with experimental results and correlations in the literature. The
improvement of the current method in accuracy and stability may significantly enhance the capability of LBM in the
simulation of complex multiphase flows under realistic conditions.


Bubbles exist ubiquitously in many natural and industrial
systems. However, their behaviour is far from been fully
understood. Gas bubbles rising in a still liquid can have
various deformed shapes and follow different patterns in
their motion as the consequence of the complex interactions
between the surface force, the buoyancy force, and the
viscous force. For a large number of bubbles in a gas-liquid
chemical reactor, the interactions between bubbles further
complicate their dynamics by affecting their velocity and
spatial distribution. The bubbles dynamics further
influences the mass transfer between the gas and liquid
phase as well as the rate of chemical reactions, and affects
the overall performance of the reactor (Fan, 1989).
Therefore, the dynamics of bubbles not only poses
intriguing scientific questions, but is also of great
importance to various engineering applications.

Direct numerical simulation is an effective technique to
study bubble dynamics, since they are able to capture the
motion and deformation of the bubble by tracking the
location of the gas-liquid interface. Among various direct
numerical simulation techniques, the Lattice Boltzmann
Method (LBM) is a unique approach that is based on the
mescoscale description of the distribution of the fluid
molecules. Compared to other more traditional
Computational Fluid Dynamics (CFD) approaches such as

Volume of Fluid (VOF), Level Set, and Front Tracking
Method, the LBM is able to model the segregation between
the gas and liquid phase using physical concepts such as
molecular interactions in the fluid, and therefore avoid the
need to explicitly track or capture the gas-liquid interface.
In addition, the LBM has efficient algorithm and is easy to
implement. These advantages make LBM a promising
technique for bubble simulation.

However, there are still some limitations of the LBM that
restrict its application in various two-phase flow problems.
One of theses issues that have received considerable
attention recently is the numerical instability of LBM at low
viscosities, even for single phase flows. In fact, a number of
techniques have been proposed to overcome this problem.
Some researchers attribute the instability to the
nonexistence of the H theorem in LBM, and accordingly
they developed the entropic LBM by introducing the
entropy function and a variable relaxation time (Ansumali
and Karlin, 2002). The entropic LBM has shown better
stability in several single phase flow simulations. However,
its extension to multiphase flow problems has not been
reported so far. Other researchers have used an implicit
formulation of the collision term, which is treated using
central difference in both space and time (Sankaranarayanan,
et al, 2002). The implicit scheme has achieved enhanced
stability at low viscosity values, and is able to simulate
bubbles with Reynolds number up to 400. The third class of

Paper No

techniques is based on the multiple-relaxation-time (MRT)
scheme, which employs multiple relaxation parameters in
the collision step, instead of the single-relaxation-time BGK
collision in traditional LBM. A number of the two-phase
models have been reformulated with the MRT algorithm,
including the color function model (Tolke et al. 2006), the
incompressible index function model (Premnath and
Abraham, 2007), and the pressure-evolution model
(Mukherjee and Abraham, 2007). The interaction potential
model for non-ideal gas has also been recently extended
with MRT algorithm, and was reported to achieve moderate
improvement in stability (Kuzmin, et al, 2008). However,
only simple flow conditions such as static droplets in
equilibrium and capillary waves were simulated.

From the application point of view, successful application of
LBM in bubble dynamics is restricted in a small parameter
window as the consequence of its numerical instability at
low viscosity values. A gas bubble rising in a viscous liquid
is often characterized by the dimensionless groups such as
Morton number (Mo) and Reynolds number (Re):
Mo =gp (1)

Re = (2)

In the above equations, g is the gravitational acceleration,
p and V are the density and kinematic viscosity of the
liquid, a is the surface tension, ub is the terminal rise
velocity of the bubble, and d, is the equivalent bubble
diameter. The most common system with an air bubble
rising in water has a Morton number of 0(1011), and
millimeter-size air bubbles rising in water usually have a
Reynolds number on the order of 102 103. Unfortunately,
with its instability for low-viscosity fluids, the LBM is often
limited to higher Mo (Mo>10 5) and lower Re (Re<102)
conditions. For example, Frank et al (2006) used the free
energy LBM approach to simulate millimeter bubbles in
99.5% glycerol, which gave a Moron number of 0(101) and
0.033 Poisson equation was used to simulated bubbles with
Mo=105, although at a high density ratio of 103 (Inamuro et
al, 2004). A phase field LBM for immiscible fluids was also
applied for bubbles with 10- (Kurtoglu and Lin, 2006). For the interaction potential LBM
approach, similar applicable range is achieved, as shown by
the coarse bifurcation studies (Mo-10-4) (Theodoropoulos et
al, 2004) and multiple bubble dynamics studies
(10-6 To the best of the authors' knowledge, the closest LBM
simulation for millimeter air bubbles in water to date was
based on the implicit LBM formulation (Sankaranarayanan,
et al, 2002), which achieved a Morton number of 10-10, and
Reynolds number up to 400. However, even that does not
match the conditions of the air-water system exactly.
Therefore, an apparent gap exists between the capability of
current two-phase LBM models and the actual flow
conditions in many low-viscosity liquids, and there is
clearly a demand to further improve the LBM model.

Another critical issue in direct simulation of bubble
dynamics is the accurate representation of the interface

7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30-June 4, 2010

between the two phases. As in many other numerical
approaches, the LBM use a diffused interface which spans
typically to a thickness of 3-5 grid spacing. On the other
hand, at least 16 grid points across the diameter are usually
required to faithfully represent the bubble shape. The
minimum grid number required increases when the bubble
undergoes large deformation. At the same time, the
simulation domain needs to be significantly larger than the
bubble size to avoid boundary effects. However, both
computation time and memory usage increase rapidly with
the increasing mesh resolution, and therefore the resolution
requirement often has to be compromised given current
computation capability. As the result, most simulations to
date consider only small bubbles with spherical or
ellipsoidal shapes, while simulations for large spherical cap
and skirted bubbles that frequently appear in engineering
applications are scarcely reported.

Adaptive Mesh Refinement (AMR) technique is a natural
candidate to meet the requirement in both resolution and
computation cost. In AMR, fine mesh resolution is used
near the bubble surface to ensure the accuracy, while coarse
mesh is applied in the bulk fluid faraway from the interface
to reduce computation cost. The mesh resolution is updated
dynamically during the simulation in accordance with the
motion of the interfaces. Although there have been a few
studies that apply AMR to traditional CFD approaches using
VOF or Front Tracking method, the application of AMR in
the LBM has not been widely reported. An important reason
is that the basic LBM algorithm, in which the particle
populations collide on lattice sites and propagate along
different lattice directions, is closely coupled with a uniform
mesh. Recently, some research has been conducted to
extend the LBM to non-uniform grids, mainly for the
modelling of single phase flows. In such simulations, fine
mesh is applied near the irregular solid boundary to enhance
the accuracy of the boundary condition, or near the
high-vorticity region to resolve the complex flow field.
Since they do not involve moving objects, the mesh
resolution is specified in the beginning of the simulation,
without being dynamically updated. Despite the promising
results obtained for single phase flows, there is no
straightforward extension of the non-uniform grid LBM
algorithm to two-phase flow problems. There seems to be
only one publication concerning the multiphase LBM-AMR
(Tolke, et al., 2006). In that study, the interface algorithm
was based on the color function model, which was similar
to the traditional interface capturing methods. Only a single
condition of bubble rising was demonstrated in the result,
without experiment comparison for broader regimes.

The motivation of the current wok is to develop a new
interaction potential model based two-phase LBM to
address the numerical instability and accuracy problem. By
using a slightly different way to incorporate the interaction
force, the interaction potential model is integrated into the
MRT algorithm in a straightforward manner. Novel
formulation of the interaction force, such as the utilization
of multi-range potential, can be readily adapted in the newly
developed model. By using a novel AMR scheme, the
simulation is able to accurately track the bubble
deformation and solve the flow field near the bubble surface.
Simulation results in 2D and 3D for bubbles under different

Paper No

conditions are compared to the experimental results and
correlations in the literature. The new LBM technique is
shown to result in significant improvement in numerical
stability and enhanced capability in capturing bubble


lattice velocity
bubble diameter
Edtvds number
particle distribution function
gravitational constant
interaction strength
Morton number
Reynolds number
lattice weight

Greek letters
A Collision matrix
p density
a surface tension
r relaxation parameter in BGK collision
y interaction potential
0 Collision integral
v kinematic viscosity

eq equilibrium

Numerical Scheme

The multiphase LBM method in this work is based on the
interaction potential model originally developed by Shan
and Chen (1993). In this model, an interaction potential is
defined as a function of the local fluid density, and the
interaction force calculated from the interaction potential is
employed to results in the segregation of the two phases.
Unlike the traditional model, the interaction potential based
LBM in this study has several advanced features. First, the
governing equation is written in a more general form and
force is incorporated in a straightforward manner. This not
only makes it easy to reformulate the equation using the
MRT scheme, but also ensures the independence of the fluid
property from the grid size. Second, the collision step is
computed in the moment space with MRT scheme to
enhance the numerical stability. And finally, the algorithm is
incorporated with the capability to detect the interface and
adjust the grid resolution accordingly. The detailed
numerical method is presented in the following sections.

* Multiphase multi-relaxation-time LBM
The LBM is based on the Boltzmann equation, which can be
written as
af + f F f (3)
at ax p av
where f(x, V, t) is the particle distribution function, and x,
v, t, p F are the spatial coordinates, particle velocity, time,
fluid density, and force, respectively. The right hand side of

7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30-June 4, 2010

Equation (3), Q, is the collision term that describes the
change of particle distribution function due to particle

In the LBM, the velocity space is discretized into a finite set
of velocities {ci corresponding to a regular lattice
structure in space, and accordingly the particle distribution
is discretized into {fi (x, t)}. The collision is described as
a relaxation process towards local equilibrium. The
equilibrium distribution feq (p, ) is a truncated
Maxwell distribution that only depends on local fluid
density p and velocity U which can be calculated
directly from the distribution functions fi In general the
collision process involves multiple physical quantities that
may relax on different time scales, and information for
those time scales can be given using a full constant
matrix A :
ni (f-ff) (4)

As the result, the LBE in the general form can be expressed
f,(x, +c,At, t + At)- f(x,,t) =-A,,(f -f,"r) + AtS, (5)

The left hand side of Equation (5), which is often referred to
as the streaming or propagation step, corresponds to the
unsteady and convective terms in the Boltzmann Equation
(3). The right hand side of Equation (5) is referred to as the
collision step. The last term Si in Equation (5) represents the

force term F f in Equation (3), and its detailed
p 3v
implementation will be discussed in more detail in later

Equation (5) provides the basic working equation for
general LBM algorithms. It is noted that traditional LBM
with BGK collision can be regarded as a special case in
which the collision matrix is a diagonal matrix with
identical diagonal elements ofl/T :

A = i (6)

Interaction potential model
For gas-liquid two-phase flow, the interaction potential
model, also called the Shan-Chen model, is one of the most
widely used LBM models (Shan and Chen, 1993). It
employs a mean-filed interaction force to mimic the
molecular interactions that cause the segregation of the gas
and liquid phases. An interaction potential If/ is defined
based on the local fluid density, and the interaction force is
calculated from the interaction potential to induce proper
phase separation. Two versions of interaction models have
been developed: the single component model for a non-ideal
fluid, and the multi-component model for a mixture of two
different fluids.

In the single component model, the interaction force is
given as the summation of the pair-wise interactions

Paper No

between particles at a given lattice site and those at
neighboring sites, and can be written as
F(x) = -G y(x)Z wi(x + ci)ci (7)
In the above equation, G is a scalar constant that represents
the strength of the interaction. The interaction potential
Sis a local quantity that depends on the fluid density.. In
this study, the following form of the potential is employed:
y(p) = 1 exp(-p) (8)
Upon Taylor expansion of Equation (7), the force can be
written as

F = -V( cGy)- c2GyVV2+o(At3), (9)
2 2
in which the first term on the right hand side corresponds to
the non-ideal part of the EOS, while the second term
contributes to the surface tension. More detailed effects of
the interaction force on the thermodynamic properties of the
gas-liquid interface are provided in a recent work by Shan
(2008) by considering the discrete lattice effect. Using the
approach developed in that study, the properties of the
gas-liquid interface can be found analytically, including the
equilibrium densities, stress tensor, density profile across
the interface, and the surface tension coefficient.

On the other hand, when two fluid components are involved,
two distributions fja (a = 1,2) are used for the two
components, and they each evolve according to the
governing LB Equation (5). The interaction force now
includes both the interaction between particles of the same
component and the particles of different components:
2 q
F, () = G (x) wy, (x + c,)c (10)
o-=1 i=l
The interaction strength parameters G11 and G22 describe the
interaction within each individual component. When
modeling a gas-liquid system, usually the first component is
assumed to be a non-ideal fluid with its EOS determined by
a potential as the one given in (8), while the second
component is an ideal gas with G22=0 and V2 (p )= p2,.
G12 and G21 have the same value, and they describe the
interaction between the two components and control the
immiscibility of the mixture.

In the original interaction potential model, only the first
layer of neighbor sites is employed in the calculation of the
interaction force. However, recent progress shows that
several favorable features can be obtained by calculating the
interaction force with the mutli-range potential, which
involves an enlarged set of grid points in more layers of
neighbor sites (Sbragaglia, et al, 2007). First, the extra
degree of freedom leads to a surface tension that can be
adjusted independently from the equation of state. Second,
the magnitude of the spurious velocity can be reduced by
either enlarging the interface thickness, or by employing a
higher-order isotropic discretization. Since it is only
involved in the formulation of the interaction force, the
mutli-range potential can be readily integrated into the MRT
model as well as the traditional interaction potential model.

7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30-June 4, 2010

Incorporation of the force term
Incorporating the force by shifting the equilibrium velocity,
which is called the Shan-Chen forcing scheme, is a simple
and unique feature of the original interaction potential
model. However, although the Shan-Chen forcing scheme is
easy to implement with the BGK collision, its extension to
the MRT framework is not straightforward. On the other
hand, the force term in the continuous Boltzmann equation
can be approximated using the gradient of the equilibrium
distribution (He, et al, 1998):
Faf_ Faf (v-u)F(11)
p av p av pRT

With a finite set of discrete velocities, the forcing term in
the LBM becomes
(ci u) F
,- = feq (12)
Then, following the derivation given by Premnath and
Abraham (2007), using a transformation:

fi = fi SiAt, the general LB Equation (5) is written

f(x+cAt,t+A t)- = f,(x,t) A, +

A1A S, At
2 )

in which lij is the components of the identity matrix.
Since Equation (13) applies to both single and multiple
relaxation time algorithms, incorporation of the interaction
force into the MRT approach becomes straightforward. The
macroscopic fluid density is calculated as in traditional

p=>/j (14)

And as the result of the transformation, the momentum of
the fluid is calculated from

pu = fc + AtF (15)
It is also noted that, when using a single relaxation time,
Equation (12) reduces to the forcing scheme proposed by
Guo et al. (2002).

* Solution of MRT-LBE
After incorporating the force term, Equation (13) becomes
the governing equation that dictates the evolution of the
distribution functions in MRT-LBM. Since the collision
matrix A is in general a full matrix, directly solving
Equation (13) involves complex matrix manipulations.
However, in practice, the solution of Equation (13) takes
advantage of a special linear transformation to diagonalize
the collision matrix. The transformation matrix T
transforms the distribution functions fi, which lie in the

velocity space, into their moments f which correspond
to macroscopic physical quantities such as density,
momentum, energy, and their fluxes that lie in the moment

Paper No

f = Tf (16)
The specific form of the transformation matrix depends on
the lattice structure. The procedure presented in this section
corresponds to the D2Q9 lattice in 2D. The 3D cases using
D3Q19 lattice has been given in the reference (Premnath
and Abraham, 2007). With D2Q9 lattice, the transformed
vector f in the moment space is given explicitly by:

f= ,e,e2,jx, q,jy,qy,pxxp (17)
where e is the energy, jx and j are the momentum in x

and y direction; qx and q are the energy fluxes; px

and py are the diagonal and off-diagonal components of
the stress tensor. The transformation matrix is expressed
explicitly as:

The equilibrium moments are obtained by applying the
transformation to the equilibrium distributions originally in
the velocity space, and can be computed directly from the
hydrodynamic quantities such as density and momentum.

2p+(j + j2)/9p
p-(jl+ jy)/9p
jx (19)
f = Tfeq = -
(j2 jy2)p
liJ y//P
In the same way, the forcing term is also transformed into
the moment space.
6(uxF, +uFy)
-6(uxF +uyFy)
F, (20)
S =TS= -F
2(uxFx -uyFy)
uxF, +uF
By multiplying the transformation matrix T, the right hand
side (RHS) of Equation (13) can be entirely transformed
into the moment space:

7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30-June 4, 2010

T(RHS)= -1, p- f ,9) +z lA j g At (21)
9 + "(,> (I 2at (21)
The transformed collision matrix = TAT- is now
diagonal in moment space:
A = diag[sl,s2, s3 s4, s, s, S,, Ss9 (22)
The diagonal elements sl through s9 are the new relaxation
parameters associated with each components of f Since
the new collision matrix is diagonal, the relaxations of
different physical quantities to fjeq are now decoupled.
Different time scales determined by sl through s9 can be
adjusted independently, although a few constraints still
apply. To be consistent with the macroscopic hydrodynamic
equation, it is required thatS1 = S4 = S6 = 1. Meanwhile,
s8 and s9 are related to the kinematic viscosity V,

S8 = Sg =1/ r, and V= T- j CAt (23)

And s2 is related to the bulk viscosity of the fluid by

S= I At (24)
S2 2
Further, symmetry requires that S5 = S As the result,
three of the relaxation parameters s2,S3, S )(= S7)remain
to be independently adjustable, and they can be used to tune
the stability of the MRT model.

In practice, the MRT-LBM algorithm is often realized in
such a way that the collision step takes place in the moment
space, while the propagation step still operates in the
velocity space. In the actual computation procedure,
transformed quantities f feq and S are first
calculated using Equations (17), (19) and (20), followed by
the collision step as described by Equation (21). Then the
RHS is transformed back to the velocity space to perform
the propagation step:
f (x + cAt, t + At) f(x, t) = T1 [T(RHS)] (25)
The hydrodynamic quantities such as density, momentum,
potential and force are calculated in the velocity space in the
usual way as discussed in previous sections.

In summary, the current MRT interaction potential method
models the two-phase flow based on the interaction force
calculated from the interaction potential, as in the classic
Shan-Chen model. However, the interaction force is directly
incorporated into the LBM without shifting the equilibrium
velocity. In the collision step, the particle distributions are
transformed into their moments which have distinct
physical meanings. Their relaxation time scales are
uncoupled and can be adjusted independently with the
parameters Sl-S9. The separation of different relaxation
scales brings improved stability to the MRT interaction
potential model. In addition, since the transformation
between velocity space and moment space is a local
operation, the transformed moments f do not need to be
stored. Therefore, the MRT algorithm only requires a small
amount of additional computation time compared to

Paper No

conventional LBM, without significant increase in memory

* Adaptive mesh refinement in LBM
When applying different mesh resolutions inside the
computation domain, it is crucial to determine which
variables change with the mesh size and which variables are
independent of the grid. In order to achieve the same lattice
speed of sound 1 Ax across the entire domain, the
1c At
1 2 3


i t I I I I
t t I I I I
tA :^ :

t^ :^ :

7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30-June 4, 2010

time step needs to be proportional to the mesh size.
Following the convention of LBM, the non-dimensional
lattice speed of sound is kept to be c throughout

the entire domain.


7 (=1)

_i A_

Collision in explode propagation Collision in propagation coalesce
fine and in fine and fine in fine
coarse coarse
Figure 1. Schematic illustration of the change of particle distributions in the cells near the refinement jump in a 2D

The fact that the time step At varies with Ax has two
consequences. Firstly, according to Equation (23), to keep
the fluid kinematic viscosity V independent of grid size,
the relaxation factor T also has to be adjusted with Ax
and must be calculated from Equation (23). Explicitly, this

r 1 = T 1- Ax (26)
S2 c 2) Ax,
Secondly, the LBM algorithm performs more steps on the
fine mesh since it requires a smaller At than on the
coarse mesh. Therefore, interpolation/average operations
need to be performed at the coarse/fine resolution interface
in both space and time.


S (1) collide (3) propagate
Coarse level --
(2) explode (6) coalesce
Fine level
(1) collide (3) propagate (4) collide (5) propagate

Figure 2. Flowchart of the computation procedure during
one time step for two refinement levels. The sequence of
sub-steps (1-6) is in accordance with those in Figure 1.

The newly developed LBM-AMR approach is based on the
multi-block structured grid. The entire computation domain

is divided into a number of blocks that have identical grid
structure but may have different grid sizes. The blocks are
categorized into different refinement levels, and grid
resolutions in each two consecutive levels differ by a factor
of 2. In order to have a smooth transition of the grid size, at
most one refinement level difference is allowed at the
boundary between two neighboring blocks. Inside each
block, uniform structured grid is used as in the traditional
LBM simulation. In this way, the original LBM algorithm
can be applied within each individual block without
modification. The classical D2Q9 lattice is employed for
the 2D computations in this study, and the D3Q19 lattice is
used for the 3D simulations.

When a refinement level jump is present between two
neighboring blocks, a special algorithm is devised to
communicate the information across block boundaries.
Similar to the volumetric formulation of Rohde et al.
(2006), the particle distributions in the current method are
located at the center of each computational cell, so that
rescaling of the non-equilibrium part of the distribution at
coarse/fine grid interface is not necessary. The inter-block
communication process is shown schematically in Figure 1.
The computational cells at the boundary between the
coarse and fine grid are referred to as the "interface" cells,
as shown in the shaded zones in Figure 1.They act as the
transition layer between two blocks, and alternate between
coarse and fine levels during a computation step through
splitting and merging. Two specially designed operations
take place in these interface cells, and they are named
"explosion" and "coalescence" after the terminology used
by Chen et al. (2006). In the 2D example illustrated in
Figure 1, a coarse interface cell is uniformly divided into


Paper No

four smaller cells during the explosion step (step 2-3 in
Figure 1), and its particle population is either evenly
distributed or linearly interpolated in these new cells. In the
coalescence step (step 6-7 in Figure 1), the four fine
interface cells again merge into a single coarse level cell,
whose particle population is the summation of those in the
four fine cells.

The original "collision-propagation" algorithm is now
complicated by the additional "explosion" and
"coalescence" operations. The flowchart of the
computation procedure on a mesh with two refinement
levels is shown in Figure 2. As noted previously, the
computation on coarse level uses dt as its time step, while
the fine grid has the time step of 0.5dt. Therefore, the
"collision-propagation" operations happen twice on the
fine level in each dt, but only once on the coarse level. The
"explosion" and "coalescence" steps each operates once
during dt. The "explosion" takes place after the first
collision step at 0.5dt, which transfers the information from
coarse level to the fine level. And the "coalescence"
happens at the end of the time step so as to pass the
information back from the fine level to the coarse level.

The change of particle distributions in each sub-step is
illustrated graphically in Figure 1, in which only one pair
of opposing lattice velocities are shown for clarity. The
first collision step operates on both fine and coarse levels
as in ordinary LBM, but with different relaxation
parameters. The interface cells, which are originally on the
coarse level, are then divided into 4 fine cells in the
"explosion". The distribution functions f, can be either
uniformly distributed or interpolated linearly in the new
fine cells. Propagation step then occurs on both coarse and
fine levels, including the interface cells. Up to this point,
the coarse level has completed a time step dt, while the fine
level only completes a half dt. Therefore, the second
collision and propagation process only occurs on the fine
level, and distribution functions only change on the fine
and interface cells while the information on the coarse
level is unaffected. Finally, the four fine interface cells
again merged into a coarse level cell, and the distribution
functions in the new coarse interface cell is just the sum of
those in the four fine cells. This finishes the entire
computation loop and both levels come back to the
pre-collision status before step 1.

The generation and refinement of the mesh is performed
with the open-source code PARAMESH (MacNeice et al,
2000), which use a tree data structure to manage the
hierarchy of block-structured grid with different refinement
levels. It has been integrated with the non-uniform grid
LBM code developed by the authors.

The computation procedure of multiphase LBM on a
non-uniform mesh is essentially the same as the single
phase algorithm, and the only additional step is
incorporation of the interaction force. Calculation of the
force involves taking the gradient of the potential, which
requires information from neighboring cells. When this
occurs at the coarse-fine grid interface, the potential on the
other side of the interface can be easily
averaged/interpolated and then communicated across the

7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30-June 4, 2010

block boundary. Once the force is obtained, the forcing
scheme itself is a local operation and can be performed at
the same time as the collision step.

The adaptive refinement of the mesh takes place every few
steps in the simulation. For multiphase flows which
involve two fluid phases with different densities, it is
natural to use the magnitude of density gradient as the
criterion for mesh refinement. The interface between two
phases is located by finding the blocks whose density
gradient is above a critical value, and these blocks are
assigned to be at the highest refinement level. Regions
away from the interface have successively larger grid size.
It is important that the highest and lowest refinement levels
are specified to avoid infinite refinement or excessive
coarsening during the simulation.

Results and Discussion

The accuracy and numerical stability of two-phase LBM
simulations often deteriorate at low viscosities due to the
growing magnitude of the unphysical spurious velocity,
which is the small circulating velocity near the interface
between two phases in the numerical results. The origin of
the spurious velocity has been recently identified to be the
insufficient isotropy in evaluating the gradient terms for
force calculation (Shan, 2008). Since the MRT scheme
involves additional freely adjustable relaxation parameters,
it is expected to benefit from the extra flexibility to reduce
the spurious velocity.


Figure 3: Magnitude of the spurious velocity as a
function of viscosity in the single component model. The
SC model with BGK collision is unstable with viscosity
less than 0.025.

The magnitude of the spurious velocity as a function of
viscosity is compared in Figure 3 for three different
versions of single component model: the single relaxation
BGK model, the MRT model, and the MRT model with
multi-range interaction potential. The spurious velocity
shown is the maximum value measured from the numerical
results for a 2D circular bubble at equilibrium. For all three
models, the spurious velocity increases as the viscosity
decreases. When the viscosity approaches 0.025, the BGK
model results in a sharp increase of the spurious velocity,
which ultimately destabilizes the simulation so that no
results could be obtained for viscosity smaller than 0.025
in the BGK model. On the other hand, the MRT model
with the same interaction force gives consistently smaller

SMRT mulirange


0.01 0.1

Paper No

spurious velocity compared to the BGK model. As a result,
fluid with viscosity as low as 0.002 can be easily simulated
using the MRT model.

Recent studies have shown that the multi-range potential
can be used to reduce the spurious velocity in the
interaction potential model by utilizing more neighboring
lattice sites during the evaluation of the interaction force
(Sbragaglia, et al, 2007). It should be noted that
adjustments in the multi-range potential may have the
effect of changing the surface tension value. This effect can
be understood from Equation (9), in which the coefficient
in front of the term yVV2V/is proportional to the square
of the surface tension. A larger value of this coefficient
corresponds to both higher surface tension, and wider
interface. These two factors actually have opposite effects
on the amplitude of the spurious velocity, but numerical
results show that their combined effect is to decrease the
spurious velocity.

The spurious velocity in the MRT scheme with multirange
potential is also plotted in Figure 3. In these simulations,
two layers of lattice neighbors (24 points) are employed in
the force calculation, and result in an interaction force that
is isotropic up to the 8th order. In comparison, the original
force evaluation using only the first layer of 8 lattice points
can only reach 4th order isotropic descretization. The
spurious velocity in the MRT-multirange model is
consistently smaller than that in the MRT model under all
viscosities, and its magnitude can often be reduced to
below half of that in the MRT model with regular force
evaluation method.

7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30-June 4, 2010

identical condition where G=-5.0 and V = 0.05 The
dilute phase inside the circle has a density of 0.11, and the
dense phase outside the circle has a density of 1.85. The
interface shown in the figures are plotted with density
contours of 0.5 and 1.0. The maximum amplitude of the
spurious velocities in BGK mode, MRT model, and MRT-
multirange model are 0.028, 0.0053, and 0.0016,
respectively. The vector plots Figure 4(a)-(c) clear
demonstrates that the combination of MRT and multi-range
potential model can effectively reduce the spurious

Due to the improved numerical stability at small viscosities,
the newly developed MRT-LBM is able to successfully
simulate millimetre-sized bubbles in a low viscosity liquid,
which are characterized by small Morton number and high
Reynolds number, and are therefore difficult for traditional
two-phase LBM. The 3-D bubble simulations are carried
out using the two-component MRT model on D3Q19
lattice. The equivalent diameter of the gas bubble is about
30 lattice units. The size of the simulation domain is
150*150*300 lattice units, which is believed to result in
negligible boundary effect on the bubble dynamics. The
simulation first runs without gravitational force for a few
thousand steps to let the system equilibrate. Once the
bubble volume reaches a constant value, the gravity is
turned on and simulation continues to run until the bubble
reaches a steady velocity. The simulation results are
compared with experimental results or empirical
correlations in the literature using the dimensionless groups
Mo, Eo and Re to validate current MRT-LBM model.

The first set of simulations are conducted at Mo=6.2*10-7,
corresponding to millimeter gas bubbles in silicone oil
DMST-05 (Zenit and Magnaudet, 2008). The equivalent
bubble diameters are in the range between 1.5mm and 3.5
mm, which give the Eitvis numbers ( Eo= pdg
between 1.0 and 5.6. Figure 5(a) shows the Reynolds
number of the bubbles at their terminal velocities. The
simulation results are in good agreement with both the
values measured in the experiment (Zenit and Magnaudet,
2008) and the empirical correlation given by Fan and
Tsuchiya (1990):

S Eo + 2C Eo1f2 (27)
u (KMo1/4 1Eo1/2 2
In the above equation, Kb, C, and n are parameters
determined by the liquid property and surface condition,
and are chosen to be 15, 1.2, and 1.0, respectively.

Figure 4: Spurious velocity at the gas-liquid interface in
2D. The spurious velocity is plotted under same
magnification in the three figures. (a) BGK model,
|u1 =0.028 ,u, =0.0053
lu.= 0.028 ; (b) MRT model, ;0.005 (c)
MRT model with 8th order isotropic force,
Ilm.I = 0.0016

The spurious velocity near a 2D circular interface is shown
in Figure 4. The three simulation are carried out under

Paper No

0 1 2 3 4 5

6 7 8

2 4 6

1.5mm 2.0mm 2.65mm



Figure 5: (a) Reynolds numbers for 1.5-3.5mm air
bubbles in silicone oil DMST-05; (b) Aspect ratios of
1.5-3.5mm air bubbles in silicone oil DMST-05; (c)
Bubble shape obtained from corresponding to the
simulation results in (a) and (b)

The bubble shape under the simulated conditions is in the
ellipsoidal regime, and can be characterized by the aspect
ratio, which is the bubble width divided by bubble height.
Various experimental correlations exist for the prediction
of aspect ratio in different liquids. For example, the
Vakhrushev-Efremov (1970) correlation gives the relation
between aspect ratio and the Tadaki number (Ta) in the
following form:
b/h= {0.81+0.2tanh[1.8(0.4 -log,, Ta]} (28)
where Ta is given by

Ta=g4 P1 du (29)

In order to compare it with the simulation results, the
terminal velocities measured from experiment (Zenit and
Magnaudet, 2008) are used in Equation (33) and (34) to
calculate the aspect ratio. The correlation prediction for the

7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30-June 4, 2010

aspect ratio is plotted as a function of the Eo number rather
than the Ta number for better comparison, as seen in Figure
5(b). In general, the simulation, experiment, and
correlation give similar bubble shapes, although the
simulated aspect ratio is slightly smaller than the
experimental results and the correlation. Particularly, the
agreement between the simulation and experiment is good
for small Eo numbers, but the difference grows as Eo
increases. For the largest bubble with Eo=5.6, the
simulation under-predicts the bubble aspect ratio by 20%.
The simulated bubble shapes for the five bubble sizes are
shown in Figure 5(c). It should be noted that the bubble
dynamics in a low-viscosity liquid is in fact an extremely
complicated phenomenon. For example, both the bubble
surface condition and the initial bubble shape are found to
have significant effect on the bubble behavior. While a
bubble in a clean liquid and with a large initial shape
disturbance will have a higher velocity and larger aspect
ratio, the bubble of same size in the same liquid with either
contaminated bubble surface or small initial distortion
usually has a lower velocity and more spherical shape (Wu
and Gharib, 2002). In the current simulations, the bubbles
rise with initially spherical shapes, which might be the
reason for their smaller aspect ratio compared to the

Equivalent diameter (mm)

Figure 6: Terminal velocity and shape of 0.5-5 mm air
bubble in water. Simulations are performed using the
two-component MRT-LBM model. The solid curves are
experimentally obtained bubble velocity in pure and
contaminated water (Clift, et al, 1978).

Figure 6 shows the simulated shape and terminal velocity
of 0.5-5 mm air bubbles in water. The corresponding Eo
numbers are in the range between 0.034 and 3.37, while the
resulted Re numbers vary from 35 to 1030. In literature,
there is considerable scattering of the experimentally data
for the terminal velocity for millimeter bubbles in water,
mostly due to the different degree of contamination on the
gas-liquid interface. As a comparison, the experimentally
observed bubble velocities in pure and contaminated water
are plotted in Figure 6 (Clift, et al, 1978). The MRT-LBM
simulation results are found to fall between the upper and
lower limits of the experimentally measured velocities. For
bubbles larger than 3mm, the simulation results are very
close to the pure water condition. The bubble near 2mm are
usually more difficult to simulate accurately, because they
are found to experience a transition in their flow pattern,
which changes from a rectilinear path to a zig-zag or spiral

Paper No

path. The transition is a complex phenomenon that is
affected by many factors, and to capture the condition for
the onset of the oscillation still requires further
investigation. In fact, no oscillation behavior is observed
for the bubbles in the current 3D simulations, although a
high Re number of 103 is reached. The aspect ratio of the
bubble is also found to be smaller than the experimental
results (Veldhuis et al, 2008), as in the previous case of
bubbles in silicone oil. The smaller aspect ratio may help to
explain the absence of oscillation in the current simulation,
as indicated by a recent experimental study (Zenit and
Magnaudet, 2008).

Figure 7: Grid arrangement for a 2-D ellipsoidal bubble
with Eo=49.82 and Mo=137.42. (a) Entire computation
domain; (b) Close-up of bubble region

When the bubble experiences large deformation, the
LBM-AMR approach is employed to balance the need for
fine grid resolution and affordable computational resource.
Such simulation usually employs 5-7 refinement levels,
and a typical grid set-up is shown in Figure 7 for a 2D
simulation. In the figure, each block represents 8*8 grids.
In this case, an oblate ellipsoidal cap bubble with Eo=49.82
and Mo=137.42 is simulated using the single component
interaction potential model. Such bubbles usually appear in
liquids with high viscosities, and have a dented bottom
region. When 7 levels of grid are used, the size of the finest
mesh located near the bubble surface is only 1/64' of that
of the largest grids in the far field. If the entire domain is
computed with the finest grids, the number of total grids
would be 35 times large than that in the AMR scheme,
which highlights the advantage of AMR in obtaining
sufficient resolution while saving tremendous
computational resources. The thin tail region of a skirted
bubble can also be captured by the AMR, as shown in
Figure 9. The flow structure with two pairs of
counter-rotating vortices in the bubble wake is also
correctly presented in the simulation.

7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30-June 4, 2010

(a) (b)
Figure 8: Skirted bubble. Eo=86.8, Mo=54.3. (a) grid
configuration (b) streamlines of relative velocity

Figure 9 shows a series of snapshots of simulation results
for an oblate ellipsoidal bubble in 3D. The simulation is
performed with the multi-component interaction potential
model using the D3Q19 lattice, and 4 levels of mesh
refinement are employed. The bubble initially has a
spherical shape, and gradually evolves into an oblate
ellipsoidal shape as it rises under the buoyancy force.
Meanwhile, the grid resolution is adjusted accordingly to
provide finest resolution near the bubble surface, as
shown in the figure. The bubble is in the oblate ellipsoidal
regime, with Eo=9.81 and Mo=1.19*10-4.

Figure 10 shows the streamlines near the ellipsoidal
bubbles calculated from the relative velocity of the liquid
phase. For the bubble with a higher Eo of 9.81, the bubble
shape has a larger deformation and higher aspect ratio. A
pair of symmetric vortices is seen in the streamlines
behind the bubble, and the Reynolds number is about 36.
In contrast, for the bubble with a lower Eo number of 2.45,
the shape is more spherical, and no vortices are found in
the bubble wake.

t=0 t=100


t=300 t=400 t=500
Figure 9: Time sequence of the rise of a 3D bubble with
Eo=9.81 and Mo=1.19*104.

Paper No

(a) (b)
Figure 10: streamlines near bubble surface (a) Eo=9.81
and Mo=1.19*10-4, (b) Eo=2.45, Mo=2.98*10-5


In this study, an advanced two-phase LBM based on the
interaction potential model is developed for the simulation
of bubble dynamics. The new LBM approach employs the
same interaction force as that in the classic Shan-Chen
model. However, the governing LB equation is expressed
in a more general form so that the interaction force is easily
incorporated. The MRT scheme is used in the collision step,
and relaxation time scales for different physical quantities
are decoupled from each other to obtain enhanced
numerical stability. For bubbles with large deformation, the
AMR technique is developed to improve the simulation
accuracy near the bubble surface while reduce the total
computational cost. To validate the newly developed MRT
and AMR techniques, 2-D and 3-D simulations are
performed to investigate the buoyant rise of gas bubbles in
viscous liquids. The simulated bubble shape and velocity
are in good agreement with experimental observations and
empirical correlations in the literature. It is shown that the
MRT scheme can successfully simulate millimeter bubbles
in low-viscosity liquids such as silicone oil DMST-05
(Mo-107) and water (Mo-10-O) due to the improved
numerical stability. For bubbles with large deformation, the
AMR scheme is able to accurately predict the bubble shape
and the flow field near the bubble surface. The improved
performance of the new LBM shows its great potential in
simulating a broad range of two-phase flow problems that
are difficult for the traditional multiphase LBM techniques.


This work was supported in part by an allocation of
computing time from the Ohio Supercomputer Center.


Ansumali, S. and Karlin, I. Single relaxation time model
for entropic lattice Boltzmann methods, Physical Review E,
65, 056312 (2002)

Chen, H., Filippova, O., Hoch, J., Molvig, K., Shock, R.,
Teixeira, C., and Zhang, R. Grid refinement in lattice
Boltzmann methods based on volumetric formulation.
Physica A. 362 158-167. (2006)

Clift, R., Grace, J.R., and Weber, M.E. Bubbles, drops and
particles, Academic Press, New York. (1978)

7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30-June 4, 2010

Fan, L.-S. Gas-liquid-solid fluidization engineering,
Butterworth Publishers, Boston, (1989)

Fan, L.-S. and Tsuchiya, K. Bubble Wake Dynamics in
Liquids and Liquid-Solid Suspensions,
Butterworth-Heinemann. Boston. (1990)

Frank, X., Funfschilling, D., Midoux, N. and Li, H.Z.
Bubbles in a viscous liquid: Lattice Boltzmann simulation
and experimental validation, Journal of Fluid Mechanics,
546, 113 (2006)

Guo, Z., Zheng, C., and Shi, B. Discrete lattice effects on
the forcing term in the lattice Boltzmann method, Physical
Review E, 65, 046308 (2002)

Gupta, A. and Kumar, R. Lattice Boltzmann simulation to
study multiple bubble dynamics, International Journal of
Heat and Mass Transfer, 51, 5192 (2008)

He, X., Shan X., and Doolen, G.D. Discrete Boltzmann
equation model for nonideal gases, Physical Review E, 57,
13 (1998)

Inamuro, T., Ogata, T., Tajima, S., and Konishi, N. A lattice
Boltzmann method for incompressible two-phase flows
with large density differences, Journal of Computational
Physics, 198, 628 (2004)

Kurtoglu, 1.0. and Lin, C.L. Lattice Boltzmann Study of
bubble dynamics, Numerical Heat transfer B, 50, 333

Kuzmin, A., Mohamad, A. A. and Succi, S.
Multi-relaxation time lattice Boltzmann model for
multiphase flows, International Journal of Modern Physics
C, 19, 875, (2008)

MacNeice, P., Olson, K. M., Mobarry, C., deFainchtein, R.,
and Packer, C. PARAMESH : A parallel adaptive mesh
refinement community toolkit. Computer Physics
Communications, 126, 330-354. (2000)

Mukherjee, S. and Abraham, J. A pressure-evolution-based
multi-relaxation-time high-density-ratio two-phase
lattice-Boltzmann model, Computers and Fluids, 36, 1149

Premnath, K.N. and Abraham, J. Three-dimensional
multi-relaxation time (MRT) lattice Boltzmann models for
multiphase flow, Journal of Computational Physics, 224,
539 (2007)

Rohde, M., Kandhai, D., Derksen, J.J., and van der Akker,
H.E.A A generic mass-conservative local grid refinement
technique for lattice Boltzmann scheme. International
journal for numerical methods in fluids. 51 439-468.

Sankaranarayanan, K., Shan, X., Kevrekidis, I.G., and
Sundaresan, S. Analysis of drag and virtual mass forces in
bubbly suspensions using an implicit formulation of the
lattice Boltzmann method, Journal of Fluid Mechanics, 452,

Paper No 7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30-June 4, 2010


Sbragaglia, M., Benzi, R., Biferale, L., Succi, S., Sugiyama,
K. and Toschi, F Generalized lattice Boltzmann method
with multirange pseudopotential, Physical Review E, 75,
0026702 (2007)

Shan, X. and Chen, H. Lattice Boltzmann model for
simulating flows with multiple phase and components,
Physical Review E 47 (1993) 1815-1819

Shan, X. Pressure tensor calculation in a class of non-ideal
gas lattice Boltzmann models. Physical Review E, 77:
066702 (2008)

Theodoropoulos, C., Sankaranarayanan, K., Sundaresan, S.
and Kevrekidis, I.G. Coarse bifurcation studies of bubble
flow lattice Boltzmann simulation,. Chemical Engineering
Science, 59, 2357 (2004)

Tolke, J., Freudiger, S., and Krafczyk, M. An adaptive
scheme using hierarchical grids for lattice Boltzmann
multiphase flow simulations. Computer and Fluids, 35, 820

Vakhrushev, I.A. and Efremov, G.I. Chem. Techn. Fuel Oils
(USSR) 5/6, 376 (1970)

Veldhuis, C., Biesheuvel, A., and van Wijngaarden, L.
Shape-oscillations on bubbles rising in clean and in tap
water. Physics of Fluids 20, 040705 (2008)

Wu, M. and Gharib, M. Experimental studies on the shape
and path of small air bubbles rising in clean water, Physics
of fluids, 14, L49 (2002)

Zenit, R. and Magnaudet, J. Path instability of rising
spheroidal air bubbles: A shape-controlled process, Physics
of Fluids, 20, 061702 (2008)

University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - Version 2.9.7 - mvs