Liquid compressibility effects during the collapse of a single

cavitating bubble
D. Fustera)
Division of Engineering and Applied Science, California Institute of Technology, Pasadena, California 91125

C. Dopazo and G. Hauke

Laboratorio de Investigacion en Tecnologias de la Combustion—Area de Mecanica de Fluidos, Centro
Politecnico Superior, Universidad de Zaragoza, C/Maria de Luna 3, 50018 Zaragoza, Spain

(Received 23 December 2009; revised 17 September 2010; accepted 22 September 2010)

The effect of liquid compressibility on the dynamics of a single, spherical cavitating bubble is studied.
While it is known that compressibility damps the amplitude of bubble rebounds, the extent to which
this effect is accurately captured by weakly compressible versions of the Rayleigh–Plesset equation is
unclear. To clarify this issue, partial differential equations governing conservation of mass, momen-
tum, and energy are numerically solved both inside the bubble and in the surrounding compressible
liquid. Radiated pressure waves originating at the unsteady bubble interface are directly captured.
Results obtained with Rayleigh–Plesset type equations accounting for compressibility effects, pro-
posed by Keller and Miksis [J. Acoust. Soc. Am. 68, 628–633 (1980)], Gilmore, and Tomita and
Shima [Bull. JSME 20, 1453–1460 (1977)], are compared with those resulting from the full model.
For strong collapses, the solution of the latter reveals that an important part of the energy concentrated
during the collapse is used to generate an outgoing pressure wave. For the examples considered in
this research, peak pressures are larger than those predicted by Rayleigh–Plesset type equations,
whereas the amplitudes of the rebounds are smaller. V C 2011 Acoustical Society of America.

[DOI: 10.1121/1.3502464]
PACS number(s): 43.35.Ei, 43.25.Yw, 43.35.Hl, 43.25.Gf [DLM] Pages: 122–131

I. INTRODUCTION strong resemblance to the waterhammer phenomenon in a

duct: As the liquid flow is halted by the abrupt closing of a
Different processes taking place during the collapse of
valve, pressure waves propagate upstream, reflect at the duct
cavitating bubbles have attracted the attention of the scientific
inlet, travel downstream to the valve, this sequence repeating
and industrial communities for many applications. At implo-
hereafter. The role of the valve is played in the present case
sion, pressures of the order of thousands of atmospheres and
by the bubble interface, which opposes the inward liquid flow.
temperatures of thousands of degrees are reached inside the
Experiments on single bubbles do often take place inside
bubble, promoting chemical and physical processes not only
spherical containers. A transducer, with a well defined freque-
inside the bubble but also in the surrounding liquid. The role
ncy, causes the vibration of the outer vessel wall. In the ab-
of the liquid on those processes has been widely studied by
sence of a bubble, standing waves are generated. A single
several authors. Besant1 found the radial dependence in an
bubble of radius R0, much smaller than the forcing wave-
infinite mass of an incompressible liquid of the instantaneous
length, can be made to oscillate under the externally imposed
pressure alteration with respect to its constant value very far
pressure field.
away from a spherical cavity as it is filled up. Rayleigh2 esti-
The Rayleigh–Plesset (RP) equation, and its weakly com-
mated that before complete collapse the pressure near the
pressible generalizations, such as the Gilmore3 and Keller–
boundary becomes very great. Even at the implosion of vapor
Miksis (KM) equations,4 have been widely used to investigate
and/or gas containing cavities, large pressures can also be
single bubble dynamics.5–11 The validity of different models
attained under some conditions within the bubble, and also in
based on the RP equations accounting for compressibility
the liquid. As a bubble collapses, its interface moves inward
effects has been discussed.12–15 Moreover, these studies are
at high velocities and this motion produces an intense com-
focused on the bubble interior behavior. Analysis of the liquid
pression of the gas/vapor mixture which, simultaneously,
dynamics during the bubble collapse in scarce.16,17 The works
increases its temperature, inducing dissociation and chemical
of Prosperetti and Lezzi,18,19 theoretically quantifying the
reactions. For example, near the implosion of a 10 lm diame-
order of accuracy of various RP-type equations as a function
ter bubble, the kinetic energy of the surrounding liquid, when
of the Mach number, are relevant examples.
the bubble compresses at 100 m/s is on the order of 108 and
Johnson and Colonius17,20 have developed numerical
107 J, which must be brought to a stop through a pressure
methods in order to capture pressure waves generated by
increment at the cavity interface between hundreds and thou-
non-spherical bubble implosion in the liquid. However, a nu-
sands of atmospheres, respectively. This situation bears a
merical study of the validity of RP-type equations to predict
strong spherical bubble collapses is missing.
a)
Author to whom correspondence should be addressed. Electronic mail: Pressure waves emitted in single bubble sonolumi-
fuster@caltech.edu nescence (SBSL) have been experimentally measured by

122 J. Acoust. Soc. Am. 129 (1), January 2011 0001-4966/2011/129(1)/122/10/$30.00 C 2011 Acoustical Society of America
V

Downloaded 14 Mar 2011 to 131.215.220.185. Redistribution subject to ASA license or copyright; see http://asadl.org/journals/doc/ASALIB-home/info/terms.jsp
  
