7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Shock Theory of a Bubbly Liquid in a Deformable Tube
K. Ando*, T. Sanadat, K. Inaba ,
J. E. Shepherd*, T. Colonius* and C. E. Brennen*
Department of Mechanical Engineering, California Institute of Technology, Pasadena, CA 91125, USA
t Department of Mechanical Engineering, Shizuoka University, Hamamatsu 4328561, Japan
t Department of Mechanical Engineering and Science, Tokyo Institute of Technology, Tokyo 1528550, Japan
kando @tcaltech.edu, tIlikId 6 ipc shi/uolli c ip inaibj mcch li cch ac ip
joseph.e.shepherd@caltech.edu, colonius@caltech.edu and brennen@caltech.edu
Keywords: Fluidstructure interaction, bubbly water hammer
Abstract
Shock propagation through a bubbly liquid filled in a deformable cylindrical tube is considered. A steady shock
theory that accounts for the compressibility of the structure, the liquid and the bubbles is reviewed. Experiments are
conducted in which a freefalling steel projectile impacts the top of an air/water mixture in a polycarbonate tube, and
stress waves in the tube material are measured. The experimental data indicate that the linear theory cannot properly
predict the propagation speeds of shock waves in mixturefilled tubes; the shock theory is found to more accurately
estimate the measured wave speeds.
Introduction
A fundamental understanding of fluidstructure interac
tion (FSI) is of great importance in industrial piping sys
tems, underwater explosions and turbomachinary (Wylie
and Streeter 1993; Cole 1948; Brennen 1994). These
flows often involve gas (or vapor) bubbles that alter
the dynamics of the fluid dramatically (Brennen 1995,
2005). Dynamic loading of fluidfilled, deformable
tubes has been extensively studied as an FSI model prob
lem (Tijsseling 1996; Ghidaoui et al. 2005). Liquid
filled tubes were first studied by Korteweg (1878) and
Joukowsky (1898), who introduced a linear wave speed
that accounts for the compressibility of both the liquid
and the structure. The KortewegJoukowsky wave speed
in the case of bubbly liquids was later validated by Ko
bori et al. (1955). A (nonlinear) shock theory that in
cludes both structural compressibility and bubbles has
been developed by Ando et al. (2010).
The goal of this paper is to examine the steady shock
theory for a bubbly liquid in a deformable cylindrical
tube. We review the derivation of the steady shock rela
tions, and report on experiments in which a freefalling
steel projectile impacts the top of a polycarbonate tube
filled with air/water mixtures with void fractions up to
one percent. Stress waves in the tube material are mea
sured and used to infer wave speeds. The measured wave
speeds are compared to the shock theory, and the model
limitations are discussed.
Nomenclature
Roman symbols
a Internal tube radius (m)
A Internal tube area (m2)
B Tensile strength of a liquid (Pa)
c Sonic speed (ms 1)
cj KortewegJoukowsky wave speed (m s 1)
E Young's modulus (Pa)
h Tube wall thickness (m)
K Bulk modulus (Pa)
Ms Shock Mach number()
n Bubble number density (m 3)
p Pressure (Pa)
R Bubble radius (m)
t Time (s)
u Velocity (ms )
Us Shock speed (ms )
x Spatial coordinate (m)
Greek symbols
a Void fraction ()
7 Stiffness of a liquid ()
co Hoop strain ()
Fluidstructure interaction parameter ()
p Density (kg m 3)
Subscripts
g Gas property
H Quantity (far) behind a shock
1 Liquid property
s Solid property
0 Initial (undisturbed) value
FSI bubbly flow model
The quasionedimensional bubbly flow equations that
account for structural deformation (Ando et al. 2010) are
reviewed. The liquid and disperse phases are treated as
a continuum in order to evaluate the average mixture dy
namics. Based on the ensembleaveraging technique of
Zhang and Prosperetti (1994), the continuum model can
formally be derived. The continuum assumption may be
effectively valid in the dilute limit. The derivation of the
ensembleaveragedbubbly flow equations and the model
limitations are further discussed in Ando (2010).
In what follows, we include the effect of FSI in the
bubbly flow equations. Let A be the internal cross
sectional area of the cylindrical tube. Assuming that the
tube inertia is negligible and the liquid pressure is only
balanced by the hoop stress, the tube area may be given
quasistatically by (Tijsseling 1996)
A Ao 1+2a (pi o) i (1)
I h Eh
where a is the internal tube radius, h is the wall thick
ness, E is Young's modulus of the tube material, pi is the
averaged liquid pressure, and the subscript 0 denotes the
initial (undisturbed) values. With a conventional control
volume analysis, the quasionedimensional mixture
averaged equations become
dpA dpuA
at dx
OpuA + a dPUA+pA
at 09X
Aoao 2
Eh1)
0, (2)
ipA
Ox(
(3)
DnA DnuA
+ 0 0, (4)
at o9x
where p is the mixture density, u is the mixture velocity,
and n is the number of bubbles per unit volume of the
mixture. For dilute mixtures, the mixture density is well
approximatedby (1 a)pl where pi is the liquid density
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
and the a is the void fraction defined as
a n j R3f(Ro)dRo. (5)
3 o
Here, R is the bubble radius, Ro is the equilibrium bub
ble radius corresponding to the ambient pressure (plo
101 kPa), and f(Ro) represents the normalized bubble
size distribution. Assuming that the liquidphase flow
is homentropic, the averaged liquid pressure will be de
scribed by the Tait equation of state (Thompson 1972),
Pl+B 1
pio + B p0 1
Pa '
where plo is the reference liquid density at plo, and 7 and
B denote the stiffness and tensile strength of the liquid,
respectively. For water, we take 7 7.15 and B
304 MPa. The term 1 in equation (3) represents pressure
fluctuations due to the phase interactions and vanishes in
the equilibrium state.
Steady shock speeds
Sonic speeds. We first derive the sonic speeds of the
bubbly liquid with and without FSI. For convenience,
we define the bulk modulus of the mixture, K, as
1 1 a a
+  (7)
K K1 K' (7)
where K1 and K. are the bulk moduli of the liquid and
gas, respectively. For Tait liquids, we have Ki = (pl +
B). If the effects of the vapor and the surface tension
are neglected, we take K pi. With the mixture bulk
modulus (7), the sonic speed of the mixture (in which
the bubbles behave quasistatically) becomes
K Ki/p (8)
P \1+af 1i
We now include the effect of the structural compress
ibility on the mixture sonic speed (8). The Korteweg
Joukowsky wave speed for the mixture is evaluated as
Q c K (9)
where and j determine the extent of fluidstructure
coupling for the cases of the mixture and the liquid
alone, respectively:
2Kao
Eh
2Kjao
Eh
Note that the structural compressibility reduces the lin
ear wave speed in the mixture (i.e. cj < c).
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Steady shock relations. The Hugoniot conditions for
equations (2) to (4) can readily be derived with the
assumption of equilibrium bubbles in front of and far
downstream of the shock. The detailed derivation of the
steady shock relations is found in Ando et al. (2010).
The steady shock speed is written as
Ig(plH) g(pio)
poAo 1 oAo
pHAH
where the subscript H denotes quantities behind the
shock, and
10 6 10 5 104 10 3 10 2
l ( A P loao Aoao 2
g(p) Ao 1 Eh Eh
The induced velocity far downstream of the shock front
is given by
ttH 1
poAo
PHAH
It is readily shown that the shock speed (10) approaches
the Joukowsky wave speed (9) if the shock strength is
infinitesimal. Consequently, the shock Mach number is
defined as
Gasphase nonlinearity. We document the steady
shock relations for the case of bubbly water with Ql =
6.14 where the value of El is computed based on the
properties of the polycarbonate tube that is used in the
experiments. For simplicity, we ignore the effects of va
por pressure and surface tension. Figure 1 demonstrates
the effects of the initial void fraction and the shock pres
sure on the shock speed and Mach number. Note that
PlH = pio indicates the linear wave cases, in which the
shock speeds reduce to the sonic speeds. It follows from
figure 1 (top) that the reduction in the shock speed due
to structural compressibility is minimized for finite val
ues of ao since the gasphase compressibility dominates
over the compressibility of the water and structure. It
is also seen that the shock speeds are greatly reduced by
even a tiny void fraction. Moreover, unless the void frac
tion is extremely small, the finite shock strength yields a
significant deviation from the linear wave speed due to
the nonlinearity associated with the gasphase compress
ibility. As a result, the shock Mach number increases as
the void fraction increases as seen in figure 1 (bottom).
We note that the shock Mach numbers are only slightly
greater than 1 for the case of water alone (ao = 0) since
the pressure perturbations up to several hundred atmo
spheres remain very weak (Thompson 1972).
106 10 5 10 4 103 102
ao0
Figure 1: Steady shock speeds (top) and shock Mach
numbers (bottom) in bubbly water. The curves are pa
rameterized by PlH/plo = 1, 2.5, 5, 10.
Waterhammer experiments
Experimental setup. Experiments were conducted in
order to examine the steady shock theory. The experi
mental apparatus depicted in figure 2 is similar to that
of Inaba and Shepherd (2010), and consists of a verti
cal polycarbonate tube (PCT0021.25, San Diego Plas
tics; E 2.13 GPa, ps 1200 kgm 3, ao 3h
19.1 mm) filled with an air/water mixture. A barrel is
mounted above the tube and a 1.50kg cylindrical steel
projectile falls under gravity, g. The freefalling pro
jectile (with drop height Hp 2 m or 0.5 m) im
pacts a 0.42kg polycarbonate buffer inserted into the
top of the tube rather than directly hitting the bubbly
liquid surface. Stress waves in the tube are measured
using six strain gauges (SR4, Vishay; denoted by gl
to g6 in figure 2) placed at intervals of 100 mm along
the tube and oriented in the hoop direction; the sig
nals are processed using a signal conditioning ampli
fier (2300 System, Vishay), and are stored in a digital
recorder (NI 6133, National Instruments; sampling rate
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
ao 0.0013 (before)
9, x
Strain
gauges
co = 0.0013 (after)
ro = 0.0081 (after)
Figure 2: Schematic of the experimental setup.
2.5 MHz).
Method of bubble generation. The bubbles are cre
ated using a bubble generator consisting of an aluminum
plate and capillary tubes (TSP020150, Polymicro Tech
nologies; inner diameter 20 pm). Up to an initial void
fraction of ao 0.0056, 91 capillary tubes are used;
for higher void fractions, the number increases to 217.
The capillary tubes are located in the drilled holes of the
plate and are fastened with epoxy. One side of the plate
is tightly covered with a chamber. The chamber is pres
surized, and the air is injected, due to the pressure head,
into the fluid column. The injected bubbles rise upward
to the column surface, and eventually escape from an air
outlet in the buffer.
Distilled water is used for the case of no air injection;
otherwise, tap water is used. Note that the number of
tiny bubbles present in tap water is negligible compared
to that of the injected air. The water temperature is kept
23 C so that the vapor pressure is much smaller than
one atmosphere.
Images of the bubbles are captured by a highspeed
video camera (Phantom v7.3, Vision Research). A white
LED lamp (Model 900445, Visual Instrumentation Cor
poration) is used for backlighting. A water jacket is at
tached outside the tube to minimize image distortion.
Images of the initial bubbles and the collapsing bub
Figure 3: Images of the injected bubbles before and af
ter the primary wave passage. The frame rate is 20000
frames per second.
bles due to the shock loading with Hp = 2 m are
shown in figure 3. The initial bubble size is found to
be broadly distributed, and the mixture is nearly homo
geneous. Moreover, for the higher void fraction case,
we observe fission of the larger bubbles. One possible
mechanism for the fission is a reentrant jet as seen in
figure 4 where the compression wave propagates down
ward and the bubbles collapse subsequently.
The initial void fraction is estimated based on the dif
ference in the column height with and without the air
injection. Uncertainty in this measurement is 0.1 mm
except for the case of the highest void fraction, ao
0.01, in which the column surface waves increase the
uncertainty to 0.5 mm. In the waterhammer experi
ments, the following void fractions were tested: ao = 0
(no air injection), 0.0013 0.0001, 0.0024 + 0.0001,
0.0056 0.0001, 0.0081 0.0001 and 0.010 + 0.001.
Buffer dynamics. In the present experiments, the liquid
pressure is unknown. Hence, the buffer velocity b (or
the piston velocity uH in the shock theory) is critical to
estimate to validate the shock theory. Three experimen
tal runs were conducted for each case of H, and ao. For
every run, the buffer position zb was recorded using the
highspeed camera with recording rate of 32000 frames
o = 0.0081 (before)
Mbb = AplAo,
(14)
where Mb is the mass of the buffer and the righthand
side represents the pressure force acting on the bottom of
the buffer. In the linear case, this pressure force may be
approximated by AplAo = pocrjbAo as can be derived
from equations (2) and (3). Integrating equation (14)
once, we get a solution of the form,
Xb bO exp  I
(15)
Figure 4: Sequence of the bubble collapse for Hp 2 m
and ao 0.0081. The frame rate is 51500 frames per
second.
where b0 is the initial buffer velocity and T is the relax
ation time for the exponential decay:
(16)
pocjAo
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
10
8 0 0
6 7.9 0.5 ms
4
2 4
2.9 + 0.2 m s
0
0 0.005 0.01
8
 linear theory
6 o H =2m
p
H =0.5m
SP
4
2 0
2 N  o 0o o
0
0
0 0.005 0.01
ao
Figure 5: The initial buffer velocity (top) and the relax
ation time (bottom).
per second, and the position history was extracted from
the movies with MATLAB image processing.
The buffer dynamics may be described by Newton's
second law (Dashpande et al. 2006; Shepherd and Inaba
2009). For simplicity, the buffer is treated as a rigid body
and wall friction is neglected. The equation of motion of
the buffer is then given by
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Thus, the evolution of the buffer position is expressed by
an exponential function.
The measured buffer positions were fitted to an expo
nential by the leastsquares method. The resulting buffer
velocity 4bo and relaxation time T are presented in fig
ure 5. We find that the buffer velocity 1b0 is only weakly
dependent on the void fraction and therefore seems to
depend only on the drop height Hp; it is 7.9 0.5 m s 1
and 2.9 0.2 ms 1 for Hp 2 m and 0.5 m, respec
tively, where the error bounds represent the standard de
viation of all the runs. It follows from figure 5 (bottom)
that the linear theory (16) is in qualitative agreement
with the experimental data and is particularly good for
the case of no air injection. Moreover, the relaxation
time decreases as the initial buffer velocity increases.
Since a stronger shock leads to a more violent bubble
collapse, the buffer momentum may decay more rapidly
due to bubbledynamic energy dissipation.
Primary wave speeds. We examine the shock (or pri
mary) waves in the experiments that are characterized
by the drop height Hp and the initial void fraction ao.
To confirm repeatability in the measurements, three ex
perimental runs were conducted for each case of Hp and
o0.
The evolution of the hoop strains for the case of
Hp 2 m and ao 0 is shown in figure 6 (top). Every
strain gauge records the primary wave following a small
amplitude precursor that propagates essentially with the
sonic speed of the tube material (Skalak 1956). It also
records a wave reflected from the tube bottom, and sig
nificant wave dispersion. Moreover, the primary wave
noticeably decays within the measurement period; thus
the decay time is comparable to the relaxation time of
the piston velocity computed in figure 5. Consequently,
the unsteadiness of the buffer dynamics cannot be ig
nored.
In figure 6 (bottom), the primary wave speeds are
computed. For comparative purposes, three thresh
old strain values (30, 40 and 50 percent of the max
imum strain measured at the strain gauge gl before
the reflected wave is observed) are used to determine
the positions of the wave fronts. The wave speed
was then obtained from the slope of a linear least
squares fit to the wave front positions. The computed
speed (524 ms 1) is in reasonable agreement with the
KortewegJoukowsky wave speed (cj 553 ms 1),
and the the dispersion due to the thresholding is very
small. This suggests that the linear theory is effectively
valid for the case of pure water, even though the wave is
dispersive and unsteady.
The bubbly water case (Hp = 2 m, ao 0.0081)
is presented in figure 7. Due to the bubble dynamics,
the structural response manifests more complex struc
I .\ *M i '
I. I ..
0 1 2
t (ms)
0.6
+ measured 1
o measured 2
x measured 3
E 0.4  fitted
a 0.2
Slope = 524 + 2 m s
0 0.2 0.4 0.6 0.8 1
t ti (ms)
Figure 6: Evolution of hoop strains (top) and the loca
tions of the primary wave fronts (bottom) for Hp 2 m
and ao 0. The dotted lines (top) denote the threshold
values used to determine the wave fronts.
tures than in the pure liquid case. The comparison of
the figures 6 and 7 (top) reveals that the bubbles reduce
the tube deformation. Moreover, the wave speed is re
duced by the bubbles. To further see the effect of the
bubbles, the case of the lower drop height (Hp 0.5 m,
ao 0.0081) is presented in figure 8. In this case,
the wave propagation is evidently unsteady. The un
steadiness may also result from the fact that the relax
ation time for the piston velocity and the measurement
time are comparable. As a result of the unsteadiness, the
threshold value becomes more critical, and the standard
deviation in the computed wave speed becomes larger.
We should note that the lower piston velocity reduces
the wave speed. This is the effect of the gasphase non
linearity as pointed out in figure 1.
We now compare the steady shock theory to the ex
perimental data. Despite the unsteady buffer dynamics,
we use the initial buffer velocity b0 = 7.9 m s 1 (or
2.9 ms 1) for Hp 2 m (or 0.5 m) to complete the
steady shock relations. Since vapor pressure at the room
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
1 2
t (ms)
0.6
S0.4
S0.2
0.5 0.1
t tgi (ms)
Figure 7: As figure 6, but with Hp 2 m and ao
0.0081.
temperature is negligible compared to the atmospheric
pressure and surface tension may not be very important
for the bubble sizes in the experiments, we neglect the
vapor pressure and the surface tension in computing the
shock relations. Figure 9 compares the measured wave
speeds to the theory including both the sonic and shock
speeds. The error bars are small (as seen in figures 6 to
8) and omitted for clarity. The measured speeds for the
cases with air injection clearly show differences from
the sonic speeds, and those for H = 2 m are larger than
those for H = 0.5 m. These indicate the effect of the
gasphase nonlinearity on the wave speeds. The shock
theory with FSI is found to more accurately capture the
trend with increasing ao than the shock theory without
FSI and the linear theory. It is therefore concluded that
both FSI and the nonlinear effect need to be considered
to accurately estimate the propagation speeds of finite
amplitude waves in mixturefilled pipes. However, the
agreement between the present theory and the experi
ments is qualitative rather than quantitative.
One of the most obvious limitations in the theory
is related to the assumption of steady wave propaga
o
S0.4
0.2
0
0.6
S0.4
H 0.2
t (ms)
+ measured 1
o measured 2 x +
x measured 3
fitted x
x slope = 200 10 m s
xO
0 1 2
t tgi (ms)
Figure 8: As figure 6, but with Hp 0.5 m and ao
0.0081.
shock (w/FSI, u =79ms1)
shock (w/o FSI, u 9 9ms 1)
Shock (w/FSI, u= 29ms 1)
1500
S_ .shock (w/oFSI, u=29m s1)
 somc (w/ FSI)
\ .. .... somc (w/o FSI)
o o o measured (H = 2 m)
S* + measured (H = 0 5 m)
S1000
r,
500
0 i. i . .
0 0.002 0.004 0.006 0.008 0.01
co
Figure 9: Comparison of the theoretical and measured
wave speeds.
tion. The strain evolution in figures 6 to 8 (top) demon
strates unsteady wave propagation in the sense that the
piston velocity decays during the measurement period.
Hence, the unsteadiness will deteriorate the model valid
ity. If one quantifies the decay rate in the relaxation pro
cess, both the structural and the bubbledynamic damp
ing need to be included in the thoery. Other limitations
are the neglect of the tube inertia and the quasione
dimensional assumption. Since the waves are dispersive
(due to the tube inertia) and of finite wave length, the
wave speed as a function of the wave length cannot be
accurately predicted by the current model. To quantify
the effect of the tube inertia, a fourequation model de
scribing both the fluid and the tube dynamics could be
calculated (Tijsseling 1996).
Conclusions
A quasionedimensional conservation law governing
dilute bubbly flows in a deformable cylindrical tube and
the steady shock relations are reviewed, and the water
hammer experiments are conducted. The present FSI
shock theory is found to be in better agreement with
the measured wave speeds than the linear theory or
the shock theory without FSI. This suggests that both
FSI and the gasphase nonlinearity need to be included
to accurately predict the propagation speeds of finite
amplitude waves in mixturefilled pipes.
Acknowledgements
The authors would like to express their thanks to
T Nishiyama for his help with the experimentation, J. S.
Damazo and R. Porowski for the bubble images, and
S. Hori for his observations about the experimental data.
This work was supported by ONR Grant No. N00014
0610730.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
C. E. Brennen. Fundamentals ofMultiphase Flow. Cam
bridge University Press, 2005.
R. H. Cole. Underwater Explosions. Princeton Univer
sity Press, 1948.
V S. Dashpande, A. Heaver, and N. A. Fleck. An under
water shock simulator. Proc. R. Soc. A, 462:10211041,
2006.
M. S. Ghidaoui, M. Zhao, D. A. Mclnnis, and D. H. Ax
worthy. A review of water hammer theory and practice.
Appl. Mech. Rev., 58:4976, 2005.
K. Inaba and J. E. Shepherd. Flexural waves in fluid
filled tubes subject to axial impact. J. Pressure Vessel
Technol., 132:021302, 2010.
N. E. Joukowsky. Memoirs of the Imperial Academy So
ciety of St. Petersburg. Proc. Amer Water Works Assoc.,
24:341424, 1898.
T. Kobori, S. Yokoyama, and H. Miyashiro. Propagation
velocity of pressure wave in pipe line. Hitachi Hyoron,
37:3337, 1955.
D. J. Korteweg. Ueber die Fortpflanzungs
geschwindigkeit des Shalles in elastischen R6hren.
Annalen der Physik und Chemie, 5:525542, 1878.
J. E. Shepherd and K. Inaba. Shock loading and fail
ure of fluidfilled tubular structures. In A. Shukla,
G. Ravichandran, and Y. D. S. Rajapakse, editors, Dy
namic Failure of Materials and Structures. Springer,
2009.
R. Skalak. An extension of the theory of water hammer.
Trans. ASME, 78:105116, 1956.
P. A. Thompson. CompressibleFluid Dynamics.
McGrawHill, 1972.
A. S. Tijsseling. Fluidstructure interaction in liquid
References filled pipe systems: A review. J. Fluids Struct., 10:109
146, 1996.
K. Ando. Effects ofpolydispersity in bubbly flows. PhD
Thesis, California Institute of Technology, 2010.
K. Ando, T. Sanada, K. Inaba, J. E. Shepherd, T Colo
nius, and C. E. Brennen. Shock propagation through a
bubbly liquid in a deformable tube. Submitted to J. Fluid
Mech., 2010.
C. E. Brennen. Hydrodynamics ofPumps. Oxford Uni
versity Press, 1994.
C. E. Brennen. Cavitation and Bubble Dynamics. Ox
ford University Press, 1995.
E. B. Wylie and V L. Streeter. Fluid Transients in Sys
tems. Prentice Hall, 1993.
Z. D. Zhang and A. Prosperetti. Ensembleaveraged
equations for bubbly flows. Phys. Fluids, 6:29562970,
1994.
