7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Unsteady Sheet Cavitation on Threedimensional Hydrofoil
A.H. Koop* and H.W.M. Hoeijmakers&
MARIN, Haagsteeg 2, 6708 PM Wageningen, The Netherlands
& Department Mechanical Engineering, University of Twente, PO Box 217, 7500 AE Enschede, The Netherlands
h.w.m.hoeijmakers@utwente .nl
Keywords: Cavitation, numerical simulation, phase transition
Abstract
An unstructured grid computational method for solving the Euler equations for inviscid compressible unsteady flow is employ
ed to predict the structure and dynamics of threedimensional unsteady sheet cavitation as it occurs on stationary hydrofoils,
placed in a steady uniform inflow. For phase transition an equilibrium cavitation model is employed, which assumes local ther
modynamic and mechanical equilibrium in the twophase flow region, which depends on the equations of state of the vapor and
(compressible) liquid only, see (Koop 2008). The method is similar to the method developed by Schnerr, Sezal & Schmidt
(2008) for structured grids. The unsteady behavior of the flow, i.e. the formation of reentrant jets, their impingement on the
cavity surface, the resulting shedding of the sheet cavity and the subsequent collapse of the bubble cloud, is inertia driven and it
is shown that a numerical method based on the Euler equations is capable of capturing the phenomena associated with unsteady
sheet cavitation. Due to the collapse of shed vapor structures strong pressure pulses are generated continually. To be able to
predict the dynamics of the pressure waves, it is essential that the fluid is considered as a compressible medium by adopting
appropriate equations of state for the liquid phase, the twophase mixture and the vapor phase of the fluid. The
threedimensional unsteady cavitating flow about model configuration (3D Twist 11 hydrofoil) is calculated for which experi
mental data is available (Foeth 2008). It is shown that the formation of reentrant flow and of a cavitating horseshoevortex is
captured by the present numerical method. The calculated results agree reasonable well with experimental observations. Fur
thermore, the collapse of the shed vapor structures and the resulting high pressure pulses is demonstrated.
Introduction
Cavitation is the evaporation of a liquid in a flow when the
pressure drops below the saturation pressure of that liquid.
The importance of understanding cavitating flows is related
to their occurrence in various technical applications, such as
pumps, turbines, ship propellers and fuel injection systems,
as well as in medical sciences such as lithotripsy treatment
and the flow through artificial heart valves. Cavitation does
not occur in water only, but in any kind of liquid such as li
quid hydrogen and oxygen in rocket pumps or the lubricant
in a bearing. The appearance and disappearance of regions
with vapor may be a major cause of noise, vibration, erosion
of surface material and efficiency loss in hydraulic machi
nery. In many technical applications high efficiency is re
quired, which results in cavitation being hardly avoidable
within the range of practical operating conditions. When
cavitation occurs it needs to be controlled. Therefore, de
tailed insight is required in the mechanisms that govern the
cavitation phenomena.
The present paper concerns the dynamics and structure of
sheet cavitation, which occurs on hydrofoils. There are a
number of closely related important aspects to sheet cavita
tion:
Shape and volume of the cavity. The topology of a
sheet cavity is strongly related to the load distribution of the
hydrofoil and thus to the pressure distribution on the surface
of the hydrofoil. Variations in volume of the sheet cavity
cause pressure fluctuations in the liquid that might lead to
impact loads on and strong vibrations of nearby structures.
Reentrant flow at the closure region of the sheet
cavity. At the back end of the cavity liquid may flow in up
stream direction, along the surface of the hydrofoil, (partly)
separating the vapor cavity from the hydrofoil surface.
These socalled reentrant and sideentrant wall jets dictate
the way in which the cavity sheet breaks up and sheds from
the hydrofoil. The shape of the closure region of the cavity
sheet governs the direction of the reentrant and sideentrant
jets.
\1,./.l,,, and collapse of vapor structures. The
breakup of a sheet cavity causes a vortical flow of vapor
clouds that convect to regions with higher pressure. Here,
these clouds collapse resulting in strong pressure pulses
leading to unsteady loads on nearby objects as well as noise
production and possibly erosion of surface material.
Since the 1990s numerical methods using the Euler or
NavierStokes equations have been developed to numeri
cally simulate cavitating flows. The main problem in the nu
merical simulation of multidimensional unsteady cavita
ting flow is the simultaneous treatment of two very different
flow regions: (nearly) incompressible flow of pure liquid in
most of the flow domain and lowvelocity highly compres
sible flow of (pure) vapor in relatively small parts of the
flow domain. In addition, the two flow regimes can often
not be distinguished that clearly, for example in the transi
tion region between vapor and liquid, i.e. the (bubbly) mix
ture region of liquid and vapor.
Furthermore, calculations of unsteady threedimensional
cavitating flows require an appropriate highresolution
mesh in the cavitating flow region, which demands sub
stantial computer resources both in terms of memory and of
speed.
In the present research a numerical method for solving the
Euler equations for 3D unsteady cavitating flow is devel
oped. The accurate prediction of the direction and momen
tum of the reentrant and sideentrant jets and their impinge
ment on the cavity surface form the indispensable basis of
an accurate prediction of the shedding of the cavity sheet.
The reentrant jets are thought to be inertia driven, so it is
expected that a mathematical model of the flow based on the
Euler equations is able to capture the major structure of
sheet cavitation.
Preliminaries
The cavitation number o is a dimensionless quantity defined
as:
p p (T)
0  (1)
where p, [Pa], p, [kgm 3] and U. [ms1] are the pressure,
density and velocity, respectively, of the freestream. Fur
thermore, pt(T) is the saturation pressure of water which is
a function of temperature T [K].
The pressure coefficient Cp is defined as the dimensionless
quantity
C P 2 (2)
2 p U2
with p the local pressure in the flow field. Neglecting skin
friction, the force on the hydrofoil can be obtained from
F =J' pdS, (3)
S
with S the surface of the hydrofoil, p the pressure on this
surface and ii the unit normal pointing into the object, i.e.
out of the flow domain. We choose the coordinate system
such that the lift force L is equal to Fz and the drag force D is
equal to Fx. The dimensionless CL and CD coefficients are
defined as
L D
C C D= pU (4)
_pU2S' U2S '
respectively, with S the projected surface area of the hydro
foil.
The void fraction a [] of the vapor within a volume V [m3]
of a fluid follows from the fluid density p:
p = ap,,t (T) + (1 a) pst (T) (5)
Here yV P a (T) (6)
V Psa (T) Psa (T)
with V, [m3] the volume of vapor within the volume Vofthe
fluid and pvsa(T) [kgm3] and i,,saT) [kgm3] the saturated
vapor and liquid density at temperature T respectively.
The dimensionless total vapor volume Vvap is defined as
1 N
Vvap (7)
C ,(
with c the chord length of the hydrofoil and N the total
number of control volumes in the computational domain.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
The Strouhal number St, based on the chord length c is de
fined as
sfC
St fc (8)
U
withf[Hz] the frequency of the periodic cavitation phenom
enon.
Mathematical Formulation
Following Saurel, Cocchi & Butler (1999) and Schmidt,
Sezal & Schnerr (2006) the equilibrium cavitation model is
adopted. This physical model is based on the assumption
that the twophase flow regime can be described as a ho
mogeneous mixture of vapor and liquid that remains in
thermodynamic and mechanical equilibrium. This implies
that in the mixture locally the vapor and liquid have identi
cal temperature, pressure and velocity. Under these as
sumptions, the flow of the mixture can be described by the
homogeneous mixture equations together with appropriate
equations of state that cover all fluid states possible: the
compressible pure liquid state, the compressible twophase
mixture state and the compressible pure vapor state. The
equations of state must be such that the hyperbolic nature of
the system of governing equations is preserved, which is ne
cessary to represent the propagation of pressure waves in
the fluid. The governing equations are the Euler equations
in conservation form for the mixture variables written here
for Cartesian coordinates as:
au aF (U) F,(U) DF,(U)
++ + =0, (9)
Dt ax ay Dz
or in integral conservation form as:
SJIJ UdQ + A F(U) idA = 0. (10)
at Q A=M
Here, U = [p, pu, pv, pw, pE]T is the vector of conserved
variables and F(U) = [F (U),F (U),F (U)] the flux
vector given as
r pu pv pw
F(U)=
pu2 +p
puv
puv
pv2 +p
puw
PVW
puw pVW pW2 + p
pHu pHv pHw
In Eq. (11) the total specific enthalpy H is equal to
H =E+ =h+ iiu (12)
P 2
with the total specific energy E defined by
I 
E=e+uu, (13)
2
and with e the internal specific energy. The normal compo
nent of the inviscid flux vector F(U) i in equation (10)
can now be found by
F(U) i=F (U) n + F(U) n + F(U) n (14)
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
which yields
pu
puu + pnx
F(U)fi= puv+ pny (15)
putw + pnz
puH
where U is the contravariant velocity component normal
to the surface A defined by
S= i ii = un + vn + wn (16)
For closure of the system of equations it is necessary to a
dopt equations of state p = p(p, T) and h = h(p, T) that des
cribe each of the three possible states: the liquid state, the
vapor state and the mixture state. In the following the liquid
phase is denoted by subscript 1, the vapor phase by subscript
v and saturation conditions by subscript sat. The speed of
sound c within each state is obtained from Koop (2008).
(ah)
p .
C2
!ap; (hC l (C p {i _O }ah I (17)
aPJ, Ir'maI,
Liquid Phase
Following Saurel, Cocchi & Butler (1999) a modification of
Tait's equation of state is used to describe the pressure in the
liquid state as a function of the density and temperature:
KP1
P (PT)= Ko [) 1 +p.la(TI) (18)
where for water K = 3.3108 Pa and N= 7.15 are taken to be
constant. A caloric equation of state that is a good approxi
mation, see (Koop 2008), is given by
e, (p,,) = e, (,)= C (T, T)+ e, (19)
where Cvi is the specific heat at constant volume for water,
To is a reference temperature and elo is the internal energy at
this reference temperature. For water these constants have
the values Cv, = 4180 JkgK1, To = 273.15 K and elo = 617.0
Jkg', respectively.
Vapor Phase
For the vapor phase the ideal gas equation of state is used:
P,(p,,e)=(yl1)pe,(p,,TO ), (20)
with y the ratio of specific heats equal to y = 1.327 for water
vapor. The corresponding caloric equation of state is
ev (p, ) = ev(L) = C,( T)+L (T,)+e, (21)
where L, represents the latent heat of vaporization and C,
the specific heat at constant volume with values equal to
L,(To) = 2.3753106 Jkg1 and C,, = 1410.8 JkgK'.
Mixture Phase
In the mixture the two phases are assumed to be in thermal
and mechanical equilibrium and the pressure in the mixture
phase is taken to be equal to the saturation pressure:
S= = P.= (T). (22)
The mixture density p and mixture internal energy e are
defined by
P = Oap t (T)+(1a) p,,, (T), (23)
pe = apsa, (T)e (T)+(1 a)p,,sa, (T)e, (T) (24)
The void fraction a of the vapor is obtained from Eq. (6).
The saturation values as a function of temperature for the
pressure psat(), the liquid density PI,,a7) and vapor density
Pv,sa(T) are given by the following analytical expressions
(Schmidt 1989) which are valid for the range T = [T, Tc]
with T, the triple point and Tc the critical point:
Inps(T) = aa, (25)
,1s(t b_ (26)
Pc 1=1
In (P 'pv (T) 7 (27)
= NcjOC. (27)
0, Pc 1=1
where 0 = 1T/Tc, withpc and p, the pressure and density at
the critical point and where the coefficients al, bl, bl,b c,
and C, are included in Table 1. Table 2 gives an overview of
values of the parameters used in the equations of state for li
quid and vapor. In Figure 1 the pv diagram for water with
the isotherm at reference temperature T = 298 K is pre
sented. The employed equations of state, i.e. Tait's modified
equation of state, the mixture state and the perfect gas
equation of state are compared with IAPWS experimental
data, see (IAPWS).
Table 1: Parameters for the saturation relations Schmidt
(2006). T, = 647.16 K,p, = 221.2105 Pa, p, = 322.0 kgm3,
T,= 273.15 K.
I a a b b c, c
1 0 0 1 0 0 0
2 7.85823 1 1.99206 1/3 2.02957 2/6
3 1.83991 3/2 1.10123 2/3 2.68781 4/6
4 11.7811 3 0.512506 5/3 5.38107 8/6
5 22.6705 7/2 1.75263 16/3 17.3151 18/6
6 15.9393 4 45.4485 43/3 44.6384 37/6
7 1.77516 15/2 675615 110/3 64.3486 71/6
Table 2: Parameters for liquid and vapor phase of water.
Liquid Vapor
Parameter Value Parameter value
N 7.15 y 1.327
Ko 3.310 Pa C,, 1410.8 Jkg'K
Ci 4180 To 273.15 K
_Jkg K 1
To 273.15 K Lv(TO) 2.75310Jkg1
e ]o 617.0 Jkg_ ejo 617.0 Jkg'
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
2
10
101
10
10
10 S
SL,
2
104
104
Modified Tait
,c,^
,'''z
o IAPWS data
\, Vap
Ss
Saturated mixture
 Liquid saturation
curve
10 0 0 0 1 0
0 n. n n. 0 0 n
T= 0 t0,0 t,y t, ,T= 0 ny
0 t2,x t2,y t2,z 0 0n
O 0 0 0 1 0 0
or saturation
curve
SPerfect
102
10
v = 1/p [m3 kg11
Figure 1: pv diagram for water with isotherm at reference
temperature T,= 298 K. Presented are IAPWS experimental
data (Schmidt 1989) (open circles), saturation curves for
liquid and vapor (dashed lines), saturation points SL and SV
for liquid and vapor, respectively, Tait's modified equation
of state for liquid, mixture state and perfect gas equation of
state (solid lines). C is critical point.
Numerical Method
The homogeneous mixture equations presented in equation
(10) are solved employing an unstructured edgebased
finitevolume method. The physical domain is divided into
control volumes V, with boundary aV,. Then equation (10)
can be considered for each control volume:
a JIJ'Ud + J'J' (U) fidS=0. (28)
t V, S=av,
Defining the control volume averages U, as
U, = _JJUd (29)
with VI the volume of control volume V, and taking into
account that S = S ,with N, the total number of
faces S, of control volume V,, equation (28) can be rewritten
as
+ F (U) idS=0. (30)
The flux F(U) ii can be evaluated by exploiting the rota
tional invariance property of the Euler equations (Koop
2008), which states that
F (U)n. + F, (U)ny + F (U)n, = T'F (TU),
T is the rotation matrix and T1 its inverse given in terms of
elements of unit normal vector i = [n, ny, n ] and unit
tangential vectors t, = [tl,, t, ] t2 = ,[t, ,,t ,z ]
Note that the unit vectors n, t,, t2 form an orthogonal sys
tem. Equation (30) can now be rewritten as:
A1 N
S JT 1F (TU)dS=O. (33)
at v J=1 S,
The intercell flux F (TU) over face S, can be approxi
mated by the numerical flux H(TUL,TUR), where the
extrapolated face values UL and UR depend on the control
volume averages U, and Uj on either side of interface S,.
Before proceeding the following notation is introduced for
the numerical flux:
T'F (TU) T1H(TU,,TUR) H(U,,UR,ii ), (34)
where the unit normal vector ni of face S, defines the ro
tation matrices T and T1. Applying equation (34) to (33)
and assuming that the numerical flux H is constant over
face S,, the semidiscretized form of the finitevolume for
mulation reads:
SNf
oU, 1
+ eH(UL,UR, S1=0. (35)
0t =1
Note that when the face S, belongs to the boundary of the
computational domain, then it is necessary to determine UL
or UR from boundary conditions. Also, when U, = U and
U = U then the finitevolume scheme is firstorder ac
curate when an approximate Riemann solver is used. In this
study the extension to higher order is achieved by consi
dering piecewise linear data reconstructions in each control
volume based on the MUSCL approach of van Leer (1979).
The extrapolated face values UL and UR are then obtained
through
UL =U, + [(VU)cg,(j xg,)], (36)
UR =Uj+I [(VU)g,j ji jg, )], (37)
where (VU),, is the gradient of the variables U at the cen
troid of control volume V, and where x and are the
location of the center of the face S, and the center of gravity
of control volume V,, respectively. Furthermore, the limiter
function T e [0,1] is applied to avoid spurious oscilla
tions at the cell interfaces. Here, the limiter function of
(31) Venkatakrishnan (1995) is used.
Defining the residual %" at time level t" for control volume
1= H (U,Uni",) Sl,
I +1 jTl
then equation (35) can be written as
U + ": =0. (39)
at
To advance the solution from time level t to time level t"+ a
standard lowstorage fourstage RungeKutta timeintegra
tion method is employed, which is defined as follows:
U) = U (40)
Uk) = U ) akAt9j ,fork=1,...,4 (41)
U7+1 = U(4), (42)
with the coefficients ak equal to [0.1084,0.2602,0.5052,1.0],
resulting in secondorder accuracy in time, see Blazek
(2005). For unsteady flow calculations the time step At is a
global time step defined as the minimum of the local time
steps of all control volumes:
At = min At,, (43)
where the local time step is defined by
At= =C' (44)
max ( ii +c, ii,)
with f, a characteristic length of control volume V, defined
as the diameter of the smallest inscribed sphere of control
volume V, and where ii and c, are the local velocity vector
and speed of sound in control volume V,, respectively. The
constant CFL number C is set to 0.8.
In the present research a number of classical flux schemes
has been investigated for the numerical flux H(U, U, fii).
For cavitating flows with its phase interfaces and strong
pressure pulses, a robust, stable and accurate scheme is re
quired. Various flux schemes have been investigated, such
as the HLLC scheme (Harten, Lax &van Leer (1983), (Toro,
Spruce & Spears 1994) and (Batten, Clarke, Lambert &
Causon 1997) and the AUSM scheme and its extensions
(Liou 2006).
The stability and accuracy of the scheme is not only influ
enced by sharp gradients of the density in the flow. The
socalled lowMach number problem is also to be engaged.
A popular technique to solve this problem has been to in
troduce suitable preconditioners, e.g. (Weiss & Smith
1995). In more recent years, adaptation of the flux schemes
by scaling or modifying the numerical dissipation in regions
with lowMach number flows have been found to be suc
cessful, e.g. (Guillard & Viozat 1999), (Liou 2006) as well
as (Schmidt, Sezal & Schnerr 2006) and (Schnerr, Sezal &
Schmidt 2008). In this research the flux scheme of Schnerr,
Sezal & Schmidt (2008) is employed, which overcomes the
inaccurate calculations of the pressure field for (lowMach
number) liquid flows.
Boundary Conditions
The treatment of the boundary conditions is based on the
conventional ghostcells approach, which implies that the
numerical flux over a boundary interface is determined with
the same numerical flux as used for internal cells. When the
face S, considered is located on the boundary of the com
putational domain the right state UR for the numerical flux
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
H(U ,U ,i ) at the cellinterface can be reconstructed
from the state Ug in a virtual control volume or "ghost" cell.
The state Ug in the ghost cell is specified by applying suit
able boundary conditions to obtain the controlvolume
averaged values of the ghost cell.
During cavitation strong pressure pulses are generated.
nonreflecting boundary conditions are required that allow
waves in the flow solution to leave the computational do
main without any reflection at the in and outflow boundary
of the computational domain. This is necessary to avoid that
spurious reflected waves interfere with the timedependent
cavitating flow solution.
For the in and outflow boundaries, the timedependent,
nonreflecting boundary conditions for the hyperbolic Euler
equations with arbitrary equations of state, which are de
rived in Koop (2008), are employed. These boundary con
ditions are based on the method of Thompson (1990), who
introduced a unified formalism for the timedependent
nonreflecting boundary conditions for the Euler equations
using the perfect gas law as the equation of state.
The central concept is that a hyperbolic system of equations
represents the propagation of waves and that at any bound
ary some of the waves are propagating into the computa
tional domain while others are propagating outwards. The
number and type of conditions at a boundary of a multidi
mensional domain are defined by the eigenvalues of the Ja
cobian of the component of the flux vector in the direction
normal to the boundary. Each of the five eigenvalues A, re
presents the characteristic velocity at which a particular
wave propagates. The behavior of outgoing waves is com
pletely described by the solution within and at the boundary
of the computational domain, while the behavior of the in
coming waves is specified by data external to the computa
tional domain.
Thus, because of the wave structure of the hyperbolic
equations, the number of boundary conditions that may be
imposed depends on the physics of the problem and may not
be specified arbitrarily, i.e. the number of boundary condi
tions which must be specified at a point on the boundary is
equal to the number of incoming waves at that point. The
number of boundary conditions required at any point on the
boundary may change in time. Also, the number of bound
ary conditions required at any time may vary with position
on the boundary. For a subsonic inflow and a subsonic out
flow four and one physical boundary conditions need to be
specified, respectively.
For the solid walls on the hydrofoil and for the tunnel walls
free slip boundary conditions are employed by using the
nonpermeability condition, i.e. u fi = 0 The values in
the ghost cells at the solid walls are obtained by applying
the Curvature Corrected Symmetry Technique for unstruct
ured grids as proposed by Wang & Sun (2003).
Geometry 3D Twist Hydrofoil
Foeth (2008) has carried out experiments on a threedi
mensional twisted hydrofoil placed in the cavitation tunnel
of Delft University of Technology. The hydrofoil is denoted
as Twist 11 hydrofoil because of its varying geometric angle
of attack from 0 at the tunnel walls to 11 angle of attack at
midsection, see Dang (2001), Koop, Hoeijmakers, Schnerr
& Foeth (2006) and Foeth (2008). The chord length of the
foil is equal to c = 0.15m. The foil spans the cavitation tun
nel from wall to wall and it is symmetric with respect to its
midspan plane.
The spanwise varying distribution of the local geometric
angle of attack a (y) is designed to avoid interaction of
the cavitation sheet with the tunnel wall. The local geo
metric angle of attack a(y) is defined by a cubic poly
nomial, such that it is a,,a at the tunnel wall, aax at
midspan and that its derivative in spanwise direction is zero
at the wall as well as at midspan:
a(y)= am, (2 y13 3y2 +1)+ wa o, (45)
where = y / c is the dimensionless coordinate in span
wise direction: y e [1,1], with y = 0 defined to be at the
midplane of the span and y = 1 at the tunnel wall at the
starboard side. Note that a,,w is the rotation angle of the
entire hydrofoil, which is always equal to the local geo
metric angle of attack at the tunnel wall. The sections of the
hydrofoil rotate around x= x/c =0.75 to reduce the
optical blocking of the midsection plane by the hydrofoil
when viewing from the sides of the hydrofoil, which is illus
trated in Figure 2(c), where the side view of the hydrofoil is
presented. The hydrofoil considered is the Twistl 1 at 2 deg
geometric angle of attack, yielding a.m = 11 deg and aw,a =
2 deg in equation (45). The chord of the foil is 0.15m and
the span is 0.3m. The hydrofoil has a NACA0009 section
with its halfthickness distribution given by
0.2
with x = x c where the coefficients are equal to ao =
0.2969, a, = 0.126, a2 = 0.3516, a3 = 0.2843, a4 = 0.1015
and with the thickness parameter t = 0.09. The hydrofoil is
presented in Figure 2, where a 3D view, top view, side view
and front view are shown.
Note that the distribution of the hydrodynamic angle of
attack will be different from the geometric angle of attack.
Due to the spanwise lift distribution there is a downwash on
the central part of the hydrofoil and an upwash near the
tunnel walls. This follows from Prandtl' s lifting line theory
applied to the twisted hydrofoil, see Koop, Hoeijmakers,
Schnerr & Foeth (2006).
(a) (b)
vx
Figure 2: 3D Twistl hydrofoil at AoA = 2 deg, flow is in
xdirection. (a) 3D view (b) top view (c) side view (d) front
view.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Computational Domain and Mesh
For numerical purposes the length of the test section is in
creased to minimize the effects of the implementation of the
inlet and outlet boundary conditions. The hydrofoil is lo
cated in a channel with height 2c, a length of 3c upstream of
the leading edge, a length of 2c downstream of the trailing
edge and a width of s = 2c. Note that for the numerical flow
simulation only the starboard half of the test section and the
hydrofoil is considered, because of its geometric symmetry
and the assumed hydrodynamic symmetry.
The computational domain is divided into tetrahedral ele
ments utilizing the software package ICEMCFD. The sur
face of the hydrofoil is divided in 7 subsurfaces, i.e. one
surface wrapping around the leading edge and three sur
faces on either side of the foil. Each surface has its own size
of surface elements to ensure a fine mesh around the nose of
the foil and a somewhat coarser mesh on the surfaces closer
to the trailing edge. Note that the tetrahedrons close to the
hydrofoil are much smaller than the tetrahedrons away from
the hydrofoil.
Following a limited grid sensitivity study (Koop, Hoeij
makers, Schnerr & Foeth 2006) it was concluded that for
singlephase water flow a total number of around 350,000
tetrahedrons is adequate for a sufficiently accurate solution
on the surface of the hydrofoil, which results in approxi
mately 69,000 control volumes in the complete computa
tional domain. For cavitating flow a tetrahedral grid with a
refinement along the suction side of the hydrofoil is consi
dered resulting in 205,000 control volumes. Both meshes
are presented in Figure 3.
Results for SinglePhase Flow
The numerical method is validated utilizing the experimen
tal data of Foeth (2008) for singlephase flow of water. The
angle of attack of the hydrofoil is 2 deg. The freestream
velocity U,, temperature T. and pressure p, are equal to U,
= 6.75 ms1, T,= 297 K and p = 0.97105 Pa, respectively
At these conditions the pressures at the surface of the hydro
foil and the lift force on the hydrofoil have been measured.
Foeth (2008) reports a lift force of 455N, which amounts to
a lift coefficient of CL = 0.46.
The calculations are carried out with the firstorder spatial
reconstruction method for the first 50,000 iteration steps,
after which the calculation is continued with the MUSCL
type secondorder reconstruction method applying the lim
iter method of Venkatakrishnan (1995) on the primitive
variables [p, u, v, w, e] The calculations have been carried
out on the coarse and fine mesh presented in Figure 3.
It is found that during the firstorder reconstruction the cal
culation converged rapidly. Switching to secondorder re
sulted in the residuals decreasing by two to three orders
only. The reason for this behavior of the convergence has
been analyzed (Koop 2008), to correspond to the very low
Mach number in combination with the used limiter method
for the secondorder spatial reconstruction on unstructured
grids. To improve the convergence behavior an additional
calculation for the freestream velocity increased from 6.75
ms'1 to of U.= 50 ms1 has been carried out. For this higher
velocity also a stagnating convergence is obtained, but at a
two order lower level.
(a)
7
Figure 3: Computational mesh for 3D Twistll hydrofoil.
AoA = deg. Flow is from left to right. (a) 357,000 tetrahe
drons, 69,000 control volumes (b) 1,096,000 tetrahedrons,
205,000 control volumes
In Table 3 the measured value for the lift coefficient is
compared with the calculated values for the freestream ve
locity U = 6.75 ms1 on both grids and for U = 50 ms1 on
the coarse grid. The lift force for U. = 6.75 ms1 is calcu
lated to be equal to 451.8N. This value corresponds to
within 1% with the measured value of 455N. Furthermore,
the experimentally obtained lift coefficient CL = 0.46 is
sufficient accurately predicted by the present numerical
method. For the coarse and the fine mesh the calculated lift
coefficient are approximately equal, which verifies that the
coarse grid is adequate for predicting these integrated
quantities for single phase flow. For inviscid flow at low
Mach number the lift coefficient for the velocities U = 6.75
ms1 and U. = 50 ms1 should be about equal, which is
verified in Table 3.
For this case of threedimensional inviscid flow over a
configuration that can be considered periodic in spanwise
direction, the drag force will not be zero. The wake down
stream of the trailing edge contains vorticity stemming from
the difference in direction of the velocity on the suction side
and the flow on the pressure side at the trailing edge result
ing in a trailing vortex sheet. This trailing vorticity induces
an upwash/downwash distribution at the foil which in
creases/decreases the local angle of attack experienced by
the hydrofoil. This results in the socalled induced drag ex
perienced by hydrofoils of finite span, see also Koop, Hoeij
makers, Schnerr & Foeth (2006) and Koop (2008). From the
good agreement between the experimentally obtained lift
force and the calculated value and the agreement between
the calculated lift coefficients for both velocities, we con
clude that the present numerical method is capable of accu
rately predicting the lift force for low speed 3D single
phase flow of water.
In Figure 4 Foeth's (2008) measured distribution of the
surface C, are compared with the numerical results for a
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
number of span positions, i.e. 50% (center), 40%, 30% and
20%. The experimental values are denoted by the open
squares, the calculated values by the closed circles.
Table 3: 3D Twistll hydrofoil at 2 deg angle of attack.
Measured lift coefficient and calculated lift and drag coef
ficients CL and CD for velocities U. = 6.75 ms1 for coarse
and fine grid and U.= 50 ms1 for coarse grid. p, = 0.97105
Pa, p. = 998.3 kgm 3, T = 297 K.
U. Numerical Experimental
[msi] mesh CD I CL CL [1
6.75 coarse 0.010 0.442 0.46
6.75 fine 0.0083 0.445 0.46
50 coarse 0.0098 0.454
Figure 4 shows that the numerical results correspond
reasonably well with the experimental data. However, the
experimental value at x/c = 0.3 at 40% span appears to devi
ate from the numerical results. According to Foeth (2008)
the value from this sensor may not be reliable. Furthermore,
an overshoot in p value at the trailing edge is observed in
the numerical results. This is caused by the relatively irregu
lar shaped control volume around the trailing edge, typical
for a nodecentered dual mesh. A solution has been found,
see Hospers (2007) by splitting the control volumes around
the trailing edge into an upper and lower control volume
allowing a discontinuous solution at the trailing edge.
In Figure 5 the distribution of the surface Cp coefficient is
presented for singlephase flow of water for U. = 50 ms1.
In the leading edge region on the suction side of the hydro
foil a clear low pressure region is visible. The design of the
hydrofoil has been such that cavitation will occur in the
midspan region of the hydrofoil and that cavitation is a
voided near the tunnel walls
(a) y/s = 0.5
e exp
S num
1 4
x/c2 []
(c) y/s = 0.3
E exp
Snum
x [.. ]
x/c[]
(b) y/s = 0.4
S exp
num
.. . ,
xe [f
(d) y/s = 0.2
E exp
* num
x/c []
Figure 4: 3D Twistl hydrofoil at AoA = 2 deg. Experi
mental (open squares) and numerical (solid circles) distrib
ution of Cp coefficient for U = 50 ms1, T. =297K, p=
0.97. 105 Pa, p = 998.2 kgm3 (a) y/s = 0.5, (b)y/s = 0.4, (c)
y/s = 0.3, (d) y/s =0.2.
Indeed, due to the so chosen spanwise varying geometric
angle of attack, in combination with the down/upwash in
duced by the trailing vortex wake, the pressure near the
walls of the cavitation tunnel is much higher resulting in a
lower C value.
C, []
28
26
24
22
18
16
14
12
08
06
04
Z 02
0
06
08
Figure 5: 3D Twistll hydrofoil at AoA= 2 deg. Distribu
tion of Cp coefficient on surface of hydrofoil. U. = 50ms 1,
T = 297K, p = 0.97 105 Pa, p = 998.2 kgm3.
Results for Twophase Flow
The unsteady cavitating flow about the 3D Twist 11 hydro
foil at 2 deg angle of attack is calculated and compared
with the results of the experiments of Foeth (2008).
Considered is the situation for steady inflow at a cavitation
number equal to a = 1.1. The geometry is symmetric with
respect to the mid span plane of the hydrofoil as it is
installed in the cavitation tunnel. Assuming the flow to be
symmetric with respect to the same plane, only the
starboard half of the configuration needs to be considered,
which saves computational time. This assumption of
hydrodynamic symmetry is supported by the experimental
findings of Foeth. To speed up the formation and periodic
shedding of the cavity sheet in the numerical simulations,
the calculations are performed at a freestream velocity of
U~ = 50 ms'. Note that in order to obtain a cavitation
number equal to = 1.1 at a temperature of T. = 297K, a
freestream pressure equal to p. = 13.75105 Pa is chosen.
To accelerate the calculations the numerical method has
been made suited for parallel computation by decomposing
the computational mesh in 16 equalsized blocks.
First, a solution employing the firstorder reconstruction
method, is calculated on the coarse mesh. The sheet cavity
is found to become steady for this mesh and order of re
construction. The calculation is continued with the second
order reconstruction method on the coarse mesh. Then, it is
found that resolution of the cavity sheet and its shedding is
not represented adequately. The coarse mesh is too coarse to
resolve the reentrant jet and the shed vapor structures
properly. The reentrant flow is captured within one com
putational control volume and the sheet cavity occasionally
sheds vortical vapor regions, which quickly dissipate.
To improve the resolution in the region with cavitation, the
finer mesh as presented in Figure 3(b) is employed which
has a refined region along the suction side of the hydrofoil
to approximately 10% chord length in normal direction.
This ensures that the sheet cavity is located in this refined
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
region. The solution obtained on the coarse grid is used as
the initial solution for the calculations on the fine grid by
employing a solutioninterpolation method.
On the fine grid the numerical time step Atfi is equal to Atfl
= 1.01108 s. This small time step is necessary to meet the
CFL condition and to be able to capture the pressure waves
generated by the collapse of vapor structures. In water these
waves travel at a velocity exceeding 1500 ms1.
 I
Figure 6: Dimensionless total vapor volume Vvap, Equation
(7), for cavitating flow about 3D Twist 11 hydrofoil at AoA
= 2 deg, ( = 1.1. Fine grid, Atcf=1.0110ls.
The transient evolution of the cavitating flow can be illus
trated through the time history of the total vapor volume
Vvap as defined in equation (7). This dimensionless quantity
as a function of dimensionless time tU~c is presented in
Figure 6. This figure clearly indicates that after a startup
time, the evolution of the total vapor volume becomes peri
odic in time: the amplitude and period of the cycles become
constant in time. The cycle illustrated by the dashed box in
Figure 6 has a dimensionless period of TU/c = 1.7, which
corresponds to a Strouhal number based on the chord length
equal to St = 0.59.
In Figure 7 the cavitation flow structures on the Twist 11
hydrofoil are presented at the time instants 16 indicated in
Figure 6. The time AT between each picture is equal to
8.5104s or AT UJc = 0.28. Two isocontours of the void
fraction are chosen to visualize the structure of the cavity
sheet as can be seen most clearly in Figure 7(6). In Figure
7(1) the total vapor volume is at its minimum. In the
midspan region of the hydrofoil the sheet cavity has dis
appeared and only a small spanwise cavitating vortex is
visible. Also visible are two side lobes of the sheet cavity at
one quarter and three quarters of the hydrofoil in spanwise
direction, respectively. From Picture (1) to just after Picture
(2) the total vapor volume increases, mainly due to an in
crease in volume of the detached vapor regions near
midspan. This detached vapor region, shaped as a horse
shoe, decreases in volume between Picture (2) and (3) until
it disappears just after Picture (3). Meanwhile, the sheet
cavity starts to grow again from the leading edge at mid
span from Picture (2) to Picture (6), which explains the in
crease of total vapor volume up to Picture (5).
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
(1) Top view
Figure 7: Twistll hydrofoil at
Formation of cavitating flow str
Presented are the isocontours of
103 and a = 0.5. Pictures (1)(6)
6 in the time history given in Fig'
This increase is counteracted due
entrant jet between Picture (4) an
region of the cavity sheet and e
the leading edge of the hydrofoil
the cycle is repeated again.
In his thesis Foeth (2008) visu
structure of the shedding of the
especially the formation of a cav
the shape of the sheet cavity with
formation of re and sideentrani
mulation the sidelobs of the sl
Figure 7(4) just after the collap
shoe vortex. It should be noted
experiments are slightly different
1.13. However, the effect of the
attack is counteracted by the sli
number and thus, the two conditi
similar. Foeth mentions that th(
sheet cavity did not change muc
tions of the numerical simulations
1.1.
In Figure 8 the streamlines on th
are presented at the timeinstan
streamlines are colored by the
(1) Side view whether it is a liquid flow or a flow of vapor. The vectors in
dicate the direction of the flow.
Sg In Figure 7(1) the vapor sheet in the center of the foil is at
L the leading edge of the foil and starts to grow. At the same
time a flow of liquid is moving upstream, denoted by the
(2) blue color presented in Figure 8(1), and has reached the
leading edge of the foil. Note the sharp transition from li
4 quid to vapor at the midspan of the hydrofoil in Figure
L 8(1). There, the flow moves upward and impinges on the
~ interface of the cavity sheet above the reentrant jet. The
fluid impinging on the interface cuts off a region of vapor as
(3) visible in Figure 7(1). Around this vapor region a circula
rP tory flow pattern is observed creating a spanwise cavitating
vortex at midspan of the hydrofoil, which is then convected
Lx with the flow as presented in Figure 8(2)(5).
The shed cavitating vortex has disappeared between Pic
(4) tures (3) and (4) and the sheet cavity grows to its longest ex
tent at midspan, see Figure 7(3)(6). At the closure line of
z S the sheet cavity a region with high vorticity has developed,
Sd which will drive the reentrant flow. As visible in Figure
8(4) the reentrant flow at the closure line of the cavity starts
(5) as a sideentrant jet. Note that during the growth of the
sheet, a reentrant flow is already moving upstream under
neath the vapor sheet; see Figure 8(3)(6), which is confirm
x ed by the experimental results of Foeth (2008). In the pre
sent numerical simulations it appears that although an up
(6) stream moving liquid flow is already present just behind the
cavity sheet, it does not have enough momentum to reach
the leading edge until the vapor cloud has collapsed just
Lz after timeinstant 3, see Figure 7(3) and Figure 8(3). More
Research, both experimentally as numerically, should be
carried out to verify this observation and to investigate the
2 deg AoA for a = 1.1. influence of the pressure pulses originating from the
uctures during one cycle, collapse of the vapor structures on the shedding mechanism
void fraction equal to a = and the formation of the reentrant jet.
correspond with points 1 We remark that in Figure 8(2) the reentrant flow appears to
ure 6. be directed away from the plane of symmetry resulting in a
flow pattern on the surface of the hydrofoil which is in very
e to the formation of a re close agreement with that observed in the experiments of
d (5) starting at the closure Foeth (2008).
extending forward towards Avery distinct feature of the shedding of the sheet cavity on
in Picture (6) after which the 3D Twist 11 hydrofoil is the formation of a cavitating
horseshoe vortex around midspan. In Figure 7(6)(3) the
lalized and explained the formation and convection of such a cavitating horseshoe
sheet cavity in detail and vortex is clearly present in the numerical simulation. In
stating horseshoe vortex, Figure 7(6) the sheet cavity has reached its maximum extent
distinct sidelobes and the in the center of the hydrofoil. Already a reentrant flow is
Sets. In the numerical si moving upstream along the surface of the hydrofoil gene
ieet are clearly visible in rating a circulation at the closure region of the cavity sheet.
se of the cavitating horse As soon as the reentrant flow reaches the leading edge, a
that the conditions in the spanwise cavitating vortex is detached from the sheet cavity
Si.e. AoA = 1 deg and a = as visible in Figure 7(1). This cavitating vortex is convected
slight increase in angle of with the flow as visible in Figure 7(2) and (3). The center of
ght increase in cavitation the vortex is convected upward, primarily by its selfin
ons are expected to be very duced velocity, giving the vortex its horseshoe shape. The
Overall shedding of the height of the horseshoe vortex is clearly visible in the side
.h compared to the condi view in Figure 7(3). It reaches up to 23 times the thickness
, i.e. AoA = 2 deg and C = of the sheet cavity. In the experiments of Foeth the height of
the shed vapor cloud is found to be an important feature.
ie surface of the hydrofoil Between Pictures (3) and (4) the shed vapor regions col
its 16 of Figure 6. The lapse near the surface of the hydrofoil creating strong pres
void fraction to indicate sure pulses that propagate over the surface of the hydrofoil.
Three timeinstants, tit3, are selected between Pictures (3)
and (4). Picture (3) corresponds with to, then t, to t3 are
equal to t =to + 1.8104 s, t2 = t+ 105 s, and t3= t2+105 s,
respectively. In Figure 9 (I)(III) the predicted pressure and
two isocontours of the void fraction are presented for the
times tit3, illustrating the origin and propagation of the
a[]
07
o5
04
02
01
a []
05
a 
09
08
07
02
01
a []
14
09
0 8
7
a []
09
S4
SO3
O2
05
08
S3
02
0 1
Figure 8: Streamlines on the surface of the 3D Twistll hy
drofoil at AoA= 2 deg, a = 1.1. The streamlines are colored
by the value of the void fraction. Pictures (1)(6) corres
pond with timeinstants 16 presented in Figure 6.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
pressure pulses above and on the surface of the hydrofoil.
In Figure 9(I) the convected vapor region reaches a region
with higher pressure. Inside the vapor region the pressure is
equal to the saturation pressure. This pressure difference
induces a local flow field directed towards the center of the
vapor region and causes the vapor region to collapse. The
inward moving liquid impacts at the center of the former
vapor region and initiates an outward propagating spherical
shock wave traveling faster than the speed of sound c in the
medium as presented in Figure 9(11). The maximum pres
sure in Figure 9(11) is 103 bar, which is 7.5 times the free
stream pressure p,.
In Figure 9(111) the shock wave is travelling radially out
ward from the region of collapse. Downstream the shock
wave runs over the surface of the hydrofoil inducing a high
loading on the surface of the hydrofoil. Upstream the shock
wave hits the vapor sheet. Due to the lower acoustic impe
dance pc of the twophase flow region compared to that of
the liquid, the shock wave is reflected as an expansion wave
from the interface of the cavity sheet along which the pres
sure is nearly constant. For the same reason the shock wave
also reflects from the remaining vapor structures above the
region of collapse.
As visible in Figure 9(111) the pressure again returns close to
the saturation pressure at the location of the former vapor
region, which causes the liquid to cavitate again. In the
present numerical simulations for threedimensional flow
the resolution of the computational mesh is not fine enough
to capture this behavior in more detail. However, in Koop
(2008)0 we have shown that for twodimensional cases this
phenomenon is the reason for the socalled rebound of
vapor clouds, which occasionally is observed during exper
iments.
From Figure 9 we conclude that with the present mathe
matical model and numerical formulation it is possible to
capture pressure pulses in the liquid flow, which are gene
rated by the collapse of vapor structures. These pressure
pulses are believed to be the origin of erosion of surface ma
terial due to cavitation. Further calculations on a finer mesh
should be carried out to focus more on the collapse phase of
the shed vapor structures.
Note that between Figure 9(11) and (III) 1000 numerical
time steps have been taken. It might be worthwhile to inves
tigate methods such as multigrid, preconditioning and/or
implicit timeintegration, which allow larger numerical time
steps to be taken. However, the larger admissible numerical
time steps should still resolve the propagation of highfre
quency pressure pulses, which we believe to have a major
influence on the selfoscillatory behavior of the sheet cavity
and its shedding.
Conclusion and Discussion
In this paper the numerical results for the cavitating flow
about the 3D Twist 11 hydrofoil at 2 deg angle of attack are
presented employing an equilibrium cavitation model in
combination with the Euler equations. The used mathemat
ical model for a compressible homogeneous watervapor
mixture at equilibrium saturation conditions offers a general
applicable model for cavitation and does not have any free
empirical parameters for phase transition.
It is shown that the numerical method is able to accurately
calculate single phase water flow about the threedimen
(I) Top vi
(II)
(III)
(I) Side vi
ewv
iew
ew
(III)
Figure 9: Top and side view of collapse of shed vapor
structures on 3D Twistl 1 hydrofoil at 2 deg AoA. Present
ed are two isocontours of the void fraction, i.e. a = 103 and
a = 0.5 and the solution for the pressure p. Denoting time to
with timeinstant 3 in Figure 6, Pictures (I)(III) correspond
with times t, to t3, where t = to + 1.8104 s, t2 = + 10 s, and
t3 = t2 + 105 s, respectively.
sional Twist 11 hydrofoil by comparison of predicted and
measured surface pressure data.
For cavitating flow we have found that the shape of the
sheet cavity and the outline of the closure region as pre
dicted by the present numerical method compare quite well
with the experimental results of Foeth (2008).
The 3D Twist 11 hydrofoil has been designed to have a clear
and controllable threedimensional sheet cavity. It has been
shown both experimentally by Foeth as numerically in this
F"T9
L
p [bar]
p [bar]
I
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
paper that the shape of the cavity and the closure line of the
cavity determine the direction of the reentrant flow and that
the reentrant flow from the sides dictates the behavior of
the shedding cycle. Therefore, the shedding of a sheet ca
vity is governed by the direction and momentum of the re
entrant and sideentrant jets and their impingement on the
cavity surface. The impingement of the reentrant flow on
the cavity interface and the detachment of a spanwise ca
vitating vortex are captured in the present numerical simula
tions based on the Euler equations. This confirms that these
phenomena appear to be inertial in nature.
The development of a reentrant flow is predicted in close a
greement with that seen in the experiments. During the
growth of the sheet, a reentrant jet is already moving up
stream underneath the vapor sheet. Note that the predicted
reentrant flow is a flow of liquid.
The formation of a cavitating horseshoe vortex and its ad
vection with the flow is captured in the present numerical si
mulations and they agree quite well with the experimental
observations.
We have shown that the present numerical method is capa
ble of predicting the collapse of shed vapor structures and
the subsequent high pressure pulses on the surface of the
hydrofoil, which is important for the prediction of erosion
and noise. In experiments it is a difficult task to capture
and/or visualize these pressure pulses and the associated
unsteady loading of the hydrofoil. In order to further vali
date the numerical method it is important to gain more
knowledge experimentally and numerically on the strength
of the pulses in combination with unsteady sheet cavitation.
Furthermore, the influence of the pressure pulses on the
shedding mechanism should be investigated both numeri
cally as experimentally, especially the influence on the for
mation of the reentrant jet.
At present the grid resolution is insufficient to be able to
capture the phenomena related to smaller scale vapor struc
tures. In the future methods need to be investigated to speed
up the present numerical method. It might be worthwhile to
investigate methods such as multigrid, preconditioning
and/or implicit timeintegration, which allow larger numer
ical time steps to be taken. However, the larger admissible
numerical time steps should still resolve the highfrequency
pressure pulses, which we believe to have a major influence
on the selfoscillatory behavior of the sheet cavity and its
shedding.
The present paper focused on sheet cavitation on a statio
nary hydrofoil located in uniform inflow conditions. Foeth
(2008) conducted experiments with a configuration con
sisting of the present hydrofoil placed behind two other,
stacked, hydrofoils with oscillating flaps generating an un
steady inflow. Huckriede (2008) and Huckriede, Koop,
Hospers & Hoeijmakers (2010) modeled these oscillating
hydrofoils employing socalled 'transpiration' boundary
conditions. The present research will be extended to nume
rically investigate the influence of unsteady inflow on the
shedding of the cavitation sheet and compare the results
with those of the experiments of Foeth (2008).
In Koop (2008) we also conducted numerical simulations
for the steady cavitating flow about a finite three dimen
sional wing generating a cavitating tip vortex. Future work
has to be carried out to investigate the interaction between
the sheet cavity and the cavitating vortex. One of the diffi
culties lies in accurately predicting the internal structure of
the cavitating vortex. Ton (2008) and Ton, Koop, Hoeij
makers & de Vries (2009) has investigated the socalled
'Vorticity Confinement' method. There it is shown that it is
indeed possible to improve the numerical results for tip
vortices. However, it is also shown that the vorticity con
finement method is not yet robust and needs to be explored
further.
Acknowledgements
This research has been funded by the Dutch Technology
Foundation STW project TSF.6170. See http:/www.stw.nl
for more details. The authors would like to thank Prof. G
Schnerr and S. Schmidt from the TU Munich, Germany for
their support during this research. Furthermore, Prof. T. van
Terwisga and E.J. Foeth from the TU Delft, the Netherlands
are acknowledged for their collaboration by conducting the
experimental part of this STW project.
References
Batten, P, Clarke, N. Lambert, C. & Causon, D.M. On the
Choice of Wavespeeds for the HLLC Riemann Solver.
SIAM J. Sci. Comput. 18:15531570 (1997)
Blazek, J. Computational Fluid Dynamics: Principles and
Applications. Elsevier, 2nd edition (2005)
Dang, J. Numerical Simulation of Unsteady Partial Cavity
Flows. PhD thesis, Delft University of Technology (2001)
Foeth, E.J. The structure of ThreeDimensional Sheet Ca
vitation, PhD Thesis, Delft University of Technology
(2008)
Guillard, H. & Viozat, C. On the Behavior of Upwind
Schemes in the Low Mach Number Limit. Computers &
Fluids, 28:6386 (1999)
Harten, A., Lax, PD. & van Leer, B. On Upstream Differ
encing and Godunovtype Schemes for Hyperbolic Conser
vation Laws. Siam Review, 25:3561 (1983)
Hospers, J.M. Development and Application of Hybrid
Grids for Finite Volume Methods. MSc Thesis, University
of Twente (2007)
Huckriede, K.W. Transpiration Boundary Condition for
Numerical Study of Unsteady 2D Sheet Cavitation on
Hydrofoils. MSc Thesis, University of Twente (2008)
Huckriede, K.W, Koop, A.H., Hospers, J.M. & Hoeij
makers, H.W.M. Finite Volume Method with Transpiration
Boundary Conditions for Flow about Oscillating Wings.
AIAA Paper 2010866 (2010)
IAPWS, The International Association for the Properties of
Water and Steam; http:/www.iapws.org/.
Jameson, A. Schmidt, W. & Turkel, E. Numerical Solution
of the Euler Equations by Finite Volume Methods Using
RungeKutta TimeStepping Schemes. AIAA Paper, 81
1259 (1981)
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Koop, A.H., Hoeijmakers, H.WM. Schnerr, G.H. & Foeth,
E.J. Design of Twisted Cavitating Hydrofoil using a Baro
tropic Flow Method. CAV2006, Wageningen (2006)
Koop, A.H. Numerical Simulation of Unsteady ThreeDi
mensional Sheet Cavitation. PhD Thesis, University of
Twente (2008)
Leer, van B. Towards the Ultimate Conservative Difference
Scheme V. A Secondorder Sequel to Godunov's Method.
Journal of Computational Physics, 32:101136 (1979)
Liou, M.S. A Sequel to AUSM, Part II: AUSM+up for all
Speeds. Journal of Computational Physics, 214:137170
(2006).
Saurel, R. Cocchi, J.P & Butler, J.B. A Numerical Study of
Cavitation in the Wake of a Hypervelocity Underwater
Profile. Journal of Propulsion and Power. 15(4):513522
(1999)
Schmidt, S.J., Sezal, I.H. & Schnerr, G.H. Compressible
Simulation of HighSpeed Hydrodynamics with Phase
Change. ECCOMAS CFD (2006)
Schmidt, E. Properties of Water and Steam in SIUnits;
0800"C, 01000 bar. SpringerVerlag. R. Oldenbourg, 4th
enlarged printing edition (1989)
Schnerr, G.H. Sezal, I.H. & Schmidt, S.J. Numerical Inves
tigation of ThreeDimensional Cloud Cavitation with Spe
cial Emphasis on Collapse Induced Shock Dynamics.
Physics of Fluids, 20:040703 (2008)
Thompson, K.W. Time Dependent Boundary Conditions for
Hyperbolic Systems II. Journal of Computational Physics,
89:439461 (1990)
Ton, T.A. Vorticity Confinement in Compressible Flow,
Implementation and Validation. MSc Thesis, University of
Twente (2008)
Ton, T.A., Koop, A.H., Hoeijmakers, H.W.M. & de Vries,
H. Investigation of Vorticity Confinement in Compressible
Flow. AIAA 092407 (2009)
Toro, E.F. Spruce, M. & Speares, W. Restoration of the
Contact Surface in the HLLRiemann Solver. Shock Waves,
4:2534 (1994)
Venkatakrishnan, V Convergence to SteadyState Solutions
of the Euler Equations on Unstructured Grids with Limiters.
Journal of Computational Physics, 118:120130 (1995)
Wang, Z.J. & Sun, Y CurvatureBased Wall Boundary
Condition for the Euler Equations on Unstructured Grids.
AIAA Journal, 41:2733 (2003)
Weiss, J.M. & Smith, W.A. Preconditioning Applied to
Variable and Constant Density Flows. AIAA Journal,
33(11):20502057 (1995)
