Group Title: 7th International Conference on Multiphase Flow - ICMF 2010 Proceedings
Title: 3.1.2 - Unsteady Sheet Cavitation on Three-dimensional Hydrofoil
ALL VOLUMES CITATION THUMBNAILS PAGE IMAGE ZOOMABLE
Full Citation
STANDARD VIEW MARC VIEW
Permanent Link: http://ufdc.ufl.edu/UF00102023/00072
 Material Information
Title: 3.1.2 - Unsteady Sheet Cavitation on Three-dimensional Hydrofoil Cavitation
Series Title: 7th International Conference on Multiphase Flow - ICMF 2010 Proceedings
Physical Description: Conference Papers
Creator: Koop, A.H.
Hoeijmakers, H.W.M.
Publisher: International Conference on Multiphase Flow (ICMF)
Publication Date: June 4, 2010
 Subjects
Subject: cavitation
numerical simulation
phase transition
 Notes
Abstract: An unstructured grid computational method for solving the Euler equations for inviscid compressible unsteady flow is employed to predict the structure and dynamics of three-dimensional 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 thermodynamic and mechanical equilibrium in the two-phase 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 two-phase mixture and the vapor phase of the fluid. The three-dimensional unsteady cavitating flow about model configuration (3D Twist11 hydrofoil) is calculated for which experimental data is available (Foeth 2008). It is shown that the formation of re-entrant flow and of a cavitating horse-shoe-vortex is captured by the present numerical method. The calculated results agree reasonable well with experimental observations. Furthermore, the collapse of the shed vapor structures and the resulting high pressure pulses is demonstrated.
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: VID00072
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: 312-Koop-ICMF2010.pdf

Full Text

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


Unsteady Sheet Cavitation on Three-dimensional 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 three-dimensional 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 two-phase 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 two-phase mixture and the vapor phase of the fluid. The
three-dimensional 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 re-entrant flow and of a cavitating horse-shoe-vortex 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.


Re-entrant 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 so-called re-entrant and side-entrant 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 re-entrant and side-entrant
jets.
\1,./.l,,,- and collapse of vapor structures. The
break-up 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
Navier-Stokes equations have been developed to numeri-
cally simulate cavitating flows. The main problem in the nu-
merical simulation of multi-dimensional 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 low-velocity 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 three-dimensional
cavitating flows require an appropriate high-resolution
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 re-entrant and side-entrant 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 re-entrant 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. [ms-1] are the pressure,
density and velocity, respectively, of the free-stream. 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) [kgm-3] and i,,saT) [kgm-3] 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 30-June 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 two-phase 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 two-phase
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+ i-iu (12)
P 2
with the total specific energy E defined by
I -
E=e+-u-u, (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 30-June 4, 2010


which yields
pu
puu + pnx
F(U)-fi= puv+ pny (15)
putw + pnz
puH
where U is the contra-variant 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.3-108 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 Jkg-K-1, 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)=(y-l1)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.3753-106 Jkg-1 and C,, = 1410.8 Jkg-K-'.


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)+(1-a) 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 = 1-T/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 p-v 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.2-105 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.753-10Jkg1
e ]o 617.0 Jkg_ ejo 617.0 Jkg'





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


2
10

101

10
10


10- S
SL,
-2
104


10-4


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


Figure 1: p-v 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 edge-based
finite-volume 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 T-1 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 inter-cell 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) T-1H(TU,,TUR) H(U,,UR,ii ), (34)

where the unit normal vector ni of face S, defines the ro-
tation matrices T and T-1. Applying equation (34) to (33)
and assuming that the numerical flux H is constant over
face S,, the semi-discretized form of the finite-volume 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 finite-volume scheme is first-order ac-
curate when an approximate Riemann solver is used. In this
study the extension to higher order is achieved by consi-
dering piece-wise 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 low-storage four-stage Runge-Kutta time-integra-
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 second-order 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
so-called low-Mach 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 low-Mach 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 (low-Mach
number) liquid flows.


Boundary Conditions

The treatment of the boundary conditions is based on the
conventional ghost-cells 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 30-June 4, 2010

