Shock Propagation in Polydisperse Bubbly
Liquids
Keita Ando, Tim Colonius, and Christopher E. Brennen
Abstract
We investigate the shock dynamics of liquid
fl
ows containing small gas
bubbles, based on a continuum bubbly
fl
ow model. Particular
attention
is devoted to
the effects on shock dynamics of distributed bubble sizes and gas-phase nonlinearity.
Ensemble-averaged conservation laws for polydisperse bubbly
fl
ows, together with a
Rayleigh–Plesset-type model for single bubble dynamics,
form the starting point for
these studies
. Numerical simulations of one-dimensional shock propagation reveal that
phase cancellations in oscillations of differ
ent-sized bubbles can lead to an apparent
damping of the averaged shock dynamics.
Experimentally we study
the propagation of
fi
nite-amplitude waves in a bubbly liquid in a de-formable tube. The
ensemble-averaged bubbly
fl
ow model is extended to quasi-one-dimensional cases and
the corresponding steady shock relati
ons are derived. These account for
the
compressibilities associated with
the
tube deformation as well as
the
bubbles and host
liquid. A comparison between the theory and
the
water-hammer experiment suggests
that the gas-phase nonlinearity plays an essential role in the propagation of strong
shocks.
1 Introduction
A number of complex
phenomena in bubbly liquids
can be
attributed to the interaction
of the
fl
ow and
the
dynamics of dispersed bubbles. Even for the case of dilute
fl
ows
(
bubble
volume fraction less than 1%), the
fl
ow structure can be signi
fi
cantly altered
by bubble dynamics. The
presence of compressible bubbles
in a stiff liquid means that
the effective compressibility of
Keita Ando Keio University, e-mail: kando@mech.keio.ac.jp
Tim Colonius California Institute of
Technology, e-mail: colonius@caltech.edu
Christopher E. Brennen California Institute of Technology, e-mail: brennen@caltech.edu
the mixture can deviate drastically from
the
liquid compressibility. If bubbles respond
in phase with
the
fl
ow excitation, the mixture becomes less stiff and the resulting wave
speed is reduced. Since
the
bubble response depends on the excitation frequency (as
well as amplitude), the mixture compressibility and the corresponding phase velocity
are regulated by the length scale of the waves. That is, waves are
dispersive
due to
bubble dynamics. Nonlinearity associated w
ith gas-phase compressibility is another
issue, particularly evident in strong shocks.
In the modeling and simulation of the dynamics of dispersed bubbly
fl
ows, it has
been a challenge to resolve wave dispersion
and nonlinear bubble dynamics. Direct nu-
merical simulations of wave interaction with individual bubbles [27, 28, 53, 70] can
shed light on the detailed
fl
ow structures but are im
practical for large-scale
fl
ow com-
putations. For practical applications, for ex
ample in hydraulic m
achinery [5, 14],
underwater explosions [22, 46], or shockwave lithotripsy [6, 49], we favor simulating
mixture-averaged dynamics (rather than solutions to individual realizations) based on
continuum models in which the bubbly mixt
ure is treated as a continuum. Averaging
approaches allow us to remove stiffness resulting from individual bubble dynamics and
thus reduce the computational effort. In such an approach, a choice of the constitutive
relation for a model closure is required. The simplest choice is a barotropic relation
[15, 17] that ignores the bubble dynamics and is therefore incapable of capturing wave
dispersion. Alternatively, to account for the unsteadiness associated with the bubbles, a
Rayleigh–Plesset-type model for the single bubble dynamics can be coupled with
mixture-averaged conservation laws [25, 57, 86]. However, to further reduce the
computational effort, monodisperse bubbles
(of the same size) are often assumed.
We consider the canonical example of shoc
k propagation through a bubbly liquid in
a tube. In the pioneering work of Campbell and Pitcher [20] and the subsequent
experiments (see for example [11, 44, 45, 60, 79]), dispersed bubbly
fl
ows were
created in a vertical tube attached to
a shock tube in order to generate
shock
propagation. Shocks in monodisperse bubbly
fl
ows exhibit an oscillatory structure that
can also be numerically predicted using continuum models coupled with
single-bubble-dynamic equations [44, 45, 58, 70, 78, 87]. In particular, the work of
Kameda et al. [45] reveals that the shock structure in the monodisperse case is
sensitive to thermal dissipation in the bubble oscillations. However, experimental
observations on shocks in the polydi
sperse case are rather limited.
The above-mentioned experiments were con
fi
ned to the case of weak shocks.
Further increase in the shock strength results in bubble
fi
ssion [16] and tube
deformation [71]. Dynamic loading of a
fl
uid-
fi
lled deformable tube is a classic
example used to study
fl
uid-structure interactions in
water-hammer events [89]. For
tubes
fi
lled with liquids (without bubbles), a linear wave speed is predicted by the
Korteweg–Joukowsky model that accounts for
the compressibility associated with the
liquid and the structure [43, 48]. This was later extended to bubbly liquid cases [47].
On the experimental side, an underwater shock simulator has been developed to study
coupled stress waves propagating in the axial direction of a
fl
uid-
fi
lled tube [26, 31,
40]. However, the study of
fi
nite-amplitude waves in a mixture-
fi
lled, deformable tube
had not been previously reported.
The aim of this chapter is to review recent progress in the modeling and simulation
of shock dynamics of dispersed bubbly
fl
ows. In particular, we
focus on the effects of
polydispersity [3] and gas-phase nonlinearity
[4] on the mixture-averaged shock dynamics. In Sect. 2, the ensemble-averaged bubbly
fl
ow model (together with single-bubble-dynamic equations) is introduced and the
acoustic properties of polydisperse bubbly li
quids are then examined. In Sect. 3, the
numerical method for shock computations is presented and one-dimensional shock
propagation in a polydisperse bubbly liqui
d is simulated to quantify the effect of
polydispersity on averaged shock dynamics. To examine the effect of gas-phase
nonlinearity as well as
fl
uid-structure interaction (FSI), shock propagation through a
bubbly liquid in a deformable tube is considered in Sect. 4. The ensemble-averaged
model is extended to quasi-one-dimensional con
fi
gurations in order to incorporate FSI;
the corresponding steady shock relations are
compared to water-hammer experiments.
Finally, we summarize the our
fi
ndings in Sect. 5.
2 Modeling of Continuum Bubbly Flows
2.1 Mixture-Averaged Conservation Laws
In deriving a continuum model for bubbly
fl
ows, there is a need to select either time,
volume, or ensemble averaging [41, 86, 91] in order to de
fi
ne mixture-averaged
quantities. If the system is ergodic [8, 12],
averaged values do not depend on a choice
of the averaging manner under appropriate scale separation [3, 56, 65]. Here, we use
the emsemble-averaging technique of Zhang and Prosperetti [90, 91, 92] to formally
derive mixture-averaged conservation laws. To be speci
fi
c, physical quantities are
statistially averaged over a large number of realizations of spherical bubbles in the
fl
ow; any scattering effects in an individual realization are discarded but the
statistically averaged dynami
cs are explored. For dilute
fl
ows, direct interactions
between neighboring bubbles are often ig
nored. In this case, the system is
two-way-coupled
in a sense that bubble interactions
are captured indirectly through the
averaged
fi
eld. An attempt to accoun
t for direct bubble/bubbl
e interactions, which is
useful in particular for high void fraction
fl
ows (see for example [13, 19, 35, 70]), is
beyond our scope.
For the two-way-coupled case, we writ
e the emsemble-averaged conservation
equations as
∂ρ
+
∇
·
(
ρ
u
)=
0
,
(1)
∂
t
∂ρ
u
+
∇
·
(
ρ
uu
)+
∇
(
p
l
−
p
̃
)=
0
,
(2)
∂
t
∂α
R
2
R
̇
+
∇
·
(
α
u
)=
3
α
,
(3)
∂
t
R
3
where
ρ
is the mixture density,
u
is the mixture velocity,
p
l
is the averaged liquid
pressure,
α
is the volume fraction of bubbles (i.e., void fraction), and
R
and
R
̇
are the
bubble radius and the bubble wall velocity, respectively. Here, the mixture density is
approximated by
ρ
=(
1
−α
)
ρ
l
where
ρ
l
is the liquid density. It is assumed that the
bubbles follow the ambient liquid motion (i.e., no-slip), for relative motion between
the phases plays a minor role in shock propagation [44, 70]. The liquid-phase pressure
is given by the Tait equation of state [81] that assumes isentropic processes in the
liquid phase:
()
m
p
l
+ B
1
ρ
= ,
(4)
p
l
0
+
B
ρ
l
m
0
1
−α
where
ρ
l
0
is the reference liquid density at the undisturbed pressure
p
l
0
, and
m
and
B
denote the liquid stiffness and tensile st
rength, respectively. For water, we use
m
=
7
.
15 and then determine the value of
B
, specifying the speed of sound in the liquid at a
given temperature.
The phase interaction term
p
̃in the momentum equation (2), which does not appear
in classic volume-averaged equations formulated with heuristic reasoning [85, 86], is
R
3
R
̇
2
R
3
p
bw
p
̃
=
α
p
l
−−ρ
,
(5)
R
3
R
3
where
p
bw
is the bubbl wall pressure in the liquid phase [15]:
4μ
l
R
̇
2
S
p
bw
=
p
b
+
−
.
(6)
RR
Here,
p
b
is the internal bubble pressure, μ
l
is the liquid viscosity, and
S
is the surface
tension. The overbar in (5) as well as (3) denotes moments with respect to the
distribution of equilibrium bubble radius
R
0
at
p
l
0
. With no-slip and monodisperse
assumptions, the phase interaction term (5) reduces to that of Biesheuvel and van
Wijngaarden [12] or Zhang and Prosperetti [91].
To account for the effect of polydispersity
, we need to evaluate the moments in
(3) and (5):
l
∞
φ
̄
(
x
,
t
)=
φ
(
x
,
t
;
R
0
)
f
(
R
0
)
d
R
0
,
(7)
0
where
φ
represents particle quantities (such as
R
3
)
and the (normalized) bubble size distribution,
f
(
R
0
)
, is assumed spatially uniform. The
equilibrium bubble radius,
R
0
, is assumed unchanged within
fl
uid-dynamic time scales,
which are typically
1
−
2
−
10 1
10101010
*
R
0
Fig. 1
Lognormal distributions (8) with varying the standard deviation
σ
. The distributions are
normalized in order to satisfy
∞
f
(
R
∗
0
)
d
R
∗
=
1.
00
much shorter than those associated with mass diffusion of dissolved gases (such as air
in water) [30, 62]. As an illustrative exam
ple, we model the disribution using a
lognormal function (Fig. 1) with the probable size
R
ref
and standard deviation
σ
:
0
1 ln
2
R
∗
0
f
(
R
∗
0
)=
√
exp
−
,
(8)
2
πσ
R
∗
0
2
σ
2
= R
0
/
R
ref
where
R
∗
0
. In the limit of
σ
→
0, the distribution (8) reduces to the Dirac
0
delta function that models monodisperse mixtures. The measured distribution of gas
bubble nuclei in water tunnel
s or seawater is reasonably
fi
t by the lognormal
distribution with
σ
≈
0
.
7 [3, 24]. For example, the void fraction is computed as
4
π
α
=
nR
3
,
(9)
3
where
n
is the number of bubbles per unit volume of the mixture. Provided that
fi
ssion and coalescence of the bubbles do not occur, the bubble number density is
conserved in time; namely
∂
n
+
∇
·
(
n
u
)=
0
.
(10)
∂
t
However, for strong shocks, there will arise bubble
fi
ssion possibly due to a Rayleigh–
Taylor-type shape instability or a re-entrant jet [16, 29, 42, 68], so that the bubble
number conservation no longer holds.
2.2 Single-Bubble-Dynamic Equations
A closure of the mixture-averaged equations (1) to (3) requires evaluating individual
bubble dynamics. Taking into account acoustic damping (to
fi
rst order) on bubble
oscillations, we use the Gilmore equation [36] that is an improvement to the incom-
pressible Rayleigh–Plesset equation [61, 69]:
() ()()()
R
̇
3
R
̇
R
̇
RH
̇
R
̇
RR
̈1
−
+
R
̇
2
1
−
=
H
1
++
1
−
,
(11)
C
23
C CCC
where the dot denotes the substantial time derivative, and
H
and
C
are the enthalpy and
the speed of sound in the liquid, respectively, at the bubble wall:
l
′
p
bw
d
p
d
p
l
H
=
l
,
C
=.
(12)
ρ
l
(
p
′
d
ρ
l
p
l
l
)
p
bw
In calculating
H
, the term
p
l
is taken from the averaged liquid pressure in the con-
tinuum model. Note that the averaged liquid pressure can be regarded, in the dilute
limit, as far-
fi
eld pressure from the bubble [18, 77].
Finally, the internal bubble pressure
p
b
in (6) needs to be determined. The simplest
choice is the polytropic gas law:
()
−
3
κ
R
p
b
=
p
v
(
T
w
)+
p
g
0
,
(13)
R
0
where
p
v
is the vapor pressure at bubble wall temperature
T
w
,
p
g
0
is the partial pressure
of the noncondensible gas in the equilibrium state, and
κ
is the polydropic index that
ranges from 1 (isothermal) to the ratio of speci
fi
c heats
γ
g
(adiabatic). Here, it is
assumed that the bubble contents have spatially uniform pressure. This homobaric
assumption is valid since the inertia of the bubble contents is neglibible compared to
the liquid inertia [66]. However, for nonlinear computations, using
(13) often gives rise to innacurate evaluations of the thermal behabior of the bubble.
For resolving the thermal effect on the dy
namics of bubbles, the full conservation
equations for both liquid and gas phases need to be solved. If the bubble contents are
perfect, the bubble energy equation can be simpli
fi
ed to [59, 66],
( )
3
γ
b
γ
b
−
1
∂
T
′′
p
̇
b
−
̇
m k
bw
(14)
=
Rp
b
+ R
v
T
w
̇
+ ,
v
R
γ
b
∂
r
w
where
γ
b
is the speci
fi
c-heat ratio of the bubble contents,
R
v
is the gas constant for
′′
the vapor,
m
̇
is the mass
fl
ux of the vapor at the bubble wall, and
k
bw
is the ther
v
mal conductivity of the bubble contents. The last term in the parenthesis stands for
heat conduction, across the bubble wall,
driven by the temperature gradient (
∂
T
/
∂
r
)
in the gas phase. Provided that the liquid is cold (well below the boiling point) and
phase change occurs instantaneously, the vapor pressure can be assumed unchanged,
except near the the end of a violent bubble collapse [34], with constant wall tem-
perature. In this case, the liquid-phase thermodynamics becomes irrelevant. In the
simulations of shocks in monodisperse bubbly
fl
ows [44, 45, 87], the full conservation
equations for the bubble contents are solved
to accurately evaluate the heat and mass
fl
uxes in (14). However, computations of the full partial differential equations are still
expensive for polydisperse cases. To further reduce the conputational effort, one may
use the reduced-order model of Preston, Colonius, and Brennen [63] which constitutes
a set of coupled ordinary differential equations describing diffusive effects on the
single bubble dynamics.
We now write the bubble-dynamic equations in conservation form, which is suit-
able for shock computations [50]. For example, the evolution of the bubble radius can
be described by
∂
nR
+
∇
·
(
nR
u
)=
nR
̇
,
(15)
∂
t
where the bubble radius is de
fi
ned as Eulerian
variable
R
(
x
,
t
;
R
0
)
, enen though it is attributed to Lagrangian particles. Now that a
large number of bubbles are assumed to exist for the mixture to be a continuum, the
bubbles may be considered to be distribute
d continuously in space [87]. The evolution
of other bubble-dynamic variables (such as
R
̇
and
p
b
) can be written in the same
fashion.
2.3 Acoustics of Polydisperse Bubbly Liquids
Before proceeding to shock computations in
Sect. 3, we examine the effect of poly-
dispersity on linear waves in a dilute bubbly liquid. By linearizing the continuum
bubbly
fl
ow equations (Sect. 2.1) coupled with the single-bubble-dynamic model
(Sect. 2.2), one may derive the dispersion relation [25], which is indentical to that from
multiple scattering theory [21, 33]:
l
11
∞
R
0
f
(
R
0
)
d
R
0
=+
4
π
n
,
(16)
c
2
2
ω
2
m
c
l
0
N
−ω
2
+
i2
βω
where
c
m
is the complex speed of sound in the mixture,
c
l
is the sonic speed in the
liquid alone,
β
is the bubble-dynamic damping constant,
ω
is the angular frequency (
ω
=
2
π
f
), and
ω
N
is the natural frequency of bubble oscillations. Since the damping
constant and the natural frequency depend on the equilibrium bubble size, their
contributions need to be integrated over the distribution,
f
(
R
0
)
, for polydisperse
mixtures. For the case of gas bubbles (without vapor), the damping constant and the
natural frequency can be written as [64]
2μ
l
p
g
0
ω
2
R
0
/
2
c
l
β
=+
ℑ
{
Υ
}
+
(17)
ρ
l
R
2
2
ρ
l
ω
R
2
1
+(
ω
R
0
/
c
l
)
2
,
00
()
p
g
0
2
S
(
ω
R
0
/
c
l
)
2
ω
2
=
ℜ
{
Υ
}−
+
ω
2
,
(18)
N
ρ
l
R
2
p
g
0
R
0
1
+(
ω
R
0
/
c
l
)
2
0
where the
fi
rst, second and third terms on the right-hand side of the damping constant
stand for viscous, thermal
and acoustic contributions,
respectively. The quantity
Υ
is a
function of the Peclet number (Pe
=
ω
R
2
0
/
α
T
;
α
T
is the thermal diffusivity of the gas):
3
γ
g
Υ
= (
√
).
(19)
√
1
−
i3
(
γ
g
−
1
)
Pe
−
1
iPecoth iPe
−
1
The effective polytropic index in (13) for th
ermal behavior of the gas is therefore given
by
κ
=
ℜ
{
Υ
}
/
3. In the limit of Pe
→
0, the thermal behavior is isothermal (
κ
→
1).
Large values of Pe, on the contrary, imply that the thermal boundary layer is thin
compared to the bubble radius and therefor
e the bubble tends to behave adiabatically
(
κ
→γ
g
). For free oscillations of such large-
sized bubbles, the frequency is well
approximated by the Minnaert frequency [54]:
13
γ
g
p
l
0
ω
M
= ,
(20)
R
0
ρ
l
where the effects of surface tension, heat tr
ansfer and acoustic radiation in (18) are
ignored.
The phase velocity
c
ph
and the attenuation
a
att
are de
fi
ned as [25]
−
1
1
ph
c
=
ℜ
,
(21)
c
m
att
1
a
=
−
20
(
log
10
e
)
ω
ℑ
,
(22)
c
m
where
a
att
has the unit of decibels per length. The acoustic theory is known to overes
-
timate attenuation under resonance [88], but the error associated with resonance can be
deemphasized by including the broad bubble size distribution since the probability of
fi
nding bubbles of certain size under resonance is low among a wide spectrum of
R
0
[32].
In Fig. 2, we examine the dispersion relation for linear waves in an air/water of void
fraction
α
=
0
.
001 at standard temperature and pressure (STP; 20
◦
C, 101kPa). The
bubble size in the mixture is assumed lognormally distributed about
R
ref
=
0 10μm
with varying the standard deviation
σ
. Figure 2(a) shows the phase velocity in
frequency space. While the phase velocity reduces to the speed of sound in the liquid
alone far above the resonant frequency (at which bubbles cannot respond to such high
frequency excitation), the sonic speed of the mixture in a low frequency limit is given
by [15]
m
(
p
l
+
B
)
p
l
ph
c
=
lim
c
= .
(23)
f
→
0
ρα
m
(
p
l
+
B
)+(
1
−α
)
p
l
Thus, in the quasistatic regime (where all the phases are in dynamic equilibrium at all
times), the sonic speed can be even smaller
than the speed of sound in air, as in
−
2
−
10 1
−
2
−
10 1
1010101010101010
f
(MHz)
f
(MHz)
Fig. 2
Dispersion relation for linear waves in an air/water mixture of
α
=
0
.
001 at STP. The bubble
size is assumed lognormally distributed about
R
ref
=
10
μm with varying
σ
. The isothermal natural
0
frequency for
R
ref
is0
.
29MHz. (
a
) Phase velocity. (
b
)
Attenuation.
0
this particular example. That is, the inclusion of a tiny amount of bubbles yields sig-
ni
fi
cant reductions in the
fl
uid compressibility, while the
fl
uid is yet dense. However,
around the resonant condition, the bubble size distribution tends to have a striking
impact on the phase velocity. If the mixture is monodisperse (i.e.,
σ
=
0) or the
distribution is narrowly peaked (say,
σ
=
0
.
1), the phase velocity is reduced as the
external frequency increases to the resonant frequency, for the amplitude of bubble
oscillations becomes larger and the mixture compressilibity is thus enhanced. If further
increasing the frequency above the resonant frequency, the phase of the oscillations
changes. Speci
fi
cally, an increase in the
ambient pressure leads to bubble expansion.
Conceptually, the mixture becomes stiffer
than the host liquid so that the phase
velocity above the resonant frequency goes beyond the speed of sound in the liquid.
On the other hand, for the case of broad size distributions, the abrupt transition across
the resonant frequency is eff
ectively eliminated due to the
fact that having bubbles of
certain size under resonance is less likely among a broud spectrum of the distribution.
Furthermore, it should be noted that the bubble size distribution increases the
attenuation below the resonant frequency (see Fig. 2(b)). To interpret this trend, we
consider linear bubble oscillations under sinusoidal excitation (without mutual
interactions between the bubbles) [64],
ε
p
l
0
χ
̈
+
2
βχ
̇
+
ω
N
2
χ
=
−
sin
(
ω
t
),
(24)
ρ
l
R
2
0
where
ε
is the in
fi
nitesimal amplitude of sinusoidally oscillating far-
fi
eld pressure (
|ε
|≪
1) and
χ
is the corresponding perturbed bubble radius:
[]
p
l
=
p
l
0
1
+
ε
sin
(
ω
t
),
(25)
R
=
R
0
(
1
+
χ
) .
(26)
0123456
t*
Fig. 3
Evolution of the
fi
rst moment (35) as a function of time normalized by the Minnaert frequency
(
ω
M
/
2
π
) for
R
ref
0
. The (inviscid) bubbles are set into linear oscillations according to a
stepwise change in the far
fi
eld pressure. The bubble size is assumed lognormally
distributed with the standard deviation
σ
.
The particular solution to the linearized problem (24) is
−ε
p
l
0
χ
p
=
�
sin
(
ω
t
+
φ
),
(27)
2
ρ
l
R
2
ω
N
2
−ω
2
+
4
β
2
ω
2
0
where the phase shift
φ
is given by
ω
N
2
−ω
2
φ
=
cos
−
1
�
.
(28)
2
ω
N
2
−ω
2
+
4
β
2
ω
2
Noting that the damping constant and the natural frequency depend on the bubble size
in the quasistatic regime, the phase shift can also vary with different values of
R
0
so
that there may arise cancellations due to the different phases among the different-sized
bubbles. Hence, (collective) volume
oscillation of a polydisperse cloud is
deemphasized due to the phase cancellations
. In other words, the phase cancellation
effect can be regarded as an
apparent
damping of the mixture-averaged dynamics and
is augmented as the distribution broadens [2]. This additional damping associated with
polydispersity may account for the observation at low frequency in Fig. 2(b).
The phase cancellation effect appears also
in free oscillations in one-way-coupled
fl
ow, which may be modeled by replacing the sinusoidal excitation term on the
right-hand side of (24) with stepwise forcing; namely
ε
p
l
0
χ
̈
+
2
βχ
̇
+
ω
N
2
χ
=
−
H
(
t
),
(29)
ρ
l
R
2
0
where
H
(
t
)
is the Heaviside step function. If the bubbles are initially stationary (i.e.,
χ
=
0 and
χ
̇
=
0 at
t
=
0), the solution to (29) is then given by
[]
ε
p
l
0
Ω
−
e
−β
t
Ω
cos
(
Ω
t
)
−β
sin
(
Ω
t
)
χ
=
−
,
(30)
ρ
l
R
2
0
Ωβ
2
+
Ω
2
where
Ω
=(
ω
N
2
−β
2
)
1
/
2
. For the inviscid case (
β
=
0), the bubbles keep oscillating
without any damping and the solution (30) is simpli
fi
ed to
[]
ε
p
l
0
cos
(
ω
N
t
)
−
1
χ
= .
(31)
ρ
l
R
2
0
ω
2
N
Moreover, if the natural frequency is
approximated by the Minnaert frequency
ω
M
in
(20), the inviscid solution further reduces to
which can be rede
fi
ned as a normalized perturbation from the new equilibrium:
()
3
γ
g
ε
t
3
γ
g
p
l
0
′
χ
=
χ
+=
cos
.
(33)
ε
3
γ
g
R
0
ρ
l
To observe the dynamics of clouds in which the bubbles are freely oscillating ac-
cording to (33), we de
fi
ne the moments of the (perturbed) bubble radius with respect to
a smooth bubble size distribution
f
(
R
0
)
:
l
∞
[]
j
′
μ
j
(
t
)=
χ
′
(
t
;
R
0
)
f
(
R
0
)
d
R
0
.
(34)
0
′
For example, μ
1
stands for the mean bubble radius
. It follows from the Riemann–
Lebesgue lemma [9] that the
fi
rst moment vanishes as
t
→ ∞
, provided that the
∞
χ
=
ε
3
γ
g
cos
t R
0
3
γ
g
p
l
0
ρ
l
−
1
,
(32)
sources) are
integral
|
f
(
R
0
)
|
d
R
0
is bounded:
0
l
∞
t
3
γ
g
p
l
0
′
μ
1
(
∞
)=
lim cos
f
(
R
0
)
d
R
0
=
0
.
(35)
t
→∞
0
R
0
ρ
l
That is, bubble oscillations eventually reach a
statistical equilibrium
at which the
different-sized bubbles oscillate totally out
of phase and the polydisperse bubble cloud
can thus be considered to be stationary, regardless of inviscid bubble oscillations. The
existence of the statistical equilibrium is numerically veri
fi
ed in Fig. 3
′
where the moment μ
1
is evaluated using the lognormal function (8). While (inviscid)
oscillations continue in the monodisperse mixture (
σ
= 0), the inclusion of bubble
size distributions eventually yields a ti
me-invariant value. Furthermore, the
broader the distribution, the more
quickly cancellation between bubbles at
different
phases of their oscillation cycles is achieved. If the distribution is suf
fi
ciently
broad, the moment can converge to a stationary state, due to apparent damping
associated with the strong cancellation effect, only within a couple of oscillation
periods for the probable size.
Nonlinear oscillations of inviscid bubbles in one-way-coupled
fl
ow can also be
shown to reach a statistical equilibrium a
nd the apparent damping associated with
polydispersity may dominate over physical dampings, in the averaged sense, if the
distribution is suf
fi
ciently broad (for details see [3, 24, 75]). In the next section, we
simulate one-dimensional shock propagation based on the continuum model and show
that the phase cancellation effect comes into
play in two-way-coupled poly-disperse
fl
ows as well.
3 Simulation of Averaged Shock Dynamics
3.1 Numerical Method
Because shocks in bubbly
fl
ows often contain oscillatory structures due to bubble
dynamics, we favor the properties of high-order accurate resolution in complex smooth
structures as well as shock capturing and robustness. Here, we select a
fi
nite-volume
(FV) weighted essentially non-oscillatory
(weighted ENO or WENO) scheme [52, 72],
which has proven to be stable and accurate for shock computations in various
examples, together with monotonicity preserving bounds [7]. ENO reconstruction [38]
is based on adaptive stencils in the sense that interporation across discontinuities is
automatically avoided. WENO schemes consist of a convex conbination of all the
ENO candidate stencils for more ef
fi
cient and accurate eval
uations. To guarantee
entropy solutions, it is safe to implement
the WENO procedure in characteristic space
[67].
The system of the governing equations in one dimension can be written as
∂
q
∂
f
(
q
)
∂
f
s
(
q
)
+= +
s
(
q
).
(36)
∂
t
∂
x
∂
x
The column vectors (conserved variables, numerical
fl
uxes, and bubble-dynamic
n
φ
̇
where
u
is the
x
-component velocity,
φ
represents the bubble-dynamic variables such
as
R
and
R
̇
, and the superscript
s
denotes sources that vanish in the equilibrium state.
Given a computational cell
[
x
i
−
1
/
2
,
x
i
+
1
/
2
]
where
i
denotes the grid index, the system
(36) is discretized in FV fashion:
χ
=
ε
3
γ
g
cos
t R
0
3
γ
g
p
l
0
ρ
l
−
1
,
(32)
sources) are
ρ
ρ
u
0
0
q
= ,
α
ρ
u n
φ
f
= ,
α
u
ρ
u
2
+
p
l
n
φ
u
f
s
= ,
0
p
̃0
s
= ,
R
3
0 3
α
R
2
R
̇
(37
)
f
s
i
+
1
/
2
−
f
s
d
q
̄
i
f
i
+
1
/
2
−
f
i
−
1
/
2
i
−
1
/
2
=
−
++
s
̄
i
,
(38)
d
t
Δ
x
i
Δ
x
i
where
Δ
x
i
is the cell width (
x
i
+
1
/
2
−
x
i
−
1
/
2
) and the overbar denotes the cell-averaged
quantities. In the FV method,
q
̄
i
is reconstructed at each side of the cell edge and the
numerical
fl
ux
f
i
+
1
/
2
is determined from a local Riemann problem. Here, we implement
the
fi
fth-order monotonicity-preserving FV-WENO scheme in the characteristic space
for reconstructing the cell-edge conserved variables and then calculate the numerical
fl
ux using the HLLC approximate Riemann so
lver [84, 83] that automatically satis
fi
es
the entropy condition. The HLLC Riemann solver restores contact waves that the HLL
Riemann solver [39] ignores, and thus gives better resolutions of the contact
discontinuities.
Once the HLLC
fl
uxes and the sources are dete
rmined, the system (38) in
semi-discrete form can be integr
ated in time. In the examples
in Sect. 3.3, a third-order
TVD Runge–Kutta scheme [37, 73] is used to integrate both the convective terms and
the bubble-dynamic sources. Such an un
split integration method works for weak
shocks, for the bubble-dynamic sources are not very stiff. It is instructive to note,
however, that if the system is stiffer in
particular for cases with strong shocks or
cavitation bubbles, one may need to apply time-step splitting techniques [35, 51, 83] or
implicit methods [23, 80] for stable time integration.
In Sect. 3.3, the FV-WENO scheme togeth
er with the HLLC Riemann solver is
used to compute the numerical
fl
uxes and the system is inet
grated in time using the
unsplit method. The computational grid is uniform with
Δ
x
=
R
ref
. For the polydis
0
perse case, the moments (7) associated with the distribution of equilibrium bubble
sizes are evaluated using Simpson’s rule with 101 quadrature points. This method has
been shown to be accurate en
ough to resolve wave dispersion in a wide range of
frequency. Moreover, a comparison to the experiment [45] validates the continuum
approach to predict mixture-averaged
dynamics for weak shocks. For further
infomation, see [3].
3.2 Steady Shock Relations
The simulation of shock propagation requires
the steady shock relations to be used as
initial conditions. While shocks in single phase
fl
ow have in
fi
nitesimal thickness (in a
continuum sense), the thickness of shocks in bubbly
fl
ow becomes
fi
nite due to
bubble-dynamic relaxation. In front of the shock, the bubbles are in equilibrium at (
R
0
,
T
0
,
p
l
0
) where
T
0
is the (undisturbed) liquid temperature. Far downstream of the shock
front, the bubble dynamics are
fi
nally damped out and the bubbles are once again in
equilibrium at (
R
H
,
T
0
,
p
lH
) where
R
H
is the new equilibrium radius corresponding to
the shock pressure
p
lH
. The speci
fi
cation of
T
0
in the
fi
nal equilibrium state implies
that the bubble temperature
eventually returns to the li
quid temperature due to heat
conduction across the bubble wall. Now, the one-dimensional conservation laws for
mass, momentum and bubble number are writte
n in a coordinate system moving with
the constant shock velocity
U
s
:
′
d
ρ
u
= 0, (39)
d
x
′
d
′
2
+ p
l
−
̃
(
ρ
up
)=
0
,
(40)
d
x
′ ′
d
nu
= 0, (41)
d
x
′
′′
where
u
is the velocity in the coordinate system with
x
=
x
−
U
s
t
. Integrating (39) to
(41) from upstream to far downstream, it transpires that, independent of the detailed
shock structure,
′
−ρ
H
u
H
=
ρ
0
U
s
,
(42)
′
2
ρ
H
u
H
+
p
lH
=
ρ
0
U
s
2
+
p
l
0
,
(43)
′
−
n
H
u
=
(44)
H
n
0
U
s
,
where the subscripts 0 and
H
denote upstream and downstream quantities, respec-
tively, and
ρ
0
=(
1
−α
0
)
ρ
l
0
and
ρ
lH
=(
1
−α
H
)
ρ
lH
.
The new equilibrium bubble radius
R
H
is related to the shock pressure
p
lH
:
()( )
−
3
2
Υ
R
H
2
Υ
p
lH
=
p
l
0
−
p
v
++
p
v
−
,
(45)
R
0
R
0
R
H
where the bubbles are assumed to behave isothermally so that the bubble temperature
fi
nally returns to
T
0
. From (42) and (44), we write the bubble number density far
behind the shock as
()
1
/
m
−
1
p
l
0
+ B
4
π
n
H
=(
1
−α
0
)+
n
0
R
3
.
(46)
p
lH
+
B
3
H
From (42) and (43), the steady shock speed becomes
p
lH
−
p
l
0
U
s
= (),
(47)
ρ
0
ρ
0
1
−
ρ
H
and the induced velocity far downstream of the shock front is then given by
()
ρ
0
′
u
H
=
u
H
+
U
s
=
1
−
U
s
.
(48)
ρ
H
It can be shown that the shock speed (47) approaches the sonic speed (23) if the shock
strength is in
fi
nitesimal.
(a)
/
p
x 10
−
2
α
0.3
σ
= 0.7
σ
= 0
0.2
02468
(c)
1.1
(b)
1
1.5 1
0.5
0.4
R
/
R
0.9
0.8
0.7
x
(mm)
Fig. 4
Spatial evolution of the unsteady shock propagation, 5
.
2μs after the steady shock relations are
imposed at
x
=
0, through an air/water mixture of
α
0
=
0
.
005. The equilibrium bubble radius is
lognormally distributed about
R
ref
=
10
μm with the standard deviation
σ
.(
a
) Averaged liquid
0
pressure. (
b
) Void fraction. (
c
) Bubble radius with different equilibrium sizes.
(a)
t
= 15 μs
t
= 36 μs
0 2 4 6 810
0 2 4 6 810
x
(mm)
Fig. 5
Spatial evolution of the shock propagation, 15μs (solid lines) and 36μs (dashed lines) after
the steady shock relations are imposed at
x
=
0, through the polydisperse mixture of
α
0
=
0
.
005,
R
ref
0
=
10
μm, and
σ
=
0
.
7. (
a
) Ensemble-averaged liquid pressure. (
b
) Bubble radius with different
equilibrium sizes.
3.3 One-Dimensional Shock Propagation
Here, we consider unsteady and steady shock propagation through a polydisperse
bubbly liquid. As a representative example, we simulate shock propagation in an
air/water mixture of void fractions below
α
0
=
0
.
005 at STP. The equilibrium bubble
size is assumed lognormally distributed with varying the standard deviation
σ
. The
shock pressure is set to
p
lH
=
2
p
l
0
. The steady shock relations in Sect. 3.2 are initially
imposed by a diaphragm at
x
=
0. Diffusive effects on the dynamics of single bubbles
(i.e., heat and mass diffusion at the bubble wall; see Sect. 2.2) are evaluated using the
reduced-order model of Preston, Colonius, and Brennen [63], which is accurate for the
case of micron-sized air bubbles in water.
We judge steadiness by observing the
fi
rst
peak of pressure oscillations due to bubble dynamics;
2
1.5
1
x
(mm)
Fig. 6
Spatial evolution of the averaged liquid pressure for steady shock propagation through an
air/water mixture of
α
0
=
0
.
005. The equilibrium bubble radius
is lognormally distributed about
R
ref
=
10μm with
σ
ranging from 0 to 0
.
7. The position where the pressure is
(
p
l
0
+
p
lH
)/
2 is set
0
at
x
=
0.
if the peak pressure is unchanged, the shock
propagation is consider
ed to be in a steady
state. Note that we limit our attention to the case of weak shocks in a dilute liquid in
order to avoid issues associated with bubble
fi
ssion [16, 29] and direct bubble
interactions [35, 70].
In Fig. 4, we examine unsteady shock propagation in an air/water mixture of
α
0
=
0
.
005 and
R
ref
=
10μm at STP; the mixture is monodisperse (
σ
=
0) or polydisperse
0
(
σ
=
0
.
7). In addition to the averaged liquid pressure and void
fraction
fi
elds, the spatial evolution of the bubble radius for different equilibrium
sizes (
R
∗
0
=
0
.
25, 0
.
5, 1, 2, 4) is plotted to visualize the individual bubble
dynamics. It follows from the pressure
fi
eld (evidently for the polydisperse case)
that high-frequency waves precede the
primary shock wave and propagate
essentially with the speed of sound in (p
ure) water. However, these precursory
waves do not perturb the void fraction
fi
eld (Fig. 4(b)), for most bubbles
(excluding very small bubbles) do not respond to such high-frequency forcing
(Fig. 4(c)). While the precursory pressure waves in the monodisperse mixture is
damped out, those in the polydisperse
mixture are still on the decay. This is
because the bubble size distribution decreas
es the attenuation
of high-frequency
acoustic waves as demonstrated in Fig. 2(b). It is also interesting to note that
oscillatory shock structure is obtained in
the monodisperse mixture but not in the
polydisperse mixture in which the different
-sized bubbles oscillate with different
phases in the neighborhood of the shock front. To see the details of the
polydisperse case, Fig. 5 presents the spatial evolution of the shock propagation at
two different times at which the larger-s
ized bubbles still show oscillations with
less effective damping. Despite the (undamped) bubble oscillations, the shock
pro
fi
le in the averaged pressure
fi
eld seems unchanged during this period; the
shock propagation can thus be said to reach a steady state. Because the
different-sized bubbles oscillate with di
fferent phases as in the case of linear
waves
p
/
p
ll
0
p
ll
0
−
4
0
0
−
2
0
0
0
2
0
0
4
0
0
(b)
2
/
p
1.5
r
e
f
=
5
μ
m
0
R
r
e
f
=
1
0
μ
m
R
0
1
R
ref
= 20
μm
0
−
400
−
200 0 200
400
r
e
f
x
/
R
0
Fig. 7
Effect of the probable bubble size
R
ref
on steady shock structure in an air/water mixture
of
0
α
0
=
0
.
005. The position where the pressure is
(
p
l
0
+
p
lH
)/
2 is set at
x
=
0. (
a
) Monodisperse case (
σ
=
0). (
b
) Polydisperse case (
σ
=
0
.
7).