Numerical Simulation of Unsteady Turbulent Cavitating Flows
Numerical Simulation of Unsteady Turbulent Cavitating Flows
Numerical Simulation of Unsteady Turbulent Cavitating Flows
Abstract
Computational Fluid Dynamics (CFD) is a very important tool for the study of complex fluid
flows and the design of hydraulic fluid flow machinery. At the same time, experimental
analysis is very difficult to perform. Thus, for a better understanding of the behaviour of such
complex flows, including turbulence, unsteadiness and cavitation, a suitable knowledge of
CFD is indispensable. Generally, the specific applications of CFD codes for solving this type
of engineering problems are not well documented and a previous work for the acquirement of
the CFD code capabilities is necessary.
This work presents numerical investigation concerning complex unsteady flows, including
turbulence and cavitation. The main objective was to acquire deeper knowledge about the
software potentials for solving this kind of flows. To reach this objective several cases have
been studied with a commercial CFD code (Fluent v6.1). The turbulence models being used
were mainly the Spalart-Allmaras, the Standard k-ε, and the k-ω model [1], [2], [3], [4]. The
utilized cavitation model was from Singhal, 2002 [5]. Based on long term considerations, this
investigation aims at the application of the acquired knowledge and experience for further
investigations relative to the cavitation phenomena in real fluid flow machines.
Several steps were necessary to understand the suitable simulation process of unsteady
turbulent cavitating flows. The case of an unsteady and turbulent, non-cavitating flow around
a 2D circular cylinder was studied as a first step using different turbulence models at
Reynolds numbers around the critical drag-crisis region. Compared with experimental data,
the results are quite divergent, but similar numerical researches (J.S.Cox et al., 1997, [6],
M.M.Zdravkovich, 1997, [7]) revealed comparable conclusions as does the present work.
Mainly 3D effects are the cause of the non accuracy of the findings (e.g. P.D.Ditlevsen, 1996)
[8].
As a second step, the cavitation phenomena has been studied in several applications. First,
the full cavitation model implemented in Fluent (Singhal, 2002) has been tested comparing
findings with corresponding experimental data. The first case was a steady, cavitating flow
through a sharp-edged orifice by Nurick, 1976 [9]. Further, the unsteady turbulent flow around
a 2D NACA 0015 hydrofoil has been simulated using the cavitation model. This work was
based on the publications by Kubota, 1992 [10] and Berntsen et al., 2001 [11]. Results
revealed that depending on the case, the cavitation model offers useful results, but only in a
qualitative way. Accurate fittings with experiments are obtained only in few cases. The
theoretical validity of the present cavitation model could be questioned. Future work consists
of the prediction of damage caused by cavitation (comparing numerical results with
experimental databases e.g. Escaler, 2001) using adequate software tools. The final goal is
to apply the knowledge obtained to damage prediction in turbomachinery.
Numerical Simulation of Unsteady Turbulent Cavitating Flows Page 3
Contents
Abstract ________________________________________________ 1
Contents ________________________________________________ 3
1. Nomenclature ________________________________________ 5
2. Introduction _________________________________________ 9
3. Physical Background_________________________________ 13
4.3. Hydrofoil............................................................................................................... 58
4.3.1. Simulation based on the case of Kubota (1992)................................... 60
4.3.1.1. Computational domain and boundary conditions .................. 61
4.3.1.2. Simulation results and discussion.......................................... 62
4.3.2. Simulation based on the case of Berntsen et al. (2001)....................... 67
4.3.2.1. Computational domain and boundary conditions .................. 68
4.3.2.2. Simulation results and discussion.......................................... 68
Conclusions ____________________________________________71
Acknowledgements ______________________________________73
Annexes________________________________________________75
Bibliography ____________________________________________89
Numerical Simulation of Unsteady Turbulent Cavitating Flows Page 5
1. Nomenclature
Symbols
ε dissipation rate of k
C chord length
η, ηb, ηc dynamic viscosity
y wall distance
κ von Karman constant
y+ normalized wall distance
λ constant (Stokes hyp.)
Yν destruction of turbulent
viscosity ν kinematic viscosity
Yk,Yω destruction of k and ω νt turbulent kinematic viscosity
σs surface tension
τw wall stresses
Ωb bubble radius
Ψ volume forces
b bubble (phase)
c liquid (phase)
e element
g non-condensable gases
i,j,k tensor notation
Mathematical tools1
⎧0 for i ≠ j
δij Î δ ij ≡ ⎨ (Kronecker delta)
⎩1 for i = j
r n
∂
∇, ∇ Î ∇=∑ (Nabla = divergence operator)
i =1 ∂xi
D ⎧ ∂ rr ∂ ∂ ⎫
Î ⎨≡ + (∇u ) = + u i ⎬
Dt ⎩ ∂t ∂t ∂xi ⎭
1
Note that the number format in tables and graphics is the European (i.e. 0,1 is “0 point 1” ; 10.000 is “ten thousand”).
Page 8 Nomenclature
∂u i n
∂u
Tensor notation: = ∑ i for n-dimensional vectors
∂xi i =1 ∂xi
∂u i n
∂u ∂u i n n
∂u
Example: uj = ∑u j i or δ ij = ∑∑ i δ ij
∂x j j =1 ∂u j ∂x j i =1 j =1 ∂x j
Numerical Simulation of Unsteady Turbulent Cavitating Flows Page 9
2. Introduction
The current tendencies in the design of hydraulic turbomachinery of every type, especially
those requiring a high grade of optimisation as well as functional guarantee at critical regions,
are strongly influenced by the CFD (Computational Fluid Dynamics) simulation techniques. In
this context there have been developed many methods and software programmes for the
elaboration and analysis of fluid mechanical problems.
CFD has been applied for about 25 years and has been spread as a very promising tool for
present and future design and development of turbo-engines. The main interest is to establish
a solid and reliable shifting system between pure theory and experimental research, providing
a tool for calculating flow phenomena and enabling engineering applications. The
development of CFD strongly depends on the progress made by the computer science. In
spite of its rapid development, experiments still remain of crucial importance. The present
work is an example of the close connection and dependence that exist between theory and
practice.
Several important works analysing cavitation phenomena have been done in the past (e.g.
Kubota, 1992 [10]; Arndt et al., 2000-2002 [12], [13]; Kunz et al. 2000 [14]; Song and Qin,
2002 [15]). In the case of the cavitating flow around a hydrofoil, the work of Kubota shows
interesting results of a numerical investigation of the unsteady structure (velocity distribution)
of cavitation around a 2D-hydrofoil. These results lead to more detailed information about
processes related to cavitation as well as further theories upon cavitation and its
classification. Later on, other works and numerical investigations found that the cloud
Page 10 Introduction
cavitation observed in experiments using hydrofoils was a large-scale vortex with many small
cavitation bubbles. The importance of the interaction between large-scale coherent vortices in
the flow field and cavitation bubbles was recognized. A bridge between cavitation,
unsteadiness and turbulence was established.
Regarding unsteady turbulent flows, much more extensive works are available. However,
theories about turbulence modelling are strongly based upon empiric data and
acknowledgement. The first chapters of this work offer several aspects of the current
tendencies of CFD and the methods used by Fluent to simulate unsteady, turbulent flows.
Further, the cavitation model that has been used for all computations will be discussed. For
details upon cavitating flows please refer to the chapters describing the cases treated within
the present work or to the listed bibliography.
Within this work the main goal was to provide advanced knowledge for two-phase fluid flow
simulations. By means of the software program FLUENT for flow simulation, an unsteady,
incompressible two-phase fluid-flow has been simulated and analysed. The results will
contribute to the main project, which aims to analyse cavitating flows in hydraulic
turbomachinery. The detection of cavitation phenomena and the understanding of the
influence of the streaming processes on the collapse of the cavitation bubbles and the
material damage caused by this phenomenon, would lead to major progress in the
development of hydraulic turbines. The results would enable the simulation and the prediction
of the position and the intensity of bubble collapse in higher pressure zones, at a given flow
set-up. Experimental research has been carried out by Escaler (2001) upon a cavitating flow
around a hydrofoil containing a cylindrical obstacle at its leading edge. In order to accelerate
the obtaining of data, numerical simulation must replace experimental research.
As a first step, steady and unsteady one-phase flows on simple 2D geometries have been
modelled (circular cylinder, sharp-edged orifice) in order to acquire experience before
simulating more difficult cases. Changing characteristic flow parameters, like the Reynolds
number, made it possible to compare the results obtained to existing experimental data.
Concerning the flow past a circular cylinder, many questions about the physical processes
that occur still exist [6], [7], [8]. Experimental and previous numerical research data was used
for judging the best method (turbulence model) to be applied.
Numerical Simulation of Unsteady Turbulent Cavitating Flows Page 11
As a further step, a two-phase flow has been analysed for a sharp-edged orifice and a
hydrofoil using the cavitation model for multiphase flows. The sharp edged orifice could be
simulated with the steady state equations, while using the Nurick correlation to validate the
results. In the case of the hydrofoil, the flow was unsteady, turbulent and depending on the
settings, composed by two phases.
Before discussing the mentioned cases that are treated within this work, an introduction into
state-of-the-art methods for turbulent flow simulations will be described from a physical point
of view. The turbulent models provided in FLUENT as well as different theories and
approaches for the simulation of unsteady, two-phase fluids will be exposed.
Numerical Simulation of Unsteady Turbulent Cavitating Flows Page 13
3. Physical Background
Turbulent flows are characterized by fluctuating velocity fields. These velocity fluctuations
also lead to fluctuations of momentum or energy. Since these fluctuations can be of small
scale and high frequency, they are too computationally expensive to be simulated directly in
practical engineering calculations. Instead, the instantaneous (exact) governing equations
can be time averaged, ensemble-averaged, or otherwise manipulated to remove the small
scales, resulting in a modified set of equations that are computationally less expensive to
solve (e.g.: Reynolds averaged). However, the modified equations contain additional
unknown variables, and turbulence models are needed to determine these variables in terms
of known quantities (closure’s problem).
We generally distinguish several main turbulence models, being the Standard k-ε model the
mostly used to use engineering tasks. The following turbulent models provided in FLUENT
(2D) have been used:
• Spalart-Allmaras model
The basic equations which lead to the Reynolds Averaged Navier-Stokes equations and the
introduction of turbulence models are the momentum equations (2nd Newton’s law) and the
continuity equation. The continuity equation postulates that the change in mass in the volume
is equal to the fluxes through the surfaces of a volume element. In tensor notation and after
dividing by the control volume, this leads to:
Page 14 Memoria
∂u i ∂u 1 ∂p ∂Ψ 1 ∂τ ji
+uj i = − ⋅ − + ⋅ (3.1)
∂t ∂x j ρ ∂xi ∂xi ρ ∂x j
∂ρ r r r r
+ ∇( ρ ⋅ u ) = 0 ; for incompressible fluids: ∇( ρ ⋅ u ) = 0 (3.2)
dt
- 3 velocity components ui
- pressure p
- density ρ
⎡τ 11 τ 12 τ 13 ⎤
-
⎢
stress tensor τ ij = τ 21 τ 22 τ 23
⎥
⎢ ⎥
⎢⎣τ 31 τ 32 τ 33 ⎥⎦
The number of unknown variables is superior to the given equations, which obliges to
assume several simplifications. The main goal is to find an acceptable relation for the stress
tensor:
1 ⎛⎜ ∂u i ∂u j ⎞
⎟ ,
With the strain rate in tensor notation eij = ⋅ +
2 ⎜⎝ ∂x j ∂xi ⎟
⎠
2) The Stokes hypothesis which assumes, that the pressure is the average normal
stress. This results in:
λ 2
=−
η 3
Introducing these assumptions into the momentum equations (3.1), and with ν = η/ρ, one
obtains the Navier-Stokes equations [3]:
∂u i ∂u ⎡ ⎛ ∂u i ∂u j ⎞ ⎤
1 ∂p ∂Ψ
+uj i = − ⋅ − +
∂
⎢ν ⎜⎜ + ⎟ − 2 ⋅ ∂u k δ ij ⎥ (3.4)
∂t ∂x j ρ ∂xi ∂xi ∂x j ⎟ 3 ⋅ ρ ∂x
⎣⎢ ⎝ j ∂xi
∂x ⎠ k ⎥⎦
Together with the continuity equation (3.2), 4 equations and 5 unknowns are obtained. By
adding a supplementary transport equation for a scalar quantity and the relation between
pressure, density and the transported variable, one obtains a closed set of equations. An
example would be the use of the transport equation for the temperature together with the
ideal gas law (for a gas flow). Another one is the simulation of two-phase flows with the
cavitation model, where the volume fraction equation brings the relation between the state
variables of the two phases, and closes the Navier-Stokes equation set [chapter 2.4].
The Navier-Stokes equation computes the exact velocity field in a flow. Simulating a turbulent
flow using the theory described above leads to a computational technique called the DNS-
method (Direct Numerical Simulation). Somehow, the present-day possibilities of such
computations quickly reach the actual computer limits, even in basic engineering problems. It
is not even of major interest to calculate the exact magnitude of fluctuations within the whole
eddies’ scale [3]. CFD-methods use models, the so called turbulence models, introducing
further equation sets in order to simplify and reduce the computation problem. The LES-
method (Large Eddy Simulation) uses the Navier-Stokes equation set for computing big
vortexes within a turbulent flow. In the smaller scales, the turbulence is simulated with the
help of turbulence models. For complex flow-fields, this method is not yet applicable in
complex industrial problems. Mostly used is the RANS-method (Reynolds Averaged Navier-
Stokes) which has been developed for the use of turbulence models into the simulation of
turbulent flows, and builds the basis of the major part of computational software for fluid
mechanics at the same time. It largely leads to best results and approaches for complex
engineering tasks. Somehow, its limits can also be quickly reached as it has been
experienced within this work.
For simplicity, the present work is restricted to incompressible flows with constant
temperature: the compressibility of the flow is treated separately within the cavitation model.
How the cavitation model treats the fluid’s compressibility will be treated in later chapters. In
Page 16 Memoria
every case for which the cavitation model was not enabled, the flow was set to be
incompressible. Therefore, the continuity equation is reduced to:
rr
∇u = 0 (3.5)
∂u i ∂u 1 ∂p ∂Ψ ∂ ⎡ ⎛ ∂u i ∂u j ⎞⎤
+uj i = − ⋅ − + ⎢ν ⎜⎜ + ⎟⎥
⎟⎥
(3.6)
∂t ∂x j ρ ∂xi ∂xi ∂x j ⎢⎣ ⎝ ∂x j ∂xi ⎠⎦
As said before, turbulence is characterized by fluctuations within the velocity-field around the
steady/unsteady fluid flow. The instantaneous velocity values u i are split into a mean value
u i = U i + u i' (3.7)
t 0 +τ t 0 +τ
1
The mean values are defined by U i = ∫ u dt ∫ u dt = 0
'
, and the fluctuations:
τ t0
i
t0
i
∂U i ∂U i 1 ∂P * ∂ ⎡ ⎛ ∂U i ∂U j ⎞ ⎤
+U j =− ⋅ + ⎢ν ⎜⎜ + ⎟ − u i' u 'j ⎥
⎟
(3.9)
∂t ∂x j ρ ∂xi ∂x j ⎢⎣ ⎝ ∂x j ∂xi ⎠ ⎥⎦
This system of equations is not closed because the correlations of the velocity fluctuations
are unknown. Since the correlations have the form of the stress tensor, they are called
Reynolds stresses or turbulent stresses. The turbulent stress tensor ρ⋅ u i' u 'j is symmetric,
Numerical Simulation of Unsteady Turbulent Cavitating Flows Page 17
and brings 6 unknown components –the turbulent stresses- into the equation set. Similar to
the momentum equations, the other transport equations have also to be averaged [3].
Here comes the role of the turbulent models, which serve as a tool for the calculation of the
turbulent stresses, and which describe the influence of turbulence on the mean flow (it does
not describe the turbulent vortexes in detail!).
The basis of a major part of all turbulence models, is the turbulent viscosity assumption of
Boussinesq, which introduces the turbulent viscosity ν t to describe the turbulent correlation
u i' u 'j analogue to the formulation for the viscous stresses τ ij (3.3):
⎛ ∂U ∂U j ⎞ 1 ' '
u i' u 'j = −ν t ⋅ ⎜ i + ⎟ + u i u i δ ij (3.10)
⎜ ∂x ∂xi ⎟ 3
⎝ j ⎠
1 ' '
where k = u i u i is the turbulent kinetic energy.
2
∂U i ∂U i 1 ∂P * ∂ ⎡ ⎛ ∂U ∂U j ⎞⎤
+U j =− ⋅ + ⎢(ν − ν t )⎜⎜ i + ⎟⎥
⎟⎥
(3.11)
∂t ∂x j ρ ∂xi ∂x j ⎢⎣ ⎝ ∂x j ∂xi ⎠⎦
2
where P* = P + ⋅ρ ⋅k
3
Beyond the turbulence models, there are those that use the turbulent viscosity assumption of
Boussinesq and those who do not [3]. The models discussed in this work, all use this
principle. Within this category, we can distinguish between turbulence models using one or
two equations.
The turbulence models treated in FLUENT are all 2-equation models except the Spalart-
Allmaras model which is a 1-equation model and the RSM-model that uses seven equations.
Furthermore, the possibility of simulating a laminar flow without using any turbulent model is
provided (LAMINAR), as well as an option for simulations of inviscid flows (INVISCID), which
can be reasonably used for high Reynolds number flows, where the friction forces are
negligible.
Page 18 Memoria
The Spalart-Allmaras model is a relatively simple one-equation model that solves a modelled
transport equation for the turbulent viscosity [3]. This embodies a class of one-equation
models in which it is not necessary to calculate a length scale related to the local shear layer
thickness. This model is conceived especially for boundary layers1 subjected to adverse
pressure gradients. In its original form, the Spalart-Allmaras model is a low-Reynolds-number
model, requiring the viscous-affected region of the boundary layer to be properly solved. It
could be important to notice, that the near-wall gradients of the transported variable in the
model are much smaller than the gradients of the transported variables in the k − ε models.
Likewise, it will be difficult to rely on the Spalart-Allmaras model to predict the decay of
homogeneous, isotropic turbulence. Furthermore, this model can be criticized for its inability
to rapidly accommodate changes in length scale, such as might be necessary when the flow
changes abruptly from a wall-bounded to a free-shear flow.
The model was developed starting with a transport equation for the turbulent viscosityν t .
Advanced versions of the Spalart-Allmaras model deal with the modified form of the kinematic
turbulent viscosity, ν~t . In this case, the transported variable is identical to the turbulent
kinematic viscosity except in the near-wall (viscous-affected) region [3]. The transport
equation for ν~t is:
∂ν~t ∂ν~t ⎧⎪ ∂ ~ ~ 2
⎫⎪
+Ui = Gν +
1 ⎡
(
ν + ~ ) ∂ν t ⎤ + C ⎛⎜ ∂ν t
ν
⎞
⎟⎟
⎨ ⎢ ⎥ b2 ⎜ ⎬ − Yν (3.12)
∂t ∂xi σν ⎪⎩ ∂xi ∂xi ⎦ ⎝ ∂xi
t
⎣ ⎠ ⎪⎭
where Gν is the production of turbulent viscosity and Yν is the destruction of turbulent viscosity
that occurs in the near-wall region due to wall blocking and viscous damping. The terms σν
and Cb2 are constants. The model-formulations for Gν,Yν can be found in respective literature
(e.g. Coussirat, [3]; Hellster [4]).
1
For more information about the boundary layer theory, please refer to corresponding literature (e.g. Batchelor)
Numerical Simulation of Unsteady Turbulent Cavitating Flows Page 19
χ3
ν t = ν~t ⋅ f v1 , where the viscous damping function fv1 is given by f v1 = and
χ 3 + C v31
ν~t
χ= , where Cv1 is a constant.
ν
The standard k-ε model is a very well tested model [3,6]. There are several known flow
situations, for which this model provides insufficient results. These are for example strongly
swirling flows, flows with strong streamline curvature, or flows with strong positive pressure
gradients; they are all cases that do not follow the Boussinesq approach. For most of these
problems, there are correction terms that can be implemented in the model. The k-ε model is
only valid for completely turbulent flows. It cannot predict flows, where viscous effects play an
important role, as it is the case for example in near wall regions [3]. The k-ε model describes
the transport and dissipation of turbulent kinetic energy. The model bases on the one-
equation model of Prandtl & Kolmogorov, where the characteristic flow velocity Ut can be
described by the kinetic energy of the fluctuation movement [3]. The characteristic length
scale Lt is proportional to Prandtl’s mixing length lm:
1
∑
2
U t = k , where k = U i'
2 i
ν t= k Lt , where Lt = c d ⋅ l m
The aim of the k-ε model is to eliminate the mixing length dependence of the turbulent
viscosity. Therefore a k-transport-equation is deduced from the Navier-Stokes equation.
Instead of using the empiric mixing length for the definition of the characteristic length scale,
the dissipation rate ε of k is introduced as follows:
k2 3
k2
ν t= C µ ⋅ , with Lt = C µ ⋅ (Cµ is a constant)
ε ε
Page 20 Memoria
k- transport equation:
∂k ∂k ∂ ⎛ ν t ∂k ⎞ ⎛ ∂U ∂U j ⎞ ∂U i
+Ui = ⎜⎜ ⋅ ⎟⎟ + ν t ⎜ i + ⎟ −ε (3.13)
∂t ∂xi ∂xi ⎝ σ k ∂xi ⎠ ⎜ ∂x ∂xi ⎟ ∂x
⎝ j ⎠ j
ε- transport equation:
∂ε ∂ε ∂ ⎛ ν t ∂ε ⎞ ε ⎛ ∂U ∂U j ⎞ ∂U i ε2
+Ui = ⎜⎜ ⋅ ⎟⎟ + C1ε ν t ⎜ i + ⎟ − C 2ε (3.14)
∂t ∂xi ∂xi ⎝ σ k ∂xi ⎠ k ⎜⎝ ∂x j ∂xi ⎟ ∂x
⎠ j k
All model constants ( C µ , C1ε , C 2ε , σ k ) have been optimised for different types of flows [3].
At this time, several notices upon the boundary layer theory have to be figured out. At the wall
(y = 0), the wall adhesion condition is valid. In very proximate wall regions, cross movements
are tending to zero. A thin zone, called viscous sublayer, or boundary layer δ (laminar) is
created. While viscous effects on the energy containing turbulent motions are negligible
throughout most of the flow, the no-slip condition at solid interfaces always ensures that, in
the immediate vicinity of a wall, viscous effects are influential [3]. Although the thickness of
this viscous-affected zone is usually some orders of magnitude smaller than the overall width
of the flow, its effects extend over the whole flow field since, typically, 50% of the velocity
change from the wall to the freestream occurs in this region.
Near-wall zones are the main source of vorticity and turbulence because in these regions
exist high gradients of variables and the momentum transport is more vigorous. Also,
turbulence mixing is suppressed by the proximate boundary, causing a high reduction of the
transport across this layer. To take into account the non-isotropic nature of turbulence in
near-wall regions, a special wall treatment has to be defined [3]. The k- ε model does not
solve the regions affected by viscous effects. Instead, empirical functions (“wall functions”)
are used. Following values are defined:
τw
• the stress velocity uε =
ρ
U
• the normalized wall velocity u+ =
Ut
Numerical Simulation of Unsteady Turbulent Cavitating Flows Page 21
y ⋅U t
• the normalized wall distance y+ =
ν
The values for u+ and y+ are normalized with the respectively highest values for u and y. The
logarithmic wall-function can be written as follows:
1
+
ln(1 + κ ⋅ y + ) + C for 10 < y + < 150
u = κ (3.15)
+ +
y for y < 10
Using the k-ε models, the logarithmic law is applied when a “standard wall-function” is
selected for the approximation of the near-wall region (when 10 < y+ < 150) without adverse
pressure gradients. Modified wall functions called non-equilibrium wall functions are
implemented in the standard model for the computation of boundary layers affected by
adverse pressure gradients.
Fig. 3.1. Nomenclature of zero and adverse pressure gradient boundary layers.
Other near-wall treatment should be used, when y+ < 10, (e.g. using in the near wall region a
solver based in a simplified Navier-Stokes equation). Note that the Spalart-Allmaras and the
k-ω models automatically applies the appropriate near-wall treatment.
Page 22 Memoria
Apart from the standard k-ε model, FLUENT provides two further variants of the k-ε model:
the RNG k-ε model and the realizable k-ε model.
The RNG k-ε model was derived using a rigorous statistical technique (called renormalization
group theory) [3]. It is similar in form to the standard k-ε model, but includes the following
refinements:
1) The RNG model has an additional term in its ε equation that significantly improves the
accuracy for rapidly strained flows.
2) The effect of swirl on turbulence is included in the RNG model, enhancing accuracy
for swirling flows.
3) The RNG theory provides an analytical formula for turbulent Prandtl numbers, while
the standard k-ε model uses user-specified, constant values.
4) While the standard k-ε model is a high Reynolds-number model, the RNG theory
provides an analytically-derived differential formula for effective viscosity that
accounts for low Reynolds-number effects. Effective use of this feature does,
however, depend on an appropriate treatment of the near-wall region [3].
The realizable k-ε model is a relatively recent development and differs from the standard k- ε
model in two important ways:
1) The realizable k- ε model contains a new formulation for the turbulent viscosity.
2) A new transport equation for the dissipation rate ε has been derived from an exact
equation for the transport of the mean-square vorticity fluctuation.
The term “realizable” means that the model satisfies certain mathematical constraints on the
Reynolds stresses, consistent with the physics of turbulent flows. Neither the standard k-ε
model nor the RNG k- ε models are realizable [3].
k2
ν t= C µ ⋅
ε
The difference between the realizable k- ε model and the standard and RNG – models is that
Cµ is no longer a constant. In Fluent it is computed from a given formula (see source [3]), as a
function of the mean strain and rotation rates, the angular velocity of the system rotation, and
the turbulence fields k and ε.
Numerical Simulation of Unsteady Turbulent Cavitating Flows Page 23
The realizable model is likely to provide superior performance for flows involving rotation,
boundary layers under strong adverse pressure gradients [4]. For flows with complex
secondary flow features, the realizable k-ε model also shows best results.
The k-ω model is a semi-empirical model based on model transport equations for the
turbulent kinetic energy (k) and the specific dissipation rate (ω), which can also be thought of
as the ratio of k and ε.
Two versions of the k-ω model are implemented in Fluent. Both, the standard and the shear
stress transport (SST) k-ω model have similar forms, with transport equations for k and ω.
The major differences are a modified turbulent viscosity formulation for the SST model and a
gradual change at the standard k-ω model in the inner region of the boundary layer to a high-
Reynolds-number version of the k-ε model in the outer part of the boundary layer [4]. Only the
standard k-ω model will be discussed here.
The turbulence kinetic energy, k, and the specific dissipation rate, ω, are obtained from the
following transport equations:
∂ ⎛ ⎞
(ρ ⋅ k ) + ∂ (ρ ⋅ k ⋅ ui ) = ∂ ⎜⎜ Γk ∂k ⎟⎟ + Gk − Yk + S k (3.16)
∂t ∂xi ∂x j ⎝ ∂x j ⎠
⎛ ⎞
and
∂
(ρ ⋅ ω ) + ∂ (ρ ⋅ ω ⋅ ui ) = ∂ ⎜⎜ Γω ∂ω ⎟⎟ + Gω − Yω + Sω (3.17)
∂t ∂xi ∂x j ⎝ ∂x j ⎠
In these equations, Gk represents the generation of turbulent kinetic energy due to mean
velocity gradients. Gω represents the generation of ω. Γk and Γω represent the effective
diffusivity of k and ω, respectively. Yk and Yω represent the dissipation of k and ω due to
turbulence. All of the above terms are calculated as described below. Sk and Sω are source
terms [4].
µt µt
Γk = µ + and Γω = µ +
σk σω
where σk and σω are the turbulent Prandtl numbers for k and ω. The turbulent viscosity, µt, is
computed by combining k and ω as follows:
Page 24 Memoria
ρ ⋅k
µt = α *
ω
⎛ α 0* + Re t / Rk ⎞ ρ ⋅k
α * = α ∞* ⎜⎜ ⎟,
⎟ where Rk and α0* are constants and Re t =
⎝ 1 + Re t / Rk ⎠ µ ⋅ω
In the high-Reynolds number form of the k-ω model, in the outer part of the boundary layer,
α * = α ∞* = 1 .
∂u j
Gk = − ρ ⋅ u i' u 'j ⋅
∂xi
Gk = µ t ⋅ S 2 ,
where S is the modulus of the mean rate-of-strain tensor, defined in the same way as for the
k-ε model.
ω
Gω = α ⋅ ⋅ Gk ,
k
α∞ ⎛ α 0* + Re t / Rω ⎞ ρ ⋅k
where α= ⎜⎜ ⎟⎟ , Rω and α0* are constants and Re t =
α* ⎝ 1 + Re t / Rω ⎠ µ ⋅ω
Yk = ρβ * f β * kω Yω = ρβ f β ω 2
The parameters used for the calculation of Yk and Yω are mainly depending on the strain rate
tensor and will not be discussed in detail here.
Numerical Simulation of Unsteady Turbulent Cavitating Flows Page 25
The SST k-ω model presents a modified turbulent viscosity definition to account for the
transport of the principal turbulent shear stress. Other modifications include the addition of a
cross-diffusion term in the ω equation and a blending function to ensure that the model
equations behave appropriately in both the near-wall and far-field zones [4].
The wall boundary conditions for the k equation in the k-ω models are treated in the same
way as the k equation is treated when enhanced wall treatments are used with the k-ε
models. This means that all boundary conditions for wall-function meshes will correspond to
the wall function approach, while for the fine meshes, the appropriate low Reynolds-number
boundary layer conditions will be applied. For further details on the wall treatment with the
standard or the SST k-ω model, especially the way they are treated in Fluent, we suggest to
use literature references and Fluent’ s User’s Guide.
Cavitation is a phenomenon that appears in low pressure zones (e.g., hydraulic machinery)
causing significant degradation in the performance, inducing reduced flow rates, lower
pressure increases in pumps, load asymmetry and vibrations, and noise. Further
consequences like erosion on the material surfaces can lead to a irreparable damage of
turbomachinery. It is exactly this phenomenon that will be dealt with, as a major goal of this
project, to streaming processes around different geometries.
Liquid fluids can lift up tensile stresses when nuclei for the inception of cavitation are missing.
Real liquid fluids always contain nuclei for the inception of cavitation. When the vaporisation
pressure pv is reached, the liquid begins to gasify. A gasification process without heating, only
caused by low pressure, is called cavitation. When cavitation bubbles arriving in higher-
pressure zones collapse, they induce near the geometry-surface microscopic erosion-
processes which lead to surface damages at a bigger scale. The vaporisation pressure of
liquid fluids strongly depends on temperature, pv = pv(T).
Page 26 Memoria
There are many levels of modelling that may be used in multi-phase computations. In
general, one may distinguish between methods that employ an Eulerian framework for
both phases and those that employ Eulerian for the gas-phase and Lagrangian for the
liquid-phase. In the Eulerian-Eulerian framework, the simplest approach is to employ a
single continuity equation for both phases, with the fluid density being described as a
continuous function varying between the vapour and liquid phases [14]. At a more detailed
level of modelling, separate continuity equations for the liquid and vapour phases are
employed along with appropriate mass transfer terms to represent the phase-change
phenomena [14]. The gas-liquid interface is, however, assumed to be in dynamic
equilibrium and, consequently, mixture momentum equations are used. This is also the
case in the full cavitation model by Singhal et al., 2002 which has been implemented in
Fluent, although the continuity equations are computed only for the mixture.
The crucial requirement of multiphase algorithms is the ability to accurately and efficiently
span both incompressible and compressible flow regimes. For single phase applications,
time-marching techniques have long been established as the method of choice for high-
speed compressible flows, while artificial compressibility or preconditioning techniques
have enabled the extension of these methods to the incompressible and low-speed
compressible regimes [14]. Preconditioning methods essentially maintain proper
conditioning of the controlling time-scales of the time-marching system by introducing
appropriate pseudo-time derivatives. Indeed, it is widely recognized that the careful
selection of these derivatives is paramount for ensuring efficiency and accuracy over a
wide range of Mach numbers, Reynolds numbers and Strouhal numbers [14], [25].
In this work, we are confronted to isothermal two-phase flows, in which the densities of the
fluids are assumed to be functions of the pressure, but not the temperature. Under this
assumption, the energy equation is not solved and only the continuity and momentum
equations are considered (see also paragraph 3.1). It should be noted that the system is
still compressible because the pressure dependence of the densities gives rise to finite
acoustic speeds [14]. The isothermal compressible model serves as a useful intermediate
Numerical Simulation of Unsteady Turbulent Cavitating Flows Page 27
step between the incompressible model and the fully compressible system. The
development of the fully compressible model is currently underway.
To fulfil the main goal of this work, it is necessary to know if the isothermal compressible
model is generally valid for the class of cavitation problems considered here. The primary
interest lies in applications where sheet- and cloud cavitation appear. These flow fields
are characterized by large density ratios between the liquid and vapour states. In addition,
the flows are fully turbulent at the Reynolds numbers considered within this work. Most
problems also exhibit large-scale unsteadiness because of cavity re-entrant jets and
cavity pinching [14] (see paragraph. 4.3.).
The cavitation model implemented in Fluent is a mixture model following the Euler-Euler
approach, as described above. It has been developed from the so called “full cavitation
model” by Singhal et al. (2002). The phases are treated as interpenetrating continua. The
model solves the continuity equation for the mixture, the momentum equation for the mixture,
and the volume fraction equation for the secondary phase. Singhal proposes an isothermal
compressible model for cavitation. The compressibility of the phases remains “artificial”. A
shortcut through Singhal’s model and the mixture model in Fluent is given in paragraph 3.3.
Furthermore, two-phase flows can be generally characterized by the Stokes number, which is
a measure of the degree of influence of the bubble-movement on the stream of the
continuous phase [1]:
τ dyn ρ b ⋅ d b2 D
Sto = ; τ dyn = ; τ str =
τ str 18 ⋅ η U
where τ dyn is the response time relative to the velocity change, τ str the characteristic term of
the stream (D: characteristic length, U: velocity). The response time τ dyn can be calculated
from the momentum equation, considering only the drag force acting on the bubbles:
r
du b r r π ρ r r
m⋅ = 3π ⋅η ⋅ d b ⋅ (u c − u b ) = c D ⋅ d b2 ⋅ c ⋅ (u c − u b ) 2 , (3.18)
dt 4 2
when
r r r
du b 24 ⋅η π ⋅ d b2 ρ c ⋅ (u c − u b ) 2 6 ⋅ 3 ⋅η π ⋅ d b3 r r m
= r r ⋅ ⋅ = ⋅ (u c − u b ) ; = ρb
dt ρ c (u c − u b ) ⋅ d b 4 2 d b2 6⋅m π 3
d
6
Page 28 Memoria
r
du b 18 ⋅η r r 1 r r
= 2 (u c − u b ) = (u c − u b ) (3.19)
dt d ⋅ ρb τ dyn
r r
24 ρ (u − u b ) ⋅ d b
with: cD = ; Re = c c
Re η
Sto << 1 the bubbles are following passively the continuous flow
Sto >> 1 the bubbles are not influenced by the continuous phase
uc ub
t t
τ dyn
Fig. 3.2. Interaction between bubble-movement and the main stream of the continuous
phase: response time and velocity change.
In the case of a cavitating flow with largely subsonic velocities and bubble diameters between
10·10-6 m and 10 mm, it can be assumed that Sto << 1 [1]. All the studied cases within this
work will be regarded as flows with very low Stokes number, for which the mixture model
(Euler-Euler approach) for a homogenous flow is recommended. This leads to several
simplifications in modelling two-phase flows. The homogeneous repartition of the bubbles
within the continuous phase can be assumed, the phases are interpenetrating.
Numerical Simulation of Unsteady Turbulent Cavitating Flows Page 29
The full cavitation model is an extension of the VOF-Model (Volume of Fluid). The cavitation
model relies on the vapour mass fraction concept which introduces the vapour mass fraction
fb. The basic approach of the full cavitation model consists in using the standard Navier-
Stokes equation set for variable fluid density and a conventional turbulence model [5]. The
fluid mixture density ρ is a function of vapour mass fraction fb, which is computed by solving a
transport equation coupled with the mass and momentum conservation equations. The
density, f b relationship of the substitute fluid is:
1 fb 1 − fb
= + (3.20)
ρ ρb ρc
ρ
αb = fb ⋅ (3.21)
ρb
Following relation can be established between the volume fractions of both phases:
αc + αb = 1 or α c = (1 − α b ) (3.22)
ρ = ρ b ⋅ α b + ρ c ⋅ (1 − α b ) (3.23)
where αc and αb are the volume fractions of the continuous and the vapour phases
respectively.
∂ ∂
( ρ ⋅ f b ) + ui ( ρ ⋅ f b ) = ∇ ( Γ ∇f b ) + Re − Rc (3.24)
∂t ∂xi
Page 30 Memoria
The source terms Re and Rc denote vapour generation (evaporation) and condensation rates,
and can be functions of flow parameters like pressure or characteristic flow velocity, or fluid
properties like the liquid and gaseous phase densities, saturation pressure and liquid-vapour
surface tension.
The present model focuses on the use of simple rational formulations for phase change rates
(Re and Rc). To do this, the bubble dynamics consideration is firstly required. As specified in
the introduction, the liquid water phase contains plenty of nuclei for the inception of cavitation.
Thus, the primary focus is on proper account of bubble growth and collapse. In a liquid with
zero velocity slip between fluid and bubbles, the bubble dynamics equation can be derived
from the generalized Rayleigh-Plesset equation [5]. Neglecting the viscous damping term and
the surface tension term, the equation can be written as:
D 2 Ω b 3 ⎛ DΩ b ⎞ ⎛ pv − p ⎞
2
Ωb + ⎜ ⎟ = ⎜⎜ ⎟⎟ (3.25)
Dt 2 2 ⎝ Dt ⎠ ⎝ ρ c ⎠
D ∂ ∂
where = + ui (3.26)
Dt ∂t ∂xi
This equation provides a physical approach to introduce the effects of bubble dynamics into
the cavitation model. Combining the Rayleigh-Plesset equation with the two-phase continuity
equations (see e.g. [5]), one finally obtains for the source terms:
1
⎡ 2 ⎛ p − p ⎞⎤
1
⎡4 ⎤ 3 ρ ⋅ρ 2
Re − Rc = ⎢ π ⋅ n ⋅ α b2 ⎥ ⋅ b c ⋅ ⎢ ⋅ ⎜⎜ d ⎟⎟⎥ (3.27)
⎣3 ⎦ ρ ⎣ ⎝
3 ρ c ⎠⎦
This is the simplified equation for vapour transport, the right side representing the generation
rate of the bubbles (source term). Though the bubble collapse process is expected to be
different from that of the bubble growth, as a first approximation, equation (3.27) is also used
to model the collapse, when p > pv , by using the absolute value of the pressure difference
and treating the right side as a sink term [5]. The local far-field pressure p is taken to be the
same as the cell centre pressure. The bubble pressure pb is equal to the saturation vapour
pressure pv in the absence of dissolved gas [5].
In equation (3.27) all terms except “n” are either constants or dependent variables. In the
absence of a general model for estimation of the number density, the phase change rate
expression can be rewritten in terms of bubble radius, which is determined by the balance
Numerical Simulation of Unsteady Turbulent Cavitating Flows Page 31
between aerodynamic drag and surface tension forces. The correlation used in the nuclear
industry can be used:
0.061 We σ s
Ωb =
2 ρ c u rel
2
For bubbly flow regime, urel is generally fairly small (5-10%) of liquid velocity. By using various
limiting arguments, for example Ωb → 0 as α → 0, and the fact that per unit volume phase
change rates should be proportional to the volume fractions of the donor phase, the following
expressions for vapour generation \ collapse rates are obtained in terms of the vapour mass
fraction f. In the bubble flow regime, the phase change rate is proportional to u2rel [5].
However, in most practical two-phase flow conditions, the dependence on velocity can be
assumed to be linear [5]. The characteristic velocity uch reflects the effect of the local relative
velocity between liquid and vapour. It can be assumed, that the local turbulent velocity
fluctuations in most turbulent flows are of the same order that ui (1-10% of the mean velocity).
Therefore, as a first pragmatic approximation, ui can be expressed as the square root of local
turbulent kinetic energy k [5]. With the empiric coefficients Ce and Cc, one arrives to the
final formulation of the full cavitation model:
1 1
k ⎡ 2 pv − p ⎤ 2 k ⎡2 p − p⎤ 2
Re = C e ρb ρc ⎢ ⎥ (1 − f ) and Rc = C s ρc ρc ⎢ v ⎥ f (3.28)
σs ⎣ 3 ρc ⎦ σ ⎣ 3 ρc ⎦
Within this formulation, the effect of non-condensable gases has not been taken into account.
For further specification, and for information about the determination of the empirical
constants Ce and Cc, please refer to the source paper.
'
pturb = 0,39 ρ k
⎛ p' ⎞
p v = ⎜⎜ p sat + turb ⎟
⎟
⎝ 2 ⎠
Unlike this approach by Singhal assuming single-phase, variable fluid density flows, the
cavitation model implemented in Fluent has a more complicated configuration due to its
capability to account for the effects of the slip velocities between phases. This capability has
not been taken into account in our simulations for the reasons that have been specified at the
beginning of this chapter. This could be done by switching off the “slip velocity” button in the
mixture model panel. Another assumption made by Fluent, this one having much more
impact upon the accuracy of the simulation results and affecting the physics of the model
basements, is the incompressibility of the liquid phase: only the secondary phase is treated
as compressible for the formulation of the density of the substitute fluid.
The presence of nuclei of inception in real fluids is taken into account by introducing the
volume fraction of “non-condensable gases” αg, which can have significant effects on both the
physical realism and the convergence characteristics of the solution, even in a very small
amount [20]. As a consequence, equation 3.22. is replaced by:
αc + αb + α g = 1 (3.29)
The volume fraction of non-condensable gases is constant and can be set in the mixture
model panel in Fluent. The dynamic viscosity of the mixture is computed from:
η = α c ⋅η c + α b ⋅η b (3.30)
Numerical Simulation of Unsteady Turbulent Cavitating Flows Page 33
4. Numerical results
Several steps were necessary for the simulation of complex flows-fields. To gain insight in the
capabilities of the models that have been used (turbulence and cavitation), several cases
were modelled. To check the capabilities of the turbulence models, first, the unsteady flow
around a circular cylinder was simulated at different Reynolds-numbers. Then the cavitation
model was used in several examples, being the initial case a circular sharp edged orifice with
steady turbulent flow, and finally the turbulent, unsteady, cavitating flow around a hydrofoil
was simulated. For all cases, there are reliable experimental investigations, as well as some
numerical research that have been taken into consideration for comparison and evaluation of
the present simulation results.
The complexity of cavitating flows grows when the phenomenon becomes unsteady. This
happens for example in a cavitating flow around a hydrofoil. For the reasons mentioned
above, and in order to acquire more knowledge about unsteady flow fields, a classical
example has been chosen for a preliminary study of a non-cavitating flow: the circular
cylinder.
The flow past a circular cylinder can be studied separately at low Reynolds numbers and at
high Reynolds numbers. A complete understanding of a flow past a circular cylinder is
particularly elusive because the transition from laminar to turbulent flow occurs in a distinct
succession over a huge range of Reynolds numbers, and each transition state is sensitive to
extremely small disturbances [6]. The simulation will be done mainly at high Reynolds
numbers. Moreover, the simulation with FLUENT shows a big spread of results, especially in
the drag crisis region, at Reynolds 5.3 x 105, where the laminar boundary layer becomes
turbulent.
The flow around a circular cylinder falls into three distinct flow regimes: subcritical,
supercritical, and transcritical. The subcritical flow indicates purely laminar boundary layer
separation. In this regime, regular vortex shedding at a Strouhal number of about 0.2 is
observed over a range of Reynolds numbers from roughly 200 to 100000. The supercritical
regime, from Reynolds numbers of roughly 100000 to 4 million, is characterized by both a
dramatic rise in the Strouhal number and a loss of organized vortex shedding altogether [6]. It
Page 34 Memoria
is somewhere in this regime that transition to a turbulent boundary layer begins to occur on
the body at or near the point of separation. In the transcritical regime, above a Reynolds
number of roughly 4 million, periodic vortex shedding re-establishes but at a higher Strouhal
number of 0.26 – 0.30. The cylinder now experiences fully turbulent boundary layer
separation.
The critical Reynolds number varies greatly with the surface roughness, the intensity of
existing fluctuations (degree of steadiness) within the outer irrotational flow [17]. For example,
the critical Reynolds number becomes lower if either the roughness of the wall surface or the
intensity of fluctuations in the free stream is increased. In this work, we will consider the
geometries having smooth surfaces, which means, that the roughness shape fluctuations are
smaller than the thickness of the boundary layer.
Neglecting the volume forces, the general Navier-Stokes equations become (eq. 3.6):
∂u i ∂u 1 ∂p ∂ ⎡ ⎛ ∂u i ∂u j ⎞⎤
+uj i = − ⋅ + ⎢ν ⎜⎜ + ⎟⎥
⎟⎥
(4.1)
∂t ∂x j ρ ∂xi ∂x j ⎢⎣ ⎝ ∂x j ∂xi ⎠⎦
Within the near-wall region (y=0), we can assume that uy = ux = 0. For an unsteady flow, the
NS-equations are reduced to:
1 ∂p ⎛ ∂ 2u ⎞
⋅ = ν ⎜⎜ 2 ⎟⎟ (4.2)
ρ ∂x ⎝ ∂y ⎠ y =0
Upstream of the highest point of the cylinder the streamlines of the outer flow converge,
resulting in an increase of the free stream velocity u ∞ and a consequent fall of pressure with
x. Downstream of the highest point the streamlines diverge, resulting in a decrease of u ∞ and
a rise in pressure. Therefore, the shape of the boundary layer will strongly depend on the
pressure gradient, as shows figure 4.1.1.
⎛ ∂ 2u ⎞ ⎛ ∂ 2u ⎞
⎜⎜ 2 ⎟⎟ < 0 (decelerating), ⎜⎜ 2 ⎟⎟ > 0 (accelerating) (4.3)
⎝ ∂y ⎠ wall ⎝ ∂y ⎠ wall
In contrast to an accelerating external flow, where the velocity profile decreases with the y-
coordinate from a positive value to zero, a point of inflection characterises the velocity profile
in a decelerating flow, where the curvature (velocity second derivatives) is zero. The shape of
the velocity profiles in the lower picture (fig.4.1.1.) suggests that a decelerating pressure
gradient tends to increase the sickness of the boundary layer.
Numerical Simulation of Unsteady Turbulent Cavitating Flows Page 35
Besides, the existence of this point of inflection implies a slowing down of the region next to
the wall, a consequence of the uphill pressure gradient. Under a strong enough adverse
pressure gradient, the flow next to the wall reverses direction, resulting in a region of
backward flow. The reversed flow meets the toward flow at a certain wall point, where
∂u
= 0 . Here the fluid near the surface is transported out into the mainstream: the flow
∂y
separates from the wall.
In general, analytical solutions of viscous flows can be found only in two limiting cases,
namely Re << 1 and Re >> 1. In the case of the circular cylinder, it is interesting to study the
flow at a Reynolds number increased beyond 40, at which point the wake behind the cylinder
becomes unstable [17]. Photographs show that the wake develops a slow oscillation in which
the velocity is periodic in time and downstream distance, with the amplitude of the oscillation
increasing downstream. The oscillating wake rolls up into two staggered rows of vortices with
opposite sense of rotation. Von Karman investigated the phenomenon as a problem of
superposition of irrotational vortices. He concluded that a non-staggered row of vortices is
unstable, and a staggered row is stable only if the ratio of lateral distance between the
vortices to their longitudinal distance is 0.28. Such a staggered row of vortices behind a blunt
body is called von Karman vortex-street. Figure 4.1.2. shows this phenomenon for different
Reynolds numbers. The vortices move downstream at a speed smaller than the upstream
velocity u ∞ . This means that the vortex pattern slowly follows the cylinder if it is pulled
through a stationary fluid. In the range 40 < Re < 80, the vortex street does not interact with
Page 36 Memoria
the pair of attached vortices. As Re is increased beyond 80 the vortex street forms closer to
the cylinder, and the attached eddies (whose downstream length has now grown to be about
twice the diameter of the cylinder) themselves begin to oscillate. Finally the attached eddies
periodically break off alternately from the two sides of the cylinder. While an eddy on one side
is shed, that on the other side forms, resulting in an unsteady flow near the cylinder. Vortices
of opposite circulation around the cylinder change sign, resulting in an oscillating “lift” or
lateral force.
Fig. 4.1.2. Streaklines in the wake behind a cylinder for different Reynolds numbers Re
(Re = 32, 55, 65, 73, 102, 161).
The passage of regular vortices causes velocity measurements in the wake to have a
dominant periodicity. The frequency f is expressed as a non-dimensional parameter known
as the Strouhal number, defined as:
f ⋅d
St =
u∞
Experiments show that for a circular cylinder the value of St remains close to 0.21 for a large
range of Reynolds numbers [17].
Numerical Simulation of Unsteady Turbulent Cavitating Flows Page 37
Below Re = 200, the vortices in the wake are laminar and continue to be so for very large
distances downstream. Above 200, the vortex street becomes unstable and irregular, and the
flow within the vortices themselves becomes chaotic. However, the flow in the wake
continues to have a strong frequency component corresponding to a Strouhal number of St
= 0.21. Above a very high Reynolds number, say 5000, the periodicity in the wake becomes
imperceptible, and the wake may be described as completely turbulent.
At high Reynolds numbers the frictional effects upstream of separation are confined near the
surface of the cylinder, and the boundary layer approximation becomes valid as far
downstream as the point of separation. For Re < 5.3x105, the boundary layer remains
laminar, although the wake may be completely turbulent. The laminar boundary layer
separates at approximately 82º from the toward stagnation point. The pressure in the wake
downstream of the point of separation is nearly constant and lower than the upstream
pressure. As the drag in this range is primarily due to the asymmetry in the pressure
distribution caused by separation, and as the point of separation remains fairly stationary in
this range, the drag coefficient also stays constant at cD ≅1.2.
Important changes take place beyond the critical Reynolds number of Recr ≅ 5.3x105. In the
range 5.3x105 < Re < 3x106, the laminar boundary layer becomes unstable and undergoes
transition to turbulence. Because of its greater energy, a turbulent boundary layer is able to
overcome a larger adverse pressure gradient. In the case of a circular cylinder, the turbulent
boundary layer separates at 125º from the toward stagnation point, resulting in a thinner wake
and a pressure distribution more similar to that of potential flows. The next figure compares
the pressure distributions around the cylinder for two values of Re, one with a laminar and the
other with a turbulent boundary layer. It is apparent that the pressures within the wake are
higher when the boundary layer is turbulent, resulting in a sudden drop in the drag coefficient
from 1.2 to 0.33 at the point of transition. For values of Re > 3x106, the separation point
slowly moves upstream as the Reynolds number is increased, resulting in an increase of the
drag coefficient. It must be noted that the critical Reynolds number at which the boundary
layer undergoes transition is strongly affected by two factors, namely the intensity of
fluctuations existing in the approaching stream and the roughness of the surface, an increase
in either of which decreases Recr.
At Reynolds numbers of roughly 200 or less, many researchers have successfully computed
the Strouhal number and mean drag over a circular cylinder [6]. At higher Reynolds numbers,
however, two-dimensional numerical methods cannot predict the lift and drag forces
accurately, due to the increasingly prominent three-dimensionality of the real flow field.
Page 38 Memoria
Additional forces acting along the cylinder-axis are the main error sources regarding the 3D-
effects. Furthermore, according to latest researches, the difficulties of three-dimensional
modelling of flows around circular cylinders could arise due to non accurate physical models
and mathematical formulations concerning the conservation of helicity, which could block the
energy cascade and alter the Kolmogorov k-5\3 scaling (e.g. see reference [8]).
Any attempt to numerically model circular cylinder flow is complicated not only by the fact that
the flow above a Reynolds number of around 180 is three-dimensional, raising doubts about
the applicability of two-dimensional simulations. Additionally, transition occurs off-body in the
wake or shear layer at Reynolds numbers between roughly 200 and the supercritical regime
[6]. Without performing very expensive direct numerical simulations (DNS), this behaviour is
not captured by numerical methods that solve the Navier-Stokes equations on typical grids
used for hydrodynamic analysis. This deficiency may or may not be important at lower
Reynolds numbers, depending on how far behind the cylinder transition occurs and what
feature of the flow is of interest [6].
But it certainly has an adverse effect at higher Reynolds numbers for which transition occurs
at or near the separation point on the cylinder. For Reynolds numbers at and above the
supercritical regime, Reynolds-averaging with the use of a turbulence model is one way to
introduce the important effect of turbulence into a numerical simulation. However, without an
accurate built-in transition model, it is difficult, if not impossible, to model the important effects
of transition, particularly when it occurs on or near the body. The turbulent models used here
do not fulfil optimal conditions to predict these complex flow characteristics. It is not
surprising, then, that most numerical studies of the flow around a circular cylinder have
focused primarily on low Reynolds number flows less than about 1000.
In this work, an investigation is made into the ability to numerically predict the hydrodynamic
properties of a circular cylinder across a range of Reynolds numbers from 3000 to 1 million,
and thus confirm previous simulation researches whose results we partially described above.
The unsteady flow field is computed with the Reynolds-averaged Navier Stokes (RANS) flow
solver Fluent v6.1, employing a variety of turbulence models, as well as the laminar Navier-
Stokes equation to prove accuracy. The computations are two-dimensional and time-
accurate. Water was used as an incompressible fluid in all computations. Body forces have
been neglected, viscous and pressure forces having been assumed to be much larger.
Second order discretisation schemes have been set for all terms to compute the flow as well
as for the unsteady solver. The SIMPLE (semi-implicit pressure linked equation) velocity -
Numerical Simulation of Unsteady Turbulent Cavitating Flows Page 39
pressure correction algorithm has been enabled for all computations. Four different hybrid
meshes have been used for comparison: a fine mesh (mesh 01, ≈ 50000 nodes), an
intermediate mesh (mesh 02, ≈ 32000 nodes), a coarser mesh (mesh 03, ≈ 14000 nodes)
and a very coarse mesh (mesh 04, ≈ 10000 nodes). The smallest cell size in the direction
from the cylinder wall is 10-5 m for the fine mesh. The computational domain is shown in
figure 4.1.3. The intermediate mesh (mesh 02) around the cylinder is shown in figure 4.1.4.
The reference pressure is the ambient pressure, Poutlet = 1.01325 bar. The inlet velocity using
a constant profile was calculated from the suited Reynolds number:
u∞ ⋅ d ⋅ ρ
Re = ,
η
where d is the diameter of the cylinder, ρ the water density, and η the dynamic viscosity of
water at ambient temperature, and u∞ the constant flow velocity.
Page 40 Memoria
The drag and the lift forces (FD and FL respectively) are the x and y-components of the
resulting force upon the cylinder. They are very important for the study and the analysis of the
results of the flow over a circular cylinder. They are calculated as follows:
1 1
FD = ⋅ cD ⋅ ρ ⋅ u∞ ⋅ A and FL = ⋅ cL ⋅ ρ ⋅ u∞ ⋅ A ,
2 2
Fig. 4.1.4. Grid distribution for the circular cylinder, mesh 02.
Numerical Simulation of Unsteady Turbulent Cavitating Flows Page 41
where A is the specific area of the submerged body, cL the lift coefficient and cD the drag
coefficient. For a circular cylinder, A can be written as:
A = d ⋅b,
In a 2D case, it can be assumed that the cylinder length is infinite, or b>>d (although, 3D
effects along the axis seem to play an important role). Therefore, for the computation of the
drag and lift coefficients and especially the plot of the drag and lift convergence, the
computation area was made by unit of length, i.e., b = 1m, with a value of the area resulting in
A = 0.02m2. A typical plot of lift and drag convergence is shown in fig.4.1.5. The convergence
criterium is the constant sinusoidal periodicity of the drag and lift fluctuations, without
neglecting of course, the drop of the residual curves.
As it can be seen on the graph (fig.4.1.5.), the drag has a double frequency compared to the
lift: whereas one vortex appears upward, two vortices simultaneously shed in the wake
behind the cylinder. The time step of the calculations defined in Fluent has been chosen to be
between 1/10 and 1/20 of the vortex shedding. Also, the lift coefficient oscillates between
negative and positive values symmetrically around zero, while the drag is always positive.
Obtaining these convergence curves, the average drag value could be read off the graphic
and compared to the experimental data. The time step of one full oscillation could be
averaged from at least ten cycles of the drag and be converted into frequency for the
calculation of the Strouhal number.
Page 42 Memoria
The most difficult task regarding the simulation of a flow past a circular cylinder was the
setting of the time step. The velocity of the flow field could be calculated for each Reynolds
number. The corresponding time step could be obtained from the Strouhal number, which
was taken from literature (e.g. Zdravkovich, 2003) and set to be equal to 0.2 for all
computations. The resulting time step (1/f) corresponds to one cycle of the lift force. In order
to visualise the oscillating phenomenon with Fluent, the time step (1/f*) to be set for
simulation was chosen to be a fraction of the time step obtained from the Strouhal number
(e.g. f* = 20f).
Nevertheless, this adjustment required some experience since finding an adequate fraction of
the frequency in order to obtain clear convergence curves as shown in figure 4.1.5. was a
difficult task. Depending on the Reynolds number and on the turbulent model, the chosen
time step could be too large and in some cases even too small! The time step definition was
not only problematic for the simulations of the flow past a circular cylinder, but as we will see,
also for the computation of cavitating flows. It did often occur that a change in amplitude and
frequency of the lift coefficient led to convergence after changing the time step.
A mesh sensitivity analysis has been done using the Spalart-Allmaras and the SST k-ω
turbulent model at Reynolds 5.0e+04 and at Reynolds 1.0e+05. Asymptotic stability has
been obtained because the mesh 01 and the mesh 02 offer similar results for the drag
coefficient in all cases, as it can be verified in the figures from appendices A.1 and A.2.
However, this asymptothicity does not guarantee good numerical results. Choosing an
adequate time step for example, also influences the accuracy of the results.
Surprisingly, it can also be observed, that the mesh 04 offers best results compared to the
experimental data. It is widely recognized that it is very difficult to obtain good results near the
drag crisis region, and a big spreading of numerical results have been reported, see e.g.
Cox, 1997) [6] (see more details in appendix A.7-A.8). The illusive accuracy of the drag
coefficient obtained with the mesh 04 is probably due to pure coincidence, since the same
computation using the mesh 03 gives very unrealistic findings.
The Spalart-Allmaras model has also been used for the mesh sensitivity control (see A.3,
A.4). It is known from previous studies (see e.g. Coussirat, 2003 [3]) that this turbulent model
is very sensitive for this kind of analysis because it revealed a very strong mesh dependence,
so that the k-ω SST model should be rather considered as reference.
Evidently, all cases simulated for refinement have been executed with identical computational
set-ups and parameter settings. The Standard k-ε model has not been used for grid
adaptation, as computations using the mesh 01 and the mesh 02 require a different
boundary-layer approach compared to the coarser meshes. It has been assumed that the
Numerical Simulation of Unsteady Turbulent Cavitating Flows Page 43
present approach is valid for all Reynolds numbers, so that the mesh 02 could guarantee
accuracy for all simulations.
At Reynolds 5.0e+04 a big difference in the results using the meshes 01 and 02 (see
appendix A.2) is due to an imperfect convergence of the drag coefficient: (the average drag
could be clearly determined, but the drag oscillations are not perfectly sinusoidal,). This
example perfectly shows the impact of the time step definition on the drag and lift
convergence and thus on the obtained Strouhal number.
The computations were done in two stages. First a steady model was used to find a steady
solution and after this one converged, the unsteady module was started with initial values
corresponding to the steady solution. Using the 2nd order scheme for the unsteady solutions
led in most cases to divergence, so that only 1st order schemes have been used.
Several computations have been done at low Reynolds number using a solver for steady
flows (without turbulent models). At Reynolds 2300, the flow becomes turbulent. As show the
annexes A.5 and A.6, the results become worse with increasing Reynolds number, especially
the drag values. This is due to the effects of turbulence that cannot be captured with this
model.
Further, near the drag crisis region, the turbulent models discussed in the third chapter of this
work have been used for simulation. The table in fig.4.1.6. gives an overview of the solved
cases. From the results, no clear conclusions can be made concerning the most adequate
turbulence model to be used at these Reynolds numbers. It can still be confirmed that all
used turbulent models show difficulties to accurately predict the important effects of transition
(refer also to chapter 3 and the preliminary discussion within this paragraph). Regarding the
different wakes generated at the different Reynolds numbers, some plots show good
accuracy with visualisation experiments and other numerical simulations. In fig.4.1.7. the
Page 44 Memoria
vorticity plot at Re = 100 can be compared with the corresponding picture in fig.4.1.2. (even
though this graphic shows the streaklines). As the computational domain is shorter than the
experimental domain from fig.4.1.2., the regular vortices that begin to form at about 30-35
cylinder diameters in the wake cannot be seen in the simulation.
In fig.4.1.8. the vorticity magnitude at Re = 1.0e+05 is plotted and compared with the results
obtained by Cox, 1997 [6].
a) b)
The plots of velocity magnitude at different Reynolds numbers in fig.4.1.9. show a good
concordance with the theoretical approach discussed at the beginning of this chapter.
Numerical Simulation of Unsteady Turbulent Cavitating Flows Page 45
a)
b)
c)
In this part of the work, cavitation characteristics were computationally determined for sharp
edged orifices. Cavitation bubbles form due to a very low static pressure that occurs in this
case near the sharp inlet corner. This low static pressure is predicted by incompressible
potential flow theory, which indicates that flow around a sharp corner, (e.g. a corner with
infinitesimally small radius of curvature), will have infinite negative pressure [18]. This
physically impossible result is a direct consequence of the constant density restriction, which
can be encountered in some cases solved in Fluent without using the cavitation model. In real
cavitating flows, the density of the fluid (water) decreases with decreasing pressure, most
likely leading to a change in phase. The sharper the corner and the higher the velocity, the
more likely cavitation is to occur. High velocities near the orifice entrance generate those
zones of very low pressure right after the constriction, which reduces the flow rate (chocking
type phenomenon, [5]) and can lead to surface damage downstream of the orifice.
In the case of a sharp inlet, where the flow separates at the corner, the flow experiences a
vena contracta. A drawing in figure 4.2.1. shows a sharp entrance flow. After the vena
contracta, the flow re-establishes and a boundary layer starts growing until exit from the
orifice.
The resulting velocity profile will depend on the Reynolds number (based on the orifice
diameter) and development distance (i.e. x-position within the orifice). If the inlet pressure is
high enough and the cavitating edge sufficiently sharp, the cavity inside the nozzle grows and
reaches its outlet. As this occurs, the downstream ambient air finds a way to flow into the
orifice, resulting in a so called hydraulic flip (e.g. Vahedi et al, 2004) [18].
Fig. 4.2.2. Cavitation characteristics in Nurick’s experiments with a lucite orifice, L/d = 5.
Nurick, 1972 [9], has published extensive experimental data for cavitation in a sharp-edged
circular orifice. In his paper, the author develops a cavitation model to establish accuracy to
justify the experimental results. Using the continuity equation, the contraction coefficient cc
can be defined as:
Ac
cc = (4.4)
A1
where Ac represents the effective flow area through the contraction and A1 represents the
geometrical area of the orifice. The value of the contraction coefficient varies with the orifice
geometry and cavitation characteristics. For a very rounded entrance, the flow will not
separate and the coefficient of contraction will be unity. For a short orifice with a sharp
entrance, as it will be used in the present case, an empirical value of cc = 0.62 was taken
[5],[9].
Another relevant integral property of the flow is the coefficient of discharge, cd. The coefficient
of discharge represents the efficiency of the orifice between inlet and outlet and thus it is a
Page 48 Memoria
measure of whatever losses occur in the orifice. Nurick defines the coefficient of discharge
as:
m&
c d = cc ⋅ σ = , (4.5)
A1 ⋅ 2 ⋅ ρ c ⋅ ( P0 − P1 )
P0 − Pv
σ= (4.6)
P0 − P1
These definitions assume the losses to occur between the vena contracta and the orifice exit,
when the Bernoulli formulation is used.
Nurick plotted the coefficient of discharge versus the cavitation number on log-log axes in
order to verify the square root dependence of cd on σ. He observed that the data from the
cavitating region lay on a straight line with a slope of one-half, where cc is the value of the Y-
intercept. At some point, the value of the cavitation number is high enough so that the orifice
no longer cavitates. The higher values of σ occur when the difference between the upstream
and downstream pressure is small. At high values of σ the coefficient of discharge stays fairly
constant. This variation occurs because the coefficient of discharge is no longer a function of
σ, but depends on the Reynolds number instead. Thus, the point of inflection of Nurick’s
graphic indicates the inception of cavitation (see also fig.4.2.5.).
Within the present work, it has been intended to reproduce Nurick’s correlation, predicting
different discharge coefficients and comparing with the experimental results. The main
objective was to check if the cavitation model used would predict inception of cavitation.
Geometrical parameters of the circular orifice are D/d=2.88 and L/d=5, where D,d,and L are
respectively the inlet diameter, the orifice diameter and the orifice length. All sizes correspond
to the lucite orifice used by Nurick to visualize his results. Experiments and simulation have
been done with a fixed exit pressure, Pb=0,95 bar. The pressure taps that Nurick placed at ¼
and ½ d downstream of the inlet of the orifice to measure the cavity pressure have also been
simulated by defining 2 points at these locations and plotting the corresponding local static
pressure. The taps have been defined at 0.095mm in radial direction from the wall, where the
diameter of the orifice is 7.62mm (in the case of the coarse mesh this corresponds to 50% of
the distance from the wall of the first mesh point). The upstream static pressure, P1, was
Numerical Simulation of Unsteady Turbulent Cavitating Flows Page 49
varied between 9.6e+04pa and 5.0e+07pa to generate the flow at different coefficients of
discharge (fig. 4.2.5.).
Pressure Taps
a) b)
Fig. 4.2.3. a) Outline of the computational domain and the boundary type settings
b) Grid closeup and taps at ¼ and ½ orifice diameter (grid with 3700 nodes).
The flow is 2-D but axisymmetric. The inlet boundary condition is the specified upstream
pressure with an inlet velocity assumed to be zero. Upper boundaries defining the contours of
the orifice are no-slip walls. At the pressure outlet boundary, the fixed pressure of Pb = 0.95
bar has been set for all cases. Figure 1 shows a close-up of the coarse grid (fig. 4.2.3.a) and
an outline of the computational domain with all boundary condition types (fig.4.2.3.b). The
smallest cell size in the near-wall region is for the coarse grid a 2.5% of the orifice diameter.
All cases have been calculated with double precision. The solving strategy used was the
steady, SIMPLE algorithm. Second order discretisation schemes have been used for the
calculation of the flow terms.
All of the described applications have been generated using the Standard k-ε model with
Standard Wall Functions. A plot of the y+ values has shown that applying this model is
accurate, since 25 < y+ < 60 along the orifice wall for the intermediate grid.
The cavitation model (Multiphase Model – Mixture) was configured with a vaporization
pressure of Pvap = 3540 Pa, a percentage of non condensable gas of 1.5e-05 and no slip
velocity between the phases. Water-vapour with a density of 0,02558 kg/m3 and a viscosity of
1.26e-06 kg/(m·s) has been defined as the secondary phase.
Page 50 Memoria
To guarantee mesh independence results, different meshes have been tested for the
simulation. Best results have been obtained with only one mesh, used by Singhal with a total
of 3700 nodes (intermediate grid). First refinement has been done with a second mesh of
7400 nodes (fine mesh). Finer meshes did not provide acceptable results (compare with
Fluent v6.2 at the end of this chapter). For comparison, a third mesh of 1900 nodes has been
successfully tested. The refinement analysis has been done using the same computational
parameter set-up and conditions for all meshes. The inlet pressure has been set at 5.0E+05
Pa.
Fig. 4.2.4. Mesh refinement analysis for Pinlet = 5.0E+05 pa using 3 different grids.
Figure 4.2.4. shows the plot of static pressure along the orifice wall and provides evidence
about the exactitude of the results obtained with the simulations of the flow through a sharp
edged orifice. The red curve represents the pressure plot obtained with the 3700 node grid.
The fine grid (7400 nodes) is represented by the blue curve, while the black curve shows the
pressure characteristics along the wall for a grid with 1900 nodes. As the graphic shows, the
differences between both curves are minimal at the region where the cavity is located (a 0.1
to 3% of error), so that the intermediate grid seems to offer a good approach. The
unpredictable hydraulic flip as well as the constant position of the reattachment zone at the
orifice wall documented with the vector velocity plot (see next paragraph) may be a result of a
too coarse grid. The error percentage at C = 0.012 where reattachment should occur at lower
Numerical Simulation of Unsteady Turbulent Cavitating Flows Page 51
inlet pressures, is about 11% and corresponds to the maximum discrepancy of the static
pressure curves between the fine and the intermediate mesh.
Varying the inlet pressure led to a drop of the y+ values near the orifice wall. This would have
implied the use of enhanced wall treatment (k-ε model). This could not be achieved because
the condition for mesh refinement which requires a strict similarity at parameter definitions
would have been infringed (this is only valid if k-e model is used, for S-A and k-ω models this
problem does not appear, see also chapters 3.1.1. – 3.1.3.). Finer meshes could not be used
for the same reason. Because of unknown reasons, the simulation of the mentioned cases
led to insufficient results when a finer mesh was used. As a consequence of our mesh
sensitivity analysis, the intermediate mesh has been used for all simulations, assuming that
further refinement would provide similar findings.
The inlet pressures and the corresponding discharge coefficients and cavitation numbers
calculated as described before (eq. 4.5 and 4.6) are listed below [Fig. 4.2.5.].
For each case, the inlet pressure has been increased according to the table in Fig. 4.2.5. The
inception of cavitation was verified by plotting the contours of static pressure and asserting
whether the minimal value exceeded or not the vapour pressure. This was reached at σ =
1.76, which fits very well with Nurick’s results (fig. 4.2.6.).
Page 52 Memoria
Nurick’s correlation can be plotted following the data of the table: from σ = 1.00 to σ = 1.76
the flow cavitates (Pstat,min < Pvap) and is described by equation 4.5. When σ < 1.76 , the cavity
pressure should remain constant and equal to the vapour pressure until hydraulic flip occurs.
For an orifice with L/d = 5, Nurick observed that hydraulic flip occurred slightly after the
inception of cavitation. It was not possible to simulate the hydraulic flip, even perfect
sharpness of the orifice being accomplished. Hydraulic flip cannot be predicted without the
use of a VOF model (see also Vahedi et al., 2004 [18]). Also, deficiencies in the software
implementation make the cavity pressure decrease much beneath the vapour pressure,
which cannot correspond to reality. Fluent programmers argue that this is a difficult numerical
task. In the latest version, the default minimal pressure corresponds to 50% of the
vaporization pressure. Fluent assures that the default value in the Fluent 6.2 version will be
85% of pvap1.
For σ > 1.76 , equation 4.5. cannot be applied anymore; the cavity pressure must substitute
the vapour pressure and the losses described by the discharge coefficient can be assumed to
be constant when σ varies. Unfortunately, Nurick plotted his results only for L/d = 10 and L/d
= 6. Since the difference between both graphics are minimal, we assumed the same for an
orifice with L/d = 5.
1
Fluent Inc. 2004, personal communication.
Numerical Simulation of Unsteady Turbulent Cavitating Flows Page 53
Cavitation occurs after the local static pressure has reached the vapour pressure of water at
ambient temperature. Nurick observed the appearance of a fuzzy region near the inlet at P1 =
1.67E5 Pa (corresponding to σ = 2.27 if we assume that Pvap = 3540 Pa also for Nurick’s
experiments) for which the cavity pressure measured with the tap did not reach the vapour
pressure yet, see Fig. 4.2.2. In the present simulation, the contours of volume fraction also
reveal the presence of vapour within the continuous water phase before the vapour pressure
was reached, see Fig. 4.2.7. This can also be considered as a good result, although the fuzzy
region has been observed at higher inlet pressures (σ = 1.82) than in Nurick’s experiments.
Thus, Fluent predicts inception of cavitation before the static pressure reaches the vapour
pressure. This phenomenon is real, and describes what Nurick calls “fuzz”, occurring when
the flow is two-phase. Besides the differences between experiments and the present
simulation, regarding the inlet pressure at which “fuzz” can be observed, another divergence
is to be noticed.
Nurick describes a separation region that lengthens to about four orifice diameters before
reattachment occurs, and assumes it to be the consequence of the fuzzy region of the orifice
inlet. Numerical results obtained at σ = 1.82 show a separated flow region but the
reattachment occurs at less than ½ orifice diameter. It seems important to say, that according
to the obtained simulation results, the position along the x-axis at which reattachment occurs
varies very few when the inlet pressure is increased. For very high inlet pressures where the
cavity length exceeds two orifice diameters, reattachment occurs already at a distance of
about ½ orifice diameter (fig. 4.2.8.). The cavity length is about two orifice diameters and is
Page 54 Memoria
much shorter than it was observed by Nurick, who predicts hydraulic flip at σ = 1.7. Figure
4.2.8. shows a plot of the vapour volume fraction at σ = 1.0 and the corresponding vector plot
of the velocity magnitude: reattachment occurs already at ½ orifice diameter. This result can
be explained by the fact that the plot shows the velocity field of the mixture.
Finally, a series of simulations have been done to reproduce the relation between the cavity
pressure and the upstream pressure. A comparison between the cavity pressure (measured
at the two taps) obtained by Nurick as a function of the upstream pressure and the same data
obtained in the present simulation is shown in fig. 4.2.9.
The linear reduction of the cavity pressure obtained by Nurick for the tap at ¼ diameter gives
evidence of the fact that the vena contracta is at or near that location. This does not happen
in the simulation, which lets us deduce that the vena contracta is much shorter in numerical
calculations than the experiments show. The curve obtained by Nurick for the tap placed at ¼
orifice diameter corresponds to the simulation plot of the static pressure at
a) b)
Fig. 4.2.8. a) Volume fraction of the vapor phase αv and b) vector velocity plot at σ = 1.02.
x = ½ orifice diameter. The plateau region in Nurick’s curve obtained for the first tap
corresponds to the appearance of fuzz, and thus what he called the inception of cavitation.
The location of the taps was shown in figure 4.2.3. Further increasing the upstream pressure,
leads to a drop of the static pressure within the cavity beyond the vapour pressure. A cavity
forms and grows as shown in figure 4.2.10. with increasing σ. At σ = 1.02 hydraulic flip is not
reached.
Numerical Simulation of Unsteady Turbulent Cavitating Flows Page 55
The solutions show robust convergence characteristics, with residual curves dropping more
than ten orders of magnitude. Figure 4.2.11. shows typical convergence plots for the orifice
cases: while an initial plateau represents initial condition errors being convected out, the
following abrupt curve drop shows the solutions converging rapidly.
The poor results obtained for high pressure differences between inlet and outlet when using
fine meshes are probably due to software deficiencies. Other meshes with finer closure have
been tested. Only the intermediate mesh offered acceptable results with good convergence
characteristics for these conditions. The pressure correction factors for the AMG solver have
been set at 0.4. Wall function options have been adapted to the obtained y+ minimum and
maximum values for each grid. Under-relaxation factors have been set very low as well as a
maximum of iterations per time step to ensure a more stable iteration process. Other turbulent
models than the standard k-ε model have also been tested with these meshes without
success (Spalart-Allmaras, RNG k-ε, SST k-ω).
At the present time, calculations have been done with the new version of Fluent. Good results
have been obtained: very fine meshes could be tested using all above mentioned turbulence
models. The earlier version was numerically very unstable, and allowed computations using
the cavitation model only coupled with the k-ε turbulence model. A mesh refinement has been
done using five different meshes. The results obtained for this mesh sensitivity analysis can
be seen in figure 4.2.12. The better quality of the results at this mesh refinement, and
particularly the numerical stability of Fluent v6.2 encourages to perform computations of
cavitating flows using this version.
The unsteady turbulent cavitating flow past a hydrofoil described in the next chapter has been
simulated using Fluent 6.1. The latest version will have to be used for computations in future
work considerations.
Numerical Simulation of Unsteady Turbulent Cavitating Flows Page 57
Fig. 4.2.12. Mesh sensitivity control using 5 different meshes (“skw_0.5” coarsest
and “skw_4.0” finest, standard k-ω model, Fluent v6.2
Page 58 Memoria
4.3. Hydrofoil
The turbulent wake generated by streamlined bodies like hydrofoils has been investigated
numerically as well as experimentally in the past. But the details of flow near the trailing edge
are still poorly understood [15]. The unsteady turbulent wake behind a cavitating hydrofoil is
even more complex. The formation of individual bubbles and subsequent development of
attached cavities, bubble clouds, etc., makes the unsteady structures in the wake of a
cavitating hydrofoil fundamentally different from those in the wake of a non-cavitating
hydrofoil. A detailed description and analysis of cavitation phenomena with different cavitation
number and angles of attack can be found in Arndt et al., 2000, and Song and Qin, 2001.
a) b)
Fig. 4.3.1. Pressure and vorticity fields of a) sheet cavitation [σ = 1.5], b) sheet/cloud cavitation
[σ = 1.3].
It is known that cavitation inception occurs at composite parameter (σ/2α) of about 8.5, where
σ is the cavitation number and α the angle of attack [15]. It has been reported by Song and
Qin, that when σ/2α is slightly lower than 8.5 (between 6 and 8), a bubble cavity will first
appear near the leading edge on the suction side when the lift coefficient is maximum, at the
location with minimum pressure [15]. This bubble, which is also an eddy, slides down along
the suction surface and the lift decreases to a minimum when it arrives at the tail of the foil.
The lift starts to increase again as the eddy moves away in the wake. A new bubble cavity
occurs when the lift increases to a maximum again and the process repeats itself. This
Numerical Simulation of Unsteady Turbulent Cavitating Flows Page 59
phenomenon is called bubble cavitation. Bubble cavitation can change to a less stable
bubble/cloud cavitation. A cloud cavity may follow a small bubble at the instant of maximum
lift. When the cloud cavity collapses near the trailing edge, the lift drops to a minimum, the
drop being this time characterized by intense fluctuations [15].
Further reducing σ/2α (between 4 and 6) results in the generation of sheet or more unstable
sheet/cloud cavitation near the leading edge of the suction side (fig.4.3.1.). This is the most
dynamic flow regime related to partial cavitation. This cavity grows and breaks off to form a
cloud cavity at about 1/3 chord length due to the re-entrant jet [15], [19]. The cloud cavity
collapses at near ¾ chord to generate strong shock wave [15]. Song and Qin estimates the
average Strouhal number to be about 1.5 accentuating the fact that this is an unsteady and
“inherently turbulent phenomenon with significant degree of randomness”. Curiously,
Berntsen calculates a Strouhal number based on the maximum foil thickness as reference
length to be about 0.2 for similar conditions.
cL
Fig. 4.3.2. Comparison between phase locked photos and numerical simulation
of void fraction for a NACA0015 hydrofoil (σ = 1.0, α = 8º) and lift plot, Kawakami, 2004.
Figure 4.3.2. shows the contours of void fraction for a NACA0015 hydrofoil at α = 8º and σ =
1.0, which are fairly close parameters to our present simulation. These results have been
taken from Kawakami (2004) and are as interesting as some works by Song and Qin [15]
taken as a reference to this work. The numerical findings of the mentioned authors based on
Navier-Stokes equations for compressible flows and a virtual single phase equation of state,
which is quite different from the assumptions taken in Fluent 6.
In most cases, according to Song and Qin, the cloud cavity breaks up and collapses into
more than one piece. When σ/2α falls below 4, the flow pattern is changed significantly. The
maximum cavity length exceeds the chord length and the cavity oscillates between partial
cavity and super cavitation [15].
Page 60 Memoria
Within the present work, we intended to predict bubble, bubble/cloud and sheet/cloud
cavitation using Fluent 6.1. Although most of the documented investigations used LES (Large
Eddy Simulation) or/and user defined cavitation models for the simulation of this type of
cavitation phenomena for a flow around a hydrofoil, we mostly tested Fluent’ s ability to solve
the case(s) described below using RANS equations and the implemented multi-phase
model. As for the case of a sharp-edged orifice and despite of a very rigorous procedure, the
results have been good only for single settings and under very restricted conditions. We first
intended to follow the investigation of Kubota (1992) and base our research on his results for
a hydrofoil under different angles of attack and a Reynolds number of 3e+05. Further, another
case has been simulated, the calculations being set up to resemble the simulations taken
from the paper of Berntsen et al. (2001). Here, a Reynolds number of 1.2e+06 based on
chord length was used; the flow was at σ/2α of 5, with an angle of attack of 8º.
Kubota already found important results concerning bubble, bubble/cloud and sheet/cloud
cavitation. Like Arndt, Song, Qin, etc. he suggested that vortex cavitation is often observed
downstream of attached cavitation. It is caused by vorticity shed into the flow field just
downstream of the cavity. Such vortex cavitation generates a large cavitation cloud under
certain conditions [10]. The vortex cavitation impinges on the body where its subsequent
collapse results in erosion. Kubota describes the cavitation phenomenon more precisely in
his formulation of bubble two-phase flow cavity model.
The cases of Kubota 1992 and Berntsen, et al. 2001 were selected to extensively check the
capabilities of the cavitation model implemented in Fluent. These benchmarks were selected
because both experimental measurements and complete information concerning the
parameters used in the numerical computations are available.
In the case of Kubota, we have restricted our investigation to the case where the flow is at σ
=1.2. This corresponds to a composite parameter σ/2α of about 4.3. A suitable mesh has
been generated for the NACA0015 hydrofoil, see Fig 4.3.3. The computation was performed
at α = 0º, 8º and 20º. The Reynolds number Re, based on the uniform flow velocity and chord
length of the hydrofoil, was 3e+05. The computation at α=0º was performed only for non-
cavitating conditions to evaluate numerical accuracy obtained in the flow pattern prediction.
Numerical Simulation of Unsteady Turbulent Cavitating Flows Page 61
The grid is a hybrid grid type with a total of about 30700 nodes. The smallest cell size in the
direction from the wall is 0.005mm, which corresponds to approximately 1/350xRe –1/2. The
results in the first grid point yield y+ < 3.
Fig. 4.3.3. Grid distribution for the hydrofoil, the Kubota case reproduction.
The distance from the leading edge to the inlet boundary is 16 to 30C, where C is the chord
length of the hydrofoil (C = 0.1m). The distance between trailing edge and the outlet boundary
of the computational domain is 20C. Figure 4.3.3. shows a close-up of the 2D grid of the
NACA 0015 – hydrofoil used for this simulation. The parabolic-shaped inlet boundary is a
velocity inlet with constant velocity profiles. The outflow boundary condition is set with
constant pressure (reference pressure). This configuration permits the flow to be treated as
unconfined. The angle of attack was simulated by defining the inlet velocity by its components
corresponding to the appropriate Reynolds number. The Reynolds number Re, pressure
coefficient cp and cavitation number σ are defined as follows:
where u ∞ is the uniform flow velocity, ρ c is the density of liquid water, η c the dynamic
viscosity of water, pref the reference pressure or operating pressure and pvap the water vapour
pressure. Using these definitions and the constant parameters mentioned above (σ, Re, C,
ρ c , η c ) the reference pressure could be calculated. The vapour pressure was 2400 Pa.
Page 62 Memoria
Fig. 4.3.4. Contours of static pressure around the NACA 0015 hydrofoil, α = 0º, Re = 3e+05.
Kubota observed a little numerical difference near the trailing edge (not plotted here) due to
the unsteadiness of the flow field. This did not result in the present calculation, which
constitutes a better approach to the experimental findings. As shown in figure 4.3.4., the
pressure field is perfectly symmetric.
The computed result for the pressure coefficient at α = 8º agrees also very well with the
experiment and is similar to the result obtained by Kubota. Kubota found for α = 0º at the
trailing edge slightly different pressure coefficient distribution for the two sides of the hydrofoil,
which is physically not correct because of the symmetric flow conditions. The present
simulation also shows better results for this purpose. Figure 4.3.5.b) shows the time-
averaged pressure distribution on the foil surface and trailing edge at α = 8º. The result
agrees fairly well with the experimental data, especially at the trailing edge of the hydrofoil,
where the computed cp by Kubota was almost 0.3 higher than the experimental one.
The lift coefficient cL computed with the present simulation was about 0.780 (cL, experiment =
0.946), while the lift coefficient predicted by Kubota was only about 58% of that obtained by
Numerical Simulation of Unsteady Turbulent Cavitating Flows Page 63
experiment (cL, Kubota = 0.5505). As a reference, Berntsen (2001) obtained cL, Berntsen = 0.740
[11]. The upper curve describes the pressure coefficient on the suction side of the hydrofoil,
where the curve drops far into negative values, whereas the lower curve corresponding to the
plot for the pressure side indicates that the local pressure is always bigger than the defined
reference pressure. The negative cp- values on the suction side near the leading edge of the
hydrofoil points out the location of the cavitation inception. Fig.4.3.6 shows a plot of the
absolute pressure around the hydrofoil. The negative pressure is physically impossible and
must be understood as the inception of cavitation. In the wake of the hydrofoil, the formation
of unsteady vortices can be observed, although their frequency and circulation characteristics
do not seem to be realistic.
a) b)
Fig. 4.3.5. Foil surface and wake pressure coefficients, NACA0015 for non-cavitating conditions:
Kubota (Re = 3 x 105); experiment (Re = 6 x 105); Hess-Smith method (Re = 3 x 105);
5
present calculation (Re = 3 x 10 ). a) α = 0º, b) α = 8º.
This is highlighted by a comparison between the plot of vorticity magnitude of the present
computation and similar experimental data or numerical results, e.g. fig 4.3.1.a) obtained
using LES. The Strouhal number based on the maximum foil thickness as reference length
was set to 0.2 as suggested by Kubota, 1992. The time step used was varied between 1.0e-
04 and 1.0e-06. This gives 50-50000 time steps per cycle. More than 106 iterations have been
performed in isolated cases to obtain lift convergence. First, a very tapered shaped trailing
edge was suspected to be a possible cause of the steady character of the flow in the wake.
A second mesh has been generated for verification, consisting of a hydrofoil having a well
rounded trailing edge. No changes within the velocity or within the pressure field have been
noticed. This second mesh has been further used for the simulation of the case of Berntsen,
Kjeldsen and Arndt, 2001, mentioned before, and represented in figure 4.3.10.
Page 64 Memoria
At this point, as cavitation conditions were manifest, the same case was tried to be solved
activating the cavitation model in Fluent. Firstly no adequate results have been obtained. All
simulations have been done using the Spalart-Allmaras turbulence model before starting the
computation using the multiphase model in order to obtain better initial values. The value of
y+ was for all cases smaller than 4, so that enhanced wall treatment could be applied while
using the k-ε turbulence model. For the same cavitation number, different flows have been
calculated varying those parameters which mostly influence the cavitation phenomenon or
facilitates its computation: operating pressure, Reynolds number, time step, discretisation
factors, angle of attack, etc.
pabs< 0 Î CAVITATION
Fig. 4.3.6. Contours of absolute static pressure around the NACA 0015 hydrofoil, α = 8º, Re = 3 x 105.
The discretisation factors have been chosen very small at the beginning of each computation.
The pressure correction parameter in the Fluent AMG solver has been set up to 0.4 as
recommended for “complex” cases. Second order schemes have been changed to first order
schemes in isolated cases to test stability. For all simulations, the double precision mode has
been used. Changing the angle of attack to α = 20º did not offer better results. As generally
observed, the cases were solved very well as long as the cavitation model was not activated.
Also in the case α is set to 20º, the static pressure plot (similar to figure 4.3.6. for α = 8º)
shows evidence of the cavitation phenomena at the leading edge of the NACA0015 hydrofoil.
The activation of the cavitation model led to divergence/ no clear convergence of the solver.
To assure that our procedure was adequate, a tutorial case of Fluent has been reproduced.
This 2nd case computes also a symmetric hydrofoil at 8º angle of attack at the same
Reynolds and the same cavitation number as in the case of Kubota. The k-ε RNG turbulent
model is suggested here, as well as a vapour pressure of 80000 Pa. These parameters for
Numerical Simulation of Unsteady Turbulent Cavitating Flows Page 65
the case of Kubota were set and the Reynolds number was increased to 6.0e+05. To obtain
the same cavitation number as Kubota (σ = 1.2), the operating pressure was set equal to the
ambient pressure. The case surprisingly converged clearly, showing characteristics of sheet
cavitation.
Fig. 4.3.7. Contours of vapour volume fraction around the NACA 0015 hydrofoil and cavity length,
σ = 1.2, α = 8º, Re = 6 x 105.
However, no oscillations of the cavity have been detected showing evidence of the cloud
cavitation that should be caused by the sheet cavity due to re-entrant jet.
a) b)
Fig. 4.3.8. Vector plot of the velocity magnitude of the NACA 0015 hydrofoil, α = 8º;
a) leading edge, b) trailing edge.
The contours of the vapour volume fraction are shown in figure 4.3.7., as well as plotted
versus chord length for the suction side (blue curve) and the pressure side (red curve). A
cavity length of about 3/5 chord length can be deduced from the graph. The vector plot of the
Page 66 Memoria
velocity magnitude at the leading edge and the trailing edge of the hydrofoil are plotted in
figure. 4.3.8. At the rear of the hydrofoil, a detachment of the flow building a vortex is to be
observed. The detachment occurs just behind the cavity, at 2/3 chord length. The circulation
of the vortex is correctly positive and represents the so called re-entrant jet.
Finally, the cavity length, l, is normalized with respect to chord length, c, and plotted versus
σ/2α. In figure 4.3.9. the numerical results found for Re = 6.0e+05 are compared to
experimental and numerical results by Berntsen et al., 2001. For the case of Kubota, σ/2α is
4.3. The result fits fairly well with the experimental data, but is not as good as the results
obtained by Berntsen et al., 2001, for similar flow conditions. In addition, Berntsen was
modelling with Fluent 5, and thus did not benefit from the physically more correct model used
here, which predicts bubble collapse.
The same case has been then reproduced following the same steps and changing only the
Reynolds number to Re = 3 x 105. No clear convergence has been obtained. A detailed
tabular listing of all the simulations done in the present work concerning the reproduction of
the case of Kubota can be found in the annexes A.9-A.11.
As it can be seen in the graphs of the annexes, the cases solved using a Reynolds number of
3.0e+05 do not converge clearly. The density of water-vapour has been set back to the
values suggested by Kubota [ρ(vap) = 0.02558 kg/m3, η(vap) = 9.0e-06 kg/(ms)], while they
have been slightly higher for the case where Re = 6.0e+05. In both cases, the volume fraction
of non-condensable gas has been set for all cases to 1.5e-05. However, following Berntsen,
Numerical Simulation of Unsteady Turbulent Cavitating Flows Page 67
increasing the nuclei density in the numerical model may increase the dynamics of the
system, thus adding computational instability. In his work, he increased the nuclei density to a
number 20 times higher than the value suggested by Kubota, 1992. In future works, it is
suggested to vary this parameter in order to find more accurate results, although it is not clear
how far water impurities affecting the cavitating flow field are taken into account in the present
cavitation model (see also e.g. Kawakami, 2004 [20]). A further suggestion is also a further
mesh refinement. Some authors having researched on similar cases report that the
computations can be started with any initial condition, so that it is not believed that varying
parameters which affect these settings would yield better solutions.
The vapour pressure should neither be an error source as long as the cavitation number is
kept constant by changing respectively the inlet / operating pressure. Reproducing the case
of Berntsen et al., 2001 all acknowledgments from the presently elaborated case have been
taken into account, and other parameters have been changed to eliminate suspicions. This
case is presented in the following.
In order to have more reference, the case of Berntsen et al. has been simulated. The
mentioned work clearly defines almost all flow parameters, as well as geometry set-up and
computational configuration, so that we could identify possible error sources more easily. The
simulation by Berntsen, et al. was done with Fluent 5.
Some changes have been operated, where the most important was the definition of the angle
of attack, which no further was given by the corresponding velocity vector components, as it
did in our previous simulations (Kubota, 1992), but with a geometrically modified chord
orientation against the x-axis (horizontal).
The hydrofoil is also a NACA0015. The grid is a hybrid C-grid type with a total of about 22000
cells. The smallest cell size in the direction from the wall is 0.01 which corresponds perfectly
to the source case. The results in the first grid point are at y+ < 5, which implies the use of the
k-ε model with enhanced wall treatment. The inlet is set to 12m/s with constant profile, the
velocity being computed from a Reynolds number equal to 1.2e+06 based on chord length.
The unsteady SIMPLE algorithm was used, as well as second order discretisation schemes.
Internally, all numbers are stored and calculated with double precision. The time step used
was 0.0001 and the water vapour pressure was set to 2400Pa. The same effect could also be
Page 68 Memoria
achieved by having a constant reference pressure and adjusting the vapour pressure. The
cavitation number was adjusted by changing the reference pressure, Pref, in the solver, as it
was done for all previous simulations. The composite parameter was σ/2α = 5 at 8 degrees
AoA. A close-up of the grid and the computational domain including boundary conditions are
drafted in figure 4.3.10.
An instantaneous static pressure contour plot and velocity vector plot at the same instant in
time as obtained by Berntsen et al. is shown in figure 4.3.11.a. In the simulation by Berntsen
et al., we can observe (fig 4.3.11.b), the cavity length is identified from both, the static
pressure plot and the velocity vector plot, this one showing clearly the re-entrant jet and the
closure of the attached cavity.
a) b)
Fig. 4.3.10. a) Grid close-up surrounding the 2D hydrofoil, and b) outline of the
computational domain with boundary condition type
The very strong vortex structures downstream of the cavity give very low pressure regions
that effectively transport the bubbles shed from the foil. The authors can well predict the flow
field, especially the cavity length and the characteristics of cloud cavitation. The dynamics of
the cavitating system is still not very accurately computed.
Unfortunately we could not obtain comparable results. As for the case of Kubota, the cases
showed good results as long as the cavitation model was not used. The Appendix A.12
shows detailed information about the methodology used for simulating this case as well as
set-up parameters. While using the cavitation model when the local static pressure was lower
than the defined vaporisation pressure, divergence of the numerical simulation has always
been observed.
Numerical Simulation of Unsteady Turbulent Cavitating Flows Page 69
This case shows clearly a probable deficiency of the cavitation model or its implementation in
Fluent 6.1. The results obtained by Berntsen et al. may not be very realistic, but at least
according to the authors, the computations converge. Concerning their unrealistic results, the
authors argue, that the cavitation model “suffers from several simplifications”, where the most
critical would be the fact that both phases should be treated as compressible, condition which
is not fulfilled for the cavitation model implemented in Fluent.
a) b)
Conclusions
Within this work unsteady turbulent cavitating flows have been simulated. Due to the complex
character of such flow fields, a preliminary study of CFD and its possibilities was
indispensable. Also, an important bibliographic work has been done concerning previous
theoretical studies of cavitating flows and the different existing cavitation models.
Several cases of steady/unsteady turbulent or/and cavitating flows with the commercial CFD
code (Fluent v6.1) have been studied. The case of an unsteady and turbulent, non-cavitating
flow around a 2D circular cylinder was simulated as a first step using different turbulence
models at Reynolds numbers around the critical drag-crisis region. A mesh refinement control
helped choosing an adequate mesh for the computations. Compared with experimental data,
the results were quite divergent, but similar to previous numerical researches of other authors
(e.g. Cox, Zdravkovich). The non-accurate prediction of the transition of the viscous
underlayer from laminar to turbulent in the drag crisis region is mainly due to 3D effects. Also,
a deeper knowledge concerning the selection of a suitable step size and its influence on the
numerical results was obtained.
Further, the cavitation phenomenon has been studied in several applications. First, the full
cavitation model implemented in Fluent has been tested comparing findings with
corresponding experimental data. The first case was a steady, cavitating flow through a
sharp-edged orifice by Nurick, 1976. A previous mesh sensitivity study has been successfully
done. The inception of cavitation has been correctly predicted at σ = 1.76.
Finally, the unsteady turbulent flow around a 2D NACA 0015 hydrofoil has been simulated
using the cavitation model. This work was based on the publications by Kubota and Berntsen
et al., 2001 Results revealed that depending on the case, the cavitation model offers useful
results, but only in a qualitative way. Accurate fittings with experiments are obtained only in
few cases. A systematic and documented simulation process argues against any user’s error.
A software problem was mainly suspected, but also some problematic presumptions within
the cavitation model.
Numerical simulations of cavitating flows are very challenging since localized large variations
of density are present within a predominantly incompressible liquid medium. In its most
general form, the gas pockets causing these density variations can be highly compressible
with large temperature variations. However, even in the simpler case where the gas bubble
may be treated as incompressible, large density changes are still present at the liquid/gas
interface.
Page 72 Memoria
In summary, an interesting step for the understanding of the cavitation phenomena and its
coupling with vortex shedding has been done. The acquired physical background and the
experience in simulating complex flows builds an important basis for the accomplishment of
this task. Nevertheless, some further progress has to be done. In order to hopefully obtain
this progress, the steps described before have to be followed.
Future work
Many of the results presented here must be further explored. Firstly, it has to be determined if
the cavitation model implementation in Fluent 6.1 is not a possible error source, because
some inconsistencies in the obtained results point out that suspicion. A new version of the
software is presently being tested, where several changes have been done regarding the
cavitation model. The case of Nurick is being used for validation. Very promising first results
have been obtained at the simulation of these cases. Different meshes have been used, as
well as different turbulent models. A very important progress has been done concerning the
numerical stability of the software. As a further procedure, the latest version of Fluent has to
be tested using the more complex cases of Kubota and Bertnsen et al., treating of an
unsteady turbulent cavitating flow around a hydrofoil. Also, other numerical codes will be tried
out (e.g. CFX-5.7), using the cases of Kubota and Berntsen et al. for simulation.
Since at the present time RANS models are the unique real possibility in complex and
turbulent flow modelling, depending on the obtained results, a question of major importance
remains open: “Are RANS models the adequate tool to be used for the study of such complex
technical problems?” Probably and in most cases related to the cavitation phenomena, at
least big vortex structures have to be mathematically computed. The work of many authors
(e.g. the aforementioned works of Arndt, Kunz, Qin, Song) yield better results using the LES
method for the simulation of turbulent cavitating flows, but LES cannot be applied in more
complex cases due to high computer requirements.
Finally, after having obtained a deeper knowledge of the selected cavitation model, the major
goal which consists in the prediction of erosive effects of cavitation in components of fluid flow
machines has to be followed. The first step to reach this goal is the code validation in 3-D
airfoils using experiments from Escaler, 2001.
Numerical Simulation of Unsteady Turbulent Cavitating Flows Page 73
Acknowledgements
The achievement of this work has been done thanks to the bilateral work between the
Departament de Mecànica de Fluids of the Universidad Politècnica de Catalunya, Bartcelona,
Spain and the Institut für Strömungsmechanik und Hydraulische Strömungsmaschinen of the
University of Stuttgart, Germany. The academic exchange programme Erasmus was the
organisation which has made this interchange possible.
First and foremost, I would like to thank my supervisor, Mr. Miguel Coussirat for his invaluable
direction and insight, his patience, and genuine care and interest in this work and my
progress. I would also like to thank my professors and advisors in Germany, Dr. Albert
Ruprecht and Prof. Eberhard Göde who offered me the possibility to achieve my master
thesis in Barcelona and gave me a strong support for my work, helpful advice and
encouraging directives.
I am also very grateful to Prof. Eduard Egusquiza who made it possible for me to accomplish
my work at the Departament de Mecànica de Fluids in Barcelona.
I would also like to express my gratitude to all members of the Departament de Mecànica de
Fluids as well as to the staff of the Institut für Strömungsmechanik und Hydraulische
Strömungsmaschinen who offered me occasional support and guidance.
Lastly, I would like to thank my parents for their dedication and continual support; they very
much encouraged and helped me in the completion of this work.
Numerical Simulation of Unsteady Turbulent Cavitating Flows Page 75
Annexes
A.5 Drag coefficient for a circular cylinder at different Reynolds numbers (50000 nodes
grid = mesh 01), experimental data and numerical results.
A.6 Strouhal number for a circular cylinder at different Reynolds numbers (50000 nodes
grid = mesh 01), experimental data and numerical results.
A.7 Drag coefficient for a circular cylinder, experimental data and numerical results
,(J.S.Cox et al., 1997).
A.8 Strouhal number for a circular cylinder, experimental data and numerical results,
(J.S.Cox et al., 1997).
A.10 Data table, systematic numerical research, case of Kubota and adaptation from the
Fluent tutorial case.
A.11 Data table, systematic numerical research, case of Kubota at lower Reynolds number.
A.1
Numerical Simulation of Unsteady Turbulent Cavitating Flows Page 77
A.2
Page 78 Memoria
A.3
Numerical Simulation of Unsteady Turbulent Cavitating Flows Page 79
A.4
Page 80 Memoria
A.5
Numerical Simulation of Unsteady Turbulent Cavitating Flows Page 81
A.6
Page 82 Memoria
A.7
Numerical Simulation of Unsteady Turbulent Cavitating Flows Page 83
A.8
Page 84 Memoria
A.9
Numerical Simulation of Unsteady Turbulent Cavitating Flows Page 85
A.10
Page 86 Memoria
A.11
Numerical Simulation of Unsteady Turbulent Cavitating Flows Page 87
A.12
Numerical Simulation of Unsteady Turbulent Cavitating Flows Page 89
Bibliography
[3] Coussirat, M.G., Theoretical /Numeric Study of Flows with Strong Streamlines
Curvature, Universitat Politècnica de Catalunya, Barcelone, Spain, 2003.
[4] Hellsten, A., New Two-Equation Turbulence Model For Aerodynamics Applications,
Espoo 2004, Finland.
[5] Singhal, K.A., Mathematical Basis and Validation of the Full Cavitation Model, Journal
of Fluids Engineering, Sept. 2002.
[6] Cox, J.S., Rumsey, C.L., Brentner, K.S., Younis, B.A., Computation Of Sound
Generated By Viscous Flow Over A Circular Cylinder, NASA Technical Memorandum
110339 , 1997
[7] Zdravkovich, M.M., Flow Around Circular Cylinders, Oxford University Press, 2003.
[8] Ditlevsen, P.D., Cascades of Energy and Helicity in the GOY Shell Model of
Turbulence, 1996, arXiv: chao-dyn/9811009v1, 4 Nov. 1998.
[9] Nurick, W.H., Orifice Cavitation and its Effect on Spray Mixing, Journal of Fluids
Engineering, Dec. 1976.
[10] Kubota, A., Kato, H., Yamagushi, H., A New Modelling of Cavitation Flows: A
Numerical Study of Unsteady Cavitation on a Hydrofoil Section, Journal of Fluid
Mechanics, 1992.
[11] Berntsen, G.S., Kjeldsen, M., Arndt, R.E.A., Numerical Modeling of Sheet and Tip
Vortex Cavitation with Fluent 5, CAV2001: Fourth International Symposium on
Cavitation, session B5.006, California institute of Technology, Pasadena, CA USA,
June 20-23, 2001.
Page 90 Memoria
[12] Arndt, R.E.A., Song, C.C.S., Kjeldsen, M., Keller, A., Instability of Partial Cavitation: A
Numerical /Experimental Approach, Proceedings Twenty-Third Symposium on Naval
Hydrodynamics, Val de Reuil, France, 2000.
[13] Arndt, R.E.A., Kjeldsen, M, Song, C.C.S., Keller, A.P., Analysis of Cavitating Wake
Flows, Proceedings of the Hydraulic Machinery and systems 21st IAHR Symposium
Sept.9-12, Lausanne, Switzerland, 2002.
[14] Kunz, R.F., Venkateswaran, H.J., Lindau, C.L., Merkle, Preconditioning Algorithms for
the Computation of Multiphase Mixture Flows, 39th Aerospace Sciences Meeting &
Exhibit, 8-11 Jan., 2001.
[15] Song, C.C:S., Qin, Q., Numerical Simulation of Unsteady Cavitating Flows, CAV2001:
Fourth International Symposium on Cavitation, California institute of Technology,
Pasadena, CA USA, June 20-23, 2001.
[16] Song, C.C.S., He., J., Numerical Simulation of Cavitating Flows by Single-phase Flow
Approach, Third International Symposium on Cavitation, Grenoble, France, 1998.
[17] Kundu, P.K., Cohen, I.M., Fluid Mechanics, Second Edition, Academic Press, 2002.
[18] Vahedi Tafreshi, H., Pourdeyhimi, B., Cavitation and Hydraulic Flip, Nonvovens
Cooperative Research Center, North Carolina State University, Fluent NEWS, spring
2004.
[19] Qin, Q., Song, C.C.S., Arndt, R.E.A., A Numerical Study of an Unsteady Turbulent
Wake Behind a Cavitating Hydrofoil, CAV2003: Fifth International Symposium on
Cavitation, Osaka, Japan, November 1-4, 2003.
[20] Kawakami, D.T., Water Quality Effects on Cavitation, University of Minnesota, Dec.
2004.
Numerical Simulation of Unsteady Turbulent Cavitating Flows Page 91
Complementary works
[22] Sato, K., Saito, Y., Unstable Cavitation Behavior in a Circular-Cylindrical Orifice Flow,
CAV2001: Fourth International Symposium on Cavitation, California institute of
Technology, Pasadena, CA USA, June 20-23, 2001.
[23] Pozrikidis, C., Introduction to Theoretical and Computational Fluid Dynamics, Oxford
University Press, 1997.
[24] Senocak, I., Shyy, W., Numerical Simulation of Turbulent Flows with Sheet Cavitation,
CAV2001: Fourth International Symposium on Cavitation, California institute of
Technology, Pasadena, CA USA, June 20-23, 2001.
[25] Lindau, J.W., Kunz, R.F., Venkateswaran, S., Boger, D.A., Application of
Preconditioned, Multiple-species, Navier-Stokes Models to Cavitating Flows,
CAV2001: Fourth International Symposium on Cavitation, session B4.005, June 20-
23, 2001.
[26] Kunz, R.F., Boger., D.A., Chyczewski, T.S., Stinebring, D.R., Gibeling, H.J.,
Govindan, T.R., Multi-phase CFD Analysis of Natural and Ventilated Cavitation about
Submerged Bodies, Third ASME/JSME Joint Fluids Engineering Conference, San
Francisco, California, July 18-23, 1999.
[27] Kjeldsen, M., Svingen, B., Arndt, R.E.A., Transient in Closed Loop Hydraulic
Systems, 10th International Meeting of the Work Group on The Behaviour of Hydraulic
Machinery under Steady Oscillatory Conditions, Trondheim, Norway, 2001.
[28] Amromin, E., Hansberger., J., Wang, H., Wosnik., M., Arndt, R.E.A., Investigation of a
Low-drag, Partially Cavitating Hydrofoil, CAV2003: Fifth International Symposium on
Cavitation, Osaka, Japan, November 1-4, 2003.
[29] Wang, M., Catalano, P., Iaccarino, G., Prediction of High Reynolds Number Flow over
a Circular Cylinder Using LES with Wall Modelling, Center for Turbulence Research,
Annual Research Briefs, 2001.