Paper No 7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
A Novel Lattice Boltzmann Method for Direct Numerical Simulation of Bubble Dynamics
Zhao Yu and LS Fan*
William G. Lowrie Department of Chemical and Biomolecular Engineering
The Ohio State University
140 West 19th Avenue, Columbus, Ohio 43210, USA
yu.284@osu.edu and fan.1 @osu.edu
Keywords: lattice Boltzmann method, bubble, adaptive mesh refinement
Abstract
The traditional lattice Boltzmann method (LBM) for twophase flow simulation is often hindered by numerical instability at
low viscosities and insufficient resolution at the gasliquid 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 multirelaxationtime (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
gridindependent 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 highReynolds 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.
Introduction
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 gasliquid
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 gasliquid 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 gasliquid 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 twophase 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 multiplerelaxationtime (MRT)
scheme, which employs multiple relaxation parameters in
the collision step, instead of the singlerelaxationtime BGK
collision in traditional LBM. A number of the twophase
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 pressureevolution model
(Mukherjee and Abraham, 2007). The interaction potential
model for nonideal 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):
34
Mo =gp (1)
ubd
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
millimetersize air bubbles rising in water usually have a
Reynolds number on the order of 102 103. Unfortunately,
with its instability for lowviscosity 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 (Mo104) (Theodoropoulos et
al, 2004) and multiple bubble dynamics studies
(106
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 1010, and
Reynolds number up to 400. However, even that does not
match the conditions of the airwater system exactly.
Therefore, an apparent gap exists between the capability of
current twophase LBM models and the actual flow
conditions in many lowviscosity 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 30June 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 35 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 nonuniform 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
highvorticity 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 nonuniform grid LBM
algorithm to twophase flow problems. There seems to be
only one publication concerning the multiphase LBMAMR
(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 twophase 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 multirange 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
deformation.
Nomenclature
lattice velocity
bubble diameter
Edtvds number
particle distribution function
force
gravitational constant
interaction strength
Morton number
Reynolds number
velocity
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
Superscripts
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 multirelaxationtime 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 30June 4, 2010
Equation (3), Q, is the collision term that describes the
change of particle distribution function due to particle
collisions.
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 (fff) (4)
As the result, the LBE in the general form can be expressed
as:
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
sections.
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 gasliquid twophase flow, the interaction potential
model, also called the ShanChen model, is one of the most
widely used LBM models (Shan and Chen, 1993). It
employs a meanfiled 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 nonideal
fluid, and the multicomponent model for a mixture of two
different fluids.
In the single component model, the interaction force is
given as the summation of the pairwise interactions
Paper No
between particles at a given lattice site and those at
neighboring sites, and can be written as
q
F(x) = G y(x)Z wi(x + ci)ci (7)
i=1
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 nonideal 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
gasliquid 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
gasliquid 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 gasliquid system, usually the first component is
assumed to be a nonideal 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 mutlirange 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
higherorder isotropic discretization. Since it is only
involved in the formulation of the interaction force, the
mutlirange 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 30June 4, 2010
Incorporation of the force term
Incorporating the force by shifting the equilibrium velocity,
which is called the ShanChen forcing scheme, is a simple
and unique feature of the original interaction potential
model. However, although the ShanChen 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 (vu)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)
PCs
Then, following the derivation given by Premnath and
Abraham (2007), using a transformation:
fi = fi SiAt, the general LB Equation (5) is written
2
f(x+cAt,t+A t) = f,(x,t) A, +
A1A S, At
2 )
,(13)
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
LBM:
p=>/j (14)
And as the result of the transformation, the momentum of
the fluid is calculated from
pu = fc + AtF (15)
2
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 MRTLBE
After incorporating the force term, Equation (13) becomes
the governing equation that dictates the evolution of the
distribution functions in MRTLBM. 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
space.
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 offdiagonal 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 = 
Jy
ly
(j2 jy2)p
liJ y//P
In the same way, the forcing term is also transformed into
the moment space.
0
6(uxF, +uFy)
6(uxF +uyFy)
F, (20)
S =TS= F
F,
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 30June 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 MRTLBM 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 twophase flow based on the interaction force
calculated from the interaction potential, as in the classic
ShanChen 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 SlS9. 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
usage.
* 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
II I
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 30June 4, 2010
time step needs to be proportional to the mesh size.
Following the convention of LBM, the nondimensional
1
lattice speed of sound is kept to be c throughout
the entire domain.
6
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
computation.
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
requires
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.
dt
0.5dt
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
substeps (16) is in accordance with those in Figure 1.
The newly developed LBMAMR approach is based on the
multiblock 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 nonequilibrium part of the distribution at
coarse/fine grid interface is not necessary. The interblock
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
c~3
Paper No
four smaller cells during the explosion step (step 23 in
Figure 1), and its particle population is either evenly
distributed or linearly interpolated in these new cells. In the
coalescence step (step 67 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 "collisionpropagation" 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
"collisionpropagation" 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 substep 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
precollision status before step 1.
The generation and refinement of the mesh is performed
with the opensource code PARAMESH (MacNeice et al,
2000), which use a tree data structure to manage the
hierarchy of blockstructured grid with different refinement
levels. It has been integrated with the nonuniform grid
LBM code developed by the authors.
The computation procedure of multiphase LBM on a
nonuniform 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 coarsefine 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 30June 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 twophase 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.
0.09
0.08
0.07
0.06
0.05
0.04
0.03
0.02
0.01
000
0.001
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
multirange 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
*BGK
*MRT
SMRT mulirange
.
A
A E
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 multirange 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 multirange 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 MRTmultirange 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 30June 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 multirange
potential model can effectively reduce the spurious
velocity.
Due to the improved numerical stability at small viscosities,
the newly developed MRTLBM is able to successfully
simulate millimetresized bubbles in a low viscosity liquid,
which are characterized by small Morton number and high
Reynolds number, and are therefore difficult for traditional
twophase LBM. The 3D bubble simulations are carried
out using the twocomponent 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 MRTLBM model.
The first set of simulations are conducted at Mo=6.2*107,
corresponding to millimeter gas bubbles in silicone oil
DMST05 (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
a
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 gasliquid 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
Eo
(a)
6 7 8
2 4 6
Eo
(b)
1.5mm 2.0mm 2.65mm
3.mm
3.0mm
3.5mm
3.5mm
(c)
Figure 5: (a) Reynolds numbers for 1.53.5mm air
bubbles in silicone oil DMST05; (b) Aspect ratios of
1.53.5mm air bubbles in silicone oil DMST05; (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
VakhrushevEfremov (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 30June 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 underpredicts 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 lowviscosity 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
experiment.
Equivalent diameter (mm)
Figure 6: Terminal velocity and shape of 0.55 mm air
bubble in water. Simulations are performed using the
twocomponent MRTLBM 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.55 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
gasliquid interface. As a comparison, the experimentally
observed bubble velocities in pure and contaminated water
are plotted in Figure 6 (Clift, et al, 1978). The MRTLBM
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 zigzag 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 2D ellipsoidal bubble
with Eo=49.82 and Mo=137.42. (a) Entire computation
domain; (b) Closeup of bubble region
When the bubble experiences large deformation, the
LBMAMR approach is employed to balance the need for
fine grid resolution and affordable computational resource.
Such simulation usually employs 57 refinement levels,
and a typical grid setup 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
counterrotating vortices in the bubble wake is also
correctly presented in the simulation.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 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 multicomponent 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*104.
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=200
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*104, (b) Eo=2.45, Mo=2.98*105
Conclusions
In this study, an advanced twophase 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 ShanChen
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, 2D and 3D 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 lowviscosity liquids such as silicone oil DMST05
(Mo107) and water (Mo10O) 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 twophase flow problems that
are difficult for the traditional multiphase LBM techniques.
Acknowledgements
This work was supported in part by an allocation of
computing time from the Ohio Supercomputer Center.
References
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 158167. (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 30June 4, 2010
Fan, L.S. Gasliquidsolid fluidization engineering,
Butterworth Publishers, Boston, (1989)
Fan, L.S. and Tsuchiya, K. Bubble Wake Dynamics in
Liquids and LiquidSolid Suspensions,
ButterworthHeinemann. 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 twophase 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
(2006)
Kuzmin, A., Mohamad, A. A. and Succi, S.
Multirelaxation 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, 330354. (2000)
Mukherjee, S. and Abraham, J. A pressureevolutionbased
multirelaxationtime highdensityratio twophase
latticeBoltzmann model, Computers and Fluids, 36, 1149
(2007)
Premnath, K.N. and Abraham, J. Threedimensional
multirelaxation 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 massconservative local grid refinement
technique for lattice Boltzmann scheme. International
journal for numerical methods in fluids. 51 439468.
(2006)
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 30June 4, 2010
61(2002)
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) 18151819
Shan, X. Pressure tensor calculation in a class of nonideal
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
(2006)
Vakhrushev, I.A. and Efremov, G.I. Chem. Techn. Fuel Oils
(USSR) 5/6, 376 (1970)
Veldhuis, C., Biesheuvel, A., and van Wijngaarden, L.
Shapeoscillations 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 shapecontrolled process, Physics
of Fluids, 20, 061702 (2008)
