7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Interface Resolving Simulation of BubbleWall Collision Dynamics
T. Omori* H. Kayama* Z. TukoviCd and T. Kajishima*
Department of Mechanical Engineering, Osaka University, 21 Yamadaoka, Suita, Osaka, Japan
t Faculty of Mechanical Engineering and Naval Architecture, Univ. of Zagreb, Ivana LuCida 5, Zagreb, Croatia
t.omori@mech.eng.osakau.ac.jp
Keywords: bubblewall collision, liquid film, interfacetracking method, immersedboundary method
Abstract
The bubblewall collision dynamics were studied by means of numerical simulation, where the maximum Reynolds
number and the maximum Weber number was 279 and 0.503 respectively. The present simulations are restricted
to twodimensional circular bubbles, but the elucidated principle collision mechanism should also apply to three
dimensional cases. It was revealed that the bubble bounce from the wall is maintained by the dimple formation on the
gasliquid interface, and the main dissipation of energy during the collision process is caused by the enhanced liquid
motion near the gasliquid interface and the dissipation accompanying the drainage of the liquid film or the vortex
formation during the receding motion plays a minor role.
Nomenclature
Roman symbols
E energy in unit depth (N)
g gravitational acceleration (m s 2)
(= 9.81, unless otherwise specified)
1,1 interface mesh spacing in the radial direction (m)
1, in the circumferential direction (in)
n unit normal (in)
p, modified pressure (N in2)
Q flow rate in unit depth (m2 /s)
r point position (m)
R bubble radius (m)
Up solid wall velocity in the noninertial frame (rn/s)
w weighting factor ()
Greek symbols
At time step width (s)
v kinematic viscosity of gas/liquid (mi2 s 1)
p density of gas/liquid (kg in3)
c interfacial tension (N/m)
Subscripts
c interface mesh control points
i interface mesh vertex points
nb neighbor cells
P current cell
I fluid in the lower region
u fluid in the upper region
L, G liquid/gas
Superscipts
m iteration counter in a PISO loop
n current time step counter
Introduction
To provide necessary boundary conditions for the pre
diction of dispersed multiphase flows, it is of great im
portance that the bubblewall collision dynamics is prop
erly modeled. The importance may be illustrated by
some experimental evidence that in gasliquid upward
circular duct flows the bubble volume fraction showed a
sharp peak at a distance from the duct wall. In the ex
periments by Valukina et al. (1979), where the bubble
diameters were between 0.5 and 1.0 mm and the mean
gas volume fraction was in the range from 0.02 to 0.1,
the peak of the bubble volume fraction was localized
around 1.7 bubble radius from the wall, which is to be
the consequence of bouncing rather than sliding motions
of bubbles, as pointed by Tsao & Koch (1997).
The bubblewall collision has also been studied to ob
tain fundamental physical description of bubbleparticle
or bubblebubble interaction problems, which are typi
cally encountered in floatation and waste water purifica
tion processes. Fujasova et al. (2007), Legendre et al.
(2005) and Tsao & Koch (1997), among others, con
ducted experimental studies on the hydrodynamics of
collision of a bubble with a horizontal wall in a liquid.
Probably due to the difference in the fluid properties,
however, there is a scatter in their results and it is hard to
deduce a unified conclusion. The limitation in the exper
imental techniques also precludes the understanding of
the bubble motion dynamics in relation to the hydrody
namics of the surrounding fluid, which is especially im
portant in the attempt to make twoway coupling mod
els.
The objective of the present work is to investigate
the bubble motion during the collision with a horizon
tal rigid wall in a welldefined condition by means of
numerical simulation. A numerical approach has a great
advantage to relate the path, the deformation and the rise
velocity of the bubble with the energy balance and the
flow field which the bubble itself is a part of. The present
simulations are restricted to twodimensional circular
bubbles, but the elucidated principle collision mecha
nism should also apply to threedimensional cases.
Numerical Methods
The NavierStokes equation and the mass conservation
equation for incompressible fluid were discretized in the
framework of Finite Volume Method. The discretization
was carried out by the central difference scheme and the
three time level scheme for space and time respectively,
which have both secondorder accuracy. To capture the
motion and deformation of a bubble, the gasliquid in
terface was expressed by the boundary of two solution
adaptive mesh domains that were allotted either for the
gas or the liquid phase, based on the interfacetracking
method proposed by Muzaferija & PeriC (1997). To
avoid mesh distortion during the simulation, the bubble
centered noninertial frame was adopted and the rigid
wall was treated as a quasi boundary condition by an
immersedboundary approach.
GasLiquid Interface
At the gasliquid interface the dynamic and the kine
matic boundary conditions are made strongly coupled by
a PISOtype iterative partitioned method. The dynamic
condition requires the stress (pressure, viscous stress and
surface tension) balance normal to the interface and the
kinematic condition has to be satisfied by moving the in
terface to assure that there is no volume flux across the
interface.
The implementation of the kinematic boundary con
dition is illustrated in Fig. 1. The possible volume flux
VI which arises in the course of a PISO loop to correct
the velocity and to determine the pressure in the liquid
phase must be compensated by the volume swept by the
interface mesh. This swept volume is written as:
2V1 At + .1i
V 1 3 (1)
3
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Area
Previous Interface
Length
h' .
New Interface
Mesh vertex points
Mesh control points
Figure 1: Interface mesh motion for the fulfillment of
the kinematic boundary condition
for the three time level scheme by virtue of SCL require
ment (see e.g. Ferziger & PeriC (2002)). The positions
of the mesh control points, which are initially at the cell
face centers, are updated as:
n+l n+ h/ n
i'6 Cr6n Cn,
where
and the unit vector normal to the cell face is used as n,.
The interface mesh positions at the new time step are
then given by the following expression:
n+ Ni (p r)
S r+ Ni n1 n'
where p is a point on the plane which is defined by the
enclosing control points
Ec w2r +1
c W2
and Ni is the unit normal vector to the plane obtained
by the leastsquares method
Nwr ) 1 (6)
The unit normal n, at a vortex point is defined by the
surrounding cell face normal vectors by the method pro
posed by Max (1999). The mesh points in the whole do
main are updated only after the convergence of the PISO
loop for efficiency, by solving the mesh vertex motion
equation with the interface vertex motion as a boundary
condition.
SolidLiquid Interface
The noslip condition at the solidliquid interface is im
posed on a grid scale by an immersed boundary ap
proach. The NavierStokes equaiton for incompressible
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
fluid is replaced by the following expression:
Preliminary Simulations
(1 7) u.+VVV2u
+ Fp,
1 
VPm ab
P
where ab represents the acceleration of the bubble cen
ter of gravity forming the noninertial simulation frame
seen from the bubble and 7 is the volume fraction of
solid phase in each mesh cell. The velocity field where
7 1 has velocity of the solid phase by the effect of the
additional term Fp written as:
3U"+1 4u" + u 1
2At (8)
for the three time level scheme adopted for the time dis
cretization.
The coupling of the NavierStokes equation and the
mass conservation equation is accomplished by the pro
cedure based on the PISO method. The NavierStokes
equation (Eq. 7) can be written in the semidiscretized
form as the following:
up = 1 H(u)
Ap I
1 )V
(1 7)1Pm ,
p )
where
H(u) = AAnbUnb + Q
and Q represents all the source terms other than the pres
sure term.
First of all, this momentum equation is solved for the
guessed velocity using the pressure from the previous
iteration m in the PISO loop
p =1 H(u) (1
AP I
1 V/",
P)V::~
The corrected pressure is obtained by solving the follow
ing Poisson equation,
V. tp,
where
Up= (1 H( +) U;+1 (13)
and consequently the velocity is updated as the follow
ing:
1
um] up 1 V/.::1 (14)
App
The above procedure, which involves the interface track
ing method explained in the previous subsection, is re
peated until the velocity divergence vanishes.
Before discussing the bubblewall collision dynamics,
the present interfacetracking method was examined for
) two fundamental cases. The examination on the two
dimensional sloshing problem showed that the gravity
driven gasliquid interface oscillation can be correctly
predicted. The simulations were also performed for a
freely rising air bubble in an initially quiescent water,
which is directly related to the present main simulation
cases and based on whose result the sufficient mesh res
olution for further simulations was determined. The res
olution issue for the solidliquid interface will be ad
dressed at the end of the present work.
TwoDimensional Sloshing
The simulations were performed for the initialvalue
problem associated with the small amplitude waves on
the interface between two stably superposed viscous flu
ids and the results were compared to the analytic so
lutions by Prosperetti (1981). It was restricted to the
cases where the two fluids have equal kinematic viscosi
ties and the analytic solutions give explicit expression
for the interface position. The parameters for the com
putations are summarized in Tab. 1. The airwater in
terface at 28 C is considered in Case I, except that the
viscosity of the water has about ten times the value of
the real value, because of the condition that the both flu
ids must have the same kinematic viscosity and Case II
deals with the case where the natural frequency is lower
and the motion is more dissipative than Case I.
As the initial interface shape, the following expression
was assigned
y(x) = 1 + ao sin{7r(0.5 x)}
(15)
and the amplitude ao was set to 0.01m. Fig. 2 shows
the initial mesh configuration. The domain size is 1 m
and 2 m horizontally and vertically, respectively. Tan
gential to the interface it has 40 cells of the same size,
and normal to the interface the mesh was divided into 20
cells in each region and clustered with the expansion ra
tio 1.088 and the finest mesh size of 0.0199 m. At the up
per boundary, constant total pressure was assumed and
for the other boundaries slipwall condition was applied.
The time history of the interface location at the left
end, namely y(0) was compared to the analytic solution
and is shown in Fig. 3. In both cases the discrepancy
between the present results and the analytic solutions are
almost indiscernible, where the discrepancy is 0.3% of
the analytical value at t = 20 s for Case I.
Freely Rising Air Bubble in Water
A freely rising clean air bubble of diameter 1.0mm in
water was simulated on three different spacial resolu
tions. The material properties were those of air/water
\7 1 Y V/", t
( APP
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Figure 2: Mesh and coordinate arrangement for 2D
sloshing problem
Table 1: Conditions for the computation of 2D sloshing
Variable Case I Case II
Pu 1.173 10.00
pi 996.31 1000.0
v 1.58 105 1.00 103
g1 9.81 1.00
a 0.072 0.10
at 28 C (Tab. 2) and the specification of the numerical
mesh is summarized in Tab. 3. The temporal resolution
was set to 0.01 ms for all the cases. The time variation
of the bubble rise velocity shown in Fig. 4 illustrates
negligible influence of the mesh resolution next to the
interface. Even reducing the overall resolution to the
half resulted in a small difference of the terminal rise
velocity (0.86 %). It can be concluded from these re
sults that Mesh I provides sufficient spacial resolution
for further simulations. The obtained terminal veloc
ity 0.25 m/s has the same order of magnitude to the ex
perimental value 0.19 m/s at 20 C by Gaudin (1957) or
0.26 m/s by Moore (1965).
Table 2: Adopted material properties for air/water
Variable Air Water
p 1.173 996.3
v 1.58 105 8.36 107
S 0.072





. ...........
 111111111
=4444
 +H44
 HH
1015 
1 01
1 005
1
0995
0 99
0
1 015
1 01
0 5 10 t s]
Figure 3: Interface oscillation at x
Lower: Case II)
15
0 (Upper: Case I,
BubbleWall Collision
The collision between a horizontal wall and a clean air
bubble of diameter 1.0mm and primarily 1.5 mm was
studied. The adopted material properties of fluids are
the same as given in Tab. 2. The wall was assumed to be
clean and chemically neutral to the air/water interface.
The initial distance between the wall and a bubble center
was set to 4.0mm for every case. The simulations were
performed on Mesh I from the previous section, which
is shown in Fig. 5.
Fig. 6 shows the variation of the vertical coordinate
and the vertical translation velocity of the bubble center.
Table 3: Grid resolutions for the rising air bubble
Variable Mesh I Mesh II Mesh III
1"1/R 0.5 102 1.0 102 1.0 102
It/(27r) 160 160 80
Total cell # 14080 14080 3520
5 10 [] 15
t [S]
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
025
02
015
0 01
005
Mesh I
Mesh II
Mesh III
0
0 005 01 t[s] 015 02
Figure 4: Rise velocity of a 2D bubble
Table 4: Dimensionless groups based on the maximum
approaching velocity to the wall
Variable D 1.0 nmm D 1.5 mm
Rem,, 191 279
Wem,, 0.354 0.503
The maximum Re and We based on the maximum rise
velocity of the bubble and the material properties of the
water are listed in Tab. 4. It can be seen from Fig. 6
that both bubbles display a rebound motion after the first
and the second impact, which supports the experimental
observation by Tsao & Koch (1997) that bubbles with
We > 0.3 underwent a measurable rebound. Because
both bubbles showed the same nature in principle, the
following discussion will be focused on the bubble of
diameter 1.5 mm.
Figs. 14 and 15 illustrate the transition of the bubble
shape and the velocity field in the gas/liquid phases dur
ing and right after the first impact to the wall. The veloc
ity vectors were measured from the inertial frame. The
selected time instances in Fig. 14 are some specific in
stances for the bubblewall collision dynamics and ex
plained in Tab. 5. At the contact time t 31.0ms,
where the distance between the bubble center and the
wall is D/2, there is still ,igmilk.iiii portion of fluid be
tween the wall and the bubble due to the deformation
of the interface. If the bubble gets closer to the wall,
the dimple of the interface is formed, as was also ob
served experimentally by Platikanov (1964). The liquid
film did not break even when the distance between the
bubble and the wall is minimum. This issue will be dis
cussed later in detail. At t 37.2 ms, where the surface
energy turns back into the maximum rebound velocity,
the distance between the bubble center and the wall is
approximately D/2 but the bubble shape and the liquid
Figure 5: Grid for the bubblewall collision simulations
(The cyan line indicates the gasliquid interface.)
Table 5: Description of the snapshots shown in Fig. 14
Time [ms] Explanation
1.00
25.6 Maximum upward translation velocity
31.0 Contact time
34.2 Minimum distance to the wall
37.2 Maximum downward translation velocity
film thickness is totally different from those at the con
tact time on approach.
The oscillations in the bubble translation velocity dur
ing the rebound (see Fig. 6), which is also visible in the
experiments by Fujasova et al. (2007), can be explained
by the oscillatory shape deformation of the bubble as
clearly seen in Fig. 15. The time instance t 43.0 ms
corresponds to the local maximum of the upward trans
lation velocity of the bubble.
The energy variation associated with the bubble mo
tion during the collision process is depicted in Fig. 7.
The energies considered here are summarized in Tab. 6.
As the surface energy, variation from the initial state is
considered. It should be noted that the kinetic energy of
the surrounding water can be written as:
i 1
f pu udV = CMpLVbU,,
2 JL 2
(16)
using the added mass coefficient CM1, the volume Vb
1Fig. 8 shows the time variation of the added mass coefficient of the
bubble. At the contact time CM is increased to 2.56 from 1 which
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
0 002 004 0 06 t [s] 0 08
005
0 002 004 0 06 t [s] 0 08
01 012
01 012
Figure 6: Bubble center coordinate and velocity in the
collision process
and the translation velocity Ub of the bubble. This en
ergy is almost identical to Ed because of the high den
sity ratio in the present study. In Fig. 7, intimate energy
exchange between the kinetic energy and the surface en
ergy is evident during the bounce and the oscillatory mo
tion (Fig. 15).
It is interesting to note that the total energy shows an
eminent decrease during the bounce period. For the first
bounce, the decrease from the one at the contact time to
the other at the time when the bubble and the wall take
Table 6: Definition of energies considered in Fig. 7
Variable Description Definition
Ed Kinetic E. 2 ]L,G plu udV
E, Surface E. o(S So)
Ep Gravit. Potential Yb g JG(pL cG)dV
ET Total Ed + Es + Ep
is the theoretical value for the cylinder.
4e05
3e05
2e05
1e05
S0
LU
1e05
2e05
3e05
4e05
5e05
6e05
Kinetic
Elastic (Surface) Potential
Gravitational Potential
Total
V/ A> _
0 002 004 006 t [s] 008
4
3.5
3
2.5 /
2 7
01 012
0 0.02 0.04 0.06 t [s]0.08 0.1
S0.2
0.15
0.1
0.05
og
0.05
0.1
0.15
0.12
Figure 7: Energy variation during the collision process
(The bottom figure replots Fig. 6 for convenience.)
the maximum distance, divided by the initially stored
total energy (Eo E,) is 0.53, which is close to the
value 0.59 claimed by Tsao & Koch (1997) in their ex
periments. The dissipation due to the film drainage as
shown in Fig. 9 was found to be way small to explain
the observed decrease in the total energy and it should
be concluded that the loss comes from the dissipation in
the global liquid motion. It also deserves attention that
the vortex formation (Fig. 10) during the receding phase
with the shape oscillation (see Figs. 15 and 7 bottom)
does not contribute much to the energy loss.
Flow in the liquid film
To discuss further the rebound phenomena in the bubble
wall collision, the present method to express solid
boundary should be assessed. When the bubble is at the
minimum distance to the wall (t 34.2 ms), the mesh
takes the form shown in Fig. 11. The simulations were
performed on the channel which represents the liquid
film at t = 34.2 ms with the present immersedboundary
35
35
25
2
15
0005 001 0015 002 0025 003 0035 004
t [s]
Figure 8: Variation of the added mass coefficient of the
bubble during the first bounce period
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
'i ,
! i, ,
I,.','." ''I,'."",.'.
Figure 10: Vortex formation in the vicinity of the inter
face during the receding motion (Upper: t = 43.0 ms,
Lower: 45.5 ms)
Figure 9: Drainage of the water in the dimple at t
34.2ms
method and the bodyfitted grid (Fig. 12). The bound
aries were assumed fixed and the gasliquid interface
was replaced by the freeslip wall. The flow was driven
by the fixed pressure difference 200Pa, which is the
value measured in the bubblewall collision simulation,
from the left to the right boundary and on these bound
aries the zerogradient boundary condition was applied
to the velocities. The number of cells for the grids and
the calculated flow rate in the steady state are summa
rized in Tab. 7. The IB I corresponds to the spacial res
olution in the present bubblewall collision simulation
(see Fig. 11). It can be seen that the present method un
derestimates the flow rate on a coarse mesh. It must be
emphasized, however, that even if the velocity in the liq
uid film is tripled, the dissipation in the liquid film is far
too small to explain the steep decrease in the total en
ergy and also that the drainage speed is far too low to set
off the film break and it does not disprove the conclusion
that the rebound should occur. Fig. 13 demonstrates that
the present method provides a correct solution with an
enhanced resolution.
Conclusions
The numerical simulations on the twodimensional
bubblewall collision problem revealed that the dimple
Figure 11: Grid near the interface at t 34.2 ms (Red:
solidliquid interface, Blue: gasliquid interface)
formation maintains the bubble bounce by preventing
the attachment of the gasliquid interface to the wall
and the concentrated loss of energy during the contact
period, which leads to the perish of the bubble mo
tion, is due to the enhanced liquid motion near the gas
liquid interface. The present numerical methods pro
vided good results for the benchmark problems and it
was also shown that the possible numerical errors do
not affect the above conclusions. Although the straight
forward comparison to the existing experimental results
is not possible, judged from facts and figures in the
present work, it can be considered that the present simu
lations on the bubblewall dynamics give reasonable re
sults even quantitatively.
Figure 12: Grid arrangement for the liquid film channel
(Upper: bodyfitted, Lower: immersedboundary)
Table 7: Summary for the liquid film channel simulation
(BF: bodyfitted, IB: immersedboundary)
BF IB I IB II IB III
Mesh 700 x 40 25 x 7 25 x 14 200 x 56
106 Q 2.84 1.01 2.09 2.65
References
Valukina N.V., Koz'menko B.K. and Kashinskii O.K.,
Characteristics of a flow of monodisperse gasliquid
mixture in a vertical tube, J. Engineering Phys. and Ther
mophys., 364, pp. 462465, 1979
Tsao H.K. and Koch D.L., Observations of high
Reynolds number bubbles interacting with a rigid wall,
Phys. Fluids, 91, pp. 4456, 1997
Fujasova M., Vejrazka J., Ruzicka M.C. and Drahos
J., Experimental study of bubblewall collision, 6th
Int. Conf Multiphase Flow, S1 WedC_38, 2007
Legendre D., Daniel C. and Guiraud P., Experimental
study of a drop bouncing on a wall in a liquid, Phys. Flu
ids, 17, 097105, 2005
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Muzaferija S. and PeriC M., Computation of Free
Surface Flows Using Finite Volume Method and Moving
Grids, Numerical Heat Transfer B, 324, pp. 369384,
1997
Ferziger J.H. and Perid M., Computational methods for
fluid dynamics, Springer, 2002
Max N., Weights for computing vertex normals from
facet normals, J. Graphics Tools, 42, pp.16, 1999
Prosperetti A., Motion of two superposed viscous fluids,
Phys. Fluids, 247, pp. 12171223, 1981
Gaudin A.M., Flotation, McGrawHill, 1957
Moore D.W., The velocity rise of distorted gas bubbles
in a liquid of small viscosity, J. Fluid Mechanics, 234,
pp. 749766
Platikanov D., Experimental investigation on the 'Dim
pling' of thin liquid film, J. Physical Chemistry, 68, pp.
36193624, 1964
004
003
E
E
S002
001
0 001 002 003 004 005 006 007 008
U [m/s]
Figure 13: Velocity profile in the liquid film channel
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Figure 15: Interface shape and velocity field in and out
side the bubble after the first rebound (From the above,
t = 40.0, 41.0, 42.0 and 43.0ms. The color indicates
the magnitude of the velocity: Red, large; Blue, small.
The apparent discontinuity inside the bubble is due to
the way of drawing and the results do not show any dis
continuities. The red box indicates the solid wall.)
Figure 14: Interface shape and velocity field in and out
side the bubble during the first rebound (From the above,
t = 1.00, 25.6, 31.0, 34.2 and 37.2 ms. For further de
scription on the figures, see the caption of Fig. 15)
