|
Items in cart [1] |
Questions? Call 888-624-8373 |
| Copyright © 2009. National Academy of Sciences. All rights reserved. Terms of Use and Privacy Statement |
Below are the first 10 and last 10 pages of uncorrected machine-read text (when available) of this chapter, followed by the top 30 algorithmically extracted key phrases from the chapter as a whole.
Intended to provide our own search engines and external engines with highly rich, chapter-representative searchable text on the opening pages of each chapter.
Because it is UNCORRECTED material, please consider the following text as a useful but insufficient proxy for the authoritative book pages.
Do not use for reproduction, copying, pasting, or reading; exclusively for search engines.
OCR for page 868
24th Symposium on Naval Hydrodynamics
Fukuoka, JAPAN, 8-13 July 2002
Prediction of Vortex Cavitation Inception Using Coupled
Spherical and Non-Spherical Models and UnRANS
Computations
C.-T. Hsiao and G. L. Chahine (DYNAFLOW, INC., USA)
ABSTRACT
A spherical and a non-spherical bubble
dynamics model were developed to study cavitation
inception, scaling, and dynamics in a vortex flow.
The spherical model is a modified Rayleigh-Plesset
model modified to account for bubble slip velocity
and for non-uniform pressures around the bubble.
The non-spherical model is embedded in an Unsteady
Reynolds-Averaged Navier-Stokes code with
appropriate free surface boundary conditions and a
moving Chimera grid scheme around the bubble. The
effect of non-spherical deformation and bubble/flow
interaction on bubble dynamics is illustrated by
comparing spherical and non-spherical models. It is
shown that non-spherical deformations and
bubble/flow interactions are important for accurate
prediction of cavitation inception. The surface-
averaged pressure modified Rayleigh-Plesset scheme
is a significantly improvement over the conventional
spherical model and is able to capture the volume
changes of the bubble during its capture. It presents a
fast scheme to study scaling.
Scaling effects on cavitation inception is
preliminarily studied using two different Reynolds
number due to two different chord lengths. The
nuclei size effect on the prediction of cavitation
inception is also studied and its important effect
highlighted.
I. Intro(luction
Prediction and scaling of vortex cavitation
on propellers from model tests is of critical
importance to naval applications since this type of
cavitation often occurs first and affects the discretion
of the propeller. In order to do so it is essential to be
able to understand and model accurately the
mechanisms that enter into play in the cavitation
inception process of a vortex flow. An analysis of
the problem indicates that, in order to accurately
predict cavitation inception, one needs to consider the
bubble dynamics, the detailed viscous flow structure
in the low pressure regions of the propeller flow
field, and their interaction. Such an approach may
lead to the conclusion that some of the effects could
be neglected, but is needed before approximations
can be made satisfactorily. Numerical predictions
(Chahine 1990, 1994, 1995) and recent experimental
observations (Arndt and Maines 2000), indicates that
the cavitation inception starts with a spherical
nucleus drawn into the vortex core from the free
stream and which grows significantly. As the bubble
approaches the vortex core, non-spherical
deformations become significant. The bubble
deforms, elongates, and later breaks up into multiple
bubbles as it travels downstream and reaches the
vortex axis.
Numerical simulations of bubble capture by
a tip vortex (Ligneul and Latorre 1989, 1993, Hsiao
and Pauley 1999, Hsiao et al. 2000, 2002) were
primarily limited to a spherical bubble model which
assumes that the bubble does not interact with the
underlying flow field and sees a uniform pressure
around it given by the value of the liquid pressure at
the location of its center (in its absence) during its
capture. Since the traveling bubble is usually very
small relative to the vortex core size before
cavitation, the modification of the vortex flow may
be neglected. The assumption of a spherical model,
however, is not always appropriate because the
bubble may grow to a size large enough such that
non-spherical deformations occur when the bubble
approaches the vortex axis. A non-spherical bubble
dynamics model based on the boundary element
method was first developed by Chahine (1990, 1994,
1995) to study bubble deformations during capture
but neglected the effect the bubble has on the
underlying vortex flow field. To further account for
the bubble/vortex interaction, Hsiao and Chahine
(2001) developed a non-spherical bubble dynamics
model with full viscous interaction between the
bubble and the vortex based on unsteady Navier-
OCR for page 869
Stokes computations and on a moving Chimera grid
scheme. The bubble dynamics model has been
successfully validated for known cases and applied to
study bubble dynamics in a Rankine line vortex.
Their results showed that non-spherical bubble
deformations and bubble and flow interaction can
significantly influence bubble dynamics in the vortex
flow below cavitation inception. However, the
unsteady Navier-Stokes computations are time-
consuming and not practical for simulating the whole
bubble dynamics process and accounting for the
nucleus motion from upstream until drawn into the
vortex core. Recognizing this fact we apply the
spherical model to simulate the bubble capture and
then turn on the non-spherical model when the
bubble size exceeds a preset limit value to best take
advantage of each model.
In the current study, our model is applied to
predict tip vortex cavitation inception for a canonical
problem in which the tip vortex flow is generated by
a finite-span elliptic hydrofoil. The flow field is first
obtained by a steady-state Navier-Stokes computation
and provides the velocity and pressure fields for the
spherical model. To study the tip vortex cavitation
inception we assume that a bubble nucleus is released
upstream. The spherical model is first used to track
the bubble during its capture by the tip vortex. When
the non-spherical model is turn on, the flow field due
to the spherical bubble motion and volume change is
superimposed to the global steady-state blade flow
field solution to provide initial conditions for the
unsteady computations.
2. Numerical Methods
2.1 Navier-Stokes Computations
To best describe the tip vortex flow field
around a finite-span hydrofoil, the Reynolds-
Averaged Navier-Stokes (RANS) equations with a
turbulence model are solved because the RANS
computations have been shown to be successful in
addressing tip vortex flow (Hsiao and Pauley 1998,
1999) and general propeller, propulsor, ship, and
free surface flows (Taylor et al. 1998, Wilson, et al.
19981. The three-dimensional unsteady Reynolds-
Averaged incompressible continuity and Navier-
Stokes equations written in non-dimensional form
and Cartesian tensor notations are given as
~=0, (1)
axi
an +ujaxj=-axi+Re axj (2)
where ui = (u, v, w) are the Cartesian components of
the velocity, xi =(x,y,z) are the Cartesian
coordinates, p is the pressure, Re = pu*L*/,u is the
Reynolds number, u* and L* are the characteristic
velocity and length, p is the liquid density, and ,u is
its dynamic viscosity. The effective stress tensor Fit
. .
is given by:
ij =~˘axj + axi a- 3 it axkk -uiuj , (3)
where by is the Kronecker delta and u~u' is the
Reynolds stress tensor resulting from the Reynolds
averaging scheme.
To numerically simulate the tip vortex flow
filed around a finite-span hydrofoil, a body-fitted
curvilinear grid is generated. Equations (1) and (2)
are transformed into a general curvilinear coordinate
system. The transformation provides a computational
domain that is better suited for applying the spatial
differencing scheme and the boundary conditions. To
solve the transformed equations, we use and modify
the three-dimensional incompressible Navier-Stokes
flow solver, DF_UNCLE, derived from the code
UNCLE developed at Mississippi State University.
UNCLE was modified at DYNAFLOW to include
bubble, cavities and free surface effects (Hsiao and
Chahine, 2001, Chahine and Hsiao, 20021. The
DF_UNCLE code is based on the artificial
compressibility method (Chorin 1967) in which a
time derivative of the pressure multiplied by an
artificial-compressibility factor is added to the
continuity equation. As a consequence, a hyperbolic
system of equations is formed and can be solved
using a time marching scheme. This method can be
marched in pseudo-time to reach a steady-state
solution. To obtain a time-dependent solution, a
Newton iterative procedure is performed at each
physical time step in order to satisfy the continuity
equation.
The numerical scheme in DF_UNCLE uses a
finite volume formulation. First-order Euler implicit
differencing is applied to the time derivatives. The
spatial differencing of the convective terms uses the
flux-difference splitting scheme based on Roe's
method (Roe, 1981) and van Leer's MUSCL method
(van Leer, 1979) for obtaining the first-order and the
third-order fluxes respectively. A second-order
central differencing is used for the viscous terms
which are simplified using the thin-layer
approximation. The flux Jacobians required in an
implicit scheme are obtained numerically. The
resulting system of algebraic equations is solved
OCR for page 870
using the Discretized Newton Relaxation method
(Vanden and Whitfield, 1993) in which symmetric
block Gauss-Seidel sub-iterations are performed
before the solution is updated at each Newton
interaction. DF_UNCLE is also accompanied by a
Baldwin-Lomax algebraic turbulence model to model
the Reynolds stresses in Equation (31.
2.2 Spherical Model
When one assumes that the bubble remains
spherical during its volume variations, the bubble
dynamics can be described by the Rayleigh-Plesset
equation (Plesset 19481. The classical Rayleigh-
Plesset equation can be derived from the continuity
and momentum equations of a ID radial viscous
flow. In the current study, this equation is modified to
account for a slip velocity between the bubble and the
host liquid, and to account for the non-uniform
pressure field along the bubble surface. The bubble is
passively convected through the known "basic" flow
field of the blade or propeller. The bubble's travel
velocity is generally different from the liquid velocity
due to the bubble dynamics and to the presence of
local flow field pressure gradients. This results in an
added pressure term due to the velocity difference
between the liquid and the bubble. This term is
obtained by adding the velocity potential due to the
flow past a sphere to the velocity potential due to a
spherically oscillating bubble and using this in the
pressure balance equation. As a result, an additional
term of slip velocity, (u-ub)214, where u is the
liquid convection velocity and ub is the bubble
travel velocity, is added to the classical Rayleigh-
Plesset equation as:
RR+3R2 1
2 p
PA + Pgo ( R ) —PenCounter
By 4,uR
R R
(u—us ~ (4)
R is the time varying bubble radius, Ro is the initial
or reference bubble radius, ~ is the surface tension
parameter, pgo is the initial or reference gas
pressure inside the bubble, and Pv is the vapor
pressure. The time varying gas pressure inside the
bubble, pp. is related to Ego through a polytropic
compression law:
pgR = P9O~ ~ (5)
with the polytropic gas constant k. k = 1 for
isothermal behavior is used in the present study.
Pencounter is the ambient pressure "felt" by
the bubble during its travel. In the classical
Rayleigh-plessset equation, Pencounter is the
pressure that exists in the liquid, in absence of the
bubble, at the location of the bubble center. This
definition is quite acceptable if the pressure field
does not vary significantly around the bubble. In the
case of bubble capture in a vortex, such a definition
used by previous investigators, leads to unlimited
bubble growth once the bubble reaches a constant
intensity line vortex axis. Hsiao et al. (2000, 2002)
introduced the Surface Averaged Pressure (SAP)
Rayleigh Plesset equation, in which Pencounter is
defined as the average of all liquid pressures at the
various parts of the bubble surface. This concept
enabled more realistic modeling of the bubble
dynamics in line vortices. We will use here the SAP
approach and compare it with the fully coupled 3D
bubble dynamics.
To describe the bubble trajectory during
capture by the tip vortex, the motion equation
described by Johnson and Hsieh (1966) is used:
dt p 4 CD (U—ub ) | u—Ub|
+ 3 (u rub ) R.
(6)
where the drag coefficient CD is given by an
empirical equation such as that of Haberman and
Morton (19531:
24 (1+0 197ROb63 +2.6xlO Reb ) ~ (7)
eb
where the bubble Reynolds number is defined as
Reb = ~ (8)
To determine the bubble motion and its volume
variation, a Runge-Kutta fourth-order scheme is used
to integrate Equations (4) and (6) through time. The
liquid velocity and pressures ~ Pencounter ~ are
obtained directly from the RANS computations. The
numerical solution of the RANS equations, however,
offers the solution directly only at the grid points. To
obtain the values for any specified location (x,y, z)
on the bubble we need to interpolate from the
background grid. To do so, an interpolation stencil
(8)
OCR for page 871
and interpolation coefficients at any specified
location need to be determined at each time step.
To determine the interpolation stencil from
the background grid, we implemented a three-
dimensional point-locating scheme based on the fact
that the coordinates (x,y,z~of the bubble location
are uniquely represented relative to the eight corner
points of the background grid stencil:
8 8 8
x = 2, NiXi, y = ~ Ni Yi ~ Z = 2, NiZi ~ (9)
i=1 i=1 i=1
where
N1 = (1 - V)~1 -~1 - By, N2 = ˘~1 -~1 -are ~ 10y
N3 =~1-~1-~), N4 =~1 Hi,
Ns = (1 - ˘~1 -yip, N6 = ˘~1 -yip,
N7 = (1-oJyr(p, N8 = YELP
˘,yr,~ are the interpolation coefficients, and
(xi,yi,zi ) are the coordinates of the eight corner
points of a grid stencil in the background grid.
Equation (9) is solved using a Newton-Raphson
method. For a bubble point to be inside the grid
stencil requires that the corresponding ˘,yt,~
satisfy 0 < ~ < 1, 0 < jr < 1, 0 < ~ < 1. Once the
interpolation stencil and interpolation coefficients are
determined, the pressure and velocities can be
obtained by using Equation (9~.
2.3 Non-Spherical Model
To fully account for the interaction between
the bubble and the flow filed and for non-spherical
deformations, a non-spherical bubble model
embedded in the UNRANS computations with
appropriate free surface boundary conditions and a
moving Chimera grid scheme has been developed by
Hsiao and Chahine (2001~.
2.3.1 Free Surface Boundary Condition
To best describe the bubble surface
behavior, a general free surface boundary condition
satisfying both kinematic and dynamic boundary
conditions is applied. The kinematic boundary
condition is the Lagrangian condition that a particle
on the surface must remain on the surface. For a free
surface of equation Ff A, t) = 0, this can be written:
DF
Dt
-=n (
The general free surface dynamic boundary
condition is the condition of zero shear stress and
balance of normal stresses at the bubble liquid
interface. Here, the stress due to the gas inside the
bubble is neglected. With the same simplifications
used by Batchelor (1967) for deriving the dynamic
boundary condition in the Cartesian coordinate
system, Hodges et al. (1996) derived a dynamic
boundary condition in a curvilinear coordinate
system by requiring the grid to be normal to the
boundary. Following their work the current study
implements the dynamic boundary condition in non-
dimensional form as
OU ~ ~~ OW ~2 OW
a,; =-g33~g a,; +g a
at _g ( 22OW+ l23W) (13)
We ~ 14)
with gij=d~i obj. gij= axi axj (15)
axk axk ask ask
(12)
P Pgv Re a,;
Off ~e
where e is the curvature, We =pu*2 L*ly is the
Weber number, ~ is the surface tension, and
Pgv = (Pg +Pv - poO)/pu*2 . (16)
To determine the gas pressure we assume that the
amount of non-condensable gas inside the bubble
remains constant and that the gas satisfies the
polytropic relation pg~k=constant~where ~ is the
bubble volume.
2.3.2 Moving Chimera Grid Scheme
The Chimera grid scheme is a grid
embedding technique which provides a conceptually
simple method for domain decomposition. In this
approach, structured sub-grids generated around
each component in the flow field or over sub-
domains of complex geometries are put together
without requiring the mesh boundaries to join with
the global grid in any special way. In the present
study, a body-fitted sub-grid is created around the
bubble and overset with a global grid which is
generated for the hydrofoil. Figure 1 illustrates the
details of a Chimera grid system in a two-
dimensional domain. The global grid consists of
three different types of points: "regular" points,
"overlap" points and "hole" points. The sub-grid
consists of two different types of points; "regular"
points and "overlap" points. The Navier-Stokes
equations are solved separately for the global grid
and for each sub-grid at all "regular" points. The
communication between these two grids is made by
OCR for page 872
interpolating flow variables at the "overlap" points.
Data from the global grid are used in the
interpolation to supply outer boundary conditions to
the sub-grid. Hence the effect of the main flow field
is properly imposed on the bubble. For the global
grid solution to account for the bubble, points in the
global grid are blanked out to form a "hole" within
some neighborhood of the bubble. On the fringes of
this blanked-out region, data from the sub-grid
solution are used for interpolation to provide interior
boundary conditions for the global grid. This way
the presence of the bubble and its effect are
transferred to the main flow field.
The implementation of the Chimera grid
scheme includes two tasks: 1) identify the "hole"
points and the "overlap" points, and 2) determine the
interpolation stencil and interpolation coefficients.
The "hole" is here defined by a "hole boundary"
which consists of the grid surface at a constant value
of (. A global grid point is considered to be a
"hole" point if it is inside the "hole boundary". A
grid point is considered to be inside the "hole" if the
dot product between the vector from the closest
point at the "hole boundary" to the grid point and the
normal vector at the "hole boundary" at the closest
point is negative or zero. The overlap points in the
global grid are then determined based on the fringes
of the "hole". Here, two layers of grid stencil are
identified as the "overlap" points in order to enable
implementation of a third order boundary condition.
The method for determining the interpolation stencil
and the interpolation coefficients from the
background grid for each overlap point is the same
as that described in section 1.2.
In addition to the two tasks described
above, a blanking technique is implemented in the
Navier-Stokes solver to account for "holes" and
interior boundaries and to account for boundary data
transferred from interpolated data sets rather than
from the usual boundary condition routines. In the
iterative solution routines, a variable ibLi,j,k) is
introduced to blank out the "hole" and "overlap"
points such that no update of the variables takes
place at these points:
ib~i k'_J1, if apointis not blanked
'J' )0, if a point is blanked
In DF_UNCLE, a system of linear equation is written
to solve for the flow variable difference, /` Q. between
iterations and to update the flow variables by
Qn+~=Qn+AQ. To ensure that AQ is zero at blanked
points, we multiply the matrix coefficients of the left
hand side and the vector of right hand by ib end
assign the diagonal terms of the matrix to be one for
Figure 1. Illustration in 2D of a Chimera grid
system. Localization of the "hole" points and of the
"overlap " points marked as · .
2.3.3 Initial conditions
It is important to specify appropriate initial
conditions for the unsteady Navier-Stokes
computations when the non-spherical model is turn
OCR for page 873
on. In the current study, the initial condition is
constructed by superimposing the local perturbation
flow due to the bubble dynamics and the global tip
vortex flow obtained by the steady-state Navier-
Stokes computations. The combined velocity field,
Hi =(u,v,w), for the initial condition is, therefore,
obtained by
R2
Ui=Ui+ 3 R(Xi-Xio) i=1,3, (17)
where r is the distance from the bubble center
x~O = (xO'yO,zO) to the field point xi. The
combined pressure field, p, is prescribed as:
R( ~ ~ RR2 [1 (R )3~
(18)
where pO is the liquid pressure at the bubble center
and
_1
PR = PgV - R R - W ' (19)
is the liquid pressure at the bubble wall. It is noted that
the variables in Equations (17~-~19) are all non-
dimensional.
2.4 Computational Domain and Grid Generation
To compute the flow around the finite-span
elliptic hydrofoil we generate an H-H type grid with
12 blocks for a computational domain which has all
far-field boundaries located six (6) chord lengths
away from the hydrofoil surface. For the unsteady
computations, we generate an O type grid on the
bubble surface (41 x 21~. Also, an O-O type 3D
domain grid ~ 41 x 21 x 25 ~ is generated to extend the
grid from the bubble surface to the outer boundary
that we locate at25R,~s. Rns is the bubble radius
when the non-spherical model is turn on.
Since the bubble moves and its surface deforms
during the numerical simulations, an efficient grid
generation scheme that can automatically generate an
appropriate grid based on the moving boundary at
each time step must be integrated with the Navier-
Stokes solver. We apply a grid generation scheme
combining both algebraic and elliptic grid generation
techniques as described by Hsiao and Pauley (19961.
This grid generation scheme is selected because the
algebraic grid generation technique is suited to create
grid clustering and boundary-orthogonal grids at the
bubble surface. This is important for resolving the
flow field near the bubble surface and for applying
appropriate free surface boundary conditions. The
elliptic grid generation technique is then applied to
smooth out the rough algebraic grid.
With the current grid generation, the first grid
spacing from the bubble surface is kept as a constant
distance (= 0.05Rns) during the computations. To
avoid uneven distribution of the surface grid points
due to bubble deformation, a re-gridding scheme is
applied every few time steps at the bubble surface to
restore the grid distribution Cording to the initial
grid-stretching factor.
3. Results and Discussion
3.1 Single Phase Flow
To illustrate the method we consider an
elliptic hydrofoil having a NACA 16020 cross-
section with an aspect ratio equal to 3 (based on
semi-span). For the steady-state computations, a 12-
block grid with a total of 2.7 million grid points was
generated. To optimize the grid resolution for the tip
vortex, the grid clustering is adjusted using repeated
computations to align the clustering around the
location of the tip vortex centerline. The final
refined grid has at least 16 grid points in the
spanwise direction and 26 grid points in the
crosswise direction within the vortex core. This
results in the smallest grid cell in the vortex center
having a size of about 2xlO~ CO, where CO is the
hydrofoil chord length at the base section. The
streamwise grid resolution near the tip region is kept
almost constant end of a spacing of 2xlO 3 CO. The
first grid in the normal direction to the foil surface is
located 1 x 10 s CO away from the solid surface.
NACA16020 oc=1 2°
~'
.
' :// —Re=1A4x10'
;' ~ Re=2.88x106
_ .'
0.1 0.15 0.2 0.25
X
Figure 2. Pressure coefficient variations along the
NACA16020 elliptic foilfor two values of the
Reynolds number.
OCR for page 874
The flow field at an angle of attack equal to
12° was computed for various Reynolds numbers.
The computed pressure coefficients along the tip
vortex centerline for the two Reynolds numbers,
Re =1.44xlO6and Re =2.88x106 are shown in
Figure 2. It is seen that the locations of the
minimum pressure for both these cases are very
close to the hydrofoil tip at x/CO=0.05andO.04
respectively. The corresponding minima m pressure
coefficients are -Cpmin =2.27 and 2.73. It is a
common, even though not always accurate (Hsiao et
al. 2000, 2002) to equate the cavitation inception
number with the negative value of Cpmin. The ratio
between the two values of Cpmin is usually
presented as the ratio of the two Reynolds numbers
to some power ~ usually taken to be 0.4. There are
recent indications that this power should itself vary
with the Reynolds number decaying at the larger
values (Shen et al. 2001a). Here we see that the ratio
corresponds to a value of a =0.27. Further
improvements including grid refinement and higher
order turbulent modeling is needed before we draw
conclusions for that matter.
3.2 Bubble Dynamics with the Spherical Model
For computations of bubble dynamics using
the spherical model, all water properties are defined
at 20° C and nuclei present in the liquid upstream are
convected in the tip vortex flow. As described in the
previous studies of Hs~ao and Pauley (1999) and
Manes and Arndt (1999), only nuclei which pass
through a restricted space volume ("window of
opportunity") can be drown into the vortex center and
grow explosively. This "window of opportunity" can
be determined by releasing nuclei upstream and
tracking their trajectory to see if they encounter the
minimum pressure in the tip vortex. In the spherical
computations shown here all bubble were released
fromx/ CO = -0.1. The release locations in the y and
z directions were chosen such that the bubbles
encounter the minimum pressure in the tip vortex.
Figure 3 shows a typical spherical bubble trajectory.
The figure also illustrates the bubble size at the
various locations along its trajectory. Figure 4 shows
the bubble behavior during its capture by the tip
vortex using the SAP spherical model. The case
shown is for Rip =50pm, c; =2.5 and
Re = 2.88x106 based on a flow speed,
up = 2.88 m/s and a chord length, CO =1 m. The
pressure encountered by the bubble, Pe72Counter' is
shown as a function of time. Also shown is the
history of the bubble radius and of the acoustic
pressure generated by the bubble at the location
x/Co=O,y/CO=Qz/CO=0.3. It is seen that the
bubble grows continuously during its capture and
reaches its maximum size after it passes the
minimum pressure location. The bubble then starts to
collapse as it sees an increasing encounter pressure
and executes strong volume oscillations which lead to
a strong acoustic signal
Definition of the cavitation inception can be
based on either the bubble size exceeding a limit
value or on the emitted acoustic pressure exceeding a
threshold. We will apply the former definition in the
present study.
~ Cavi~ng bubble:
i: ~::~: ~ :~
... . .... ~
Figure 3. Bubble trajectory during its capture and
resulting bubble size from actual computation
0.002s
0.002
,~.0015 ~ 5\
0.001
0.0005
Re=2.88x1 o6 RO=50pm c,=2~
0.04 0.06 0.08 0.1
Time (see)
13000
12000
1 1 000
1 0 0002
_
9000 ~
a)
8000 cut
7000 Is
6000
5000
4000 U]
3000
2000
Figure 4. Example computation of bubble dynamics
for encountered pressure, bubble radius and emitted
acoustic pressure versus time during its capture in
the tip vortex of a NACA 16020 foil using the SAP
spherical model.
OCR for page 875
3.3 Bubble Dynamics for Non-Spherical Model
To include non-spherical bubble
deformations and bubble/vortex interaction effects on
bubble dynamics, the non-spherical bubble model is
turn on when the bubble size exceeds a preset value.
For consistency with the numerical scheme, this limit
size is set to correspond to the smallest grid size in
the tip vortex core region (Rns=2xlO 4Co). To
reduce the unsteady Navier-Stokes computational
burden, a computational sub-domain which contains
51x41x51 grid points is extracted from the overall
computational domain as shown in Figure ~ This
sub-domain is chosen large enough so that the flow
quantities at the boundary of the domain are not
influenced by the bubble dynamics.
Figure 5. A view of the sub-domain used for the
unsteady RANS computations with bubble
deformation.
An so type grid is then created for the
bubble overset grid for the unsteady Navier-Stokes
computations. To start these computations, the
velocity and pressure field obtained from Equations
(17) and (18) are applied as initial conditions. The
initial values at all boundaries of the sub-domain are
maintained over time during unsteady computations
except for the solid boundary at which no-slip
condition is enforced. Figure 6 shows the initial
pressure contours on three selected grid planes (two
of them for the global grid aid one for the overset
grid) with the overlap points blanked out. It is seen
that the flow field is only locally altered by the
presence of the bubble.
Figure 6. Init~alpressure contours on three
selected grid planes (two of them for the global gad
and one for the overset grid) with the overlap points
blanked out.
Figure 7 illustrates the bubble shapes
several time steps obtained for Ro = 50pm, cr = 2.5
and Re =2.88x106 while the bubble is traveling
along the tip vortex. It is seen that the bubble
elongates in the axial direction, then as in the
spherical model, the bubble starts to collapse after
reaching its maximum size. The unsteady 3D
computations, however, fail so far to continue once
strong deformations develop over the bubble surface
during the collapse. Comparisons of the bubble
radius versus time for the spherical models (the
conventional and the SAP model) and the non-
spherical model are shown in Figure 8. For the non-
spherical model we have used the equivalent bubble
radius based on the bubble volume. It is seen that the
conventional spherical model very significantly over-
OCR for page 876
predicts the maximum bubble size since the pressure,
P. applied in that spherical model is the pressure
existing at the bubble center. Such a model obviously
does not account for pressure variations around the
bubble surface. Our surface averaged pressure (SAP)
spherical modelenables a more realistic evaluation of
the bubble dynamics. Here, the average is obtained
using the six polar points at the bubble surface, the
solution of our spherical model agrees very well with
the nonspherical model as shown in Figure 8.
1
tax
my::
Figure 7. Bubble shapes at several time steps
obtainedfor Ro = 50 pm, ~ = 2.5 arid
Re =2.88x106.
R=2.88xl o6 RO=50.um cs=2.5
n nn7
n nnF
0.005
0.004
0.003
0.002
0.001
Sphertcal Bubble NoSAP
- S phedcal Bubble With SAP
Non-Sphelical Model
I .... I.... ~ .A .
0.0 4 0. 05 0 .06 0.07 0.0 8 0 .09 0.1
Tune (see)
Figure 8. Comparisons of the bubble radius versus
time for the spherical models (the conventional and
the SAP model) and the 3D UnRANS computations.
3.4 Scaling Effect on Cavitation Inception
To study scaling effects on cavitation inception, the
bubble dynamics models were applied to predict
cavitation inception in two different scales,C0=lm
and 0.5m. All computations were conducted with the
same free stream velocity, u<= = 2.88 m/s, and initial
bubble radius, Ro = 50pm .
As mentioned in Section 3.2, the bubble size
is used for defining cavitation inception. Here, the
bubble size is represented by the bubble radius for the
spherical model and by the equivalent radius for the
non-spherical model. A large series of computations
was conducted with the spherical model (since it is
very fast to run) for a range of cavitation number and
is shown in Figures 9 and 10. For selected cavitation
numbers computations were also conducted with the
3D UnRANS non-spherical model. Both
conventional and SAP spherical models were used.
Figure 9 and 10 shows the maximum size the bubble
reaches at each cavitation number studied for
Re = 1 .44x 1 o6 and Re = 2.88 x 1 o6 respectively. It
is seen that in each case the bubble experiences an
explosive growth when the cavitation number is
smaller than a limit cavitation inception value. The
predictions of the 3D non-spherical model (at least in
terms of the equivalent radius) are very close to those
obtained by the SAP spherical model. Averaging the
encounter pressures around the bubble surface
appears to be a justified approach. Close to the
cavitation inception number, a sudden change in the
maximum bubble size is seen, in which case direct
punctual comparison between the three models could
give the impression of very different results.
To determine the cavitation inception from
the current result, we must define the cavitation
inception event first. If the cavitation inception event
is defined based on the sudden change in the
maximum bubble size, the cavitation inception
numbers for both scales are read as 2.13 and 2.58
respectively. This indicates that the cavitation
inception number varies between two scales as the
ratio of the Reynolds numbers to the same power as
the ratio of the two Cpmin. In other words,
~2/~1 ~ CPmin2/CPmin1 = (Re2/Re1 )0 28 is similar
to the prediction of single phase flow. However, if
the cavitation inception event is based on when the
maximum bubble size reaches some value, say 2mm,
then the prediction of the cavitation inception number
may be quite different, e.g. Hi ~ Re036 if the
criterion for cavitation inception event is defined for
RmaX > 2mm.
OCR for page 877
3.5 Nuclei Size Effect on Cavitation Inception
The nuclei size is already known to have a
profound effect on the prediction of tip vortex
cavitation inception (Hsiao and Pauley, 1999 and
Shen et al., 2001b. Hsiao et al. 2000, 2002~. Previous
numerical studies, however, relied on a spherical
model. Here, nuclei size effect is also studied using
our non-spher~cal model. 3D computations were
conducted for two different nuclei sizes,
Ro = 20pm arid 50,um, at the same Reynolds
number, Re = 2.~8x106 .
Comparison between the bubble behavior
with these two initial nuclei sizes are also shown
using the curves maximum bubble size versus
cavitation number. It is seen from Figure 11 that the
smaller nuclei requires a much smaller value of the
cavitation number for cavitation inception. Then, for
values of cavitation number lower than the inception
limit, a more abrupt change in the RmaX versus
s curve is seen. However, as already known, we find
again that the maximum bubble size reached by the
bubble is independent of the initial nuclei size for
cavitation number much smaller than the limit
cavitation inception value.
As we did for Rip =50,um, we compare in
Figure 11 the SAP spherical model computations and
the 3D non-spherical computations for Rid = 20pm. It
is seen that the SAP spherical model over-predicts a
little the maximum bubble growth size at low
cavitation inception numbers giving larger bubble
sizes. This may reflect that using only 6 polar points
may not be very accurate for large deformations and
bubble sizes.
4. Conclusions
A bubble dynamics model combining a
spherical model and a non-spherical model embedded
in UnRANS computations was developed to predict
cavitation inception in complex flow configurations
and was applied here for tip vortex cavitation
inception on a finite-span hydrofoil.
Comparisons between the classical spherical
model and our current model showed that the non-
spherical deformation and bubble/flow interaction
can be very important. A Surface Averaged Pressure
(SAP) scheme can be applied to improve
significantly the spherical model when predicting the
cavitation inception.
Scaling effects were demonstrated by
comparing the results for two different scales. It is
again found as in Hsiao et al. (2000, 2002) that the
definition of cavitation inception or the means for
detecting inception, can significantly affect the
predicted cavitation inception number.
Comparison between two different nuclei
sizes showed that smaller nuclei sizes result in
smaller inception values, and that the maximum size
at low 6 is independent of the initial nuclei size.
0.01 ~
nmn
0.008
~ ~7
o06
.005
.004
0.003
ao02
0.001
Re=1 A4x106 RO=5Olim CO=O~m
~ \ Spherical M odel N 0 SAP
- \ — Spherical M odel Wile SAP
~ Non-Spheri:al M odel
O -I ~ ~ , I . ~ ~ ~ I A___. . . 1
.9 2 2.1 2.2 23 2.4
Cavitation no.
Figure 9. Comparison of the maximum bubble
radius size versus cavitation number between
spherical (conventional and SAP) and non-
spherical models for the small scale,
CO = 0.5m, Re = 1.44 xl 06, Rg = 50,um.
0.0125
0.01
E 0.0075
c:
0.005
0.0025
Re=2.88:~106 RO=50,um CO=lm
\ ~ -- Spherical Mod b No SAP
\— Sphedcal Mod b Wffl SAP
~~ Non Spherical Model
2.3 2.4 2.5 2.6
Cavitation no.
27 2.8
Figure 10. Comparison of the maximum bubble
radius size versus cavitation number between
spherical (conventional and SAP) and non-
spherical models for the large scale,
CO = I m, Re = 2.8 8x Jo6 ~ Ro = 50,um.
OCR for page 878
Re=2. 88x1 o6 Cn=1 m
0.005~'
0.004
X 0.003
a:
0.002
0.001
~ - - R. =20 Lun SDhertal ModelWlthSAP
_
-
_ ~
_ .— ~ I · · __~
2.2 2.3 2.4 2.5 2.6 2.7
Cavitation No.
u
Ro-50 rum Spherical Model With S AP
Ro_20 run Nm-SphericalMod~
\ ~ Ro=SO run Na~-Spherical Modd
'\ \
Figure 11. Comparison of the maximum bubble
radius size versus cavitation number between
spherical (conventional and SAP) and non-
spherical models for two different nuclei sizes,
Ro =20 arm 50pm at Re = 2.88x 106.
ACKNOWLEDGMNETS
This work was conducted at DYNAFLOW,
INC. (www.dvnaflow-inc.com) and has been
supported by the Office of Naval Research under
contract No. N0014-99-C-0369 monitored by Dr. Ki-
Han Kim. This support is greatly appreciated.
REFERENCES
t1] Arndt, R.E and Maines, "Nucleation and Bubble
Dynamics in Vortex Flows," ASME Journal of
Fluids Engineering, Vol. 122, 2000, pp. 48
493.
[2] Batchelor, G. K., Fluid Dynamics, Cambridge
University Press, 1967.
[3] Chahine, G.L., "Nonspherical Bubble Dynamics
in a Line Vortex," ASME Cavitation and
Multiphase Flow Forun~ Toronto, Canada, Vol.
98,Junel990,pp.121-126.
[4] Chahine, G.L, "Bubble Dynamics and
Cavitation Inception in Non-Uniform Flow
Fields" 20th Svmposium on Naval
Hvdrodvnamics, Santa Barbara, California,
August 1996, pp. 290-311.
t5] Chahine, G. L., "Bubble Interaction with
Vortices," Vortex Flow, Chapter 18, Ed. S.
Green, Kluwer Academic, 1995.
[6] Chahine, G.L., Hsiao, C.-T., " Numerical
Investigation of Sheet and Cloud Cavitation
Inception and Dynamics on Propeller Blades,"
Rep. 98003- 1, 2002, Dynaflow, Inc., Jessup,
MD.
[7J Chorin, A. J., "A Numerical Method for Solving
Incompressible Viscous Flow Problems,"
Journal of Computational Physics, Vol. 2, 1967,
pp. 12-26.
t8] Haberman, W.L., Morton, R.K., "An
Experimental Investigation of the Drag and
Shape of Air Bubbles Rising in Various
Liquids," Report 802, 1953, DTMB.
[9] Hodges, B., Street, R. Zang, Y., "A Method for
Simulation of Viscous, Nonlinear, Free-Surface
Flows," 20th Svmposium on Naval
Hydrodynamics, 1996,pp. 791-809.
t103 Hsiao, C.-T., "Numerical Study of the Tip
Vortex Flow Over a Finite-Span Hydrofoil, "
Ph.D. Thesis, Department of Mechanical
Engineering, The Pennsylvania State
University, Adviser L.L. Pauley, 1996
[11] Hsiao, C.-T., Pauley, L.L., "Numerical Study of
the Steady-State Tip Vortex Flow over a Finite-
Span Hydrofoil," ASME Journal of Fluid
Engineering, Vol. 120, 1998, pp. 345-349.
[12] Hsiao, C.-T., Pauley, L.L., "Numerical
calculation of Tip Vortex Flow Generated by a
Marine Propeller," ASME Journal of Fluid
Engineering, Vol. 121, 1999a, pp. 63~645.
t13] Hsiao, C.-T., Pauley, L.L., "Study of Tip Vortex
Cavitation Inception Using Navier-Stokes
Computation and Bubble Dynamics model,"
ASME Journal of Fluid Engineering, Vol. 121,
l999b, pp. 198-204.
[14] Hsiao, C.-T., Chahine, G.L., Liu, H.L., "Scaling
Effects on Bubble Dynamics in a Tip Vortex
Flow: Prediction of Cavitation Inception and
Noise," Rep. 98007-lNSWC, 2000, Dynaflow,
Inc., Jessup, MD.
[1SJ Hsiao, C.-T., Chahine, G.L., Liu, H.L., "Scaling
Effects on Prediction of Cavitation Inception in
a Line Vortex Flow," to appear in ASME
Journal of Fluids Engineering, 2002.
[16] Hsiao, C.-T., Chahine, G.L., "Numerical
Simulation of Bubble Dynamics in a Vortex
Flow Using Navier-Stokes Computations and
Moving Chimera Grid Scheme," 4t International
Svmposium on Cavitation, Pasadena, CA, June
20-23, 2001.
t17] Johnson, V.E., Hsieh, T., "The Influence of the
Trajectories of Gas Nuclei on Cavitation
Inception," Sixth Svmposium on Naval
Hydrodynamics, 1966, pp. 163-179.
OCR for page 879
[18] Ligneul, P., Latorre, R., "Study on the Capture
and Noise of Spherical Nuclei in the Presence of
the Tip Vortex of Hydrofoils and Propellers,"
Acustica. Vol. 68, 1989, pp. 1-14.
El91 Ligneul, P., Latorre, R., "Study of Nuclei
Distribution and Vortex Diffustion Influence on
Nuclei Capture by a Tip Vortex and Nuclei
Capture Noise," ASME Journal of Fluid
Engineering, Vol. 115, 1993, pp. 50~507.
[20] Plesset, M. S., "Dynamics of Cavitation
Bubbles," Journal of Applied Mechanics, Vol.
16, 1948, pp. 228-23 1.
t21] Roe, P. L., "Approximate Riemann Solvers,
Parameter Vectors, and Difference Schemes,:
Journal of Computational Physics, Vol. 43, 1981,
pp. 357-372.
t22] Shen, Y.T., Jessup, S.D., Cowing, S., "Tip
Vortex Inception Logarithmic Scaling,"
NSWCCD-50-TR-2001/060, October 2001a,
Naval Surface Warfare Center, Carderock
Division, West Bethesda, MD.
[231 Shen, Y.T., Chahine, G.L., Hsiao, C.-T., Jessup,
S.D., "Effects of Model Size and Free Stream
Nuclei on Tip Vortex Cavitation Inception
Scaling," 4th International Symposium on
Cavitation, Pasadena, CA, June 2~23, 2001b.
[24] Taylor, L.K., Pankajakshan, R., Jiang, M.,
Sheng, C., Briley, W.R., Whitfield. D.L.,
Davoudzadeh, F., Boger, D.A., Gibeling, H.J.,
Gorski, J., Haussling, H., Coleman, R., and
Buley G., "Large-Scale Simulations for
Maneuvering Submarines and Propulsors,",
AIAA Paper 98-2930, 1998.
[25]van Leer, B., "Towards the Ultimate
Conservative Difference Scheme. V. A Second
Order Sequel to Godunov's Method," Journal of
Computational Physics, Vol. 32, 1979, pp. 101-
136.
[26jVanden, K., Whitfield, D. L., "Direct and
Iterative Algorithms for the Three-Dimensional
Euler Equations," AIAA-93-3378, 1993.
[27] Wilson, R., Paterson, E., and Stern, F.
"Unsteady RANS CFD Method for Naval
Combatant in Waves", 22n~ Svmposium on
Naval Hvdrodvnamics, Washington, D.C.,
1998, pp. 183-197.
OCR for page 880
DISCUSSION
S. Cordier
Bassin d'Essais des Carenes, France
The results shown of a novel modelization of
bubble dynamics are very interesting The
importance of nuclei size and density in
cavitation (surface and vortex) is well put in
evidence. This effect has been demonstrated in
experiments in the GTH and reported in the
literature (Fruman and Briancon). Do you have
any plans to validate your results against this
type of data?
AUTHORS' REPLY
We are well aware of the work you mentioned.
The results shown in our paper are focused on
single bubble dynamics. Validation of our
current scheme requires a well-controlled
experiment in which the initial nucleus
conditions are known and the nucleus is
followed in its motion. To study the effect of
nuclei size and density on the scaling problem, a
statistical type model which allows one to
randomly distribute the nuclei according to a
known nuclei size distribution is required. Since
this will require a large number of simulations
with a significant nuclei population, the SAP
spherical model, described in the paper, seems to
be a good choice for such a study that we are
presently conducting. We will be presenting such
results in the future and will try to validate our
results with the experimental data you suggested
assuming you could provide us with the nuclei
size distribution in those tests.
DISCUSSION
R. Gornstein
Navatek Ship Ltd., USA
Can your method be extended to surface flows
(sheet cavitation) as well as vortex flows?
AUTHORS' REPLY
Actually, our numerical scheme was initially
developed to solve the sheet cavitation problem
(Chahine and Hsiao 2000~. We then extended it
by incorporating a moving Chimera grid scheme
to enable the simulation of the dynamics of
moving bubbles in a tip vortex flow. The same
scheme, of course, can be used to further
improve the simulation of sheet cavitation.
However, the current scheme is limited when the
interface becomes multi-connected. We are
working on further improving this aspect.
Representative terms from entire chapter:
tip vortex
