A boundary integral numerical method for the dynamics of nonspherical cavitation bubbles in inviscid incompressible liquids is described. Only surface values of the velocity potential and its first derivatives are involved. The problem of solving the Laplace equation in the entire domain occupied by the liquid is thus avoided. The collapse of a bubble in the vicinity of a solid wall and the collapse of three bubbles with collinear centers are considered

Lucca, G.

Prosperetti, A.

English

NASA Technical Reports Server

A aurnerical method for the dynamics of non-spherical cavitation bubbles 
L. Guerri 
Istituto di Hatematica, Universita di Pavia, 27100 Pavia, Italy 
G. Lucca anci A. Prosperetti 
Istituto di Fisica, Via Celoria 16, 20133 Hilano, Italy 
Abstract 
In this paper a boundary integral numerical method for the dynamics of non-spherical cav- 
itation bubbles in inviscid incompressible liquids is Cescribed. A very efficient feature of 
this approach is that ~t involves only surface values of the velocity potential and its first 
derivatives, thus avoiding the problem of solving the Laplace equation in the entire domain 
occupied by the liquid. To illustrate the performance of the method the collapse of a bubble 
in the vicinity of a solid wall and the collapse of three bubbles with cc'linear centers are 
considered. 
Introduction 
The development of the study of the dynamics of non-spherical bubbles has been seriously 
hindered in the past by the lack of computational tools applicable to severe deformations of 
the spherica'. shape. Existing methods I - '  are cumbersome to implement and require long exe- 
cution times. These problems have prevented their widespread usage and common acceptance. In 
this study we shall present a much more efficient method based on the boundary inteq-a1 ap- 
proach which can be implemented in a relatively small number of program instructions (500- 
600). It is flexible, precise, and requires short execution times (typically 3-5 minutes on 
a Univac 1100). The method is applicable to irrotational, incompressible, inviscid flow and 
we present an axisymmetric version of it. To illustrate the results obtainable in this way 
we consider the Rayleigh problem, one of the cases studied by Plesset and Chapman ', and the 
collapse of three bubbles with collinear centers. 
Mathematical formulation 
We consider the flow induced by one or more bubbles in a liquid occupying a domain R 
bounded by the surface of the bubbles, S , rigid boundaries, S , and a "surface at infinity" 
S .  In the hypotheses of inviscid irrotakional flow the velocity field g can be derived from 
a potential, 2 = I$, and the mathematical problem can be put in the form' 
3 = o  o n s  
an r p  
where p is the pressure, p is the density, a/an is the normal derivative, and d/dt=a/at+P@* 
is the convective derivative.The pressure in the bubble,pb,and at infinity, p-, are taken 
constant in this study, although this is by no means required by the method. Surface tension 
effects have been disregarded, but they can be included with relative ease. 
The strategy of solution is in principle quite straightforward. Suppose that the position 
of the bubble surfaces Sb and the value of @ on Sb are known at time t. Then Eq.(l) can be 
solved using the houndary conditions (2),(5), and the known 4 on Sb. From a knowledge of the 
potential it is then possible to evaluate the right-hand sides of (3) and (41, which allows 
https://ntrs.nasa.gov/search.jsp?R=19820015567 2020-03-21T16:51:53+00:00Z
the computation of the new configuration of Sb and of the potential on it at time t+At. The 
cycle is then repeated. 
The crucial step in this procedure is the solution of the potential problem. In the past 
this has been done by finite differences ' r 2 r '  or by singularity ;;lethods '@'. The first ap- 
proach is wasteful because it requires tha: the potential be knorn over the entire domain n, 
although only free surface values are needed . In addition, the position of Sb does not co- 
incide with grid nodes in general, which leads to well-known problems of accuracy, coding. 
and others. The second method avoids the explicit solution of (1) because the assumed form 
for 4 already is a solution of the Laplace equation. However a certain amount of guessing is 
necessary, and it is not clear whether very extreme bubble deformations can be described 
sattsfactorily dth this approach. Our technique, based on the boundary integral, or bound- 
ary element,method5, retains the advantages of the singularity method wnile avoiding its dis- 
advantages since it is based on an exact relation, Green's identity. This identity for a 
point !! belonging to Sb or S has the formS 
r 
In conditions of axial symmetry, which we assume, this relation can be simplified by carrying 
out the integrstion over the azimuthal angle explicitly with the result 
In this equation x and I' belong to a meridian (halflplane passing through the axis of sym- 
metry and rb and rr are the traces on this plane of Sb and Sr respectively: s' is the arc 
length on these curves r. In the meridian plane we choose a Cartesian system of coordinates 
(r,z) , of which the z-axis coincides with the axis of symmetry. Then G,H have the form 
where A = (r+r')'+ (z-z')', m=4rr1/A, and K denotes the complete elliptic integral of the 
first kind. 
In the applications to be described in this paper the existence of the riqid boundary is 
accounted for by introducing a system of "image" bubbles. Therefore condition (2) will be 
implicitly satisfied and it is no longer necessary to indicate explicitly the rigid boundary 
Sr or rr. Accordingly, we shall not carry these symbols along in the following equations and 
we shall denote rb simply by r . 
Numerical method 
Since 4 is known on r from the previous time step, Eq.(7) can be regarded as an integral 
equation for am/anS. To solve it, we approximate the line I' by a polygonal line consisting 
of n segments 1. on each of which 4 and3@/an1 are taken constants.(rJe have also made some 
numerical experiments assuming a 1 inear variation of 4 obtaining virtual1 y identical results) . 
In this way we obtain from (7) a linear system 
n 
where w = a@(yj)/ an' can be interpretecl as the value oi: the normal velocity at the mid- 
point of Tj and 
a = G r , z r , z 0  d 8  , bi= ((4) + I , I ~(r~,z~;r'.z') $(r',z8) ds' . (10)  
i j 
r j=l 
Here r ,ai are to be interpreted as the coordinates of the midpoint yi of r The integrals i' in (lot are performed numerically by a six-point Gaussian formula except for those over rj-ri 
for which the integrand has an (i-tegrable) logarithmic singularity. A series expansion in 
the neighborhood of this singularity coupled with Simpson's formula over the rest of the 
interval is used in this case.The elliptic integrals are computed using the simple approxim- 
ate formulae of Ref. 8. The tangential velocity along the bubble surface is computed to the 
same accuracy as the normal velocity as follows 
where ri+1 and zi are the extrema of the segment Ti. The normal and tangential velocities 
computed in this way are then referred to the (r,z) system of coordinates and transported 
by linear interpolation from the midpoints to the extrema of the segments ri. 
At this point Eqs. (3) and (4) can be used to compute the advanced-time values of the 
positions of the vertices of the polygonal line and of the associated velocity potentials. 
For this purpose the following second-order formula is used 
where T= At / At - and Atnut -t . A similar relation is used to compute the new values 
fie of thg posftlons of xi. ~ ~ : ~ p o ~ s i b i l i t ~  of using higher order time integration formu- lae ot this type is a distinctive advantage of the present method over finite-differences 
ones. This advantage stems from the fact that here it is possible to follow the individual 
trajectories o: +he vertices of the polygonal line during the entire calculation. The time 
step used in (12) is adjusted during the calculation so as to prevent excessive variations 
of ri, zi during a single step. 
Results 
We shall present the results in terms of lengths and times made dimensionless with re- 
spect to %, the initial bubble radius, and { p/ (p,-pb) ) ' '' . The pressure inside the 
bubble and at infinity have been kept constant in all the exam:)les discussed. 
To test the reliability of the method and of the code we have fimt of all computed the 
collapse 3f a single spherical bubble in an unbounded liquid, the well-known Rayleigh prob- 
lem . The calculation started to-develop instabilities in the velocity for a bubble radius 
R-1.411 x lo-', radial velocity Rr483.7, at time tz0.9162. It is well known that this spher- 
ical collapse is unstable . Therefore, the ability of the code to compute a decrease in rad- 
ius by nearly two orders of magnitude with no smoothing techniques applied and within a max- 
imum 1.8% deviation from sphericity is an indication of an excellent performance. The ana- 
iytical result for the radial velocity for the value of the radius indicated above is 
R~487.4, again in excellent agreement with the numerical result. Finally, the analytical 
value of the total collapse time is tc=0.91468, in very good agreement with the computed 
value, which practically corresponds to total collapse. 
As a second test we have considered one of the cases studied by Plesset and chapman1 .Here 
the bubble collapses in the presence of a plane rigid wall, from which it is separated by a 
distance equal to half the initial radius at t=O. In the calculation the plane was simulated 
by introducing an image bubble equal and symmetrically located with respect to the real bub- 
ble. By taking advantage of the symmetry of the problem it is possible to inaintain the num- 
ber of unknowns in the system (9) equal to the number of segments of the real bubble. We 
show in Fig.1 the bubble shapes at selected instants (heavy lines) and the trajectories of 
some of the points on the bubble surface (light lines). It is interesting to observe the very 
strong convergence of streamlines in the jet which evolves in the Inter stuqes of the collapse. 
The shapes of Fig.1 agree with those of Ref.1, but a more stringent comparison is afforded by 
the values of the velocity of the "north pole* of the bubble as computed by us (Fig.2, upper 
line) and as given in Ref.1 (Fig.2, open circles). Although the comparison is good, some minor 
discrepancies exist the origin of which is not clear at the present ti-. It may be noted how- 
ever from the shapes published by Plesset and chapman' that some imprecision affects their 
method near the *north pale* of the bubble, where the tangent fails to be zero as it should. 
The lower line of Fig. 2 shows the time development of the velocity of the *aouth pole* of 
the bubble. 
As a final example we consider the collapse of a system of three equal and equally spaced 
bubbles, with centers on the axis of symmetry. The evolution of the process is shewn in Fig. 
3 for initial bubble spacings of 0.5 (left), 0.2 (center), and 0.1 (right). Because of the 
symmetry of this situation only the upper bubble and the upper half of the middle bubble are 
shown in the figure. The collapse of the central bubble is strongly inhibited by that of the 
other two and its "poles* move very little. This effect of course increases with the proxim- 
ity oC the bubbles. 
Finally, to give an idea of the computati.ona1 requirements we present in Table 1 some in- 
formation on the number of segments used (for all the examples this is the number of segments 
in the first quadrant of the (r,z) plane) . the number of time steps, and the CPU needed for 
the calculation on a UNIVAC 1100/8. Note that we have not tried yet to optimize the method 
with respect to time step size and number of segments used. For instance, it is believed that 
the latter quantity, in the cases of Fig. 3, could be reduced substantially without appreci- 
able loss of accuracy. 
Table 1. Some details on the computation 
Case Number of segments Number of time steps CPU time(mins and secs) 
Rayleigh problem 16 105 1 7 
Figure 1 3 2 79 3 15 
Figure 3, left 4 8 91 8 30 
Figure 3, center 4 8 9 9 9 12 
Figure 3, riqht 4 8 100 9 7 
Conclusions 
The boundary integral method described in the present study has proven to be a very effi- 
cient and accurate tool for the study of non-spherical bubble dynamics. It makes possible 
extensive investigations of non-spherical bubble behavior at a fraction of the effort and 
cost required by other methods. In addition to its value for specific problems, we think that 
this feature is particularly important in the present state of research in this field because 
in our opinion only numerical experiments can help develop the intuition necessary for fur- 
ther progress. 
Acknowledgment 
This study has been supported in part by MinisLero della Pubblica Istruzione of Italy. 
References 
1. Plesset,M.S., and Chapman,R.B., "Collapse of an initially spherical vapour cavity in 
the neighborhhod of a solid boundary," J. Fluid Mech., Vol. 47, pp. 283-290. 1971. 
2. Chapman,R.B., and Plesset,M.S., "Nonlinear effects in the collapse of a nearly spher- 
ical cavity in a liquid," J. Basic En%., Vol. 94, pp. 142-145. 1972. 
3. Iva~iy,R.D., and Hamnitt,F.G.,"Cavitation bubble collapse in viscous compressible liq- 
uids-Numerical analysis," J. Basic Enq., Vol. 87, pp. 977-985. 1965. 
4.  Prosperetti,A., "On the dynamics of non-spherical bubbles," in Lauterborn,W. (Ed.) 
Cavitation and Inhomogeneities in Underwater Acoustics, Springer-Verleg 1980, pp. 13-22. 
5. Jaswon,M.A., and Sym,G.I., Integral Equation Methods in Potential Theory and Elasto- 
statics, Academic Press 1977. 
6. Bevir,M.K., and Pielding,P.J.,"Nunerical solution of incompressible bubble collapse 
with jetting," in Ockendon,J.R., and Hodgkins,W.R. (Eds.) Moving Boundary Problems in Heat 
Flow and Diffusion, Clarendon Press 1975. 
7. Blake,J.R., ~ n d  Gibson,D.C.,*Growth and collapse of a vapour cavity near a free sur- 
face," J. Fluid Mech., in preen. 
8. Hasting8,C.. Approxii~tione for Jiqital Computers, Princeton University Press 1955. 
t 
Figurc :. Collapse of an initially spherical bubble in tl e neighborhood of a solid plane 
:- l fro- which it is separated by 0.5 times the initial radius. The shapes shown (heavy lines) 
P 3.e at t=O, ta0.950, t11.010, and t=l.C33. The collapse was completed during the time step 
3 after the last one shown. Light lines are particle trajectories. 
t 
F
ig
ur
e 
2.
 
Co
mp
ar
is
on
 b
et
we
en
 
t
he
 v
e
lo
ci
ty
 o
f 
t
he
 "
n
o
r
t
h 
po
le
" 
o
f 
t
he
 b
ub
bl
e 
o
f 
F
ig
. 
1 
(u
pp
er
 l
in
e)
 a
n
d 
t
he
 v
a
lu
es
 g
iv
en
 b
y 
P
le
ss
et
 
a
n
d 
c
ha
pm
an
l(
op
en
 c
ir
cl
es
).
 T
he
 
lo
we
r 
li
ne
 s
ho
ws
 t
he
 
v
e
lo
ci
ty
 o
f
 t
he
 "
s
o
u
t
h 
po
le
".
 
F
ig
u
re
 3
. 
C
o
ll
ap
se
 o
f 
th
re
e 
e
q
u
al
 a
n
d 
e
q
u
al
ly
 s
pa
ce
d 
bu
bb
le
s 
w
it
h
 c
o
ll
in
e
a
r 
c
e
n
te
rs
. 
T
he
 p
ro
ce
ss
 
is
 s
yn
an
et
ri
c 
a
bo
ut
 t
h
e 
z=
O
 p
la
n
e.
 
T
he
 i
n
it
ia
l 
bu
bb
le
 s
p
ac
in
g
s 
a
r
e
 0
.5
 
(l
ef
t)
, 
0.
2 
(c
en
te
r)
, 
a
n
d 
0.
1 
(r
ig
ht
) 
ti
m
es
 t
h
e
 i
n
it
ia
l 
r
a
d
iu
s.
 
T
he
 s
h
ap
es
 s
ho
w
n 
a
r
e
 
fo
r:
 A
, 
t=
O;
 
B
, 
t=
0.
95
0;
 
C
, 
tz
1
.0
63
; 
D
,
 
t=
1
 .I
0
5 
(l
ef
t 
a
n
d 
r
ig
h
t)
, 
t=
1
 .I
03
 
(c
en
te
r)
; 
E,
 
t=
1.
13
0;
 
F,
 
t=
1.
13
5.
 


A numerical method for the dynamics of non-spherical cavitation bubbles

Abstract

Similar works

Full text

Available Versions

NASA Technical Reports Server