7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Cavitating flow simulation with a sharpinterface model
E. Lauer*, X.Y. Hu*, S. Hickel* and N.A. Adams*
Institute of Aerodynamics, Technische Universitat Miinchen, Munich, Germany
Eric.Lauer@aer.mw.tum.de
Keywords: cavitation, nonequilibrium, sharp interface, levelset, conservative
Abstract
We present a sharpinterface multifluid flow approach including an efficient evaporation/condensation model. The
model assumes that phase change is in thermal nonequilibrium, and is much slower than the hydrodynamic interaction
described by the Riemann problem. This extension gives us the opportunity to look at a wide range of problems with
phase change, including cavitating flows. Several multidimensional numerical examples are investigated. First we
present results which show the influence of different phase change rates. Then we demonstrate the capabilities of our
model for the nonspherical cavitation bubble collapse close to a solid wall.
Introduction
Cavitation is known as one of the major reasons for dam
ages in hydraulic applications. The collapse of cavita
tion bubbles is the main cause of erosion on marine pro
pellers or on blades in turbomachinery. On the other
hand, the destructive effects of cavitation can be ex
ploited in medical application such as the shock wave
lihliips
Two fundamental mechanism are observed in the cav
itation bubble collapse process. First Rayleigh (1917)
described the release of shock waves upon bubble col
lapse, second Korfeld & Suvorov (1944) observed the
development of a reentrant jet in nonspherical bubble
collapse like it occurs close to solid walls. Both, the im
pact of shock waves and of the reentrant jet on a surface
might lead to erosion.
Progress in highspeed photography allowed more de
tailed experimental investigations. Lindau & Lauter
bor (2003) visualized the propagation of shocks and
the generation of the reentrant jet for the bubble col
lapse close to a solid wall. Furthermore they show the
formation of a counter jet during the rebound, indicating
that there is still a lack in understanding of the mech
anisms in cavitationbubble collapse. However, more
detailed experimental measurements are difficult due to
small spatial and temporal scales. Thus numerical simu
lation might help to get a better understanding, but cav
itating flow is still a challenging problem for modem
computational fluid dynamics.
In this paper we present a sharpinterface method,
which is based on an earlier approach of Hu et al. (2006).
The material interface is accurately resolved on arbitrary
Cartesian grids and therefore a good insight into inter
face movement and deformation is provided.
We suppose that compressible effects such as shock
wave propagation are important in vapor and water at
least at the final stages of collapse. Viscosity and sur
face tension on the other hand, are neglected in this pa
per, since Johnsen & Colonius (2009) show that their
influence is small over the major part of a vapor bubble
collapse.
We apply the Euler equation for both fluids separately
using a standard finitevolume scheme, which is modi
fied by considering computational cells being cut by the
interface. The discretized governing equations are up
dated conservatively. The method solves the difficulty of
conservation for fronttracking methods and treats topo
logical changes naturally by combining the interface de
scription and geometric operations with a levelset tech
nique.
Interactions between both fluids at the interface are
taken into account by solving the Riemann problem at
the interface and introducing a conservative exchange
term. While the original model deals only with mo
mentum and energy exchange over the interface, the new
approach takes into account the effect of phasechange
phenomena such as evaporation and condensation.
After first results obtained on onedimensional
equidistant Cartesian grids presented by Hu et al. (2008),
we now present the fully threedimensional version of
our method working on arbitrary Cartesian grids.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Governing equations
We consider the threedimensional Euler equations in
conservative form:
au
S+V.U=0
with the solution vector U [p, pul, pu, pus, E].
Density p, momentum pu and energy E are the con
served variables. The flux vector is given as
Fi (U)
uip
uipul + Silp
uipu2 + 6i2p
UiPU3 + 6i3p
Ui (E + p)
where u is the velocity vector.
Figure 1: Two dimensional schematic of conservative
discretization for a cut cell.
Ideal gas
In case of an ideal gas, the pressure is determined by
the definition of total energy
E = p + pu2 (3)
7 1 2
The ideal gas equation of state in a dimensionless form
closes the problem
T = nMa21 (4)
P
with the Mach number Ma and the ratio of specific
heats 7. We model vapor as an ideal gas and we assume
7 1.335.
Waterlike fluids
For waterlike fluids we apply Tait's equation of state.
The temperature T is assumed to be constant, and pres
sure p and total energy E become functions of density
P
p B P) B+A
E1Po A
E (p + B A) + B
71
The additional coefficients B, A and po are taken as con
stants. For water we set B 3310 bar, A 1 bar,
Po 1 kg/m3 and 7 7.15.
Since the total energy is only a function of density and
velocity, energy conservation in the Euler equations (Eq.
(1) and (2)) is omitted.
Overview of the method
Mathematical basis
The computation domain 2 is divided by the interface
F(t) into two different subdomains Ql(t) and Q2(t).
The volume Q 1(t) accounts for the region occupied by
a first fluid and Q2(t) refers to the region occupied by
a second fluid. We solve Eq. (1) for the fluid occu
pying the subdomain Qm(t) (m = 1, 2) on a three
dimensional grid with spacings Ax(i) Ay(j) Az(k). In
tegrating Eq. (1) over the spacetime volume V n Qm (t)
of a computational cell (i, j, k) we obtain the following
finite volume discretization
n+1
*In
dt f n (t)
Jvonm(t)
dx dy dz i +V F(U,)
8I
where U, is the the solution vector of the considered
fluid. The fluid volume V n m,(t) can be expressed
through am(t)AxAyAz where am(t) is the time de
pendent fluid volume fraction of the considered cell and
1 fluid m with al(t) + a2(t) 1 Applying Gauss theo
A+ 2pu (6) rem we obtain
n+l
Sn +
+ I
di
dtvnn(t)
dt (
a(vnom
ogUr
dx dy dz 
at
dx dy dz F (U,) n
(t))
A/F
/I
1,1,:1 n
\. f 1 ::
l i,' ....... ..... ...#1.l1..
/a
cT/l
I i,
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
where d(V n Q, (t)) are the cell faces. Defining the in
terface through F(t) = V n dOQ(t), the term d(V n
R, (t)) can be expressed as the combination of two ele
ments.
The first element is the sum of the surfaces ob
tained by intersecting the interface with the cell
faces, which can be expressed in the following
form AIAyAz, A)AyAz, AAxAz, AAAxAz,
A~AxAy, A"AxAy where A' represents the so
called cell face aperture described in Fig. 1 and with
A4 + ,i 1 1. The second element represents the por
tion of the interface AF(t) within the cell (i, j, k). After
introducing the volume averaged solution
S 1 (t)V ()
am (t)V Jam(t)V
Figure 2: Local reference system for a cut cell.
U dx dy dz (9) as a local modification of the standard finite volume
scheme close to the interface.
Eq. (8) can be rewritten in the following form
V (c4 1u4 1
fn+1
fjn 1
+ I
dt AyAz r. 1I/)Fm
dt AxAz .92 iF22
dt AxAy [A (t) F
dt X, (AF(t))
jA '()F
A1 (t)F
where cmUm and Ur are the vector of the conse
quantities in the cut cell and the vector of volume
aged conservative variables respectively. FP is th
erage flux across a cell face and Xm is the average n
momentum and energy exchange across the inter
For full cells, which are not cut by the interface, vol
fractions and cell face apertures become unity, the
responding interface surface fraction vanishes and
(10) simplifies to
Geometry of the interface
Cell face apertures APm(t), volume fractions am(t)
(10) and interface area AF(t) are determined from al levelset
field b, i.e., a signed distance function from a local
1] point of the domain to the immersed interface. The
zerolevelset contour ( 0) describes the interface
21] between the two fluids. Cells that are completely inside
S the other fluid are blocked through a zero value for the
cell face apertures and consequently they will not have
31] any contributions during the computation. The levelset
field allows for the treatment of arbitrary geometries
and can deal with moving and strongly deforming
interfaces, where the levelset field is updated following
Fedkiw et al. (1999) and the volume fractions and face
*rved apertures are reconstructed at each time step.
aver
e av Interface interaction
aass, The interface interaction term Xm (AF) contains two
face. elements accounting for the contributions of pressure
lume and mass transfer across the interface.
cor
Eq.
X,(Ar) = XP + X1
V(U 1 U )
vm Um)
n+1
rn+1
nn+1
dt AyAz [F)
dt AxAz [F7
dt AxAy [Fm
recovering a standard finite volume scheme on a tl
dimensional Cartesian grid.
Being effective only in the cells that are cut by
interface, this immersed interface method can be
F21
' I
F"(
In the presence of viscosity, interface friction can be
taken into account by introducing a noslip correction
on the interface. In case of high temperature differences
between both fluids or large time scales, an additional
heat transfer term would be needed as well, but both
terms are neglected in this paper.
I Pressure term XP
(11) Following an approach proposed by Hu et al. (2006),
the evolution of the interface is determined by the inter
ree face condition given by a twomaterial Riemann prob
lem
/ the
seen
on r (13)
an Tr
m M)
J
91(U1, U2) 0
To obtain the momentum and energy exchanges
across the interface, the Riemann problem is solved
in every cut cell in interface normal direction. For
problems with only weak or moderate shocks and ex
pansion waves, simple noniterative approximate Rie
mann solvers such as those of Toro (1997) are sufficient
for reasonable results, otherwise more accurate iterative
solvers are needed. In this paper, the interface condition
is obtained by the interface interaction method of Hu &
Khoo (2004).
After solving the Riemann problem for the interface
pressure p, and the interface normal velocity us, the
contributions of the term XPm to momentum and energy
equation are computed as follows
0
piAF (nm i)
piAF (nm 3)
pIAr (nm uk)
PIAr (nm uj)
(14)
The overall conservation property is satisfied since
nl = n2, i.e., the transfer quantities for the two fluids
have the same magnitude but opposite sign.
Mass transfer term Xm
In case of vaporliquid interactions with phase
change, we apply a mass transfer term Xm at the in
terface. It covers not only the mass transfer across the
interface but also the momentum and energy exchange
due to evaporation or condensation. The subscript for
differentiating between the two different fluids at an in
terface m, is set to m v, li to identify vapor and liquid
respectively.
To model the phase change at the interface, we ap
ply the thermal nonequilibrium assumption proposed
by Fujikawa & Akamatsu (1980). During the phase
change process the pressure is in equilibrium at the va
por pressure and the temperature has a discontinuity at
the phase interface. An expression for the rate of evap
oration and condensation m can be found in the book of
Schrage (1953). We use the dimensionless form
m = ')Ma
v'Z2T ( Vv'
P
x V )
(15)
Here, 7 and Ma, are the ratio of specific heats and the
Mach number in the vapor phase; a is the accommo
dation coefficient for evaporation or condensation (as
sumed constant), equal to the ratio of vapor molecules
sticking to the phase interface to those impinging on it;
T, and Ti are the temperatures of the vapor and the liq
uid at the phase interface, respectively; p, is the actual
vapor pressure at the interface, and ps (Tij) is the equilib
rium (saturation) vapor pressure at temperature Ti and
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
is obtained from the ClausiusClapeyron equation in a
dimensionless form
Ps(Tui) poexp V[Ma LV ( 1
L V\(o
a)]
where Lv is the latent heat of vaporization, and po and
To are given by a reference saturation state. With Eq.
15, the phasechange induced interface velocity is
Aq*
where pli is the density of liquid at the interface. In cav
itating flow, Aq* is usually small and much smaller than
the advection velocity of the interface (Brennen (1995)),
which suggests that the change rates of pressure and
temperature caused by the latent heat are iiniginilik.il
Therefore, it is straightforward to assume that the phase
change is much slower than the hydrodynamic interface
interaction. If the effect of heat conduction in the bulk
are also neglected, one can simply obtain pV, put, T1i and
Tv from the solution of the Riemann problem (13).
In case of evaporation the mass transfer term of the
vapor domain Xm is finally expressed as
xM7
r AF
r AF (Ui,. i)
r4 Ar (Ui U)
rm Ar (Ui t k)
m AF (e, + 1 ui 2) + i Aq* AF
2F
In case of condensation
/ m A
r AF (U i)
Xm 4 Ar (u j)
m AF (uv k)
,j, Ar (,v + 1 I v 12)
Km AF e + u Q+ pi Aq* AF
(19)
where ul, and uv are the velocity vectors of the liquid
and the vapor in a considered cut cell, respectively; e, is
the internal energy of the vapor in this cell.
The contribution to the mass transfer term of the liq
uid domain Xm follows from conservation
XiL Xm (20)
Note that the energy exchange is usually of no im
portance for the liquid since the introduced temperature
change is negligible.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
] Fluid 1 regi_
SFluid 2 regi
Si+1, j ktargetcell
Interface
n
ijk
Al I k
Figure 3: 2D overview of mixing target cells.
Mixing Procedure
The intersection between the immersed interface and
the Cartesian grid can give rise to small cells, for which
a stable fluid state may not be reached based on the time
step calculated according to the full cell size. On the
other hand, accounting for the real size of these small
cells while computing the time step according to the
CFL condition can lead to extremely small time step
sizes.
In order to ensure numerical stability without exces
sively decreasing the time step size, a conservative mix
ing procedure is introduced, where the conserved quan
tities of small cells are mixed with larger neighbor cells.
The procedure follows the approach proposed by Hu et
al. (2006). Modifications are introduced in order to ac
count for non uniform grid spacing. The number of tar
get cells adopted for the mixing is also increased and
now includes all the possible neighbor cells available in
a three dimensional field.
The mixing procedure is applied only to cut cells with
a volume fraction below a certain threshold, which in
our computations is set to am 0.6. Increasing this
threshold leads to a better numerical stability but at the
same time makes the method more diffusive and vice
versa. Each fluid will be considered separately (m =
1,2).
Seven mixing target cells are determined from the lo
cal normal vector on the interface. For each mixing tar
get cell a weight 3m is defined
3 = Im 2 cmiXuj,k (21)
S2 .
/% Ili3^i ai,mixi,,f
m = Jnm 0nm tk
2
=% nm k i
^ = (l(nm. ) (nm IJ) amiXi,mixf,k
JY (n .).) (nm ) k 2
P/ .(ix z) (nm .) (n ( i) 2
*mixi,mij ,miXk
To avoid double mixing between two cells, the
weighting factor is set to zero whenever the cell's vol
ume fraction is larger than that of the target cell. The
weights are subsequently normalized to ensure consis
tency
Z 9 1 (22)
trg
where trg stands for trg {x, y, z, xy, xz, yz, xyz}.
The mixing flux 1'.' is calculated for each mixing
direction in the following way
"I.
am trg
Oim VPm Qtrg,m Vtrg
(* . VyUtry,m) amV
(amVUm) atry,mVtry
Conservativity is maintained by a flux formulation,
consistent with nonuniform Cartesian grids
Um
Utrg,m
* O 'I
* _1 1
amV
atrg,mVtrg
where (Um) and (Ut~g,m) are the volume averaged
conserved variables of the cut cell and of the target cell
before mixing. This procedure is applied accordingly
for all cut cells and mixing target cells at the end of each
RungeKutta substep.
It is interesting to see how the present mixing
procedure deals automatically with vanished and newly
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Solid wall
vaporwater interaction
created empty cells. In the first case, the residual
conservative quantities are all transported to the target
cells, whereas for the second case the conservative
quantities in a newly created small cell are derived
directly from its target cells.
Extending procedure
In order to allow for using unmodified, central in
terpolations stencils in the finite volume reconstruction
scheme close to the interface and to provide physical
conditions for newly created cells, the volume averaged
variables are extended across the interface using
BU,
f n VU, 0 (26)
Otf
Following Hu et al. (2006), Eq. (26) is solved in the
pseudotime tf until a steady state solution in the near
interface region is reached.
Numerical results
The following numerical examples are provided to il
lustrate the potential of our conservative sharpinterface
method with phasechange model. We use a fifthorder
WENO scheme of Jiang & Shu (1996) and a thirdorder
TVD RungeKutta scheme of Shu & Osher (1988) to
discretize the Euler equations. Computations are carried
out with the CFL number of 0.6.
All test cases are twofluid problems containing vapor
and water. We always model vapor as ideal gas and use
Tait's equation of state for water. Material properties are
given above in the correspondent material law section.
1 D Vaporwater interaction
We reconsider the 1D vaporwater interaction case
presented in Hu et al. (2008). The test problem is the
condensation of oversaturated water vapor in the pres
ence of an acoustic wave. Initially, as shown in Fig. 4,
both phases have a length of 2 cm. Vapor is on the left
side and water on the right side, with the phase interface
at the origin. Both phases are at rest and have a common
0.25 0.5
Time (ms)
0.75
Figure 5: 1D Vaporwater interaction: Evolution of the
phase interface with different accommodation
coefficients.
temperature of 293 K at the beginning. The initial va
por pressure is 93 mbar, whose saturation temperature is
343 K. The water pressure is 193 mbar. The ends of the
domain are reflecting boundaries. The number of grid
points is 200.
Figure 5 gives the evolution of the phase interface
with different accommodation coefficients: a = 0, 0.075,
0.25, 0.5 and 0.75. One can find, that if there is no phase
change (a=0), the water first expands due to its higher
pressure and then shrinks and expands periodically as
the pressure wave propagates back and forth in the liq
uid. When the phase change is included in the process,
the liquid volume still oscillates, but the oscillation is su
perposed upon an expansion due to condensation as the
vapor is oversaturated. Figure 5 also gives the depen
dence of the evolution of the interface on the accommo
dation coefficient. It can be noticed that for large a the
expansion slows down at later times due to faster con
suming of the oversaturated vapor.
Figure 6 plots the velocity profile at t = 1 ms. One
can find that the discontinuity of velocity in the center
gives the mass flux across the phase interface. With
phase change considered, as the consuming of over
saturated vapor, the pressures of both vapor and water
decrease considerably (see Fig. 7).
These observations suggest, that phase change plays
an important role in case of a deviation from the phase
equilibrium.
3D bubble oscillation
Next test case is the growth and collapse of a three
dimensional high pressure gas bubble of radius 16 cm in
water. The setup shown in Fig. 8 is based on the work
Solid wall
(p4, v, Pv,Tv)
Oversaturated
vapor
Figure 4: Sketch of the 1D
problem.
(pl,qi,Pl,T1)
Water
a=0.0
VAa= 0.075
a 0.25
S0.5
c= 0.75
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
solid wall
0.6
0.2
o&
2 1.5 1 0.5 0 0.5 1 1.5 2
x (cm)
Figure 6: 1D Vaporwater interaction: Velocity profiles
at t = 1 ms for with and without phase change.
0.16
' 0.12
, 0.08
0.04
2 1.5 1 0.5 0 0.5
x (cm)
Figure 7: 1D Vaporwater interaction: P
at t = 1 ms for with and without
water
vapor
16 cm
Figure 8: 2D sketch of the 3D bubble oscillation prob
lem. Only one eighth of the bubble is consid
ered.
of Wardlaw (1998). We take advantage of the symmetry
of the full bubble problem and consider only oneeighth
of the gas bubble in a domain with a length of 10 m in all
directions. Reflecting boundary conditions are imposed
at the domain boundaries.
In contrast to the original setup we use vapor as in
ternal gas in the bubble, since we are interested in seeing
the influence of phase change on the oscillation ampli
tude and frequency. The initial material states for this
problem are:
vapor: p = 78039 bar, p = 1.63 g/cm3
water: p = lObar, T = 293 K
Initial velocity is zero everywhere. Initial vapor temper
ature and water density are given by the corresponding
equations of state.
We investigate three different levels of grid refine
ment corresponding to 40 (coarse), 80 (medium) and
160 (fine) cells in each direction. For all refinement lev
els we use the same hyperbolic stretching. Fig. 9 gives
an impression of the coarse grid and shows the initial
1 .5 2 size of the vapor bubble.
Results for the evolution of the bubble radius on dif
ferent grids without phase change (a 0.0) are given
pressure profiles in Fig. 10. We see that the vapor bubble starts to grow
t phase change, due to large internal pressure until a maximum size of
approximately 1.8 m on the fine grid is reached. At this
point the water pressure is already higher then the in
ternal vapor pressure leading to a collapse process. An
oscillation develops since the collapse is stopped when
the vapor pressure reaches a critical value. Due to mas
, =0.075
I
I
ct 0.075
I' '
1I~
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Figure 9: 2D view of the coarse grid for the 3D bubble
oscillation problem. Initial bubble radius is
shown in the bottom left corer.
2
g 1.5
1
3
0 30 60
Time (ms)
Figure 10: 3D bubble oscillation: Evolution of the bub
ble radius without phase change (a 0.0)
on different grids.
Time (ms)
Figure 11: 3D bubble oscillation: Evolution of the bub
ble radius for different accommodation coef
ficients on medium grid (80x80x80).
sive compression towards the end of the collapse, the
bubbleradius growth rate is much higher close to the
bubblesize minimum than at its maximum. As we use
reflecting boundaries, there are deviations in each oscil
lation cycle, since the considered problem is no longer
axisymmetric after the initial shock has been reflected
at the boundaries of the domain. The vapor bubble is
not perfectly spherical anymore. Thus we calculate the
bubble radius based on the bubble volume, although the
deviation from a perfect sphere is small.
The results fit qualitatively with results by Kadioglu
& Sussman (2008) who used a twodimensional coupled
level set and volumeoffluid method assuming the bub
ble growth and collapse to be axisymmetric. A quantita
tive comparison is not possible as they don't use vapor as
internal gas. This has very strong influence on the oscil
lation amplitude and frequency due to different material
properties.
Comparing the results of different levels of grid re
finement in Fig. 10, we notice that the oscillation ampli
tude increases while the the time for one complete cycle
decreases with finer mesh resolution due to two differ
ent reasons: Less dissipation in the numerical scheme
and a sharper representation of the initial bubble which
leads to a higher initial vapor mass on fine grids. We
also observe that the relative error for increasing grid
resolution decreases strongly which indicates highorder
convergence for our method.
Results for the influence of phase change on the bub
ble behavior can be found in Fig. 11. Increasing the ac
commodation coefficient leads to a faster phase change
and therefore to a lower oscillation amplitude as vapor
fine
= 00 medium
coarse
,\ //
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Y
10 mm
symmetry
outlet
10 mm
Figure 12: 3D bubble oscillation: Evolution of the bub
ble radius with phase change (a 0.1) on
different grids.
pressure is reduced due to condensation. Additionally,
lower vapor pressures slightly decrease the time period
for one complete oscillation cycle.
Since our method ensures conservation of mass for
both phases when phase change is switched off (a
0.0), the bubble radius oscillates around a constant av
erage value. The amplitude is damped only by numer
ical dissipation which decreases with higher grid res
olution (Fig. 10). In contrast, the average value de
creases continuously as soon as the phase change model
is used (a / 0.0). For high accommodation numbers
(a = 0.1, 0.05) the bubble already reaches the grid reso
lution limit in the first three oscillation cycles. The vapor
is completely condensed and the bubble disappears (Fig.
11).
Fig. 12 provides an additional resolution study, this
time on the problem with phase change (a = 0.1).
Again, convergence for our method can be clearly seen,
as the relative error between successive grid resolutions
decreases. We can further conclude that the final bubble
collapse without rebound, already seen in Fig. 11, is
physically sound since total condensation is achieved
for all grid resolutions in the same oscillation cycle.
3D vapor bubble collapse close to a solid wall
The collapse of a threedimensional vapor bubble
close to a solid wall is considered next. The initial bub
ble radius is 400 pm, while the offset of the bubble cen
ter from the wall is 416 pm in ydirection (see Fig. 13).
We take advantage of symmetries and calculate only
one quarter of the full test problem. Like in the work
of Niederhaus et al. (2008), outlet boundary conditions
are enforced on the farfield bounds of the computational
Figure 13: 2D sketch of the 3D vapor bubble collapse
close to a solid wall problem. Only one quar
ter of the bubble is considered.
domain. The outlet condition applies a zerothorder ex
trapolation to the boundary, so that gradients across the
boundary are zero. Thus shock reflections are mini
mized but not eliminated completely. We try to avoid an
influence of reflections on the solution by keeping the
boundary far from the region of interest. A Cartesian
grid with an equidistant spacing A 4 pm in all axis
directions is used in the region close to the vapor bubble
while hyperbolic stretching is used in the farfield.
The initial material states are:
vapor: T = 293.0 K,
water: T = 293.0 K.
0.02340 bar
100.0 bar
where the common temperature is the saturation temper
ature corresponding to the initial vapor pressure. Both
fluids are initially at rest. The accommodation coeffi
cient is taken as a = 0.01.
Results for five different times can be found in Fig.
14ae. Only a small section of the total domain is shown,
although the bottom end of the pictures indeed indicates
the solid wall. As only one quarter of the full problem
has been computed, results where mirrored on the (X
Y) and (YZ)plane for visualization.
The initial situation is shown in Fig. 14a in order to
have an impression of the size of the cavitation bubble
at later times.
After 3.95 ps the bubble volume has shrunk to ap
proximately 1/17 of the initial size. The bubble's stand
off distance has not changed much, since there is a re
gion of low water pressure in the gap between bubble
and wall. A cavity starts to develop at the top, driven by
a region of high water pressure above the bubble.
S1.5
0.5
f3
"3
0 1
Time (ms)
water
vapor
solid wall
416 m
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
t 3.95 ps
K
t = 4.28 ps
p [bar]
2500.
2000.
1500.
1000.
V
V
t 4.24 ps
LI
K>A
Figure 14: a)e) Pressure fields at five different time instances for the 3D vapor bubble collapse close to a solid wall
problem. f) Bubble surface after 4.24 ps. Data is mirrored over symmetryplanes.
t = 0.0 s
t 4.24 ps
t = 4.41 ps
I
At time t = 4.24 pus, the reentrant jet is already fully
developed. While the water pressure is very high in the
region above the vapor bubble and at the bottom of the
cavity, we observe even negative pressures in the main
field of the reentrant jet caused by the high velocity in
this area. A 3D impression of the bubble surface at this
time is given in Fig. 14f, where the crossing point of the
black lines indicates the base point of the bubble on the
wall.
When the reentrantjet finally breaks through the bub
ble (Fig. 14d) the vapor bubble starts to collapse com
pletely. Strong shock waves are emitted and impinge on
the wall below. The wall pressure reaches values of more
than 2 104 bar corresponding to more then two hundred
times the initial water pressure.
After the complete bubble collapse and shock reflec
tions at the wall, we observe regions of very low pres
sure in the water. In these regions one would expect a
rebound of the vapor bubble. However, a rebound would
need a nucleation model which is subject of ongoing
work.
Conclusions
We have presented a conservative sharpinterface model
for multifluid flows, working on arbitrary Cartesian
grids. Introducing a phase change model, based on a
nonequilibrium assumption at the phase interface, we
treat evaporation and condensation. Computational re
sults show the influence of the accommodation coeffi
cient which corresponds to the rate of evaporation and
condensation respectively.
The model is currently limited to problems with ini
tially existing interfaces. Further studies will be focused
on modeling the nucleation process.
Acknowledgements
Financial support has been provided by the German
Research Council (Deutsche Forschungsgemeinschaft 
DFG).
References
Brennen C.E., Cavitation and bubble dynamics. Oxford
University Press, 1995
Fedkiw R., Aslam T., Merriman B., Osher S., A non
oscillatory Eulerian approach to interfaces in multimate
rial flows (the ghost fluid method). J. Comp. Phys., Vol.
152, 457492, 1999
Fujikawa S., Akamatsu T., Effects of the non
equilibrium condensation of vapour on the pressure
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
wave produced by the collapse of a bubble in a liquid.
J. Fluid Mech., Vol. 97, 481512, 1980
Hu X.Y., Adams N.A., Johnson E., laccarino G., Model
ing fullMachrange cavitating flow with sharp interface
model. Proceedings of the 2008 Summer Program, CTR,
Stanford, CA, July 6 August 1 2008
Hu X.Y., Khoo B.C., An interface interaction method
for compressible multifluids. J. Comp. Phys., Vol. 198,
3564, 2004
Hu X.Y., Khoo B.C., Adams N.A., Huang F.L., A con
servative interface method for compressible flows. J.
Comp. Phys., Vol. 219, 553578, 2006
Jiang G.S., Shu C.W., Efficient implementation of
weighted ENO schemes. J. Comp. Phys., Vol. 126, 202
228, 1996
Johnson E., Colonius T., Numerical simulations of non
spherical bubble collapse. J. Fluid Mech., Vol. 629, 231
262, 2009
Kadioglu S. Y., Sussman M., Adaptive solution tech
niques for simulating underwater explosions and implo
sions. J. Comp. Phys., Vol. 227, 20832104, 2008
Kornfeld M., Suvorov L., On the destructive action of
cavitation. J. Appl. Phys., Vol. 15, 495506, 1944
Lindau O., Lauterborn W., Cinematographic observation
of the collapse and rebound of a laserproduced cavita
tion bubble near a wall. J. Fluid Mech., Vol. 479, 327
348, 2003
Niederhaus J.H.J., Greenough J.A., Oakley J.G., Ran
jan D., Anderson M.H., Bonazza R., A computational
parameter study for the threedimensional shockbubble
interaction. J. Fluid Mech., Vol. 594, 85124, 2008
Rayleigh L., On the pressure developed in a liquid dur
ing the collapse of a spherical cavity. Phil. Mag., Vol.
34, 9498, 1917
Schrage R.W., A Theoretical Study of Interphase Mass
Transfer. Columbia University Press, 1953
Shu C.W., Osher S., Efficient implementation of essen
tially nonoscillatory shock capturing schemes. J. Comp.
Phys., Vol. 77, 439471, 1988
Toro E.F, Riemann Solvers and Numerical Methods
for Fluid Dynamics: A Practical Introduction. Springer,
Berlin, New York, 1997
Wardlaw A., Underwater explosion test cases. Technical
Report IHTR 2069, ADB238684, Office of Naval Re
search, 1998