Wang et al.21 and Pecha and Gompf,22 but the lack of realis- DTb Dpb 1 @ @Tb
tic models prevents establishing a correlation between the qb cp;b ¼ þ 2 r 2 kb
Dt Dt r @r @r
pressure values far away from the bubble and those encoun-  2
4 @ub ub
tered inside the bubble during the implosion. Moreover, þ lb  ; (6)
in multi-bubble systems pressure waves could interact with 3 @r r
the bubble cloud;23 although they should not significantly
alter the dynamics of the individual surrounding bubbles, R0
qb ¼ qb Tb ; (7)
this interaction may, under some conditions, generate in the W
liquid strong shock waves propagating inward to the cloud
center.24,25 where the subscript b denotes bubble variables and properties,
In order to quantify the accuracy of widely used RP-type R0 is the universal perfect gas constant, and W is the averaged
equations, a new model is formulated, improving the predic- molecular weight of the gas. Equation (7) is applied to both
tions of the dynamic behavior of the compressible liquid condensable and non-condensable gases, and, thus, the water
around an externally forced single bubble. In Sec. II the gov- vapor is considered as one of the components.
erning equations, boundary, and initial conditions are pre- The following boundary conditions are used. At the
sented; the numerical method used to integrate the system is bubble center (r ¼ 0), under spherical symmetry, derivatives
also described. Section III discusses the numerical results of of the gas variables are set to zero,
the present model and compares them with those obtained via
RP-type equations. Conclusions are detailed in Sec. IV. @qb @Tb @pb
¼ ¼ ¼ 0: (8)
@r @r @r
II. PROPOSED MODEL
For a compressible Newtonian pure liquid the governing The velocity at r ¼ 0 must also vanish
equations for continuity, momentum and energy for spheri-
cally symmetric motions are ub ¼ 0: (9)

At infinity, r ¼ R1  R0, the liquid temperature far away

1 Dpl DTl 1 @r2 ul
 bl ¼  ; (1) from the bubble is taken as a constant
ql c2l Dt Dt r 2 @r
  Tl ¼ T1 ðconstantÞ: (10)
Dul @pl 4 @ 2 @ul
ql ¼ þ r ll
Dt @r 3r 2 @r @r
For the forcing acoustic wave, either oscillating pressure or
8 ll ul 4 ul @ll
  ; (2) velocity far away from the bubble can be specified. In this
3 r2 3 r @r work a sinusoidal pressure is imposed
 
DTl Dpl 1 @ @Tl
ql cp ¼ bl Tl þ 2 r 2 kl pl ðr ¼ R1 ; tÞ ¼ pl;1 þ Dp1 sinð2pftÞ; (11)
Dt Dt r @r @r
 2
4 @ul ul where Dp1 is the wave amplitude induced, for example, by
þ ll  ; (3)
3 @r r the pressure transducer and pl,1 is the hydrostatic pressure.
At the interface [r ¼ R(t)], mass transfer is neglected,
where the subscript l denotes the liquid phase, Dv=Dt repre- limiting the strict model applicability to forcing pressure
sents the substantial derivative of any arbitrary variable v, wave frequencies above 100 kHz, where RP models have
cl is the speed of sound, q is the density, pl is the pressure, bl shown that the bubble dynamics is unaffected by evapora-
is the thermal expansion coefficient, Tl is the temperature, ul tion/condensation.26–29 Viscous stresses and radiation heat
is the radial velocity, ll is the dynamic liquid viscosity coef- transfer are also neglected. Continuity, momentum, and
ficient, cp is the specific heat at constant pressure, jl is the energy transport equations applied to an infinitesimal shell
thermal conductivity of the liquid, and r and t stand for the about the interface yield the following relations, coupling
radial coordinate and time, respectively. bubble and liquid variables,
The transport equations for a mixture of perfect non-
reacting gases, obeying the Navier–Poisson and Fourier law R_ ¼ ul ¼ ub ; (12)
constitutive relations, and neglecting radiative heat transfer, are 2r
pb ðr ¼ RÞ ¼ pl ðr ¼ RÞ þ ; (13)
RðtÞ
Dqb qb @r 2 ub
þ 2 ¼ 0; (4)
Dt r @r @Tb @T1
  kb ¼ kl ; (14)
Dub @pb 4 @ @ub @r @r
qb ¼ þ 2 r 2 lb
Dt @r 3r @r @r where R_ is the interface velocity, r is the surface tension,
8 lb ub 4 ub @lb and jb and kl are the gas and liquid thermal conductivity,
  ; (5)
3 r2 3 r @r respectively.

J. Acoust. Soc. Am., Vol. 129, No. 1, January 2011 Fuster et al.: Liquid compressibility in cavitation 123

Downloaded 14 Mar 2011 to 131.215.220.185. Redistribution subject to ASA license or copyright; see http://asadl.org/journals/doc/ASALIB-home/info/terms.jsp
 Continuity of temperature profiles at the interface is also (2) The continuity, momentum, and energy equations in both
