Group Title: 7th International Conference on Multiphase Flow - ICMF 2010 Proceedings
Title: P1.41 - Cavitating flow simulation with a sharp interface model
Full Citation
Permanent Link:
 Material Information
Title: P1.41 - Cavitating flow simulation with a sharp interface model Cavitation
Series Title: 7th International Conference on Multiphase Flow - ICMF 2010 Proceedings
Physical Description: Conference Papers
Creator: Lauer, E.
Hu, X.Y.
Hickel, S.
Adams, N.A.
Publisher: International Conference on Multiphase Flow (ICMF)
Publication Date: June 4, 2010
Subject: cavitation
sharp interface
level set method
Abstract: We present a sharp-interface multi-fluid flow approach including an efficient evaporation/condensation model. The model assumes that phase change is in thermal non-equilibrium, 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 multi-dimensional 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 non-spherical cavitation bubble collapse close to a solid wall.
General Note: The International Conference on Multiphase Flow (ICMF) first was held in Tsukuba, Japan in 1991 and the second ICMF took place in Kyoto, Japan in 1995. During this conference, it was decided to establish an International Governing Board which oversees the major aspects of the conference and makes decisions about future conference locations. Due to the great importance of the field, it was furthermore decided to hold the conference every three years successively in Asia including Australia, Europe including Africa, Russia and the Near East and America. Hence, ICMF 1998 was held in Lyon, France, ICMF 2001 in New Orleans, USA, ICMF 2004 in Yokohama, Japan, and ICMF 2007 in Leipzig, Germany. ICMF-2010 is devoted to all aspects of Multiphase Flow. Researchers from all over the world gathered in order to introduce their recent advances in the field and thereby promote the exchange of new ideas, results and techniques. The conference is a key event in Multiphase Flow and supports the advancement of science in this very important field. The major research topics relevant for the conference are as follows: Bio-Fluid Dynamics; Boiling; Bubbly Flows; Cavitation; Colloidal and Suspension Dynamics; Collision, Agglomeration and Breakup; Computational Techniques for Multiphase Flows; Droplet Flows; Environmental and Geophysical Flows; Experimental Methods for Multiphase Flows; Fluidized and Circulating Fluidized Beds; Fluid Structure Interactions; Granular Media; Industrial Applications; Instabilities; Interfacial Flows; Micro and Nano-Scale Multiphase Flows; Microgravity in Two-Phase Flow; Multiphase Flows with Heat and Mass Transfer; Non-Newtonian Multiphase Flows; Particle-Laden Flows; Particle, Bubble and Drop Dynamics; Reactive Multiphase Flows
 Record Information
Bibliographic ID: UF00102023
Volume ID: VID00449
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: P141-Lauer-ICMF2010.pdf

Full Text

7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 -June 4, 2010

Cavitating flow simulation with a sharp-interface model

E. Lauer*, X.Y. Hu*, S. Hickel* and N.A. Adams*

Institute of Aerodynamics, Technische Universitat Miinchen, Munich, Germany
Keywords: cavitation, non-equilibrium, sharp interface, level-set, conservative


We present a sharp-interface multi-fluid flow approach including an efficient evaporation/condensation model. The
model assumes that phase change is in thermal non-equilibrium, 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 multi-dimensional 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 non-spherical cavitation bubble collapse close to a solid wall.


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 turbo-machinery. On the other
hand, the destructive effects of cavitation can be ex-
ploited in medical application such as the shock wave
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 re-entrant jet in non-spherical bubble
collapse like it occurs close to solid walls. Both, the im-
pact of shock waves and of the re-entrant jet on a surface
might lead to erosion.
Progress in high-speed photography allowed more de-
tailed experimental investigations. Lindau & Lauter-
bor (2003) visualized the propagation of shocks and
the generation of the re-entrant 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 cavitation-bubble 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 sharp-interface 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
We apply the Euler equation for both fluids separately
using a standard finite-volume 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 front-tracking methods and treats topo-
logical changes naturally by combining the interface de-
scription and geometric operations with a level-set tech-
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 phase-change
phenomena such as evaporation and condensation.
After first results obtained on one-dimensional
equidistant Cartesian grids presented by Hu et al. (2008),
we now present the fully three-dimensional 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 three-dimensional Euler equations in
conservative form:


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)

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)
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.

Water-like fluids
For water-like 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 B -P) -B+A
E1Po A

E (p + B A) + B

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 sub-domain Qm(t) (m = 1, 2) on a three-
dimensional grid with spacings Ax(i) Ay(j) Az(k). In-
tegrating Eq. (1) over the space-time volume V n Qm (t)
of a computational cell (i, j, k) we obtain the following
finite volume discretization


dt f n (t)

dx dy dz i +V F(U,)

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


Sn +
+ I