H(U ,U ,i ) at the cell-interface 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 control-volume
averaged values of the ghost cell.
During cavitation strong pressure pulses are generated.
non-reflecting 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 time-dependent
cavitating flow solution.
For the in- and outflow boundaries, the time-dependent,
non-reflecting 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 time-dependent
non-reflecting 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 multi-di-
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
non-permeability 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 three-di-
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
mid-section, 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
mid-span 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
mid-span and that its derivative in spanwise direction is zero
at the wall as well as at mid-span:

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
mid-plane 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 mid-section 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 half-thickness 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
x-direction. (a) 3D view (b) top view (c) side view (d) front
view.


7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30-June 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 ICEM-CFD. The sur-
face of the hydrofoil is divided in 7 sub-surfaces, 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
single-phase 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 Single-Phase Flow

The numerical method is validated utilizing the experimen-
tal data of Foeth (2008) for single-phase flow of water. The
angle of attack of the hydrofoil is -2 deg. The free-stream
velocity U,, temperature T. and pressure p, are equal to U,
= 6.75 ms-1, T,= 297 K and p = 0.97-105 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 first-order spatial
reconstruction method for the first 50,000 iteration steps,
after which the calculation is continued with the MUSCL
type second-order 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 first-order reconstruction the cal-
culation converged rapidly. Switching to second-order 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 second-order spatial reconstruction on unstructured
grids. To improve the convergence behavior an additional
calculation for the free-stream velocity increased from 6.75
ms'1 to of U.= 50 ms-1 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 free-stream ve-
locity U = 6.75 ms-1 on both grids and for U = 50 ms-1 on
the coarse grid. The lift force for U. = 6.75 ms-1 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
ms-1 and U. = 50 ms-1 should be about equal, which is
verified in Table 3.
For this case of three-dimensional 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 so-called 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 30-June 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 ms-1 for coarse
and fine grid and U.= 50 ms-1 for coarse grid. p, = 0.97-105
Pa, p. = 998.3 kgm 3, T = 297 K.
U. Numerical Experimental
[ms-i] 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 node-centered 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 single-phase 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
mid-span 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 ms-1, T. =297K, p=
0.97. 105 Pa, p = 998.2 kgm-3 (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 span-wise 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 Two-phase 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 free-stream 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
free-stream pressure equal to p. = 13.75-105 Pa is chosen.
To accelerate the calculations the numerical method has
been made suited for parallel computation by decomposing
the computational mesh in 16 equal-sized blocks.
First, a solution employing the first-order 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 re-entrant jet and the shed vapor structures
properly. The re-entrant 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 30-June 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 solution-interpolation method.
On the fine grid the numerical time step Atfi is equal to Atfl
= 1.01-10-8 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.01-10l-s.

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 start-up
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 1-6 indicated in
Figure 6. The time AT between each picture is equal to
8.5-10-4s or AT UJc = 0.28. Two iso-contours 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
mid-span 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 span-wise
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
mid-span. 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 30-June 4, 2010


(1) Top view


Figure 7: Twistll hydrofoil at
Formation of cavitating flow str
Presented are the iso-contours of
10-3 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 side-entrani
mulation the side-lobs 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 time-instan
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 mid-span of the hydrofoil in Figure
L 8(1). There, the flow moves upward and impinges on the
~ interface of the cavity sheet above the re-entrant 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 mid-span 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 mid-span, 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 re-entrant flow. As visible in Figure
8(4) the re-entrant flow at the closure line of the cavity starts
(5) as a side-entrant jet. Note that during the growth of the
sheet, a re-entrant 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 time-instant 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 re-entrant jet.
correspond with points 1 We remark that in Figure 8(2) the re-entrant 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
horse-shoe vortex around mid-span. In Figure 7(6)-(3) the
lalized and explained the formation and convection of such a cavitating horse-shoe
sheet cavity in detail and vortex is clearly present in the numerical simulation. In
stating horse-shoe vortex, Figure 7(6) the sheet cavity has reached its maximum extent
distinct side-lobes and the in the center of the hydrofoil. Already a re-entrant 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 re-entrant 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 self-in-
ons are expected to be very duced velocity, giving the vortex its horse-shoe shape. The
Overall shedding of the height of the horse-shoe vortex is clearly visible in the side
.h compared to the condi- view in Figure 7(3). It reaches up to 2-3 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 1-6 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 time-instants, ti-t3, are selected between Pictures (3)










and (4). Picture (3) corresponds with to, then t, to t3 are
equal to t =to + 1.8-10-4 s, t2 = t+ 105 s, and t3= t2+105 s,
respectively. In Figure 9 (I)-(III) the predicted pressure and
two iso-contours of the void fraction are presented for the
times ti-t3, 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 time-instants 1-6 presented in Figure 6.


7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30-June 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 two-phase 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 three-dimensional 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 two-dimensional cases this
phenomenon is the reason for the so-called 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 multi-grid, preconditioning and/or
implicit time-integration, which allow larger numerical time
steps to be taken. However, the larger admissible numerical
time steps should still resolve the propagation of high-fre-
quency pressure pulses, which we believe to have a major
influence on the self-oscillatory 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 water-vapor
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 three-dimen








(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 iso-contours of the void fraction, i.e. a = 10-3 and
a = 0.5 and the solution for the pressure p. Denoting time to
with time-instant 3 in Figure 6, Pictures (I)-(III) correspond
with times t, to t3, where t = to + 1.810-4 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 three-dimensional 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 30-June 4, 2010

paper that the shape of the cavity and the closure line of the
cavity determine the direction of the re-entrant flow and that
the re-entrant 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 side-entrant jets and their impingement on the
cavity surface. The impingement of the re-entrant flow on
the cavity interface and the detachment of a span-wise 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 re-entrant flow is predicted in close a-
greement with that seen in the experiments. During the
growth of the sheet, a re-entrant jet is already moving up-
stream underneath the vapor sheet. Note that the predicted
re-entrant flow is a flow of liquid.
The formation of a cavitating horse-shoe 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 re-entrant 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 multi-grid, preconditioning
and/or implicit time-integration, which allow larger numer-
ical time steps to be taken. However, the larger admissible
numerical time steps should still resolve the high-frequency
pressure pulses, which we believe to have a major influence
on the self-oscillatory 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 so-called '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 so-called
'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:1553-1570 (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 Three-Dimensional 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:63-86 (1999)

Harten, A., Lax, PD. & van Leer, B. On Upstream Differ-
encing and Godunov-type Schemes for Hyperbolic Conser-
vation Laws. Siam Review, 25:35-61 (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 2010-866 (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
Runge-Kutta Time-Stepping Schemes. AIAA Paper, 81-
1259 (1981)


7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30-June 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 Three-Di-
mensional Sheet Cavitation. PhD Thesis, University of
Twente (2008)

Leer, van B. Towards the Ultimate Conservative Difference
Scheme V. A Second-order Sequel to Godunov's Method.
Journal of Computational Physics, 32:101-136 (1979)

Liou, M.S. A Sequel to AUSM, Part II: AUSM+-up for all
Speeds. Journal of Computational Physics, 214:137-170
(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):513-522
(1999)

Schmidt, S.J., Sezal, I.H. & Schnerr, G.H. Compressible
Simulation of High-Speed Hydrodynamics with Phase
Change. ECCOMAS CFD (2006)

Schmidt, E. Properties of Water and Steam in SI-Units;
0-800"C, 0-1000 bar. Springer-Verlag. R. Oldenbourg, 4th
enlarged printing edition (1989)

Schnerr, G.H. Sezal, I.H. & Schmidt, S.J. Numerical Inves-
tigation of Three-Dimensional 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:439-461 (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 09-2407 (2009)

Toro, E.F. Spruce, M. & Speares, W. Restoration of the
Contact Surface in the HLL-Riemann Solver. Shock Waves,
4:25-34 (1994)

Venkatakrishnan, V Convergence to Steady-State Solutions
of the Euler Equations on Unstructured Grids with Limiters.
Journal of Computational Physics, 118:120-130 (1995)

Wang, Z.J. & Sun, Y Curvature-Based Wall Boundary
Condition for the Euler Equations on Unstructured Grids.
AIAA Journal, 41:27-33 (2003)

Weiss, J.M. & Smith, W.A. Preconditioning Applied to
Variable and Constant Density Flows. AIAA Journal,
33(11):2050-2057 (1995)




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