7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Boundary Element Analysis for Bubble Collapse
near Moving or Deformable Boundaries
Y. Jinbo and H. Takahira
Department of Mechanical Engineering, Osaka Prefecture University
11 Gakuencho, Nakaku, Sakai, Osaka 5998531, Japan
jinbo06@fluid.me.osakafuu.ac.jp and takahira@me.osakafuu.ac.jp
Keywords: bubble, boundary element method, Kelvin impulse, floating body, elastic tube
Abstract
The Boundary Element Method (BEM) is utilized to predict the bubble collapse near moving or deformable bound
aries. First, the BEM is applied to the growth and collapse of a bubble under a floating body. The threedimensional
deformation of the bubble, the translation and rotation of the floating body, and the motion of water surface are
taken into account in the simulation. The bubble motion is analyzed using the Kelvin impulse. It is shown that the
horizontal translational motion of the bubble is much dependent on the rotational motion of a floating body: when
the moment of inertia of the floating body is small, the largest horizontal translation is realized between the axis of
flotation and the edge of the floating body. Second, a fast simulation method is proposed for underwater explosion.
In the method, the motion of the bubble generated by the underwater explosion is modeled with a point source and a
doublet. The bubble motions, coupled with the motions of free surface and a floating body, are solved by using the
BEM. The results show that the present approximation results are in agreement with those obtained from the BEM
in which the interfacial motion of a bubble is solved directly. Finally, the bubble collapse in an elastic tube is simu
lated by taking the deformation of the tube into account. The side wall of the tube is constricted by the bubble collapse.
Introduction
As is well known, when a bubble collapses in a liquid,
extremely high pressure fields are generated around the
bubbles (Brennen, 1995). The collapse of the bubble
follows the formation of liquid jets and shock waves
leading to material damage. Recently, the phenomena
related to the bubble collapse apper in various fields
of engineering and medicine. For example, the bub
ble collapse induces the cavitation erosion of the sur
face of propellers. The bubble collapse is also impor
tant in underwater explosion under a ship: the bottom
of a floating body is destroyed due to the collapse of
bubbles (Keil, 1961, Klaseboer et al., 2005, Imakita &
Yasuda, 2007). The damage of the body is primarily
caused by liquid jets or shock waves from the bubble
formed by the explosion rather than the shock waves
by the explosion itself. Thus, understanding the bub
ble collapse under a floating body is needed to evaluate
the damage and safety of a ship. Also, in Extracorporeal
Shock Wave Lithotripter (ESWL), the shock waves gen
erated outside of human body are utilized to crush the
calculi in patients' body. When the calculi are crushed,
however, cavitation bubbles are formed near the calculi.
The collapse of cavitation bubbles causes tissue dam
ages. On the other hand, the bubble collapse in an acous
tic field can be utilized effectively to crush the calculi
(Yoshizawa et al., 2009). A cell permeabilization tech
nique with shock waves or ultrasound induced by the
bubble collapse is also being developed in gene therapy
and anticancer drug delivery (Ohl et al., 2006, Kodama
et al., 2009). Thus the bubble collapse near moving or
deformable boundaries, e.g., ships or tissues, is an im
portant issue to understand the phenomena.
To simulate the bubble collapse, a lot of numerical
methods have been developed. There are two typical
approaches for bubble collapse. One is the Lagrangian
approach in which the interface location is tracked with
the fluid velocity at the interface in a Lagrangian man
ner. The other is the Eulerian approach in which the
interface is captured with a function (e.g., distance func
tion) containing interface information in an Eulerian
grid. The typical method with the Lagrangian approach
is the boundary element method (BEM). The BEM has
been widely used in the bubble dynamics (Guerri et al.,
1982, Blake et al., 1986, 1987, Blake & Gibson, 1987,
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Chahine & Duraiswami, 1992). Since the boundary in
tegral equations are solved without the grids for the am
bient liquid, the BEM is not expensive for 3D analysis
although a difficulty arises in the topology change of
the interface; special treatment is needed after the liq
uid jet threads the bubble interface. On the other hand,
the Eulerian method has an advantage that the topolog
ical change of the bubble shape can be considered eas
ily. However, the accurate threedimensional analysis is
very expensive because the resolution of the interface is
much dependent on the grid size.
In the present works, the BEM is employed to achieve
fast numerical prediction for ship damage due to under
water explosion in designing the ship body. First, the
BEM is applied to the growth and collapse of a bub
ble under a floating body. The threedimensional defor
mation of the bubble, the translation and rotation of the
floating body, and the motion of water surface are taken
into account in the simulation. The Kelvin impulse is
used to analyze the characteristic of the translational mo
tion. Second, a fast simulation method is proposed for
underwater explosion. In the method, the bubble mo
tions are modeled with a point source and a doublet;
the source and doublet represent the radial and transla
tional motions of the bubble, respectively. The bubble
motions, coupled with the motions of the free surface
and the floating body, are solved by using the BEM. Fi
nally, the bubble collapse in an elastic tube is simulated
by taking the deformation of the tube into account.
Numerical method for the bubble motion under
a floating body
Schematic of analysis The underwater explosion
is modeled with a growing and collapsing bubble be
cause the bubble collapse is a primary cause of ship dam
age. The schematic of arrangement of the bubble and the
floating body is shown in Fig. 1. Initially, the body is at
rest with the draft of Lf on the water surface. The cen
ter of gravity of the body is taken to be the origin of the
fixed Cartesian coordinate system; the initial location of
the water surface is assumed to be z 0 The direc
tion of gravity is the positive direction of the z axis. An
initially spherical bubble with high internal gas pressure
is at rest at the position (x, y, z) (Xo, Yo, Lo). The
internal gas is an ideal gas, and is assumed to obey the
adiabatic change. In the analysis, since the characteristic
radius of the bubble that is almost identical to the max
imum radius is about 0.4 m, the heat transfer through
the bubble wall and the surface tension are neglected.
Also, we assume that the liquid surrounding the bubble
is incompressible and invicid, and the flow field is irro
tational. Therefore, the velocity potential D which sat
isfies the Laplace's equation is defined in the analytical
Figure 1: Schematic of analytical model for the bubble
motion under a floating body.
domain.
Method I First, the direct 3D BEM is applied to the
bubble motion under a floating body. We call the fol
lowing method "Method I." Using Green's theorem, the
following boundary integral equation is derived:
(x) (x)
47a
S Gr (, x')
where
G2(x, x')(X').i, (1)
1
G1 G2
4xx' x I
0G1
a)n'
e is the solid angle at the observation point, x the posi
tion vector, S Sb + Sw + Ss, and a/On the outward
normal derivative. Sb, Sw, and S, denote the bubble
surface, the water surface, and the wall of the floating
body, respectively. The boundary conditions on the bub
ble and water surfaces for the boundary integral equa
tion are given by the velocity potential obtained from
the unsteady Bernoulli equation as shown in Eqs. (2)
and (3). Also, the boundary condition for the floating
body is given by the continuity of the normal velocity at
the body wall as in Eq. (4).
DxIV(X)12 + Po '. + gz (x e Sb), (2)
Dt 2 P
D(x) 1I )2 +
DL 2
an(x
9n
n(x) v (x)
(x E S,), (3)
(x E S,), (4)
where t is the time, po the atmospheric pressure, p the
liquid density, g the gravitational acceleration, n the out
ward unit normal, and D/Dt the Lagrangian derivative.
v is velocity of fluid and is calculated as v V4. The
internal gas pressure /, is obtained from the following
equation:
Vo (5)
where / ,, is the initial gas pressure, V the volume of the
bubble, V% the initial volume of the bubble, and 7 the
adiabatic index. The velocity at the wall of the floating
body, vs, is represented by the translational velocity at
the center of gravity, vc, and the angular velocity, wG,
as
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
floating body is given by (Tanizawa, 1995)
n (ia9(x)N
n 9at )
n(x), dvG dwg
l dt +t
+Wc x (w x x)+ 2we x (V$(x) s(x))
Vs V G +WG X Xr,
where x, is position vector in the local coordinate sys
tem whose origin is the center of gravity of the float
ing body. The position vector on the bubble surface and
water surface is tracked with the Lagrangian manner as
Dx/Dt = VP.
The translational and rotational motions of the rigid
floating body is determined by the following equations:
MdV j pndS, (7)
s is
d(I wG) px, x ndS, (8)
at i
where M is the mass of the floating body, I the tensor
of inertia. p is the pressure difference from the initial
pressure, and is defined as
P= P p ~V 2 + pg(z ), (9)
at 2
where zQ is the z value of the initial position vector at
the wall of the floating body. The partial derivative of
the velocity potential with time, 00l/t, is determined
by solving the following boundary integral equation:
On \21 x
(x e S),
(13)
where k, is the normal curvature on the wall. It should
be noted that dvc/dt and dwc/dt in Eq. (13) involve
the unknown parameter Od/lt implicitly. The numer
ical procedure to determine the motion of the floating
body is summarized as follows. First, a0/0t on the
floating body is obtained by solving Eq. (10). Next,
the pressure on the floating body is obtained from Eq.
(9). Using the pressure, the motion of the floating body
is determined by Eqs. (7) and (8).
The linear triangle elements are used in the 3D BEM
analysis. Since the shape and size of boundary elements
become nonuniform as the interface is tracked in the La
grangian manner, the local reconstruction of the bound
ary elements is utilized to diminish numerical instability.
An element whose area is larger than 1.35 times of the
mean area of the boundary elements is divided into small
elements, and an element smaller than 0.7 times of the
mean area is combined with neighboring elements (Un
verdi & Tryggvason, 1992). Exchange of sides of ele
ments is also introduced to maintain the elements uni
form. Smoothing of boundary nodes is also utilized to
diminish the numerical instability (Takahira & Murase,
2006). The Heun's method is used for the time integra
tion. The increment of time At is given by
4e() (x))
47 ( at )
SGi(x, x') O ( (x'))
s O' at
C
max (DF/Dt)'
G(, X') 1. (10)
The unsteady Bernoulli equation is used for the bound
ary condition in Eq. (10) on the bubble surface and the
water surface, respectively, as
IV(x) lV(X)1 I+ + g ( E Sb),
at 2 P
aTi(x) 1'
tIV (X)12 + 2Z
at 2
(x E S,). (12)
The boundary condition for Eq. (10) on the wall of the
where F is a variable to be updated, and C is typically
taken to be C=0.01.
When a liquid jet threads a bubble surface, the bubble
becomes toroidal, and the analytical domain becomes
doubly connected. We use the same numerical proce
dure as Zhang et al. (2001) to consider a toroidal bubble.
In the method, a vortex ring is introduced to consider the
discontinuity of velocity potential after the impact of the
liquid jet on the bubble surface.
Numerical method II A fast prediction method for the
bubble motion under the floating body is also developed
by using a point source and a doublet that represent the
bubble motion. We call the method "Method II" after
this. Introducing the point source and the doublet, the
velocity potential D in the analytical domain satisfies the
(14)
k\ IV, (x) V,(X)12
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
(ii) t/to= 0.828
(iii) t/to= 1.594
(iv) t/to= 1.638
F
z
(v) t/to 1.656 (vi) t/to = 1.660 (vii) t/to = 1.677 (viii) t/to= 1.692
I^^SS^S^^^' It~S~ii~ii~x11ii I ij I" si^^^
PT1_WR ^^i'l X'^^^^^ EWxOHx^lE'si^^ ^x>o
RxPQX1D1SKLX11K l EMxPQXSSXSX' MANNIAIS ~ E XSS SN ^IN
(a) Floating body and bubble
(viii)
Fx
(b) Bubble
Figure 2: Successive bubble shapes when Lf/Rc=1.0 and X/oRf=0.8.
following boundary integral equation:
4w
SGf (x, x') G2(x, x'))x')n.
is On' f
qo qn (XB x)
XB X IXB X13(15)
where XB is the position vector for the bubble centroid
and S S, + S,. The point source qo/xB x and
the doublet q, (XB x)/ XB x13 in Eq. (15) rep
resent the radial and translational motions of the bubble,
respectively, whose strength is given by (Takahira, 1997,
Takahira, 2004)
qo = R2R, (16)
R3
q1 = (VB UB), (17)
where R is the bubble radius, R the radial velocity at the
bubble wall, vB the translational velocity at the bubble
centroid, and UB the fluid velocity at the bubble centroid
induced by the boundaries except for the bubble. Defin
ing 4 as a velocity potential induced by the boundaries
except for the bubble, 4 is given by
sI(XB X') aIOn'
GC2(x, x')(X')1 .i', (18)
f
where UB VB (D in which VB is the gradient in terms
of aB.
The bubble motions are determined by the following
equations (Takahira, 1997, Takahira, 2004):
RRR+ '2 PB 1 V B
2 p 4
dl 1
it 2p
UB)} pV dU+ p
I dt
(19)
0. (20)
Equation (19) is the RayleighPlesset equation including
the influence of the translation, and Eq. (20) is the equa
tion for the translational motion. R is the acceleration
of the bubble wall. pB in Eq. (19) is the liquid pressure
at the bubble centroid determined by the flow fields in
terms of 4, and is written as
04 1 ,
PB = po pat PBI 12 +pgZB
duB/dt in Eq. (20) is the acceleration of the fluid at the
bubble centroid and is given by
dUB VB ( a) 1 2 ,B
dt ( b) + 1
dt t2
(22)
The dynamical motion of a bubble can be determined
by coupling Eqs. (19)(22) with the boundary integral
equations. Substituting uB expressed by the boundary
integral form into q1, the following boundary integral
(i) t/to= 0.000
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
(ii) t/to= 0.817
(iii) t/to= 1.565
(iv) t/to= 1.600
z
(v) t/to = 1.615 (vi) t/to = 1.630 (vii) t/to = 1.654 (viii) t /to = 1.713
a~g~gN09999=0 0999 9999wsi8s ~ 8~r~s
(a) Floating body and bubble
mx
(b) Bubble
Figure 3: Successive bubble shapes when Lf/R,=0.5 and Xo/Rf =1.0.
equation is derived:
0 ( ) D j 1 ,X')] O '
4 x7 ( 2 n'
[G2(x, x') 2Gq 12X, X'), x ),'l 
qo 1R3(xB X)
+ I+1 ( X v) (23)
\XB x 2 XB X 13
where
R3(XB X)VBG(xB, ),
G,11= VBGi(X, '),
a IB x3
R3(XB X) a
Gq12 X13B VB G1(XB, x').
XB O3 On'
The boundary conditions for Eq. (23) are given by Eqs.
(3) and (4), respectively.
The translational and rotational motions of the float
ing body are determined by Eqs. (7) and (8). 0d/lt,
which is needed to calculate the pressure acting on the
body, is calculated by using the following boundary in
tegral equation:
e(a) d9(x)
47 at
\XB X1
(XB
q1X
G2(x, x') '
(XB X)
 q *ovB X 3
x)+ (VB VB) 9 ( ) (24)
X3 (XV 13
The above equation is rewritten by expressing qo and q4
with 9A/Ot as
(x) OGD(x) ),) Gq0[(x, x')
GCo2 (X, x') + G12(X, x')]O .s'
x {2R + .
1 1
P0 + IV4( )12 + V1 IV
2 4
R
XB U2
 UB 2
gZB} + (XB {R3VB[ IVB (X )2 gzB]
Qv vBv (XB x)
qoVB + (VB VB) q ,
XB X3
(25)
where Gqoi and G0o2 are
R
Gqoi = G1(xB, x'),
XB X
R a
Gq02  G1(XB, X').
aXB X1 On'
The boundary conditions of Eq. (25) are given by Eqs.
(12) and (13). It should be noted that dO/lt is also
needed to calculate dvc/dt and dwo/dt, which requires
a similar expansion used in the derivation of Eq. (25).
All boundary integral equations are solved by using the
BEM with linear elements.
(i) t/to 0.000
fimmm
"P
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
0.01
0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
0.09
t/to
 Method I
 Method II
0 0.5 1 1.5 2
t/to
(b)
Method I
. Method II
0.5 1 1.5 2
t/to
t/to
Figure 4: Comparison between Methods I and II (Lf/Rc=l.O, Xo/Rf=0.8): (a) bubble radius, (b) internal pressure,
(c) displacement of bubble centroid in the x direction and (d) displacement of bubble centroid in the z direction.
Numerical results for the bubble motion under
a floating body
Computational conditions Computations have been
done under the following conditions. The atmospheric
pressure: po = 100 kPa. the liquid density: p = 998.2
kg/m3, the gravitational acceleration: g = 9.806 m/s2,
and the polytropic index: 7 = 1.4. The floating body
is a cylinder with the radius of Rf = 0.85 m. The ini
tial bubble radius is Ro = 6.316x10 2 m. The initial
internal pressure is /i = 11.362 MPa, Since the initial
gas pressure is higher than the atmospheric pressure, the
bubble grows in the beginning. After the bubble reaches
the maximum volume, it starts to collapse. Since the
maximum radius is about 0.4 m in the present condi
tions, the characteristic length Rc is taken to be Rc = 0.4
m. The length, time and pressure are normalized with
R6, to(= R/ po/p), and po, respectively. However,
the horizontal offset for the bubble position (the distance
between the axis of floatation (z axis) and the bubble
centroid) is normalized with the body radius. The verti
cal distance between the bubble centroid and the bottom
of the floating body is Lb/R = 1.25.
Deformation of a bubble under a floating body Fig
ure 2 shows the motion of a floating body and successive
bubble shapes when Lf/R6c 1.0 and Xo/Rf 0.8.
The overall view from the y direction and the enlarge
ment of bubble shapes at time (v) (viii) are shown in
(a) and (b), respectively.
First, the bubble grows rapidly because of its high in
ternal gas pressure, and its volume reaches the maximum
at (ii) t/to = 0.828. At this stage, the floating body ro
tates counterclockwise due to the fast expansion of the
bubble. When the bubble collapses, it becomes a egg
like shape and translates toward the floating body be
cause the inertia of the fluid surrounding the bubble is
asymmetric. The bubble volume becomes the minimum
at t/to = 1.660 in Fig. 2(vi). The liquid jet developed
in Fig. 2(v) impacts on the upper surface of the bub
ble at the t/to = 1.662. The liquid jet becomes narrow
with increasing the bubble volume during the rebound
ing stage.
Figure 3 shows the results in the case of Lf /Rc 0.5
and Xo/Rf 1.0 where the mass of the floating body
is half of that in Fig. 2. At (ii) t/to 0.817, the bubble
volume reaches the maximum. Since the mass of the
floating body in Fig. 3 is lighter than that in Fig. 2, the
rotation of the floating body is larger. Also, the water
surface near the bubble rises as the bubble grows. The
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
(b)
Method I
S  Method II
(c)
Method I
. Method II
0.45
0.4
0.35
0.3
0.25
0.2
0.15
0.1
0.05
0
0.05
0 0.5 1 1.5 2
t/to
0.5 1 1.5 2
t/to
(d)
Method I
 Method II
0 0.5 1 1.5 2
t/to
Figure 5: Comparison between Methods I and II (Lf/Rc=0.5, Xo/IRf=.0): (a) bubble radius, (b) internal pressure,
(c) displacement of bubble centroid in the x direction and (d) displacement of bubble centroid in the z direction.
liquid jet occurs in the opposite direction to the floating
body.
The time histories of (a) the equivalent bubble radius
(Req = (3V/47)1/3), (b) the internal gas pressure, (c)
the displacement of the bubble centroid in the x direc
tion, and (d) the displacement of the bubble centroid in
the x direction are shown in Figs. 4 and 5. Figures 4
and 5 correspond to the results of Figs. 2 and 3, respec
tively. The displacements of the bubble centroid in the
x and z directions are defined as bx = XB Xo and
bz ZB Lo, respectively, where (XB,ZB) is the posi
tion of the bubble centroid. The bubble does not move in
the y direction because of the symmetric configuration.
The solid and dotted lines show the results for Method I
and Method II, respectively.
The comparison of the bubble radius and the in
ternal gas pressure in Fig. 4 shows that the predic
tion with Method II is in good agreement with that of
Method I; the minimum radius at the collapsing stage
is Req/Rc=0.167 for Method I and Req/Rc=0.163 for
Method II; the maximum pressure is /' /1',,=v 7 for
Method I and / /1 ,,="11 7 for Method II. On the other
hand, Method II overestimates the displacement of the
bubble centroid in the final stage of the collapse and
the second rebounding stage although the agreement be
tween Method I and II matches before the final stage of
the collapse. One reason for the overestimation is that
neglecting the deformation of the bubble in Method II
underestimates the virtual mass of the bubble, and con
sequently, tends to overestimate the translational accel
eration of the bubble (Jinbo et al., 2010). The same ten
dency stated above is found in Fig. 5. Since the com
putational condition in Fig. 5 is close to the condition
of the neutral collapse, the difference of the displace
ment of the bubble centroid between Methods I and II
is relatively large. It is found that the computation with
Method II is about 6 times as fast as that with Method I.
The trajectory of the bubble centroid Figure 6 is the
trajectories of the bubble centroid for the simulation of
Fig. 3 (Lf/Rc 0.5). The abscissa and coordinate de
note bx/Rc and bz/Rc, respectively. The lefthand side
and righthand side of the figures correspond to the re
sults for (a) Method I and (b) Method (II), respectively.
The trajectory of the centroid is plotted for the various
initial horizontal bubble location Xo/Rf. The solid cir
cle, the cross, and the open triangle indicate the initial
bubble location, the location when the bubble takes the
maximum volume, and the location when the bubble
C 0.
0.
S0.
n
0.01
0
0.01
0.02
0.03
0.04
0.05
0.06
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
0.1
0.05
0
0.05
0.
0.1
0.15
0.2
0.25
0.3
0.08
0 0.02
0.08 0.06 0.04 0.02
b / Rc
(b) Method II
Figure 6: Trajectories of bubble centroid.
0.004
0.008
0.012
0.016
0.02
0.024
0.028
1.0
1.0
.6"',,,0.6
 Method I ,
..... Method II, "
0.2
Xo/Rf= 0.35 "0 35
0 0.5 1 1.5 2
t /to
Figure 7: The Kelvin impulse in the vertical direction.
takes the minimum volume, respectively. Thus, the du
ration of the bubble movement from the solid circle to
the open triangle corresponds to one cycle of the bub
ble motion. It is well known that the bubble is strongly
affected by the inertia of the boundaries surrounding
the bubble. In general, the bubble translates toward the
rigid boundary when it collapses near the rigid boundary,
while it does away from the boundary when it collapses
near the free surface. Therefore, when the floating body
is heavy, the bubble translates toward the floating body
during the collapsing period. However, when the initial
horizontal bubble location is close to the water surface
Figure 8: The Kelvin impulse in the horizontal direc
tion.
even if the floating body is heavy, the bubble is influ
enced more strongly by the water surface than the float
ing body. In Fig. 6, since the body is not so heavy, the
bubble translates in the opposite direction to the floating
body (i.e., b, is positive) when the initial horizontal posi
tion is larger than Xo/Rf 0.2. On the other hand, the
displacement in the x direction does not change mono
tonically with increasing the horizontal offset: the dis
placement increases when the bubble is initially located
at Xo/Rf = 0.35 and 1.2. The largest horizontal trans
lation is realized between the axis of flotation and the
edge of the floating body. Although the prediction with
Xo/Rf= 0.0
0.06 0.04 0.02
bx / R
(a) Method I
0 0.02
t / to
I
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Method II overestimates the translation, Method II cap
tures basically the same tendency as Method I for the
trajectories.
Kelvin impulse for the bubble The characteristic
of the bubble translation will be discussed using the
Kelvin impulse defined by the following equation (Blake
& Gibson, 1987, Blake, 1988):
IK l= P ndS. (26)
j Sb
The definition of the Kelvin impulse is related to the vir
tual inertia of the bubble. The value of the Kelvin im
pulse at the bubble collapse is an indicator to judge the
direction of the bubble translation and the liquid jet. The
time histories of the Kelvin impulses for Lf/Rc 0.5
in the z and x directions are shown in Figs. 7 and 8,
respectively. The Kelvin impulse is plotted for various
Xo/Rf in the figures. The Kelvin impulse is normal
ized with the R /ppo. The solid and dotted lines show
the results for Method I and Method II, respectively. If
the Kelvin impulse in the z direction is negative during
the collapsing stage, the bubble is expected to translate
in the negative z direction, while it does in the posi
tive z direction if the Kelvin impulse is positive. There
fore, the neutral collapse where the bubble almost stays
in the initial position may be realized when the Kelvin
impulse is nearly zero (Blake & Gibson, 1987, Duncan
& Zhang, 1991, Brujan et al., 2005). In the case of
Xo/Rf = 0.2 in Fig. 7, the Kelvin impulse oscillates
around zero value, which shows the condition is close
to the neutral collapse. Since the mass of the body is
light and the influence of the water surface is relatively
strong, the Kelvin impulse becomes positive during the
collapsing stage when Xo/Rf is 0.6 and 1.0 in Fig. 7;
the bubble translates away from the water surface. On
the other hand, the Kelvin impulse in the horizontal di
rection in Fig. 8 is always negative, which suggests that
the bubble translates toward the flotation axis when the
bubble collapses. The absolute value of the Kelvin im
pulse is largest at Xo/Rf=0.35 and decreases in order
of 0.2, 0.6, and 1.0. This result agrees with the result in
Fig. 6 that the largest horizontal translation of the bubble
centroid is realized when Xo/Rf=0.35. The Kelvin im
pulse in the vertical direction is predicted well using the
Method II although the Kelvin impulse in the horizontal
direction is underestimated for Method II.
In the Fig. 9, the horizontal Kelvin impulse for Lf /R
= 0.5 at different times are plotted versus the initial hori
zontal offset of the bubble where TRmax is the time when
the bubble reaches the maximum volume. The lines and
symbols are corresponding to the results for Method I
and Method II, respectively. As evident in Fig. 9, the
results with Method II agree with those with Method
I. The horizontal Kelvin impulse has two local minimal
0.005
0.01
0.015
0.02
0.025
(
22
0 X
60
a '
CW y
i''A 7 '^
S 0.4 0.8 1.2
Xo/ R
t/ TRmax 0.5
Method I
o Method II
t/Tmax 1.0
Method I
o Method II
t/ TRmax 1.5
Method I
A Method II
t/TRmax 1.9
 Method I
0 Method II
1.6 2
Figure 9: Influence of horizontal offset on the Kelvin
impulse in the x direction.
0.2 .,\
S0.4
7 0.6
0.8
S1.0
1.2
S1.4
1.6
0 0.5 1
Xo/ Rf
1.5 2
Figure 10: Influence of horizontal offset on the initial
Kelvin impulse in the x direction (Method I).
values at Xo/Rf 0.35 and Xo/Rf 1.2 at any stage
of growth and collapse. As mentioned in Fig. 6, the
horizontal translation does not change monotonically in
terms of the initial horizontal offset; the largest horizon
tal translation is found in Xo/Rf 0.35. This transla
tional characteristic under the floating body is caused by
the fact that the Kelvin impulse takes two local minima
in terms of the initial horizontal offset.
Figure 10 shows the relationship between the horizon
tal offset and the initial horizontal Kelvin impulse at one
time step (t A t) for three kinds of draft of the floating
body (Lf/Rc =0.5, 1.0, and 2.0). It should be noted
that the mass of the floating body is proportional to the
draft. Even in the very early stage of the bubble mo
tion, two local minima are observed in the case of the
Lf/Rc =0.5 and 1.0 although the absolute value of the
Kelvin impulse is small. This suggests that it is possible
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
0
0.2
0.4
4 0.6
0.8
L 1.0
S1.2
1.4
1.6
0 0.5 1
Xo / Rf
I /Iref 0.5
I/Iref 2.0
1.5 2
Figure 11: Influence of the moment of inertia on the
initial Kelvin impulse in the x direction (Method I).
to discuss the relationship between the horizontal off
set and the horizontal Kelvin impulse using the results
in the early stage of the bubble motion. As evident in
Fig. 10, the Kelvin impulse takes only one local mini
mum when the floating body is heavy (Lf/Rc = 2.0),
which suggests that the characteristic of the horizontal
translation is determined not only by the initial position
but also by the mass and moment of inertia of the float
ing body. The influence of the moment of inertia on the
horizontal Kelvin impulse is shown in Fig. 11. The mo
ment of inertia is taken as I/Iref = 0.5 and 2.0 where
Iref is the moment of inertia for the floating body with
Lf/Rc 1.0. The dimensionless moment of inertia
of the reference body, Iref/(pR5) is 16 in terms of the
principal axis in the radial direction of the floating body.
When the moment of inertia is small (I/Iref = 0.5), the
absolute value of the first local minimum of the horizon
tal Kelvin impulse near Xo/Rf = 0.35 becomes promi
nent and that near Xo/Rf 1.0 decreases. A local min
imum is found near Xo/Rf 1.0 for the floating body
with the large moment of inertia (I/Iref = 2.0). Thus the
characteristic of the horizontal translation of the bubble
is much dependent on the moment of inertia.
The bubble motion in an elastic tube
Numerical model The numerical model of an elas
tic tube is shown in Fig. 12. The radius and length of
the tube are Rt and Lt, respectively. The fluid in the
tube is initially at rest. The compressibility and viscos
ity of the fluid are disregarded. The lefthand side (inlet)
and righthand side (outlet) boundaries of the tube are
fixed to a constant pressure, and the fluid is allowed to
flow in or out according to the bubble motion. In the
present computation, the inlet pressure is equal to the
outlet pressure, and is taken to be the atmospheric pres
Figure 12: Schematic of an elastic tube.
sure po. The gas bubble with the high internal pressure
is also initially at rest at (Xo, Yo, Zo). The origin of
the Cartesian coordinate is taken to be the center of the
tube. The side wall of the tube deforms only in the ra
dial direction. The elasticity of the tube is modeled with
springmass systems (Duncan & Zhang, 1991, Takahira
et al., 2002). The boundary integral equation for the
problem is basically the same as Eq. (1). The bound
ary conditions at the inlet and outlet are given by
DtF(x) 1 V(x)12 + po I.
Dt 2 p
8, ( X) 1 iv ( )
at 2
(X) 1
at 2
(x Sb), (27)
(x E Sinlet), (28)
(x E Soutlet), (29)
where Sb is the bubble wall, Sinlet the inlet boundary,
Soutet the outlet boundary. Letting r be the displace
ment of the side wall in the radial direction, the bound
ary condition on the side wall (Sside) is written as
OK(x) =n(x). (Or(x) e)(x) (x Sside), (30)
an at
where e, is the unit normal in the radial direction The
motion of the side wall is determined by the following
equations:
a2r k pw po
at2 m m
(31)
where m is the mass per unit area, and k the spring stiff
ness per unit area. p, is the pressure at the side wall,
and is obtained from
S0i 1
Pw P PV I2 +po.
(32)
a0/at is calculated by using the similar boundary inte
gral equation to Eq. (10) whose boundary conditions on
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
(i) t/tRt = 0.0000
(ii) : =0.5397
(iii) t/tR =1.0597
(iv) t/tR = 1.8797
IZ ..
(v) t/tRt = 2.0397
(vi): =2.0931
(vii) t/tt = 2.1204
(viii) t/tt = 2.1421
(v)
"
l i l ,1 1
. : , T
'f^\
j I i'i ri li3 ++";
^..::' .
., ij' , ' *,)*;?
q $ 'tt ; i1
(viii)
y ? ,y i.
tUz t
(a) elastic tube and bubble (b) bubble
Figure 13: Successive shapes of bubbles and tubes when It/Rt 1.0 (Case A).
Sb, Sinlet, and Soutlet are
8a(X) 1V4(x) + Po .
at 2 p
at 2
8+(X) 1 V+(+)
a( 2)
a~t 3()2
S(x eSb), (33)
(x E Sinlet), (34)
SE Soutlet). (35)
The boundary condition on the side wall of the tube is
a9 (0 ) e, (x) nx(x)
an+ at at2 e n( x
at at Qt an
a(x) {V+(x)} n(x) (x E Sside). (36)
at an
Numerical results In the following computations, the
dimensionless parameters are introduced using Rt, po
and p. The dimensionless length of the tube and the ini
tial bubble radius are Lt/Rt 45 and Ro/Rt 0.05,
respectively. Xo/Rt and Zo/Rt for the initial bub
ble position are taken to be zero. The dimensionless
mass and the dimensionless spring stiffness are taken to
be m/(pRt) = 4.0 and kRt/po = 4.0, respectively.
The dimensionless initial gas pressure inside the bub
ble is /i ,,/po 113.6. The time is normalized with
tRt = RIP0/ p
Figures 13 (Case A) and 14 (Case B) are successive
shapes of bubbles and tubes. The overall views from the
x direction are shown in the figures. Also, the enlarge
ments of the 3D bubble shape are shown at the specified
time. Since the bubble is initially located in the axis of
the tube in Case A (Yo/Rt = 0), the axisymmetric mo
tion is expected in Fig. 13. On the other hand, the initial
location of the bubble in the y direction in Case B is
Yo/Rt 0.3. Letting the distance between the initial
bubble centroid and the side wall be It, It/Rt for Cases
A and B are 1.0 and 0.7, respectively. Figure 15 shows
the time histories of the displacement of the side wall at
the points N, S, and E in Fig. 12. Only one curve is plot
ted for Case A because of the symmetric arrangement.
Because of the symmetrical arrangement, no transla
tional motion occurs for the bubble in Case A. The bub
ble is elongated in the y direction at (v) t/tRt 2.0397.
Then, the instability of the bubble surface develops in
the y direction leading to the fission of the bubble. As
shown in Figs. 13 and 15, the side wall of the tube ex
pands in the growing stage of the bubble, and shrinks in
(vii)
a'^ i^ i:
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
(i) t/tRt = 0.0000
. . . . .
(ii).' =0.5381
(iii) :' =1.0581
(iv):' =1.8681
z
(v) t/tRt= 2.0381
(vi):' =2.0870
(vii):. =2.1166
(viii):. =2.1370
::::
" ;'77:,
;j'i
 I i
:* 
'3,
(viii)
y
e 2j'
L)z
t"x
(a) elastic tube and bubble (b) bubble
Figure 14: Successive shapes of bubbles and tubes when ItRt = 0.7 (Case B).
the collapsing stage. On the other hand, in Case B, the
asymmetric arrangement of the bubble causes the bub
ble translation toward the y direction and the liquid jet
is formed on the bubble surface; the liquid jet also de
velops toward the side wall (the y direction). Since the
distance between the initial bubble centroid and the side
wall is shorter in Case B, the displacement at the point
N is large in Fig. 15. The sink flow induced by the bub
ble collapse results in the further constriction of the side
wall in the vicinity of the bubble. As evident in the fig
ures, the bubble motion is significantly affected by the
side wall.
Conclusions
The growth and collapse of a bubble under a floating
body were simulated by using the threedimensional
boundary element method with linear element in which
the translation and rotation of a floating body and the
free motion of a water surface were taken into account.
Two kinds of simulation methods were developed in the
present study: the direct three dimensional dynamics of
bubbles is considered in one method (Method I); the
bubble dynamics is modeled with a source and a dou
blet that represent the radial and translational motions
in the other method (Method II). It was shown that the
0.02
0
0.02
0.04
0.06
0.08
0.1
0.12
0.14
0.16
0.18
Case A
 Case B: N \
 Case B: S
 Case B: E
0 0.5 1 1.5 2 2.5
t/ tRt
Figure 15: Displacement of the side wall of an elastic
tube.
faster simulation was done with Method II. Although the
prediction using Method II overestimated the translation
of the bubble, it was in qualitative agreement with that
using Method I. Also, the bubble collapse in an elas
tic tube was simulated by taking the deformation of the
tube. The results are summarized as follows.
(i) The bubble behavior is significantly affected by the
interactions among the bubble, the floating body, and
the water surface; the directions of the bubble transla
tion and the liquid jet are much dependent on the initial
horizontal position of the bubble.
(ii) The characteristic of the translational motion of
the bubble is evaluated with the Kelvin impulse.
(iii) The horizontal translational motion of the bubble
is much dependent on the rotational motion of a floating
body; when the moment of inertia of the floating body
is small, the largest horizontal translation is realized be
tween the axis of flotation and the edge of the floating
body.
(iv) When a bubble collapses near the side wall of the
elastic tube, the liquid jet develops toward the side wall
near the bubble. The side wall is constricted due to the
bubble collapse.
References
Blake, J. R., Taib, B. B. and Doherty, G., Transient Cav
ities near Boundaries. Part 1. Rigid Boundary, J. Fluid
Mech., Vol. 170, pp. 479497, 1986.
Blake, J. R., Taib, B. B. and Doherty, G., Transient Cav
ities near Boundaries. Part 2. Free Surface, J. Fluid
Mech., Vol. 181, pp. 197212, 1987.
Blake, J. R. and Gibson, D. C., Cavitation Bubbles near
Boundaries, Ann. Rev. Fluid Mech., Vol. 19, pp. 99123,
1987.
Blake, J. R., The Kelvin Impulse: Application to Cav
itation Bubble Dynamics, J. Aust. Math. Soc., B30,
pp.127146, 1988.
Brennen, C.E., Cavitation and Bubble Dynamics, Ox
ford Univ. Press, 1995.
Brujan, E. A., Pearson, A. and Blake, J. R., Pulsating,
Buoyant Bubbles close to a Rigid Boundary and near
the Null Final Kelvin Impulse State, Int. J. Multiphase
Flow, Vol. 31, pp.302317,2005.
Chahine, G. L. and Duraiswami, R., Dynamical Interac
tions in a MultiBubble Cloud, Trans. ASME, J. Fluids
Eng., Vol. 114, pp. 680686, 1992.
Duncan, J. H. and Zhang, S., On the Interaction of a Col
lapsing Cavity and a Compliant Wall, J. Fluid Mech.,
Vol. 226, pp. 401423, 1991.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Guerri, L., Lucca, G. and Prosperetti, A., A Numerical
Method for the Dynamics of NonSpherical Cavitation
Bubbles, Proc. 2nd Int. Colloq. on Drops and Bubbles,
pp. 175181, 1982.
Imakita, A. and Yasuda, A., Response of Cylindri
cal Floating Structure Subjected to an Underwater Ex
plosion, J. Jpn. Soc. Nay. Archit. Ocean Eng., Vol.6,
pp. 379386, 2007 (in Japanese).
Jinbo, Y., Takahira, H., Kobayashi, K. and Yasuda, A.,
Boundary Element Analysis of Bubble Dynamics under
a Floating Body (2nd Report, Approximation Method
with Spherical Bubble Model), Trans. Jpn. Soc. Mech.
Eng., Vol. 76, B, pp. 230238, 2010 (in Japanese).
Keil, A. H., The Response of Ships to Underwater Ex
plosions, Trans. Soc. Nay. Arch. Mar. Eng., Vol. 69,
pp. 366410, 1961.
Klaseboer, E., Khoo, B.C. and Hung, K.C., Dynam
ics of an Oscillating Bubble near a Floating Structure,
J. Fluids Struct., Vol. 21, pp. 395412, 2005.
Kodama, T., Tomita, Y., Watanabe, Y., Koshiyama, K.,
Yano., T. and Fujikawa, S., Cavitation Bubbles Mediated
Molecular Delivery During Sonoporation, J. Biomech.
Sci. Eng., Vol4, pp. 124140, 2009.
Ohl, C.D., Arora, M., Ikink, R., Jong, N., Versluis, M.
and Delius, M., Sonoporation from Jetting Cavitation
Bubbles, Biophys. J., Vol. 91, pp. 42854295, 2006.
Takahira, H., Growth of Traveling Bubbles Near an
Axisymmetric Body in a Potential Flow, JSME Int. J.,
Vol. 40, B, pp. 240249, 1997.
Takahira, H., Yasuda, A. and Uemura, A., Effects of
Heat Transfer of Internal Gas on the Bubble Collapse
Near a Compliant Wall, Proc. 5th JSMEKSME Fluids
Eng. Conf., CDROM, Total 6 pages, 2002.
Takahira, H., A Remark on the Pressure Terms in the
RayleighPlesset Equation for Cavitating Flows Trans.
Jpn. Soc. Mech. Eng., Vol.70, B, pp. 617622, 2004 (in
Japanese).
Takahira, H., Murase, Y. and Kamiirisa, H., Numerical
Simulation of Threedimensional Bubble Collapse un
der a Floating Body, Proc. 6th Int. Symp. on Cavitation,
Total 8 pages, 2006.
Tanizawa, K., A Nonlinear Simulation Method of 3
D Body Motions in Waves, J. Soc. Nay. Archit. Jpn.,
Vol. 178, pp. 179191, 1995.
Unverdi, S.O. and Tryggvason, G., A FrontTracking
Method for Viscous, Incompressible, Multifluid Flows,
J. Comput. Phys., Vol. 100, pp. 2537, 1992.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Yoshizawa, S., Ikeda, T., Ito, A., Ota, R., Takagi, S.
and Matsumoto, Y., High intensity focused ultrasound
liii IIips with cavitating microbubbles, Med. Biol.
Eng. Comput., Vol.47, pp. 851860, 2009.
Zhang, Y.L., Yeo, K.S., Khoo, B.C. and Wang, C.,
3D Jet Impact and Toroidal Bubbles, J. Comput. Phys.,
Vol. 166, pp. 336360, 2001.