assumed phases [Eqs. (1)–(6)] are solved using an atomic layer
epitaxy arbitrarian Eulerian Lagrangian (ALE) stabilized
Tb ¼ Tl : (15) finite element method (FEM).31 Let the spatial domain be
denoted by X and its boundary by C. The domain is sub-
At the initial time, t ¼ t0, the existence in the liquid of divided into nel elements Xe. The variational formulation
cavitation nuclei of a given radius is assumed. These gas is defined for the set of variables Y, according to Hauke
nuclei are considered to be in equilibrium with the surround- and Hughes.32,33 Given a trial solution space SY, a
ing liquid. Thus, the initial bubble radius, R0, is known, and weighting function space V, and the solution at time t,
the initial velocity is zero. The initial temperature of the bub- this formulation yields the solution Y [ SY at time t þ h
ble is equal to that of the liquid (thermal equilibrium), such that for all W [ V,
ð
Tl ðt ¼ t0 Þ ¼ Tb ðt ¼ t0 Þ: (16)  
W  U;t ðYÞW;i Fadv
i ðYÞþW;i Kij Y;j W  S dX
X
The initial bubble pressure is obtained in terms of the initial nel ð
X  
liquid pressure from the continuity of the radial stress þ LT W  sðLY  SÞdX
(mechanical equilibrium) at the interface, using the Young– e¼1 Xe
ð
Laplace equation  
¼ W  Fadv diff
i ðY Þ þ Fi ðY Þ ni dC; (19)
C
2r
pb ðr; t0 Þ ¼ pl ðr ¼ R; t0 Þ þ : (17) where the summation convention is used, Fadv is the ith con-
R0 i
vective flux, Fdiff
i the ith diffusive flux, ni is the ith component
The initial density can be obtained from the equation of state of the normal vector, and S the source term. Kij denotes the
for a mixture of perfect gases given the temperature, pres- diffusion matrices. The first and last integrals constitute the
sure, and composition. Galerkin terms expressed as a function of the variables Y,
written in conservative form to ensure that the weak solution
III. NUMERICAL METHOD is bestowed with the correct Rankine-Hugoniot relations.
The least-squares contribution is written in terms of the
To solve the system of Eqs. (4)–(17) the following set of
differential operators L and LT , which, respectively, are
unknown variables is defined:
given by
   
@ @ @ @
Y ¼ pg ; ug ; Tg ; pl ; ul ; Tl : (18) L ¼ A0 þ Ai  Kij (20)
@t @xi @xi @xj
The wide range of time scales present in the bubble expan- and
sion and implosion processes requires the use of adaptive  
@ @ @ @
stepsize control. In particular, in the numerical code a fifth LT ¼ AT0 þ ATi  KTij ; (21)
@t @xi @xi @xj
order Cash–Karp Runge–Kutta method with local truncation
error monitoring has been utilized to adjust the stepsize and where A0i s are the Euler Jacobians. The stabilizing matrix is
ensure accuracy.30 denoted by s and definitions can be found in Refs. 34 and
The initial boundary conditions determine the initial 35. Here, the diagonal version36–38 has been employed.
boundary value problem. In particular, when variables (Yt) Because transient processes are taken into account, a
are known at some arbitrary time t, the temporal integration temporal integration method is required. The finite element
supplies the solution (Ytþh) at t þ h, where h is the time step. method can be expressed as
The stepsize control is achieved using the solution of a
fifth order Runge–Kutta method and that of a fourth order MY;t þ KY ¼ F: (22)
Runge–Kutta method with the same intermediate points as
those used in the fifth order method. As a result, no extra eval- For practical reasons, an implicit method is used
uations are required. Estimating the error using both solutions,
which theoretically scales as h5, the stepsize is controlled. MY;t þ KtþDt YtþDt ¼ FtþDt : (23)
The Runge–Kutta method is based on the evaluation of
This method provides a greater robustness than explicit
the derivatives in different intermediate points of the tempo-
ones, which translates into larger temporal steps and a sub-
ral stepsize. For each intermediate point, the resolution strat-
stantial saving of central processing unit (CPU) time. The
egy is as follows:
temporal derivative is approximated by
(1) Properties are evaluated at the time in which the value of
Y,t is required. The calculation of the properties includes YtþDt ¼ Yt þ DtY;t (24)
the gas density, which is computed from the obtained
values for the gas pressure and temperature (note that and Eq. (23) can then be rearranged as
instead of the density, the pressure is used as a primitive
variable for the gas). Y;t ¼ ½M þ KDt1 ½F  KYt : (25)

124 J. Acoust. Soc. Am., Vol. 129, No. 1, January 2011 Fuster et al.: Liquid compressibility in cavitation

Downloaded 14 Mar 2011 to 131.215.220.185. Redistribution subject to ASA license or copyright; see http://asadl.org/journals/doc/ASALIB-home/info/terms.jsp
 The main disadvantage of this method is that the matrices M
and K and the vector F should be evaluated at t þ Dt. The val-
ues at t þ Dt are extrapolated with the values of the derivatives
obtained in the previous step (for further details see Ref. 28).
The system of equations represented by Eq. (25) is
decoupled as follows:
(a) The coupled liquid and gas energy equations, Eqs. (3)
and (6), are first solved inside the bubble and the sur-
rounding liquid.
(b) The continuity and momentum equations, Eqs. (1), (2),
(4), and (5), inside the bubble and the surrounding liquid
are also solved using an ALE stabilized FEM method.
As the liquid is considered to be a low compressibility
substance, it is not desirable to segregate the solution
method. The strong coupling among the equations in
both phases forces their solutions to be obtained simulta-
neously, which increases the CPU time. FIG 1. Comparison of the numerical and analytical solution of the ampli-
tude of the pressure wave at r ¼ 0 as a function of R1. Dp1 ¼ 1 atm and
It has been found that the best strategy consists of using the
f ¼ 10 000 Hz.
pressure as the primitive variable in both phases. When pres-
sure is selected in the gas as a primitive variable and the con-
tinuity and momentum equations are solved, a non-linear @ul @pl
term appears in the gas continuity equation ql ¼ : (28)
@t @r
1 Dpg DTg 1 @r2 ug
 bT ¼ 2 ; (26) The velocity is zero at the origin and the liquid pressure is
pg Dt Dt r @r
imposed at a distance R1
where bT ¼ 1/Tg for a perfect gas. In order to avoid this non-
linearity, the term 1/pg has been lagged, being evaluated at r¼0 ul ð0; tÞ ¼ 0;
time t instead of considering it as an unknown variable. r ¼ R1 pl ðR1 ; tÞ ¼ pl;1 þ Dp1 sinðxtÞ: (29)
In this case, there are two unknown variables, the pres-
sure and the velocity, per equation. In the equation of the The standing wave solution is40
node at the interface, liquid and gas velocities are equal
 
[Eq. (12)]. However a pressure jump at the interface appears, Dp1 1 R1 sinðkrÞ
and the liquid pressure there as a function of the gas pressure ul ðr; tÞ ¼ cosðkrÞ 
ql c sinðkR1 Þ r kr
must be obtained from the local force balance [Eq. (13)], in
 cosðxtÞ; (30)
which the viscous stresses have been neglected. The result is
a heptagonal system of equations from which the velocity
1 R1
and pressure temporal derivatives in the liquid and in the gas pl ðr; tÞ ¼ pl;1 þ Dp1 sinðkr Þ
can be worked out. Even though the solution of this type of sinðkR1 Þ r
system is more expensive than that of tridiagonal systems,  sinðxtÞ; (31)
the computational cost is still affordable using appropriate
strategies for band diagonal systems included in Ref. 39,
where x ¼ 2pf and the wavenumber k ¼ x/cl. At r ¼ 0 the
resulting in a much more robust technique.
amplitude of the standing pressure wave, Dpl,0, is
(3) When the complete equations are solved inside and out-
side the bubble the interface velocity is supplied by the Dpl;0 kR1
¼ : (32)
velocity field obtained in the previous step. Dp1 sinðkR1 Þ

IV. NUMERICAL RESULTS Equation (32) implies that Dpl,0/Dp1 ! 1 for

ð2f =cl ÞR1 ! 0 and predicts resonance for (2f/cl)R1 ¼ 1.
A. Standing waves in a pure liquid
The analytical and the numerical solutions for the
In the absence of a bubble, it is possible to derive an an- dimensionless standing pressure wave amplitude at r ¼ 0
alytical solution for standing spherical pressure waves in a in water forced with a frequency of 10 kHz as a function of
spherical flask. The inviscid Euler equation with negligible (2f/cl)R1 are depicted in Fig. 1 (the values of the parameters
inertia terms and the appropriate mass conservation can be are contained in Table I). The oscillating pressure amplitude
written as40 at the boundary has been set equal to 1 atm. Both solutions
are practically identical and this test serves as a validation of
1 @pl 1 @ ðul r2 Þ the code, at least, regarding the propagation of pressure
2
¼ 2 ; (27)
ql c @t r @r waves in the liquid.

J. Acoust. Soc. Am., Vol. 129, No. 1, January 2011 Fuster et al.: Liquid compressibility in cavitation 125

Downloaded 14 Mar 2011 to 131.215.220.185. Redistribution subject to ASA license or copyright; see http://asadl.org/journals/doc/ASALIB-home/info/terms.jsp
 TABLE I. Parameters of the numerical simulations.

ll ql cl f

8  104 kg/(m s) 1000 kg/m3 1480 m/s 10 000 Hz

The standing pressure wave amplitude at the flask center

depends on R1. For that reason, should a bubble at r ¼ 0 be
present, conditions for the most intense implosions depend
not only on the bubble resonance frequency but also on that
of the vessel.

B. Weak bubble oscillations FIG. 2. Maximum bubble radius as a function of the excitation frequency for
Dp1 ¼ 0.1 and 0.5 atm. Results for the complete model and the model using
When a single bubble in a compressible liquid is consid- the RP equation.
ered, no analytical solutions are available and one must
resort to numerical integration of the conservation equations.
sffiffiffiffiffiffiffiffiffiffiffi
Solutions provided by RP models should, in principle,
ql cl R_
converge to those of the present model for low amplitudes of R1  R0 ; (35)
the forcing pressure wave, where liquid compressibility Dp1
effects are expected to be negligible. It is thus interesting to
verify the coincidence of the present results with those one can expect a good agreement between RP models and
predicted by RP models, typically governed by the equation the present one.
At low amplitudes of the forcing pressure wave this
3 1  condition is easily satisfied when R1 is chosen to be at a
RR€ þ R_ ¼ pl ðr ¼ R; tÞ  pl;1 ðtÞ ; (33) pressure antinode. In this case, Dp1 is much larger than the
2 ql
pressure disturbance induced by the bubble at R1 and it is
where pl,1 (t) is the pressure at a radial distance far away possible to find out a priori the amplitude of the pressure
from the bubble, liquid compressibility effects being negligi- wave applying Eq. (32). Thus, attention is first focused on
ble. Due to the low amplitudes tested in this section, all small excitation cases. Figure 2 depicts the values of Rmax/R0
RP-type equations should tend to the same solution. predicted by both the standard RP model and the present
For a comparison among the different models, pl,1 (t) in one, for an air bubble with an initial radius of 10 lm
Eq. (33) and other RP-type equations can, in principle, be immersed in water forced ultrasonically with Dp1 equal to
approximated by the solution of Eq. (31) evaluated at r ¼ 0. 0.1 and 0.5 atm, respectively. The natural frequency in the
Except for very high frequencies, the bubble radius is much linear regime,41
smaller than a wavelength and therefore, the pressure at  

1=2
r ¼ R1, required as a boundary condition in the new model, 3kp p1 2 3kp  1 r 8vl
can be obtained from Dpl,0 using Eq. (32). In other words, xp ¼ 2pfres ¼ þ  4 ; (36)
ql R20 ql R30 R0
given the amplitude of the pressure at the flask center, one
can use Eq. (32) to impose the appropriate boundary condi- has been used to make frequencies dimensionless in Fig. 2,
tion to compare models. where vl is the kinematic viscosity of the liquid. This
In the current model, the boundary condition must be expression depends on the polytropic coefficient, kp, which,
imposed as far as possible from the bubble in order to avoid for air, may vary from 1 to 1.4. Thus, the natural frequency
any influence of R1 on the solution. To estimate the mini- (fres) ranges between 290 000 and 340 000 Hz. A numerical
mum radius, R1, for imposing the boundary condition with calculation of the amplitude of the bubble oscillation as
no effect on the bubble radius evolution, one can take Uc,1 a function of the forcing frequency for an amplitude of
¼ Dp1/(qlcl) as a characteristic velocity induced in the liq- Dp1 ¼ 0.1 atm reveals that the natural frequency is around
uid by a forcing pressure wave. Should the liquid in the pres- 3  105 Hz.
ence of an oscillating bubble be assumed incompressible, its As expected, good agreement between the two models is
velocity at R1 would be found for the range of frequencies tested here. The errors of
 2 the RP model always remain below 5% for small amplitudes,
_ Rb which validates the code and also shows the convergence of
Uc;2 ¼R : (34)
R1 the present model to the solution in the incompressible limit.
The new model also allows calculation of the pressure
As the velocity due to the bubble oscillation becomes much field in the liquid surrounding the bubble. As an example,
smaller than the forcing pressure wave induced velocity, the oscillation of a 10 lm air bubble in water forced ultra-
Uc;2  Uc;1 , the influence of the bubble in the modification sonically at the natural frequency of the bubble (3  105 Hz).
of the local velocity field is negligible. This condition The amplitude of the pressure imposed at the boundary of
implies that for the domain is set to generate a standing wave of 1.5 atm

126 J. Acoust. Soc. Am., Vol. 129, No. 1, January 2011 Fuster et al.: Liquid compressibility in cavitation

Downloaded 14 Mar 2011 to 131.215.220.185. Redistribution subject to ASA license or copyright; see http://asadl.org/journals/doc/ASALIB-home/info/terms.jsp
 the radial pressure distribution in the liquid with and without
a bubble for two cycles. The solution without a bubble is
directly computed using Eq. (31). The pressure fields are
obtained after 200 cycles, when standing waves are clearly
established in the flask and the bubble response is periodic.
The two pressure fields contained in Fig. 3 show how the nu-
merical solution obtained with bubble recovers the standing
waves patterns expected from the analytical solution given
by Eq. (31). Due to the spherical symmetry of the problem,
the amplitude of the waves decays with the distance r. The
presence of the bubble significantly modifies the pressure
distributions in the liquid in a distance of the order of the
bubble radius. During the bubble collapse, the pressure in
the liquid is significantly influenced by the presence of the
bubble up to distances of the order of 50 times the bubble ra-
dius. The collapse induces a pressure wave in the liquid
propagating outward from the bubble and quickly decaying
due to geometrical effects.
As the amplitude of the forcing pressure wave becomes
larger, obtaining a periodic solution for the bubble response
becomes more difficult. This is due to the appearance of
long-lived transients in the numerical solutions, compounded
by the interference of the driving pressure with reflections at
the edge of the computational domain of pressure waves
radiated by the bubble itself. This fact prevents a rigorous
comparison of the model results with typical experimental
setups like those of Wang et al.42

C. Validity of approaches based on the RP equation

during strong collapses
The analysis of strong implosions, where the validity of
the RP equation is questionable, is considered in this section.
Pressures in the liquid region near to the bubble can reach
values of the order of thousands of atmospheres for intense
bubble collapses and then compressibility effects in the liq-
uid become important.17,22,42 In those situations, peak pres-
sures can be significantly influenced by other physical
mechanisms like, for example, mass transfer, chemical reac-
tions, and ionization. A rigorous comparison with experi-
mental data would demand the inclusion of all possible
relevant effects in the model. In this section, the accuracy of
the RP-type equations is quantified under the assumption
that some of these effects are negligible. This fact restricts
the applicability of the current model and impedes a direct
comparison with experiments. However, it still allows a
quantitative estimation of compressibility effects as a func-
FIG. 3. Temporal evolution of the bubble radius evolution (up) and radial tion of the peak pressure at collapse.
pressure distribution for a 10 lm air bubble in water excited with a fre- In particular, the following equations are tested.
quency of 3  105 Hz in an spherical flask with bubble (middle) and with-
out bubble (bottom). The forcing pressure is set to generate a wave
The RP equation with the modification proposed by
amplitude of 1.5 atm at the flask center when there is no bubble. The cur- KM,4
rent model allows capturing pressure waves emitted during the implosion
as well as the pressure disturbances induced by the bubble in the surround-    
ing liquid. R_ 3 R_ _ 2
1 RR€ þ 1 R
cl 2 3cl
 
1 R_  
at the flask center when there is no bubble. The value of ¼ 1þ pl ðr ¼ RÞ  pl;1 ðtÞ
ql cl
the perturbation induced at the flask wall is obtained from 
RðtÞ @pl ðr; tÞ
Eq. (32). The radius of the flask is R1 ¼ 0.01 m. Figure 3 þ ; (37)
depicts the temporal evolution of the bubble radius and of ql cl @t r¼RðtÞ

J. Acoust. Soc. Am., Vol. 129, No. 1, January 2011 Fuster et al.: Liquid compressibility in cavitation 127

Downloaded 14 Mar 2011 to 131.215.220.185. Redistribution subject to ASA license or copyright; see http://asadl.org/journals/doc/ASALIB-home/info/terms.jsp
 the Gilmore equation In order to compare the different models under conditions
    where compressibility effects become relevant, the following
R_ € 3 R_ _ 2 test is considered: A bubble with a given initial pressure,
1 RR þ 1 R
cl 2 3cl pb,0, at the reference temperature and zero velocity is ini-
     tially assumed to be present in a liquid whose pressure at in-
1 R_ R_ RðtÞ @H 
¼ 1þ Hþ 1 ; (38) finity is maintained constant and equal to
q l c l c q c @t 
l l l r¼RðtÞ
pl;1 ¼ pl;0 þ Dp1 : (41)
where After a certain time, the bubble compresses to a new equilib-
ð pl ðr¼RÞ rium pressure. This case can be an approximation of the con-
dp
H¼ ; (39) ditions when the bubble attains its maximum radius in
pl;1 q
SBSL; it can also be representative of the interaction
and the expression derived by Tomita and Shima,43 proven between a strong pressure wave and a gas bubble. In practi-
by Lezzi and Prosperetti19 to be second order accurate, cal applications, when violent collapses are produced
Dp1  pb;0 and, therefore, one can take pl;1  Dp1 .
 2 !  ! The initial pressure difference between the bubble and
R_ 23 R_ 3 4 R_ 7 R_ 2 2
12 þ RR€ þ 1 þ R_ the liquid at infinity produces an acceleration of the interface
cl 10 cl 2 3 cl 5 cl
 which, assuming an incompressible liquid, can be obtained
pl ðr ¼ RÞ  pl;1 ðtÞ 3 pl ðr ¼ RÞ  pl;1 ðtÞ from Eq. (33)
¼ 1
ql 2 ql c2l 1  
 2 !    R€ðt ¼ 0Þ ¼ pl;0 ðr ¼ R0 Þ  pl;1
R0 ql
1 R_ R_ RðtÞ @pl ðr; tÞ
þ þ 12 (40) 1
2 cl cl ql cl @t r¼RðtÞ ¼ Dp1 : (42)
R0 ql

FIG. 4. (Color online) Bubble radius

and bubble pressure (at r ¼ 0) evolu-
tions as functions of dimensionless
time for the complete model, the
standard RP equation and the RP-type
equations with the compressibility
corrections suggested by KM, Gil-
more and Tomita and Shima. The
ratios between the liquid pressure at
infinity and the initial bubble pressure,
pb,0 are 10 (top) and 100 (center). At
the bottom, a zoom of the implosion
for a ratio equal to 100 is shown. The
model proposed by Tomita and Shima
displays the more accurate behavior
before the collapse.

128 J. Acoust. Soc. Am., Vol. 129, No. 1, January 2011 Fuster et al.: Liquid compressibility in cavitation

Downloaded 14 Mar 2011 to 131.215.220.185. Redistribution subject to ASA license or copyright; see http://asadl.org/journals/doc/ASALIB-home/info/terms.jsp
 Integrating the momentum equation in the liquid from collapse and during a very short time, the corrections intro-
the bubble interface to some arbitrary distance r, the follow- duced by Tomita and Shima do not properly capture the bub-
ing expression for the initial pressure field is readily obtained: ble radius evolution. The reason of this anomalous behavior
can be attributed to the fact that the accuracy of the approxi-
  
mations introduced by the corrections are only proved to
pl ðr; t ¼ 0Þ ¼ pb;0 þ
2r € 0 1  R0
 ql RR : (43)
R0 r work for low Mach numbers. However, during the quick
expansion produced after the intense collapse, the Mach
The Rayleigh collapse time,2 defined by41 number becomes of the order of 1 and the analysis presented
by Lezzi et al. does not apply. For that reason, although the
rffiffiffiffiffiffiffiffi equation proposed by Tomita and Shima provides, overall,
ql
tc ¼ 0:915 R0 ; (44) the most accurate bubble radius evolution, appreciable errors
pl;1
are encountered in the R–t curves during the first instants of
can be taken as a convenient characteristic implosion time of the expansion produced after the collapse. Compared to the
a bubble in an incompressible liquid. corrected formulas, the present model predicts larger pres-
To test the validity of the models, an air bubble of sures during the collapse and lower amplitude of the
10 lm at pb,0 ¼ 105 Pa is introduced in water. The value of rebound. That is, the corrections introduced in the RP equa-
pl,1, is used to control the bubble implosion intensity (Dp1 tion seem to overestimate the energy dissipation before the
¼ pl,1  pb,0). collapse, whereas they tend to underestimate it after the bub-
The bubble radius and pressure (at r ¼ 0) evolutions as ble reaches its minimum radius. For intense collapses, the
functions of the dimensionless time t/tc are shown in Fig. 4 large pressures reached at the interface compress the liquid,
for different bubble implosion intensities. Figure 5 depicts absorbing in that process a significant fraction of the energy
the peak pressures at the bubble center reached during the concentrated during the implosion. When the bubble starts to
bubble collapse for each of the models tested in this work. expand, part of that energy is invested in generating pressure
The larger pl,1, the stronger the implosion. As can be seen, waves in the liquid which are emitted from the bubble inter-
the Rayleigh time correctly scales the implosion occurrence face outward at every rebound. This wave is depicted in
for the different cases. From Figs. 4 and 5, one might con- Fig. 6 for the case of Dpl,1/pb,0 ¼ 100. In general, it can be
clude that the RP-type models capture relatively well the concluded that the energy absorbed by the outgoing wave is
bubble radius and pressure evolutions for implosions with not correctly captured by the first and second order models
peak pressures below 1000 atm. For more intense collapses, tested here. Pressure values of the order of those attained
the corrected RP-type equations only predict the order of inside the bubbles are reached at distances of the order of the
magnitude of the pressures reached during the collapse. As bubble radius; the effect of these pressure waves is still appa-
expected, the correction suggested by Tomita and Shima rent up to characteristic distances of the order of several
turns out to be the most accurate approximation,19 although times Rb.
the relative error for the most intense collapse tested here Thus, liquid compressibility effects are crucial in order
is still of the order of 100%. Only right after the bubble to get accurate estimations of the internal pressures reached
inside the bubble, as well as the amplitude of the rebounds
produced thereafter. During intense bubble collapses, the
new model yields pressure values up to two orders of magni-
tude lower than those predicted under the incompressible hy-
pothesis (standard RP equation). The modifications proposed

FIG. 5. (Color online) Peak pressures at the bubble center at the collapse
obtained with the five models compared in this work for different ratios of
the initial bubble pressure to the pressure at infinity. Models based on the
RP equation accurately predict the peak pressures for values of the peak FIG. 6. Liquid pressure profiles after the bubble collapse every 4 ns. For
pressures below 1000 atm. For extremely violent implosions, RP models large pressures, the wave dissipates energy and the amplitude decay is larger
dramatically fail predicting the peak pressures. The modified RP-type mod- than 1/r. As the wave propagates outward from the bubble, its amplitude
els studied in this work significantly improve predictions, especially that decreases, dissipation becomes less important and the shock wave amplitude
proposed by Tomita and Shima, but errors are still of the order of 100%. is inversely proportional to the distance.

J. Acoust. Soc. Am., Vol. 129, No. 1, January 2011 Fuster et al.: Liquid compressibility in cavitation 129

Downloaded 14 Mar 2011 to 131.215.220.185. Redistribution subject to ASA license or copyright; see http://asadl.org/journals/doc/ASALIB-home/info/terms.jsp
 by KM, Gilmore, and Tomita and Shima, to account for liq- a previous version of this manuscript. This work is a part of
uid compressibility effects, significantly improve the stand- the Ph.D. dissertation of D. Fuster28 and has been partially
ard RP predictions; however, errors are still large for supported by Ministerio de Ciencia y Tecnologia of Spain,
extremely violent implosions. Thus, the complete solution of under Grant No. CTM2004-061 84-C02-02.
the liquid continuity and momentum equations in the liquid
is compulsory if accurate values are required. As a first
1
approximation however, the modified versions of the RP W. Besant, A Treatise on Hydrostatics and Hydrodynamics (Deighton,
Bell, London, 1859), p. 228.
equation can be used to obtain an approximate value of pres- 2
L. Rayleigh, “On the pressure developed in a liquid during the collapse of
sures reached during the collapse. a spherical cavity,” Philos. Mag. 34, 94 (1917).
3
F. R. Gilmore, “The growth of collapse of a spherical bubble in viscous
V. CONCLUSIONS compressible liquid,” Technical Report No. 26-4, California Institute of
Technology (1952).
4
A spherically symmetric one-dimensional model of sin- J. Keller and M. Miksis, “Bubble oscillations of large amplitude,”
J. Acoust. Soc. Am. 68, 628–633 (1980).
gle bubble dynamics has been presented. The complete 5
M. Plesset and A. Prosperetti, “Bubble dynamics and cavitation,” Annu.
model solves the mass, momentum, and energy equations Rev. Fluid Mech. 9, 145 (1977).
6
both inside the bubble and in the surrounding liquid. The A. Prosperetti, “Bubble phenomena in sound fields: Part one,” Ultrasonics
model directly incorporates liquid compressibility effects, 22, 69–77 (1984).
7
A. Prosperetti, “A new mechanism for sonoluminescence,” J. Acoust. Soc.
allowing the study of physical phenomena where small var- Am. 101, 2003–2007 (1997).
iations of the liquid density may play an important role. The 8
R. Toegel and D. Lohse, “Phase diagrams for sonoluminescing bubbles:
new model has been numerically validated by comparison A comparison between experiment and theory,” J. Chem. Phys. 118,
with the analytical solution of an ultrasonic field generating 1863–1875 (2003).
9
Y. An and C. Ying, “Model of single bubble sonoluminescence,” Phys.
standing pressure waves within a spherical flask containing a Rev. E 71, 1–12 (2005).
compressible liquid without a bubble. 10
K. Yasui, T. Tuziuti, M. Sivakumar, and Y. Iida, “Theoretical study of
Predictions of the bubble radius amplitude as a function single-bubble sonochemistry,” J. Chem. Phys. 122, 1–12 (2005).
11
G. Hauke, D. Fuster, and C. Dopazo, “Dynamics of a single cavitating and
of frequency, obtained with both the standard RP model and
reacting bubble,” Phys. Rev. E 75, 1–14 (2007).
the present one, have also been compared for low amplitude 12
W. Moss, “Understanding the periodic driving pressure in the Rayleigh–
oscillations; a good agreement is apparent below, above, and Plesset equation,” J. Acoust. Soc. Am. 101, 1187–1190 (1997).
13
at the bubble natural frequency, showing that liquid com- S. Hilgenfeldt, M. Brenner, S. Grossmann, and D. Lohse, “Analysis of the
Rayleigh–Plesset dynamics for sonoluminescing bubbles,” J. Fluid Mech.
pressibility effects can be neglected in order to accurately 365, 171–204 (1998).
predict the bubble dynamics. 14
H. Lin, B. D. Storey, and A. J. Szeri, “Inertially driven in homogeneities
For very intense collapses, however, where pressures at in violently collapsing bubbles: The validity of the Rayleigh–Plesset equa-
the bubble interface reach values above 1000 atm, significant tion,” J. Fluid Mech. 452, 145–162 (2002).
15
A. Preston, T. Colonius, and C. Brennen, “A reduced order model of diffu-
differences are observed between the predictions of the com- sive effects on the dynamics of bubbles,” Phys. Fluids 19, 1–19 (2007).
plete model and those based on the incompressible RP equa- 16
R. Nigmatulin, I. Akhatov, N. Vakhitova, and R. Lahey, “On the forced
tion. Predicted peak pressure values are up to two orders of oscillations of a small gas bubble in a spherical liquid-filled flask,” J. Fluid
magnitude smaller than those calculated with the RP equa- Mech. 414, 47–73 (2000).
17
E. Johnson and T. Colonius, “Compressible multicomponent flow calcula-
tion. The corrections suggested by KM, Gilmore, or Tomita tions and shock-wave interactions,” in Sixth International Symposium on
and Shima significantly reduce the error to values around Cavitation (2006).
100%. As suggested by Lezzi and Prosperetti,19 the most 18
A. Prosperetti and A. Lezzi, “Bubble dynamics in a compressible liquid.
accurate predictions are provided by the correction proposed Part 1. First-order theory,” J. Fluid Mech. 168, 457–478 (1986).
19
A. Lezzi and A. Prosperetti, “Bubble dynamics in a compressible liquid.
by Tomita and Shima.43 Part 2. Second-order theory,” J. Fluid Mech. 185, 289–321 (1987).
The rigorous solution of the continuity and momentum 20
E. Johnsen and T. Colonius, “Implementation of WENO schemes in
equations in the liquid yields peak pressures larger than those compressible multi-component flow problems,” J. Comput. Phys. 219,
715–732 (2006).
obtained by RP-type models corrected to account for com- 21
Z. Wang, R. Pecha, B. Gompf, and W. Eisenmenger, “Single bubble sono-
pressibility; although pressures are higher, the expansion ra- luminescence: Investigations of the emitted pressure wave with a fiber
dius reached during the first rebound is clearly smaller. Large optic probe hydrophone,” Phys. Rev. E 59, 1777 (1999).
22
pressures attained inside the bubble during the collapse gener- R. Pecha and B. Gompf, “Microimplosions: Cavitation collapse and shock wave
emission on a nanosecond time scale,” Phys. Rev. Lett. 84, 1328–1330 (2000).
ate pressure waves in the liquid which propagate outward. 23
W. Lauterborn, T. Kurz, R. Geisler, D. Schanz, and O. Lindau, “Acoustic
The effect of these propagation waves does not seem to be cavitation, bubble dynamics and sonoluminiescence,” Ultrason. Sono-
correctly captured by the first and second order corrections. chem. 14, 484–491 (2007).
24
The inclusion of models accounting for effects missing G. Reisman, Y. Wang, and C. Brennen, “Observations of shock waves in
cloud cavitation,” J. Fluid. Mech. 344, 255 (1998).
in the current one, like mass transfer or chemical reactions, 25
Y. Wang, “Shock waves in bubbly cavitating flows,” Ph.D. thesis, Califor-
in future work should provide a more realistic correlation nia Institute of Technology, Pasadena, CA, 1996.
26
between pressures measured far from the bubble and those A. Colussi, L. Weavers, and M. Hoffmann, “Chemical bubble dynamics and
reached inside the bubble in SBSL experiments. quantitative sonochemistry,” J. Phys. Chem. A 102, 6927–6934 (1998).
27
G. Puente and F. Bonetto, “Proposed method to estimate the liquid-vapor
accommodation coefficient based on experimental sonoluminescence
ACKNOWLEDGMENTS data,” Phys. Rev. E 71 (2005).
28
D. Fuster, “Modeling and numerical simulation of the dynamics of liquid-
The authors would like to acknowledge the comments bubble cavitating systems with chemical reaction processes,” Ph.D. thesis,
and suggestions of Jean-Louis Thomas and Tim Colonius on Centro Politecnico Superior, Zaragoza, Spain, 2007.

130 J. Acoust. Soc. Am., Vol. 129, No. 1, January 2011 Fuster et al.: Liquid compressibility in cavitation

Downloaded 14 Mar 2011 to 131.215.220.185. Redistribution subject to ASA license or copyright; see http://asadl.org/journals/doc/ASALIB-home/info/terms.jsp
 29 36
D. Fuster, G. Hauke, and C. Dopazo, “Influence of accommodation coeffi- G. Hauke, “Simple stabilizing matrices for the computation of compressi-
cient on nonlinear bubble oscillations,” J. Acoust. Soc. Am. 128, 5–10 ble flows in primitive variables,” Comput. Methods Appl. Mech. Eng.
(2010). 190, 6881–6893 (2001).
30 37
W. Press, S. Teulkolsky, W. Vetterling, and B. Flannery, Numerical G. Hauke, A. Landaberea, I. Garmendia, and J. Canales, “A segregated
Recipes in Fortran 77 (Press Syndicate of the University of Cambridge, method for compressible flow computation part I: Isothermal compressible
The Pitt Building, Trumpington Street, Cambridge CB2 1RP, 1992), flows,” Int. J. Numer. Methods Fluids 190, 271–323 (2004).
38
Chap. 16. G. Hauke, A. Landaberea, I. Garmendia, and J. Canales, “A segregated
31
T. Hughes, W. Liu, and T. Zimmerman, “Langrangian-eulerian finite ele- method for compressible flow computation. Part II: General divariant
ment formulation for incompressible viscous flows,” Comput. Methods compressible flows,” Int. J. Numer. Methods Fluids 49, 183–209 (2005).
39
Appl. Mech. Eng. 29, 329–349 (1981). W. Press, S. Teulkolsky, W. Vetterling, and B. Flannery, Numerical
32
G. Hauke and T. Hughes, “A unified approach to compressible and incom- Recipes in For tran 77 (Press Syndicate of the University of Cambridge,
pressible flows,” Comput. Methods Appl. Mech. Eng. 113, 389–395 The Pitt Building, Trumpington Street, Cambridge CB2 1RP, 1992),
(1994). Chap. 2.
33 40
G. Hauke and T. Hughes, “A comparative study of different sets of varia- P. Morse and K. Ingard, Theoretical Acoustics, 927 (Princenton University
bles for solving compressible and incompressible flows,” Comput. Meth- Press, Princeton, NJ, 1986).
41
ods Appl. Mech. Eng. 153, 1–44 (1998). C. Brennen, Cavitation and Bubble Dynamics (Oxford University Press,
34
T. Hughes and M. Mallet, “A new finite element formulation for computa- New York, 1995), p. 254.
42
tional fluid dynamics: III. the generalized streamline operator for multidimen- Z. Wang, R. Pecha, B. Gompf, and W. Eisenmenger, “Single bubble sono-
sional advection-diffusion systems,” Comput. Methods Appl. Mech. Eng. 58, luminescence: Investigations of the emitted pressure wave with a fiber
305–328 (1986). optic probe hydrophone,” Phys. Rev. E 59, 1777–1780 (1999).
35 43
F. Shakib, T. Hughes, and Z. Johan, “A new finite element formulation Y. Tomita and A. Shima, “On the behavior of a spherical bubble and the
for fluid dynamics formulation: X. the compressible Euler and Navier– impulse pressure in a viscous compressible liquid,” Bull. JSME 20, 1453–
Stokes equations,” Comput. Methods Appl. Mech. Eng. 89, 141–219 (1991). 1460 (1977).

J. Acoust. Soc. Am., Vol. 129, No. 1, January 2011 Fuster et al.: Liquid compressibility in cavitation 131

Downloaded 14 Mar 2011 to 131.215.220.185. Redistribution subject to ASA license or copyright; see http://asadl.org/journals/doc/ASALIB-home/info/terms.jsp