dt (

dx dy dz --

dx dy dz F (U,) n



1-,1,:1 n

\. f- 1 ::

l- i,' ....... ..... ...#1.l-1..

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-
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


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
zero-levelset 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.
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.


X,(Ar) = XP + X1

V(U 1 U )
vm Um)




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

' I


In the presence of viscosity, interface friction can be
taken into account by introducing a no-slip 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 two-material Riemann prob-

/ the

on r (13)

an Tr
m M)


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 non-iterative 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

piAF (nm i)
piAF (nm 3)

pIAr (nm uk)
PIAr (nm uj)


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 vapor-liquid 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
To model the phase change at the interface, we ap-
ply the thermal non-equilibrium 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'Z2-T ( Vv'--

x V )


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 Clausius-Clapeyron equation in a
dimensionless form

Ps(Tui) poexp V[Ma LV ( 1
L V\(o


where Lv is the latent heat of vaporization, and po and
To are given by a reference saturation state. With Eq.
15, the phase-change induced interface velocity is


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
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


r AF
r AF (Ui,. i)
r4 Ar (Ui U-)
rm Ar (Ui t k)
m AF (e, + 1 ui 2) + i Aq* AF

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
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



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

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 =

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
=% 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-

Z 9 1 (22)
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


am trg
Oim VPm Qtrg,m Vtrg

(* . VyUtry,m) amV

(amVUm) atry,mVtry

Conservativity is maintained by a flux formulation,
consistent with non-uniform Cartesian grids



* -O 'I
* _1 1


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
Runge-Kutta 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

vapor-water 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

f n VU, -0 (26)

Following Hu et al. (2006), Eq. (26) is solved in the
pseudo-time 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 sharp-interface
method with phase-change model. We use a fifth-order
WENO scheme of Jiang & Shu (1996) and a third-order
TVD Runge-Kutta 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 two-fluid 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 Vapor-water interaction
We reconsider the 1D vapor-water interaction case
presented in Hu et al. (2008). The test problem is the
condensation of over-saturated 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)


Figure 5: 1D Vapor-water interaction: Evolution of the
phase interface with different accommodation

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 over-saturated. 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 over-saturated 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

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 set-up shown in Fig. 8 is based on the work

Solid wall

(p4, v, Pv,Tv)

Figure 4: Sketch of the 1D



VAa= 0.075

a 0.25


c= 0.75

7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 -June 4, 2010

solid wall



2 1.5 1 0.5 0 0.5 1 1.5 2
x (cm)

Figure 6: 1D Vapor-water interaction: Velocity profiles
at t = 1 ms for with and without phase change.


' 0.12

, 0.08


2 1.5 1 0.5 0 0.5
x (cm)

Figure 7: 1D Vapor-water interaction: P
at t = 1 ms for with and without


16 cm

Figure 8: 2D sketch of the 3D bubble oscillation prob-
lem. Only one eighth of the bubble is consid-

of Wardlaw (1998). We take advantage of the symmetry
of the full bubble problem and consider only one-eighth
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 set-up 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



ct 0.075-

I' '


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.


g 1.5


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
bubble-radius growth rate is much higher close to the
bubble-size 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 two-dimensional coupled
level set and volume-of-fluid 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
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 high-order
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

= 00 medium

,\ //

7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 -June 4, 2010

10 mm



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.
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 three-dimensional 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 y-direction (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 far-field 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 zeroth-order 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 far-field.
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.
14a-e. 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 (Y-Z)-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.



0 1

Time (ms)



solid wall

-416 m

7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 -June 4, 2010

t 3.95 ps


t = 4.28 ps

p [bar]



t 4.24 ps



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 symmetry-planes.

t = 0.0 s

t 4.24 ps

t = 4.41 ps


At time t = 4.24 pus, the re-entrant 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 re-entrant 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
When the re-entrantjet 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


We have presented a conservative sharp-interface model
for multi-fluid flows, working on arbitrary Cartesian
grids. Introducing a phase change model, based on a
non-equilibrium 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.


Financial support has been provided by the German
Research Council (Deutsche Forschungsgemeinschaft -


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, 457-492, 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, 481-512, 1980

Hu X.Y., Adams N.A., Johnson E., laccarino G., Model-
ing full-Mach-range 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,
35-64, 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, 553-578, 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, 2083-2104, 2008

Kornfeld M., Suvorov L., On the destructive action of
cavitation. J. Appl. Phys., Vol. 15, 495-506, 1944

Lindau O., Lauterborn W., Cinematographic observation
of the collapse and rebound of a laser-produced 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 three-dimensional shock-bubble
interaction. J. Fluid Mech., Vol. 594, 85-124, 2008

Rayleigh L., On the pressure developed in a liquid dur-
ing the collapse of a spherical cavity. Phil. Mag., Vol.
34, 94-98, 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 non-oscillatory shock capturing schemes. J. Comp.
Phys., Vol. 77, 439-471, 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

University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - Version 2.9.7 - mvs