10
FLOWS WITH BUBBLE DYNAMICS
10.1 INTRODUCTION
In the last chapter, the analyses were predicated on the existence of an effec-
tive barotropic relation for the homogeneous mixture. Indeed, the construc-
tion of the sonic speed in sections 9.3.1 and 9.3.3 assumes that all the phases
are in dynamic equilibrium at all times. For example, in the case of
bubbles
in liquids, it is assumed that the response of the bubbles to the change in
pressure,
δp
, is an essentially instantaneous change in their volume. In prac-
tice this would only be the case if the typical frequencies experienced by the
bubbles in the flow are very much smaller than the natural frequencies of
the bubbles themselves (see section 4.4.1). Under these circumstances the
bubbles would behave quasistatically and the mixture would be barotropic.
However, there are a number of important contexts in which the bubbles are
not in equilibrium and in which the non-equilibrium effects have important
consequences. One example is the response of a bubbly multiphase mixture
to high frequency excitation. Another is a bubbly cavitating flow where the
non-equilibrium
bubble dynamics lead to shock waves with substantial noise
and damage potential.
Inthischapterwethereforeexaminesomeflowsinwhichthedynamics
of the individual bubbles play an important role. These effects are included
by incorporating the Rayleigh-Plesset equation (Rayleigh 1917, Knapp
et
al.
1970, Brennen 1995) into the global conservation equations for the mul-
tiphase flow. Consequently the mixture no longer behaves barotropically.
Viewing these flows from a different perspective, we note that analyses of
cavitating flows often consist of using a single-phase liquid pressure distri-
bution as input to the Rayleigh-Plesset equation. The result is the history of
the size of individual cavitating bubbles as they progress along a streamline
in the otherwise purely liquid flow. Such an approach entirely neglects the
246
interactive effects that the cavitating bubbles have on themselves and on the
pressure and velocity of the liquid flow. The analysis that follows incorpo-
rates these interactions using the equations for nonbarotropic homogeneous
flow.
10.2 BASIC EQUATIONS
In this chapter it is assumed that the ratio of liquid to vapor density is
sufficiently large so that the volume of liquid evaporated or condensed is
negligible. It is also assumed that bubbles are neither created or destroyed.
Then the appropriate continuity equation is
∂u
i
∂x
i
=
η
(1 +
ηv
)
Dv
Dt
(10.1)
where
η
is the population or number of bubbles per unit volume of liquid
and
v
(
x
i
,t
) is the volume of individual bubbles. The above form of the
continuity equation assumes that
η
is uniform; such would be the case if
the flow originated from a uniform stream of uniform population and if
there were no relative motion between the bubbles and the liquid. Note also
that
α
=
ηv/
(1 +
ηv
) and the mixture density,
ρ
≈
ρ
L
(1
−
α
)=
ρ
L
/
(1 +
ηv
).
This last relation can be used to write the momentum equation 9.2 in terms
of
v
rather than
ρ
:
ρ
L
Du
i
Dt
=
−
(1 +
ηv
)
∂p
∂x
i
(10.2)
The hydrostatic pressure gradient due to gravity has been omitted for sim-
plicity.
Finally the Rayleigh-Plesset equation 4.25 relates the pressure
p
and the
bubble volume,
v
=
4
3
πR
3
:
R
D
2
R
Dt
2
+
3
2
DR
Dt
2
=
p
V
−
p
ρ
L
+
p
Go
ρ
L
R
o
R
3
k
−
2
S
ρ
L
R
−
4
ν
L
R
DR
Dt
(10.3)
where it is assumed that the mass of gas in the bubble remains constant,
p
V
is the vapor pressure,
p
Go
is the partial pressure of non-condensable gas at
some reference moment in time when
R
=
R
o
and
k
is the polytropic index
representing the behavior of the gas.
Equations 10.1, 10.2, and 10.3 can, in theory, be solved to find the un-
knowns
p
(
x
i
,t
),
u
i
(
x
i
,t
), and
v
(
x
i
,t
)(or
R
(
x
i
,t
)) for any bubbly cavitating
flow. In practice the nonlinearities in the Rayleigh-Plesset equation and in
the Lagrangian derivative,
D/Dt
=
∂/∂t
+
u
i
∂/∂x
i
, present serious difficul-
247
ties for all flows except those of the simplest geometry. In the following
sections several such flows are examined in order to illustrate the interactive
effects of bubbles in cavitating flows and the role played by bubble dynamics
in homogeneous flows.
10.3 ACOUSTICS OF BUBBLY MIXTURES
10.3.1 Analysis
One class of phenomena in which bubble dynamics can play an important
role is the acoustics of bubble/liquid mixtures. When the acoustic excitation
frequency approaches the natural frequency of the bubbles, the latter no
longer respond in the quasistatic manner assumed in chapter 9, and both
the propagation speed and the acoustic attenuation are significantly altered.
A review of this subject is given by van Wijngaarden (1972) and we will
include here only a summary of the key results. This class of problems has
the advantage that the magnitude of the perturbations is small so that the
equations of the preceding section can be greatly simplified by linearization.
Hence the pressure,
p
, will be represented by the following sum:
p
= ̄
p
+
Re
$
̃
pe
iωt
%
(10.4)
where ̄
p
is the mean pressure,
ω
is the frequency, and ̃
p
is the small amplitude
pressure perturbation. The response of a bubble will be similarly represented
by a perturbation,
φ
, to its mean radius,
R
o
, such that
R
=
R
o
1+
Re
$
φe
iωt
%
(10.5)
and the linearization will neglect all terms of order
φ
2
or higher.
The literature on the acoustics of dilute bubbly mixtures contains two
complementary analytical approaches. Foldy (1945) and Carstensen and
Foldy (1947) applied the classical acoustical approach and treated the prob-
lem of multiple scattering by randomly distributed point scatterers repre-
senting the bubbles. The medium is assumed to be very dilute (
α
1). The
multiple scattering produces both coherent and incoherent contributions.
The incoherent part is beyond the scope of this text. The coherent part,
which can be represented by equation 10.4, was found to satisfy a wave
equation and yields a dispersion relation for the wavenumber,
κ
,ofplane
waves, that implies a phase velocity,
c
κ
=
ω/κ
, given by (see van Wijngaar-
den 1972)
1
c
2
κ
=
κ
2
ω
2
=
1
c
2
L
+
1
c
2
o
1
−
iδ
d
ω
ω
n
−
ω
2
ω
2
n
−
1
(10.6)
248
Here
c
L
is the sonic speed in the liquid,
c
o
is the sonic speed arising from
equation 9.14 when
αρ
G
(1
−
α
)
ρ
L
,
c
2
o
=
k
̄
p
/
ρ
L
α
(1
−
α
)
(10.7)
ω
n
is the natural frequency of a bubble in an infinite liquid (section 4.4.1),
and
δ
d
is a dissipation coefficient that will be discussed shortly. It follows
from equation 10.6 that scattering from the bubbles makes the wave prop-
agation dispersive since
c
κ
is a function of the frequency,
ω
.
As described by van Wijngaarden (1972) an alternative approach is to
linearize the fluid mechanical equations 10.1, 10.2, and 10.3, neglecting any
terms of order
φ
2
or higher. In the case of plane wave propagation in the
direction
x
(velocity
u
) in a frame of reference relative to the mixture (so
that the mean velocity is zero), the convective terms in the Lagrangian
derivatives,
D/Dt
, are of order
φ
2
and the three governing equations become
∂u
∂x
=
η
(1 +
ηv
)
∂v
∂t
(10.8)
ρ
L
∂u
∂t
=
−
(1 +
ηv
)
∂p
∂x
(10.9)
R
∂
2
R
∂t
2
+
3
2
∂R
∂t
2
=
1
ρ
L
p
V
+
p
Go
R
o
R
3
k
−
p
−
2
S
ρ
L
R
−
4
ν
L
R
∂R
∂t
(10.10)
Assuming for simplicity that the liquid is incompressible (
ρ
L
= constant)
and eliminating two of the three unknown functions from these relations,
one obtains the following equation for any one of the three perturbation
quantities (
Q
=
φ
, ̃
p
,or ̃
u
, the velocity perturbation):
3
α
o
(1
−
α
o
)
∂
2
Q
∂t
2
=
3
kp
Go
ρ
L
−
2
S
ρ
L
R
o
∂
2
Q
∂x
2
+
R
2
o
∂
4
Q
∂x
2
∂t
2
+4
ν
L
∂
3
Q
∂x
2
∂t
(10.11)
where
α
o
is the mean void fraction given by
α
o
=
ηv
o
/
(1 +
ηv
o
). This equa-
tion governing the acoustic perturbations is given by van Wijngaarden,
though we have added the surface tension term. Since the mean state must
be in equilibrium, the mean liquid pressure, ̄
p
, is related to
p
Go
by
̄
p
=
p
V
+
p
Go
−
2
S
R
o
(10.12)
and hence the term in square brackets in equation 10.11 may be written in
249
the alternate forms
3
kp
Go
ρ
L
−
2
S
ρ
L
R
o
=
3
k
ρ
L
( ̄
p
−
p
V
)+
2
S
ρ
L
R
o
(3
k
−
1) =
R
2
o
ω
2
n
(10.13)
This identifies
ω
n
, the natural frequency of a single bubble in an infinite
liquid (see section 4.4.1).
Results for the propagation of a plane wave in the positive
x
direction
are obtained by substituting
q
=
e
−
iκx
in equation 10.11 to produce the
following dispersion relation:
c
2
κ
=
ω
2
κ
2
=
3
k
ρ
L
( ̄
p
−
p
V
)+
2
S
ρ
L
R
o
(3
k
−
1)
+4
iων
L
−
ω
2
R
2
o
3
α
o
(1
−
α
o
)
(10.14)
Note that at the low frequencies for which one would expect quasistatic
bubble behavior (
ω
ω
n
) and in the absence of vapor (
p
V
=0) and sur-
face tension, this reduces to the sonic velocity given by equation 9.14 when
ρ
G
α
ρ
L
(1
−
α
). Furthermore, equation 10.14 may be written as
c
2
κ
=
ω
2
κ
2
=
R
2
o
ω
2
n
3
α
o
(1
−
α
o
)
1+
i
δ
d
ω
ω
n
−
ω
2
ω
2
n
(10.15)
where
δ
d
=4
ν
L
/ω
n
R
2
o
. For the incompressible liquid assumed here this is
identical to equation 10.6 obtained using the Foldy multiple scattering ap-
proach (the difference in sign for the damping term results from using
i
(
ωt
−
κx
) rather than
i
(
κx
−
ωt
) and is inconsequential).
In the above derivation, the only damping mechanism that was explicitly
included was that due to viscous effects on the radial motion of the bubbles.
As Chapman and Plesset (1971) have shown, other damping mechanisms
can affect the volume oscillations of the
bubble; these include the damping
due to temperature gradients caused by evaporation and condensation at
the bubble surface and the radiation of acoustic energy due to compress-
ibility of the liquid. However, Chapman and Plesset (1971) and others have
demonstrated that, to a first approximation, all of these damping contribu-
tions can be included by defining an
effective
damping,
δ
d
, or, equivalently,
an effective liquid viscosity,
μ
e
=
ω
n
R
2
o
δ
d
/
4.
10.3.2 Comparison with experiments
The real and imaginary parts of
κ
as defined by equation 10.15 lead respec-
tively to a sound speed and an attenuation that are both functions of the
frequency of the perturbations. A number of experimental investigations
have been carried out (primarily at very small
α
) to measure the sound
250
Figure 10.1.
Sonic speed for water with air bubbles of mean radius,
R
o
=
0
.
12
mm
, and a void fraction,
α
=0
.
0002, plotted against frequency. The
experimental data of Fox, Curley, and Larson (1955) is plotted along with
the theoretical curve for a mixture with identical
R
o
=0
.
11
mm
bubbles
(dotted line) and with the experimental distribution of sizes (solid line).
These lines use
δ
d
=0
.
5.
Figure 10.2.
Values for the attenuation of sound waves corresponding to
the sonic speed data of figure 10.1. The attenuation in
dB/cm
is given by
8
.
69
Im
{
κ
}
where
κ
is in
cm
−
1
.
251
speed and attenuation in bubbly gas/liquid mixtures. This data is reviewed
by van Wijngaarden (1972) who concentrated on the experiments of Fox,
Curley, and Lawson (1955), Macpherson (1957), and Silberman (1957), in
which the bubble size distribution was more accurately measured and con-
trolled. In general, the comparison between the experimental and theoretical
propagation speeds is good, as illustrated by figure 10.1. One of the primary
experimental difficulties illustrated in both figures 10.1 and 10.2 is that the
results are quite sensitive to the distribution of bubble sizes present in the
mixture. This is caused by the fact that the bubble natural frequency is quite
sensitive to the mean radius (see equation 10.13). Hence a distribution in
the size of the bubbles yields broadening of the peaks in the data of figures
10.1 and 10.2.
Though the propagation speed is fairly well predicted by the theory, the
same cannot be said of the attenuation, and there remain a number of unan-
swered questions in this regard. Using equation 10.15 the theoretical esti-
mate of the damping coefficient,
δ
d
, pertinent to the experiments of Fox,
Curley, and Lawson (1955) is 0.093. But a much greater value of
δ
d
=0
.
5
had to be used in order to produce an analytical line close to the experi-
mental data on attenuation; it is important to note that the empirical value,
δ
d
=0
.
5, has been used for the theoretical results in figure 10.2. On the
other hand, Macpherson (1957) found good agreement between a measured
attenuation corresponding to
δ
d
≈
0
.
08 and the estimated analytical value of
0.079 relevant to his experiments. Similar good agreement was obtained for
both the propagation and attenuation by Silberman (1957). Consequently,
there appear to be some unresolved issues insofar as the attenuation is con-
cerned. Among the effects that were omitted in the above analysis and that
might contribute to the attenuation is the effect of the relative motion of the
bubbles. However, Batchelor (
1969) has concluded that the viscous effects
of translational motion would make a negligible contribution to the total
damping.
Finally, it is important to emphasize that virtually all of the reported
data on attenuation is confined to very small void fractions of the order of
0
.
0005 or less. The reason for this is clear when one evaluates the imaginary
part of
κ
from equation 10.15. At these small void fractions the damping
is proportional to
α
. Consequently, at large void fraction of the order, say,
of 0
.
05, the damping is 100 times greater and therefore more difficult to
measure accurately.
252
Figure 10.3.
Schematic of the flow relative to a bubbly shock wave.
10.4 SHOCK WAVES IN BUBBLY FLOWS
10.4.1 Normal shock wave analysis
The propagation and structure of shock waves in
bubbly cavitating flows
represent a rare circumstance in which fully nonlinear solutions of the gov-
erning equations can be obtained. Shock wave analyses of this kind were in-
vestigated by Campbell and Pitcher (1958), Crespo (1969), Noordzij (1973),
and Noordzij and van Wijngaarden (1974), among others, and for more de-
tail the reader should consult these works. Since this chapter is confined to
flows without significant relative motion, this section will not cover some
of the important effects of relative motion on the structural evolution of
shocks in bubbly liquids. For this the reader is referred to Noordzij and van
Wijngaarden (1974).
Consider a normal shock wave in a coordinate system moving with the
shock so that the flow is steady and the shock stationary (figure 10.3). If
x
and
u
represent a coordinate and the fluid velocity normal to the shock,
then continuity requires
ρu
= constant =
ρ
1
u
1
(10.16)
where
ρ
1
and
u
1
will refer to the mixture density and velocity far upstream of
the shock. Hence
u
1
is also the velocity of propagation of a shock into a mix-
ture with conditions identical to those upstream of the shock. It is assumed
that
ρ
1
≈
ρ
L
(1
−
α
1
)=
ρ
L
/
(1 +
ηv
1
) where the liquid density is considered
constant and
α
1
,
v
1
=
4
3
πR
3
1
,and
η
are the void fraction, individual bubble
volume, and population of the mixture far upstream.
Substituting for
ρ
in the equation of motion and integrating, one also
253
obtains
p
+
ρ
2
1
u
2
1
ρ
= constant =
p
1
+
ρ
1
u
2
1
(10.17)
This expression for the pressure,
p
, may be substituted into the Rayleigh-
Plesset equation using the observation that, for this steady flow,
DR
Dt
=
u
dR
dx
=
u
1
(1 +
ηv
)
(1 +
ηv
1
)
dR
dx
(10.18)
D
2
R
Dt
2
=
u
2
1
(1 +
ηv
)
(1 +
ηv
1
)
2
(1 +
ηv
)
d
2
R
dx
2
+4
πR
2
η
dR
dx
2
(10.19)
where
v
=
4
3
πR
3
has been used for clarity. It follows that the structure of
the flow is determined by solving the following equation for
R
(
x
):
u
2
1
(1 +
ηv
)
2
(1 +
ηv
1
)
2
R
d
2
R
dx
2
+
3
2
u
2
1
(1 + 3
ηv
)(1 +
ηv
)
(1 +
ηv
1
)
2
dR
dx
2
(10.20)
+
2
S
ρ
L
R
+
u
1
(1 +
ηv
)
(1 +
ηv
1
)
4
ν
L
R
dR
dx
=
(
p
B
−
p
1
)
ρ
L
+
η
(
v
−
v
1
)
(1 +
ηv
1
)
2
u
2
1
It will be found that dissipation effects in the bubble dynamics strongly
influence the structure of the shock. Only one dissipative effect, namely that
due to viscous effects (last term on the left-hand side) has been explicitly
included in equation 10.20. However, as discussed in the last section, other
dissipative effects may be incorporated approximately by regarding
ν
L
as a
total
effective
viscosity.
The pressure within the bubble is given by
p
B
=
p
V
+
p
G
1
(
v
1
/
v
)
k
(10.21)
and the equilibrium state far upstream must satisfy
p
V
−
p
1
+
p
G
1
=2
S
/
R
1
(10.22)
Furthermore, if there exists an equilibrium state far downstream of the shock
(this existence will be explored shortly), then it follows from equations 10.20
and 10.21 that the velocity,
u
1
, must be related to the ratio,
R
2
/
R
1
(where
R
2
is the bubble size downstream of the shock), by
u
2
1
=
(1
−
α
2
)
(1
−
α
1
)(
α
1
−
α
2
)
(
p
1
−
p
V
)
ρ
L
R
1
R
2
3
k
−
1
254
Figure 10.4.
Shock speed,
u
1
, as a function of the upstream and down-
stream void fractions,
α
1
and
α
2
, for the particular case (
p
1
−
p
V
)
/ρ
L
=
100
m
2
/sec
2
,2
S/ρ
L
R
1
=0
.
1
m
2
/sec
2
,and
k
=1
.
4. Also shown by the
dotted line is the sonic velocity,
c
1
, under the same upstream conditions.
+
2
S
ρ
L
R
1
R
1
R
2
3
k
−
R
1
R
2
(10.23)
where
α
2
is the void fraction far downstream of the shock and
R
2
R
1
3
=
α
2
(1
−
α
1
)
α
1
(1
−
α
2
)
(10.24)
Hence the
shock velocity,
u
1
, is given by the upstream flow parameters
α
1
,
(
p
1
−
p
V
)
/ρ
L
,and2
S/ρ
L
R
1
, the polytropic index,
k
, and the downstream
void fraction,
α
2
. An example of the dependence of
u
1
on
α
1
and
α
2
is shown
in figure 10.4 for selected values of (
p
1
−
p
V
)
/ρ
L
= 100
m
2
/sec
2
,2
S/ρ
L
R
1
=
0
.
1
m
2
/sec
2
,and
k
=1
.
4. Also displayed by the dotted line in this figure is
the sonic velocity of the mixture (at zero frequency),
c
1
, under the upstream
conditions; it is readily shown that
c
1
is given by
c
2
1
=
1
α
1
(1
−
α
1
)
k
(
p
1
−
p
V
)
ρ
L
+
k
−
1
3
2
S
ρ
L
R
1
(10.25)
Alternatively, the presentation conventional in gas dynamics can be
adopted. Then the upstream Mach number,
u
1
/c
1
, is plotted as a function
of
α
1
and
α
2
. The resulting graphs are functions only of two parameters,
the polytropic index,
k
, and the parameter,
R
1
(
p
1
−
p
V
)
/S
. An example is
255
Figure 10.5.
The upstream Mach number,
u
1
/c
1
, as a function of the
upstream and downstream void fractions,
α
1
and
α
2
,for
k
=1
.
4and
R
1
(
p
1
−
p
V
)
/S
= 200.
included as figure 10.5 in which
k
=1
.
4and
R
1
(
p
1
−
p
V
)
/S
= 200. It should
be noted that a real shock velocity and a real sonic speed can exist even
when the upstream mixture is under tension (
p
1
<p
V
). However, the nu-
merical value of the tension,
p
V
−
p
1
, for which the values are real is limited
to values of the parameter
R
1
(
p
1
−
p
V
)
/
2
S>
−
(1
−
1
/
3
k
)or
−
0
.
762 for
k
=1
.
4. Also note that figure 10.5 does not change much with the parame-
ter,
R
1
(
p
1
−
p
V
)
/S
.
10.4.2 Shock wave structure
Bubble dynamics do not affect the results presented thus far since the speed,
u
1
, depends only on the equilibrium conditions upstream and downstream.
However, the existence and structure of the shock depend on the bubble
dynamic terms in equation 10.20. That equation is more conveniently written
in terms of a radius ratio,
r
=
R/R
1
, and a dimensionless coordinate,
z
=
x/R
1
:
1
−
α
1
+
α
1
r
3
2
r
d
2
r
dz
2
+
3
2
1
−
α
1
+
α
1
r
3
1
−
α
1
+3
α
1
r
3
dr
dz
2
+
1
−
α
1
+
α
1
r
3
4
ν
L
u
1
R
1
1
r
dr
dz
+
α
1
(1
−
α
1
)
1
−
r
3
256
Figure 10.6.
The typical structure of a shock wave in a bubbly mixture
is illustrated by these examples for
α
1
=0
.
3,
k
=1
.
4,
R
1
(
p
1
−
p
V
)
/S
1,
and
u
1
R
1
/ν
L
= 100.
=
1
u
2
1
(
p
1
−
p
V
)
ρ
L
r
−
3
k
−
1
+
2
S
ρ
L
R
1
r
−
3
k
−
r
−
1
(10.26)
It could also be written in terms of the void fraction,
α
,since
r
3
=
α
(1
−
α
)
(1
−
α
1
)
α
1
(10.27)
When examined in conjunction with the expression in equation 10.23 for
u
1
, it is clear that the solution,
r
(
z
)or
α
(
z
), for the structure of the shock
is a function only of
α
1
,
α
2
,
k
,
R
1
(
p
1
−
p
V
)
/S
, and the effective Reynolds
number,
u
1
R
1
/ν
L
, where, as previously mentioned,
ν
L
should incorporate
the various forms of bubble damping.
Equation 10.26 can be readily integrated numerically and typical solu-
tions are presented in figure 10.6 for
α
1
=0
.
3,
k
=1
.
4,
R
1
(
p
1
−
p
V
)
/S
1,
u
1
R
1
/ν
L
= 100, and two downstream volume fractions,
α
2
=0
.
1 and 0.05.
These examples illustrate several important features of the structure of these
shocks. First, the initial collapse is followed by many rebounds and subse-
quent collapses. The decay of these nonlinear oscillations is determined by
the damping or
u
1
R
1
/ν
L
. Though
u
1
R
1
/ν
L
includes an effective kinematic
viscosity to incorporate other contributions to the bubble damping, the value
of
u
1
R
1
/ν
L
chosen for this example is probably smaller than would be rel-
evant in many practical applications, in which we might expect the decay
to be even smaller. It is also valuable to identify the nature of the solution
as the damping is eliminated (
u
1
R
1
/ν
L
→∞
). In this limit the distance be-
tween collapses increases without bound until the structure consists of one
257
Figure 10.7.
The ratio of the ring frequency downstream of a bubbly
mixture shock to the natural frequency of the bubbles far downstream as
a function of the effective damping parameter,
ν
L
/u
1
R
1
,for
α
1
=0
.
3and
various downstream void fractions as indicated.
collapse followed by a downstream asymptotic approach to a void fraction
of
α
1
(
not
α
2
). In other words, no solution in which
α
→
α
2
exists in the
absence of damping.
Another important feature in the structure of these shocks is the typical
interval between the downstream oscillations. This
ringing
will, in practice,
result in acoustic radiation at frequencies corresponding to this interval, and
it is of importance to identify the relationship between this ring frequency
and the natural frequency of the bubbles downstream of the shock. A char-
acteristic ring frequency,
ω
r
, for the shock oscillations can be defined as
ω
r
=2
πu
1
/Δ
x
(10.28)
where Δ
x
is the distance between the first and second bubble collapses. The
natural frequency of the bubbles far downstream of the shock,
ω
2
,isgiven
by (see equation 10.13)
ω
2
2
=
3
k
(
p
2
−
p
V
)
ρ
L
R
2
2
+(3
k
−
1)
2
S
ρ
L
R
3
2
(10.29)
and typical values for the ratio
ω
r
/ω
2
are presented in figure 10.7 for
α
1
=
0
.
3,
k
=1
.
4,
R
1
(
p
1
−
p
V
)
/S
1, and various values of
α
2
. Similar results
were obtained for quite a wide range of values of
α
1
. Therefore note that
the frequency ratio is primarily a function of the damping and that ring
frequencies up to a factor of 10 less than the natural frequency are to be
258
Figure 10.8.
Supersonic bubbly flow past a 20
◦
half-angle wedge at
a Mach number of 4. Flow is from left to right. Photograph taken in
supersonic bubbly flow tunnel (Eddi
ngton 1967) and reproduced with
permission.
expected with typical values of the damping in water. This reduction in the
typical frequency associated with the collective behavior of bubbles presages
the natural frequencies of bubble clouds, that are discussed in the next
section.
10.4.3 Oblique shock waves
While the focus in the preceding two sections has been on normal shock
waves, the analysis can be generalized to cover oblique shocks. Figure 10.8
is a photograph taken in a supersonic bubbly tunnel (Eddington
1967) and
shows a Mach 4 flow past a 20
◦
half-angle wedge. The oblique bow shock
waves are clearly evident and one can also detect some of the structure of
the shocks.
10.5 FINITE BUBBLE CLOUDS
10.5.1 Natural modes of a spherical cloud of bubbles
A second illustrative example of the effect of
bubble dynamics on the be-
havior of a homogeneous bubbly mixture is the study of the dynamics of
a finite cloud of bubbles. One of the earliest investigations of the collective
259
Figure 10.9.
Notation for the analysis of a spherical cloud of bubbles.
dynamics of bubble clouds was the work of van Wijn
gaarden (1964) on the
oscillations of a layer of
bubbles near a wall. Later d’Agostino and Brennen
(1983) investigated the dynamics of a spherical cloud (see also d’Agostino
and Brennen 1989, Omta 1987), and we will choose the latter as a example
of that class of problems with one space dimension in which analytical so-
lutions may be obtained but only after linearization of the Rayleigh-Plesset
equation 10.3.
The geometry of the spherical cloud is shown in figure 10.9. Within the
cloud of radius,
A
(
t
), the population of bubbles per unit
liquid
volume,
η
,is
assumed constant and uniform. The linearization assumes small perturba-
tions of the bubbles from an equ
ilibrium radius,
R
o
:
R
(
r, t
)=
R
o
[1 +
φ
(
r, t
)]
,
|
φ
|
1
(10.30)
We will seek the response of the cloud to a correspondingly small perturba-
tion in the pressure at infinity,
p
∞
(
t
), that is represented by
p
∞
(
t
)=
p
(
∞
,t
)= ̄
p
+
Re
$
̃
pe
iωt
%
(10.31)
where ̄
p
is the mean, uniform pressure and ̃
p
and
ω
are the perturbation
amplitude and frequency, respectively. The solution will relate the pressure,
p
(
r, t
), radial velocity,
u
(
r, t
), void fraction,
α
(
r, t
), and bubble perturbation,
φ
(
r, t
), to ̃
p
. Since the analysis is linear, the response to excitation involving
multiple frequencies can be obtained by Fourier synthesis.
One further restriction is necessary in order to linearize the governing
equations 10.1, 10.2, and 10.3. It is assumed that the mean void fraction in
the cloud,
α
o
, is small so that the term (1 +
ηv
) in equations 10.1 and 10.2
260
is approximately unity. Then these equations become
1
r
2
∂
∂r
r
2
u
=
η
Dv
Dt
(10.32)
Du
Dt
=
∂u
∂t
+
u
∂u
∂r
=
−
1
ρ
∂p
∂r
(10.33)
It is readily shown that the velocity
u
is of order
φ
and hence the convective
component of the material derivative is of order
φ
2
; thus the linearization
implies replacing
D/Dt
by
∂/∂t
.Thentoorder
φ
the Rayleigh-Plesset equa-
tion yields
p
(
r, t
)= ̄
p
−
ρR
2
o
∂
2
φ
∂t
2
+
ω
2
n
φ
;
r<A
(
t
)
(10.34)
where
ω
n
is the natural frequency of an individual bubble if it were alone in
an infinite fluid (equation 10.13). It must be assumed that the bubbles are
in stable equilibrium in the mean state so that
ω
n
is real.
Upon substitution of equations 10.30 and 10.34 into 10.32 and 10.33 and
elimination of
u
(
r, t
) one obtains the following equation for
φ
(
r, t
)inthe
domain
r<A
(
t
):
1
r
2
∂
∂r
r
2
∂
∂r
∂
2
φ
∂t
2
+
ω
2
n
φ
−
4
πηR
o
∂
2
φ
∂t
2
= 0
(10.35)
The incompressible liquid flow outside the cloud,
r
≥
A
(
t
), must have the
standard solution of the form:
u
(
r, t
)=
C
(
t
)
r
2
;
r
≥
A
(
t
)
(10.36)
p
(
r, t
)=
p
∞
(
t
)+
ρ
r
dC
(
t
)
dt
−
ρC
2
2
r
4
;
r
≥
A
(
t
)
(10.37)
where
C
(
t
) is of perturbation order. It follows that, to the first order in
φ
(
r, t
), the continuity of
u
(
r, t
)and
p
(
r, t
) at the interface between the cloud
and the pure liquid leads to the following boundary condition for
φ
(
r, t
):
1+
A
o
∂
∂r
∂
2
φ
∂t
2
+
ω
2
n
φ
r
=
A
o
=
̄
p
−
p
∞
(
t
)
ρR
2
o
(10.38)
The solution of equation 10.35 under the above boundary condition is
φ
(
r, t
)=
−
1
ρR
2
o
Re
̃
p
ω
2
n
−
ω
2
e
iωt
cos
λA
o
sin
λr
λr
;
r<A
o
(10.39)
261
where:
λ
2
=4
πηR
o
ω
2
ω
2
n
−
ω
2
(10.40)
Another possible solution involving (cos
λr
)
/λr
has been eliminated since
φ
(
r, t
) must clearly be finite as
r
→
0. Therefore in the domain
r<A
o
:
R
(
r, t
)=
R
o
−
1
ρR
o
Re
̃
p
ω
2
n
−
ω
2
e
iωt
cos
λA
o
sin
λr
λr
(10.41)
u
(
r, t
)=
1
ρ
Re
i
̃
p
ω
1
r
sin
λr
λr
−
cos
λr
e
iωt
cos
λA
o
(10.42)
p
(
r, t
)= ̄
p
−
Re
̃
p
sin
λr
λr
e
iωt
cos
λA
o
(10.43)
The entire flow has thus been determined in terms of the prescribed quan-
tities
A
o
,R
o
,η,ω
,and ̃
p
.
Note first that the cloud has a number of natural frequencies and modes
of oscillation. From equation 10.39 it follows that, if ̃
p
were zero, oscillations
would only occur if
ω
=
ω
n
or
λA
o
=(2
m
−
1)
π
2
,m
=0
,
±
2
...
(10.44)
and, therefore, using equation 10.40 for
λ
, the natural frequencies,
ω
m
,of
the cloud are found to be:
1.
ω
∞
=
ω
n
, the natural frequency of an individual bubble in an infinite liquid, and
2.
ω
m
=
ω
n
1+16
ηR
o
A
2
o
/π
(2
m
−
1)
2
1
2
;
m
=1
,
2
,...
, which is an infinite series
of frequencies of which
ω
1
is the lowest. The higher frequencies approach
ω
n
as
m
tends to infinity.
The lowest natural frequency,
ω
1
, can be written in terms of the mean void
fraction,
α
o
=
ηv
o
/
(1 +
ηv
o
), as
ω
1
=
ω
n
1+
4
3
π
2
A
2
o
R
2
o
α
o
1
−
α
o
−
1
2
(10.45)
Hence, the natural frequencies of the cloud will extend to frequencies much
smaller than the individual bubble frequency,
ω
n
, if the initial void frac-
tion,
α
o
, is much larger than the square of the ratio of bubble size to cloud
size (
α
o
R
2
o
/A
2
o
). If the reverse is the case (
α
o
R
2
o
/A
2
o
), all the natural
frequencies of the cloud are contained in a small range just below
ω
n
.
Typical natural modes of oscillation of the cloud are depicted in figure
262
Figure 10.10.
Natural mode shapes as a function of the normalized radial
position,
r
&
A
o
, in the cloud for various orders
m
= 1 (solid line), 2 (dash-
dotted line), 3 (dotted line), 4 ( broken line). The arbitrary vertical scale
represents the amplitude of the normalized undamped osc
illations of the
bubble radius, the pressure, and the bubble concentration per unit liquid
volume. The oscillation of the velocity is proportional to the slope of these
curves.
Figure 10.11.
The distribution of bubble radius osc
illation amplitudes,
|
φ
|
, within a cloud subjected to forced excitation at various frequencies,
ω
,
as indicated (for the case of
α
o
(1
−
α
o
)
A
2
o
/R
2
o
=0
.
822). From d’Agostino
and Brennen (1989).
263
Figure 10.12.
The amplitude of the bubble radius osc
illation at the
cloud surface,
|
φ
(
A
o
,t
)
|
, as a function of frequency (for the case of
α
o
(1
−
α
o
)
A
2
o
/R
2
o
=0
.
822). Solid line is without damping; broken line in-
cludes damping. From d’Agostino and Brennen (1989).
10.10, where normalized amplitudes of the bubble radius and pressure fluc-
tuations are shown as functions of position,
r/A
o
, within the cloud. The
amplitude of the radial velocity oscillation is proportional to the slope of
these curves. Since each bubble is supposed to react to a uniform far field
pressure, the validity of the model is limited to wave numbers,
m
, such that
m
A
o
/R
o
. Note that the first mode involves almost uniform oscillations
of the bubbles at all radial positions within the cloud. Higher modes involve
amplitudes of oscillation near the center of the cloud, that become larger and
larger relative to the amplitudes in the rest of the cloud. In effect, an outer
shell of bubbles essentially shields the exterior fluid from the osc
illations of
the bubbles in the central core, with the result that the pressure osc
illations
in the exterior fluid are of smaller amplitude for the higher modes.
10.5.2 Response of a spherical bubble cloud
The corresponding shielding effects during forced excitation are illustrated
in figure 10.11, which shows the distribution of the amplitude of bubble
radius oscillation,
|
φ
|
, within the cloud at various excitation frequencies,
ω
.
Note that, while the entire cloud responds in a fairly uniform manner for
ω<ω
n
, only a surface layer of bubbles exhibits significant response when
ω>ω
n
. In the latter case the entire core of the cloud is essentially shielded
by the outer layer.
264
Figure 10.13.
The amplitude of the bubble radius osc
illation at the cloud
surface,
|
φ
(
A
o
,t
)
|
, as a function of frequency for damped oscillations at
three values of
α
o
(1
−
α
o
)
A
2
o
/R
2
o
equal to 0.822 (solid line), 0.411 (dot-
dash line), and 1.65 (dashed line). From d’Agostino and Brennen (1989).
The variations in the response at different frequencies are shown in more
detail in figure 10.12, in which the amplitude at the cloud surface,
|
φ
(
A
o
,t
)
|
,
is presented as a function of
ω
. The solid line corresponds to the above
analysis, that did not include any bubble damping. Consequently, there are
asymptotes to infinity at each of the cloud natural frequencies; for clarity
we have omitted the numerous asymptotes that occur just below the bub-
ble natural frequency,
ω
n
. Also shown in this figure are the corresponding
results when a reasonable estimate of the damping is included in the analy-
sis (d’Agostino and Brennen 1989). The attenuation due to the damping is
much greater at the higher frequencies so that, when damping is included
(figure 10.12), the dominant feature of the response is the lowest natural fre-
quency of the cloud. The response at the bubble natural frequency becomes
much less significant.
The effect of varying the parameter,
α
o
(1
−
α
o
)
A
2
o
/R
2
o
, is shown in figure
10.13. Note that increasing the void fraction causes a reduction in both the
amplitude and frequency of the dominant response at the lowest natural
frequency of the cloud. d’Agostino and Brennen (1988) have also calculated
the acoustical absorption and scattering cross-sections of the cloud that this
analysis implies. Not surprisingly, the dominant peaks in the cross-sections
occur at the lowest cloud natural frequency.
It is important to emphasize that the analysis presented above is purely
linear and that there are likely to be very significant nonlinear effects that
265