Group Title: 7th International Conference on Multiphase Flow - ICMF 2010 Proceedings
Title: 7.2.2 - Boundary Element Analysis for Bubble Collapse near Moving or Deformable Boundaries
ALL VOLUMES CITATION THUMBNAILS PAGE IMAGE ZOOMABLE
Full Citation
STANDARD VIEW MARC VIEW
Permanent Link: http://ufdc.ufl.edu/UF00102023/00176
 Material Information
Title: 7.2.2 - Boundary Element Analysis for Bubble Collapse near Moving or Deformable Boundaries Particle Bubble and Drop Dynamics
Series Title: 7th International Conference on Multiphase Flow - ICMF 2010 Proceedings
Physical Description: Conference Papers
Creator: Jinbo, Y.
Takahira, H.
Publisher: International Conference on Multiphase Flow (ICMF)
Publication Date: June 4, 2010
 Subjects
Subject: bubble
boundary element method
Kelvin impulse
floating body
elastic tube
 Notes
Abstract: The Boundary Element Method (BEM) is utilized to predict the bubble collapse near moving or deformable boundaries. First, the BEM is applied to the growth and collapse of a bubble under a floating body. The three-dimensional 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 simulated by taking the deformation of the tube into account. The side wall of the tube is constricted by the bubble collapse.
General Note: The International Conference on Multiphase Flow (ICMF) first was held in Tsukuba, Japan in 1991 and the second ICMF took place in Kyoto, Japan in 1995. During this conference, it was decided to establish an International Governing Board which oversees the major aspects of the conference and makes decisions about future conference locations. Due to the great importance of the field, it was furthermore decided to hold the conference every three years successively in Asia including Australia, Europe including Africa, Russia and the Near East and America. Hence, ICMF 1998 was held in Lyon, France, ICMF 2001 in New Orleans, USA, ICMF 2004 in Yokohama, Japan, and ICMF 2007 in Leipzig, Germany. ICMF-2010 is devoted to all aspects of Multiphase Flow. Researchers from all over the world gathered in order to introduce their recent advances in the field and thereby promote the exchange of new ideas, results and techniques. The conference is a key event in Multiphase Flow and supports the advancement of science in this very important field. The major research topics relevant for the conference are as follows: Bio-Fluid Dynamics; Boiling; Bubbly Flows; Cavitation; Colloidal and Suspension Dynamics; Collision, Agglomeration and Breakup; Computational Techniques for Multiphase Flows; Droplet Flows; Environmental and Geophysical Flows; Experimental Methods for Multiphase Flows; Fluidized and Circulating Fluidized Beds; Fluid Structure Interactions; Granular Media; Industrial Applications; Instabilities; Interfacial Flows; Micro and Nano-Scale Multiphase Flows; Microgravity in Two-Phase Flow; Multiphase Flows with Heat and Mass Transfer; Non-Newtonian Multiphase Flows; Particle-Laden Flows; Particle, Bubble and Drop Dynamics; Reactive Multiphase Flows
 Record Information
Bibliographic ID: UF00102023
Volume ID: VID00176
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: 722-Jinbo-ICMF2010.pdf

Full Text



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
1-1 Gakuen-cho, Naka-ku, Sakai, Osaka 599-8531, Japan
jinbo06@fluid.me.osakafu-u.ac.jp and takahira@me.osakafu-u.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 three-dimensional
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 three-dimensional 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 three-dimensional 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
4-xx' 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).

D|xIV(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 Rayleigh-Plesset 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 9999wsi-8s ~ 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 x-3
R3(XB X) a
Gq12 X13B VB G1(XB, x').
XB O-3 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)
X|3 (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 X|3


(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 counter-clockwise 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 left-hand side
and right-hand 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 .,\
S-0.4
7 -0.6
-0.8
S-1.0
-1.2
S-1.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
S-1.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 left-hand side (inlet)
and right-hand 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
spring-mass 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) 1|V4(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 axi-symmetric 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 three-dimensional
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. 479-497, 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. 197-212, 1987.

Blake, J. R. and Gibson, D. C., Cavitation Bubbles near
Boundaries, Ann. Rev. Fluid Mech., Vol. 19, pp. 99-123,
1987.

Blake, J. R., The Kelvin Impulse: Application to Cav-
itation Bubble Dynamics, J. Aust. Math. Soc., B30,
pp.127-146, 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.302-317,2005.

Chahine, G. L. and Duraiswami, R., Dynamical Interac-
tions in a Multi-Bubble Cloud, Trans. ASME, J. Fluids
Eng., Vol. 114, pp. 680-686, 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. 401-423, 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 Non-Spherical Cavitation
Bubbles, Proc. 2nd Int. Colloq. on Drops and Bubbles,
pp. 175-181, 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. 379-386, 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. 230-238, 2010 (in Japanese).

Keil, A. H., The Response of Ships to Underwater Ex-
plosions, Trans. Soc. Nay. Arch. Mar. Eng., Vol. 69,
pp. 366-410, 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. 395-412, 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. 124-140, 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. 4285-4295, 2006.

Takahira, H., Growth of Traveling Bubbles Near an
Axisymmetric Body in a Potential Flow, JSME Int. J.,
Vol. 40, B, pp. 240-249, 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 JSME-KSME Fluids
Eng. Conf., CD-ROM, Total 6 pages, 2002.

Takahira, H., A Remark on the Pressure Terms in the
Rayleigh-Plesset Equation for Cavitating Flows Trans.
Jpn. Soc. Mech. Eng., Vol.70, B, pp. 617-622, 2004 (in
Japanese).

Takahira, H., Murase, Y. and Kamiirisa, H., Numerical
Simulation of Three-dimensional 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. 179-191, 1995.

Unverdi, S.O. and Tryggvason, G., A Front-Tracking
Method for Viscous, Incompressible, Multi-fluid Flows,
J. Comput. Phys., Vol. 100, pp. 25-37, 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. 851-860, 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. 336-360, 2001.




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

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