|
Items in cart [0] |
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 616
An Unsteady Three-Dimensional Euler Solver Coupled
with a Cavitating Propeller Analysis Method
].-K. Choi, S. Kinnas (The University of Texas at Austin, U.S.A.)
ABSTRACT
A fully three-dimensional unsteady Euler solver,
based on a finite volume scheme and the pressure
correction method, is developed and applied to the
prediction of the effective wake for propellers sub-
ject to non-axisymmetric inflows. The propeller is
modeled via unsteady body force terms in the Euler
equations. First, the Euler solver is validated against
the analytical solution of actuator disk theory, and
the solution in the case of uniform inflow. Then, the
unsteady effective wake is predicted for a three-cycle
axial inflow wake, and it is found that the unsteadi-
ness in the effective wake is small in this case. Finally,
the predicted cavity shape and the unsteady velocity
field are compared with the experimental data mea-
sured in cavitating propeller experiments which have
been performed at the MIT water tunnel.
1 INTRODUCTION
The effective wake is a very crucial issue during the
design and assessment of propeller performance, es-
pecially in the presence of blade cavitation, because
it is well known that the predicted cavity extent and
volume, as well as the magnitude of the predicted
pressure pulses depend strongly on the estimated ef-
fective wake inflow.
There has been experimental (Huang et al 1976,
Huang & Cox 1977) and theoretical work (Huang &
Groves 1980, Shih 1988) for the prediction of the ef-
fective wake by employing steady axisymmetric Eu-
ler equations. Later, effective wake prediction meth-
ods using Reynolds Averaged Navier-Stokes(RANS)
equations were developed (Stern et al 1988b, Stern
et al 1988a, Kerwin et al 1994, Kerwin et al 1997) for
axisymmetric flow applications and (Stern et al 1994)
for non-axisymmetric applications. In the above
methods, the propeller is represented by the body
force term in RANS equations. The method using
the body force to represent the blades was also used
in turbomachinery applications (Damle et al 1997),
where the body force was determined by iteratively
adjusting it until the velocity tangency condition on
the rotor and stator blades was satisfied. A different
approach which does not need the propeller poten-
tial solution is to solve the RANS equations directly
with the no-slip boundary condition on the propeller
blade surface. (Chen et al 1994, Stanier 1998)
In a first effort toward the better prediction of
the effective wake harmonics, the authors have al-
ready developed a steady three-dimensional Euler
solver and successfully applied it to predict the steady
non-axisymmetric effective wake for given propeller
geometry and inflow nominal wake (Choi & Kinnas
1998a, Kinnas et al 19994. In that method, the pro-
peller effect, which appears through the body force
terms in the Euler equations, is averaged with time.
Nevertheless, the propeller loading is allowed to vary
in space with strength which depends on the loading
at each blade angle.
In the present work, a fully three-dimensional
unsteady Euler solver, based on a finite volume ap-
proach and the pressure correction method, is de-
veloped and applied to the prediction of the un-
steady effective wake for propellers subject to non-
axisymmetric inflows. The proposed unsteady anal-
ysis, which is the unsteady extension of our previous
work (Choi & Kinnas 1998a), captures the truly un-
steady behavior of the interaction between the vortic-
ity in the inflow and the propeller action by adopting
the body force distribution which varies in both space
and time.
2 PRESENT METHOD
The unsteady flow field around a propeller can be de-
composed into two parts; one rotational and one ir-
rotational field. In the present method, the propeller
induced velocity field, qp, is the irrotational velocity
field. The propeller induced velocity, qp, can be ex-
pressed in terms of the perturbation potential, Up, by
the following relation.
qp(~c, t) = V>p(~c, t)
The propeller vortex lattice method (VLM), named
MPUF-3A, is used to solve for the perturbation po-
tential on the propeller as shown on the left side
of Figure 1 (Kinnas et al 1998, Kinnas et al 19994.
OCR for page 617
Blade pressure p(x,t)
or Circulation F(r,O,t)
Propeller VLM
Vop(x,t) is determined from the propeller vortex lattice
method solver (with respect to rotating system)
1 ~
;;
Vortical Flow FVM
q~(i,t) is determined from the finite volume
method Euler equation solver
~ ~ ~ r I ~ ~ ~ I ~ 1 or or - i _
Dqt / Dt = -Vp/p ~ f
Figure 1: Interaction between the potential flow vor-
tex lattice method (VLM) and the Euler equation
based finite volume method (FVM).
The perturbation potential as defined by equation
(1) with respect to a coordinate system which ro-
tates with the blade at a constant angular velocity,
Q. should satisfy the Laplace equation. Note that
the propeller rotation, Q. has only an axial compo-
nent which is negative for right-handed propellers.
The boundary condition on the blade mean camber
r
surface IS
`g = — (qe + Q x A) ri (2)
where, r' is the unit vector normal to the mean cam-
ber surface, and qe is the effective wake as defined in
the next paragraph. In addition, a Kutta condition
must be applied at the trailing edge of the blades
with the assumed trailing wake geometry.
Then, the rotational part can be expressed by
subtracting this irrotational field, qp, from the total
velocity field, A. The rotational part determined in
this way is defined as the effective wake velocity, qe,
given as follows.
qe my, t) — At my, t) —qp (0c, t)
qt(0C't)—V>p(0C't) (3)
The effective wake, qe, is used as the inflow, when
solving for the propeller perturbation potential, 8p.
The total velocity, A, iS evaluated from the Euler
solver, as described in the next paragraphs. Once the
potential he is found, two things can be calculated;
(a) the propeller loading as a function of space and
time, and (b) the propeller induced velocity qp, or
equivalently V¢p. The propeller loading is needed
to determine the body force distribution in the Euler
solver, and the propeller induced velocity is needed
when the effective wake is calculated. If unsteady
sheet cavitation exists on the blade, then the extent
and the thickness of the cavities are also determined
as part of the solution.
On the other hand, the global flow field is ro-
tational and is solved by the finite volume method
Euler solver, as shown on the right side of Figure 1.
The propeller loading obtained from the VLM is con-
verted into the body force f, and this body force is
distributed over the cells corresponding to the blade
location. When the Euler equations,
Dqt _Vp + if (4)
are solved, the total velocity qt iS obtained. The effec-
tive wake can then be computed by equation (3) US-
ing the propeller induced velocity, qp, already known
from the propeller VLM. The overall process is iter-
ative, between the potential flow solution of the pro-
peller and the Euler equation solution of the global
flow field, as shown in Figure 1.
The main assumption behind the current ap-
proach is that viscosity plays a less significant role
in the interaction between the rotational inflow and
the propeller action, even though it plays a major
role in developing the nominal wake. This is a valid
assumption in the vicinity of the propeller where the
propeller-inflow interaction is dominated mainly by
inviscid vorticity dynamics. In order to predict the
nominal wake, a high Reynolds number viscous flow
around the ship hull must be solved, and this some-
times has been found to be a very difficult task (Lars-
son 19974. Instead, one can use the measured nom-
inal wake so that the "real fluid effects" (including
the effects of viscosity) are included in the nominal
wake data. The interaction between the vorticity in
the nominal wake and the action of the propeller can
then be captured by solving the Euler equations, and
the effective wake can be calculated from the differ-
ence between the total flow field and the flow field
induced by the propeller.
In the current algorithm, all components of the
inflow must be specified at the upstream boundary of
the domain in which the Euler equations are solved.
This boundary must be placed at a sufficient dis-
tance upstream of the propeller plane so that the
specified inflow is not affected by the presence of the
propeller. Ideally, the nominal wake inflow (in the
absence of the propeller) must be measured at the
upstream boundary of the calculation. However, of-
ten the nominal wake is measured at the propeller
plane, and thus the inflow at the upstream bound-
ary needs to be deduced from the measured. This
situation was also encountered in the effective wake
analysis of (Huang & Groves 19804. In some cases,
OCR for page 618
Ship
~ propeller plane
Figure 2: Ship-fixed Cartesian coordinate system.
a reverse flow calculation is attempted in order to
find the flow at the upstream boundary which will
produce the measured nominal wake at the propeller
plane. It is concluded that placing the nominal wake
at the upstream boundary is the easiest approach,
especially when the vertical component of the inflow
is small compared to the axial component. One more
variation in the current method is that the effective
wake is evaluated not at the propeller plane but at
a plane just upstream of the propeller, due to the
difficulty in computing the induced velocity near the
propeller panels. Luckily, the location where the ef-
fective wake is evaluated has been found to have in-
significant effect on the prediction (Choi & Kinnas
1998a, Choi 20004.
3 FORMULATION
A ship-fixed Cartesian coordinate system, as shown
in Figure 2, is used in the present formulation. The
origin of the coordinate system is at the center of the
propeller disk, with the -axis pointing downstream
along the propeller shaft axis. In the coordinate sys-
tem, the positive y points vertically upward, and the
positive z-axis port side (left) of the ship. Also the
angle ~ is measured around the -axis starting from
the y-axis, as shown in Figure 2.
The variables are made non-dimensional by a
reference length which is chosen as the propeller ra-
dius R. and a reference velocity equal to the velocity
at infinity UOO (ship speed).
, A A Ai
ha, v, w
v_ (n,v,w) = U
00
A A A
a_ (fx, fy'fz)
f off ~ fY' fZ) U2 /R
(i,j,k+l)/ T
51 ~
i--` ~
W' ~
,' S
,'
,
1. a' ,
(i,j,k) B
(i,j+l,k+l) N (i+l,j+l,k+l)
Q. '''\\ ~ 17
6 /
E ,3
/(i+l,j+l,k)
(i+l,j,k)
Figure 3: Cell (i,j,Ic) on which the finite volume
method (FVM) is applied.
t=
A
t
R/Uoo'
p
A
= P 2 (8)
where, the hat(^) means a dimensional variable, and
A
f is the body force per unit mass.
3.] Steacly Euler Solver
In the steady Euler solver, named WAKEFF-3D, the
method of artificial compressibility (Chorin 1967) is
applied. In the artificial compressibility method,
pseudo-unsteady terms are introduced on the left
hand side of the steady incompressible Euler equa-
tion to mimic the solution procedure of the unsteady
compressible Euler equations.
After introducing the pseudo-unsteady terms,
the Euler equations can be written in the follow-
ing three-dimensional Cartesian form using non-
dimensional variables.
bU OF dG OH
0~* + 0~ + fly + 'pz = Q (9)
where, t* is a pseudo-time which can be regarded just
as an iteration parameter. The terms U. F. G. H.
and Q are defined as follows.
U=
13
v , F= ov , G=
w ow
w/)
H= vw , Q
w2 +p
vl)
ov
v2+p ~
vw
o
fx
fy
fz
(10)
(5)
where, ~ is the artificial compressibility factor which
<6' is a controllable constant. If steady state is reached,
the first term in equation (9) vanishes and the incom-
pressible steady flow solution is obtained. For the
(7) rest of this section, the superscript, *, on the time
OCR for page 619
"/~ G
~dex(:
Figure 4: Arrangement of the nodes and cells in the
present method.
variable, t, will be omitted with the understanding
that it is the pseudo-time.
According to the finite volume method (FVM),
the Euler equations (9) are integrated over a cell of
volume V. Then, with the use of the divergence
theorem, the integral form of the Euler equations
is applied to each cells as shown in Figure 3. The
semi-discrete equation in which only the space is dis-
cretized can be written as follows for each cell.
Ouija
fit Vile +
f=N,W,S,r,T,~
+Gf (Adolf + Hf (Adolf ~ = Fiji Vim
where, Vim is the volume of the cell (i,j,~. The
areas (Ax,Ay,Az~f are the projections of the area
of the face, f, along the (~,y,z) directions, respec-
tively. Note that the vertex based scheme, in which
all variables in U are stored at the nodes, is used in
the present method.
For the time discretization, Ni's Lax-Wendroff
method (Ni 1982) is applied. That is, the variable U
at node, "5" (shown in Figure 4), at the next time
(pseudo-time) step, r' + 1, is approximated by the
following second order formula.
U5+l ~ U5 + ( 0t ) At+ ( 0~2 ) ~ 2) (12)
where, At is the time step size, and the superscript
r' represents the value of the current time step. The
first order term, bU/0t, and the second order term,
02U/0t2, are approximated by the average of their
values at the neighboring cells, where the values at
each cell are expressed using the discretized govern-
ing equation (11) itself. Then, the resultant equation
becomes the distribution formula, in which the con-
tributions from each neighboring cell are grouped to-
gether. The contributions from the surrounding cells
to the node "5", (5U5)c' are given in (Choi 2000)
with more details of the formulation.
The artificial dissipation (or viscosity) is applied
to improve the stability of the numerical method (An-
/
z / i-index
Figure 5: Finite volume method applied to the cylin-
drical grid in the Cartesian coordinate system.
derson 19954. The second and fourth order dissipa-
tions, respectively p2 and p4, are scaled by At and
added to the distribution formula.
U'l+l = U5 + ~ (5U5)c (13)
c=A,B,C, ~ AH
+ At (~2—p4)
with
iFi(Ax if aim p2 = ~2 (5iiU5 +5jjU5 +~kU5) (~4)
~4 = ~4 (5iiiiU5 + ~jjjj U5 + 5~p U5 ~ (15)
where, bid, etc. are the finite central difference op-
erators. For example, the fourth order i-directional
difference operator, biiii, is defined as follows.
~iiii = ~ ji-2 - 4( jigs + 6( ji—4( hi+ + ~ ji+2 (16)
The artificial dissipation coefficients, ~2 and ~4, are
user specified constant parameters that control the
amount of the dissipation. For most of the cases pre-
sented in this paper, the second order dissipation co-
efficient, ~2, iS zero. The time step size should be
determined so that the CFL condition (Courant et al
1967) is satisfied. The iteration in time (pseudo-time)
is continued until the maximum change of variables
is less than a certain tolerance.
Since a cylindrical grid system fits the propeller
problem naturally, a cylindrical grid system is used
in discretizing the fluid domain, as shown in Figure
5, while a Cartesian coordinate system is used to de-
scribe the space variables and the velocity vectors, as
well as the Euler equations. The k-index of the grid
is increasing in the B-direction and the surface cor-
responding to the last k-index and that of the first
k-index are the same surface. On this surface, a pe-
riodic boundary condition is applied.
3.2 Unsteacly Euler Solver
In the unsteady Euler solver, named WAKEFF-3U,
the finite volume method is applied only for the dis-
OCR for page 620
cretization of the three unsteady momentum equa-
tions, while the continuity equation is handled by the
pressure correction method (Patankar 1980, Rhie &
Chow 1983).
3.2.1 Finite Volume Discretization for Mo-
mentum Equations
In the unsteady formulation, the following dimen-
sionless incompressible Euler (momentum) equations
in differential form are used.
bU + OF + FIG + OH = Q (17)
fit 0~ by Liz
where, the column matrices U. F. G. H and Q. each
of which has only three elements, are defined as fol-
lows.
_ ,u2 'led
= v , F= ov , G= v2 ,
w ow vw
ow -UP/ + fx
H = vw , Q = —~P/0Y + fy
w2 -0P/0Z + fz
Note that the pressure gradient terms are grouped
in Q with the body force terms because the pressure
correction method is adopted in the unsteady for-
mulation. From this point, the formulation and the
numerical implementation of equation (17) are very
similar to those of the steady case, described in the
previous section, and are not repeated here.
3.2.2 Pressure Correction Method for Con-
tinuity Equation
In the present method, the idea of the SIMPLE
method (Patankar 1980) is followed. As opposed to
the staggered grid system of the original SIMPLE
method, the collocation scheme, in which all vari-
ables are defined and stored at nodes, is applied.
To achieve this, the pressure weighted interpolation
(Rhie & Chow 1983) is applied when the pressure
gradient is discretized.
Using equation (13) and the definition of U by
equation (18), the velocity at the next time step at
node "5" can be expressed as follows.
The intermediate velocity field denoted by the
superscript, *, is defined as the one based on an inter-
mediate (or guessed at first) pressure field, p*. Since
the method is iterative, to find the correct (within
a certain tolerance) velocity and pressure fields for
time step r' + 1, the intermediate fields are updated
repeatedly. Thus, the intermediate velocity and pres-
sure are not marked with the time step index, under-
standing that they are between the time steps r' and
{1~ + 1.
U5 = U5 + ~ (5U5)c +~\t (pl2 - pl4) (20)
c=A,B,C, ~ AH
where, (5U5)C is evaluated by using the intermediate
pressure field, p*. Also, the artificial dissipations,
[l2 and [l4, are evaluated based on the intermediate
velocity field, U*.
The pressure and velocity corrections, which are
represented by prime ('), are defined as follows.
P5 = P5 +P5 (21)
U'5l+l = U5 + U5 (22)
where, p5 and U5 are the corrections to the pressure
and the velocity fields, respectively, so that U50+1
satisfies the continuity equation. According to the
derivations given in (Choi 2000), the velocity correc-
tion, U5, is expressed by the pressure correction and
the pressure correction p5 can be found by solving a
discretized equation for the pressure correction. The
pressure correction equation can be solved iteratively
using the following explicit expression for p5 with the
under-relaxation, CIpp.
(P5)rlew = (1—CIpp) (pools (23)
+ orpp —<[90 - (9lP1 + 96P6 + 98Ps
+gsEg + 9loPlo + 9llPll))
The optimum value of the under-relaxation, asp, is
found to be 0.5 in the present method. In summary,
the following operations are performed within each
time step of the present unsteady method.
Sequence of Operations
- ~u~+1 - 1. Guess the pressure field p*.
U50+1 _ v50+1 = US + ~ (5U5)c~ (19)
wr7~+1 c=A,B,C, ,H 2. Solve the momentum equation by equation (20),
+~\t (~2—~4 ~ to obtain intermediate velocities u*, v* and w*.
Note that the artificial dissipation is based on the
velocities at time step no, and that the terms (5U5)c~
include the residual, Rc, and the source, Qc, of the
surrounding cells.
3. Solve for the pressure correction, p', by itera-
tively applying the equation (234. The optimum
number of this inner iteration is found later to
be IppIter=2.
OCR for page 621
4. Calculate the (correct) pressure p0+i using equa-
tion (21), and the (correct) velocities u0+i, v0+i
and w0+i using equations (224.
5. Treat the corrected pressure p0+i as a new in-
termediate pressure, p*, return to step 2, and
repeat until either the residual of the continuity
equation, QM, is smaller than a given tolerance
or a prescribed number of iterations, IpcMax, is
reached.
Later, the optimum number of pressure correction
per each time step is determined to be IpcMax=10.
3.3 Boundary Conclitions
The six boundary conditions applied in the present
methods are summarized next.
· Upstream boundary:
(AVOW) = (~~V,w~giverl, girl = 0 (24)
· Downstream boundary:
Aid, v, w, p) = 0 (25)
· Axisymmetric shaft or hub boundary:
The free slip condition is applied by using the
tangential velocity component, uc.
—u' Or — —~ ~ — = — (26)
or Ilr or rH
Age, uc) ~ _ nx
where, the unit normal vector, Into, D~r)' pointing
into the fluid region is defined to be positive.
· Center line boundary (j=1 line in the grid):
i Nk ~ (~' v, w, P)i=2 (27)
~=l,Nk
where, No is the number of nodes in the circum-
ferential direction.
· Outer boundary (far field):
~ = 1, ~ ~ = 0, p = 0 (28)
· Periodic boundary:
(11, V, W. p)~=l = (0, V, W. p)~=Nk (29)
where, the k-index is shown in Figure 5.
By
W
~ body surface
~ .
E
Figure 6: A cell of the finite volume method (FVM)
mesh that contains the body surface.
4 BODY FORCE
To find an appropriate distribution of body force
which represents the propeller (or hydrofoil in two-
dimensions) in the Euler solver, two-dimensional lift-
ing body problems are discussed first.
4. ] Two-Dimensional Hotly Force
For simplicity a steady incompressible flow around a
two-dimensional hydrofoil is considered for the mo-
ment. Figure 6 shows the part of the body surface
intersected by a cell of the FVM mesh. Pressure p is
acting on the part of the body surface of length A, and
the cell has a size of /\~ x /\y. From the dimensionless
Euler equations for two-dimensional, incompressible,
steady flows, the following expression for the body
force can be derived for the cell shown in Figure 6.
f _ HE HW + ANON —USES PE—PW
~— A~ BY
, 2 ,,2,, HEVE —HWVW + P ~ PS (31)
cent
where, the subscripts N. S. W. and E stand for the
four edges of the cell at north, south, west, and east,
respectively. Note that the first two terms of the
discretized equations (30) and (31) are the forces due
to the momentum convection.
The last terms of the discretized equations above
can be interpreted that the body force at a cell should
be equivalent to the integrated pressure force acting
in the normal direction over the part of the body sur-
face contained in the same cell. When the cell shown
in Figure 6 is assumed to be sufficiently small, the
body force can be expressed as follows by equating
the two forces; the total force due to the pressure,
pAri and the integration of the body force over the
cell, /\~\y f.
f = A,,, ~ (32)
If the body surface is vertical, the body force simply
becomes fit = pa/\, and if the body surface is hor-
izontal, the body force becomes fy = p/~\y. These
OCR for page 622
1
no
- 0 .5 ~
-1
NACA 0015, a=0.8, f/c=0.05, a=5°
~ ~ ~ . .
. ................ _ ~ I .............. 1
~ Potential solution ~ ~ ~ ~ ~ ~ ~
............ ~ ~ 2 ~ ~ ~ ~ 2 .~. >' ), ~ ~
. ~ ~ ~ '= = =~ .~ ~ ~ ~ y,
NACA 0015, a=0.8, f/c=0.05, a=5°; p=0.01, 62=0.0, 64=0.005, At=0.1 CFL, 200x200 cells
— ,
....... ,, ,,, ,~ \ ................... .~
I ~ ~ , I ' ~ ~ I ~ , '! 1 1 I:::::: 1 /~i~l~
-0.5 0 0.5 1 -1 -0.5 0 0.5
X X
Figure 7: Stream lines and pressure contours of the
potential flow solution; NACA0015 thickness, NACA
a=0.8 camberline, 5% maximum camber, 5° angle of
attack.
0.6
0.5
0.4
0.3
0.2
n1
n
-0 1
-0.3
-0.4
X
Figure 8: Body force cells and the body force vectors
that represent the foil according to the surface distri-
bution model. (Arrows represent the force acting on
the fluid at each of 292 body force cells.)
expressions for the In, y)-directional body forces are
consistent with the last terms of the discretized equa-
tions (30) and (31) under the condition that the pres-
sure inside the body is assumed to be zero.
4.1.1 Surface Distribution Model
In the surface distribution model, the body force is
applied to the cells that contain a part of the bound-
ary element method (BEM) body surface element.
The first and second terms of the equations (30) and
(31) are calculated by using the local FVM velocity
(n, v), so that they can exactly cancel the convective
terms on the left hand side of the momentum equa-
tions, when the body forces are substituted into the
...........
pi
_ 035
_ 025
_ 015
~~b
~ 055
Figure 9: Stream lines and pressure contours of
the Euler solution (surface distribution model);
NACA0015 thickness, NACA a=0.8 camberline, 5%
maximum camber, 5° angle of attack.
momentum equations. On the other hand, the last
terms are computed from the surface pressure known
from the BEM solution by applying equation (32)
and by setting zero pressure at the nodes inside the
body.
~ fin ~
+
HE—u2W ANON —USES
~~ LAY
AN —AS + SIEVE —WOW
by ~~ FVM
Pate
fry
PLAY
. ~~ my REM
(33)
When this body force model is applied, the momen-
tum equations provide a direct balance between the
surface pressure force and the cell body force. In
other words, this method is equivalent to specifying
the pressure of the cell as the value obtained from
the boundary element method (BEM) solution at the
same location.
The surface distribution model is applied to a
two-dimensional foil made of a NACA0015 thickness
and a NACA a=0.8 camberline with 5% maximum
camber. The geometry of the foil is represented by
100 boundary elements in the BEM solution. In Fig-
ure 7, the potential flow solution around the foil is
shown for uniform inflow with an angle of attack of
5°. The potential flow solution is obtained from the
two-dimensional perturbation potential method de-
scribed and fully verified in (Choi & Kinnas 1998b).
The flow field very near the surface of the foil is not
accurate because the evaluation of the induced veloc-
ity at locations which are closer to the BEM element
OCR for page 623
than the size of the element, is inaccurate due to the
constant strength singularity approximation. How-
ever, note that the pressures on the surface control
points which are used in the calculation of the FVM
body force are accurate because they are obtained di-
rectly, by differentiating the solution from the BEM
along the surface of the body.
In Figure 8, the 292 body force cells and the
body force vectors at each cell, that represent the
foil, are shown. It is observed in the figure that the
largest part of the suction side is pulling the fluid
while the largest part of the pressure side is pushing
the fluid.
In Figure 9, the solution of the Euler equations
is shown, as obtained with the body force distri-
bution shown in Figure 8. The formulation of the
two-dimensional Euler solver is similar to that of the
three-dimensional problem explained earlier, when
the z coordinate is omitted. The definition of the
artificial dissipation is slightly different only for this
specific case. That is, the last term in equation (13)
should be replaced with /\t(,u2 - ,u44/Sij, where Sij
is the average of the areas of the cells neighboring
the node (i, j). For this computation, the square do-
main extending —8 < ~ < 8 and —8 < y < 8 is
used with the boundary conditions specifying (n, v, p)
at the upstream, top, and bottom boundaries and
0~n,v,p)/~ = 0 at the downstream boundary. The
pressure field and the streamlines obtained from the
Euler solver are very close to those of the potential
flow solution shown in Figure 7. For the most part
of the field, excluding regions very near the foil sur-
face and the wake surface, the difference between the
two velocity fields is observed to be less than 1% of
the inflow velocity. It is an important fact that the
flow field, including the effect of the foil thickness, is
well simulated with the surface distribution model.
Also, note that the FVM computation required about
5 hours of CPU time on a DEC Alpha Station 600
5/266i, while the BEM computation took only 0.1
second on the same machine.
4.1.2 Camberline Pressure Model
It is natural that very thin foils should be represented
by the body force cells along the mean camber surface
of the foil. This can be achieved easily by just relocat-
ing the pressure gradient part of the body force deter-
mined by equation (32) for the cells on the BEM sur-
face to the cells on the mean camber surface, resulting
into two cell thickness distribution of body forces. A
simpler version of this method is to have only one
cell thickness distribution along the camberline by
adding the two force vectors from the upper and the
iDEC Alpha Station 600 5/266 is approximately equivalent
to a 500 MHz Pentium PC.
lower surface of the foil. Numerical tests showed that
both body force cell arrangements produced practi-
cally the same velocity fields with only small local
differences (Choi 20004. The one cell thickness dis-
tribution has the advantage that there is no need to
set zero pressure inside the foil because the addition
of the upper and the lower pressures automatically
implies the assumption that the two pressures are
based on a certain common reference value, say zero,
and there is no need to set this explicitly. It can also
be shown that the finite volume discretization of a
cell that contains both surfaces of the foil arrives at
the same conclusion.
The above mentioned body force model can be
made even simpler by neglecting the convective terms
which used to be calculated from the local FVM so-
lution, and keeping only the pressure gradient term,
as shown next.
(Pyrex )+ + (Pyrex ~
_~ 2\y ,~ (34)
~~ my BEM
ifxl=) .
where, the superscript + and - mean the values from
the upper surface and the lower surface of the BEM
boundary, respectively. This body force model which
is simpler to apply is named as the camberline pres-
sore model. It has been confirmed by numerical tests
that the original camberline distribution model (with
both FVM and BEM terms) and the simpler camber-
line pressure model (with only BEM terms) do not
make any significant difference in the predicted ve-
locity at the locations corresponding to the effective
wake evaluation plane in propeller applications.
To test the camberline pressure model, a rela-
tively thin hydrofoil geometry which would be sim-
ilar to a section of a propeller at outer radii is se-
lected. A two-dimensional foil section constructed
from NACA0005 thickness and NACA a=0.8 cam-
berline with 2% maximum camber is used in this test.
The potential flow BEM solution around the foil,
which is represented by 100 elements, is shown in Fig-
ure 10 when uniform inflow with 5° angle of attack
is incident. The potential flow solution is obtained
again from the two-dimensional perturbation poten-
tial method described in (Choi & Kinnas 1998b).
Again, the flow field very near the surface of the foil is
not accurate because of the inaccuracy of the induced
velocity calculation routine.
The body force distribution is shown in Figure
11, where only 10 body force cells are used to rep-
resent the foil. Note that the starting point of each
force vector is the lower left corner of the cell which
the body force belongs to. The body force has smaller
magnitude than that from Figure 8 because what
OCR for page 624
NACA 0005, a=0.8, f/c=0.02, a=5°; p=0.08, O2=2.0, O4=1 .0, At=0.5 CFL, 1 50x150 cells
1
0.5 t
-n .h
............. ,.,;
_ ~ Potential solution ~ -o t
.,.,.,.~ ~
.,., ~ ................................
............... _''
................ ,.,.,.,~,,,,,,, ~~ ~ . """""'l"
\ in_
_ 0 35 0.5
~025
, ~ _ 015
._ A' _ .
..............................................................................................
025 ~ 0
_ '-055 -0.5
~1' ~ ~ .~F, , ..
-0.5 0 0.5 1
X
Figure 10: Stream lines and pressure contours of the
potential flow solution; NACA0005 thickness, NACA
a=0.8 camberline, 2% maximum camber, 5° angle of
attack.
our
0.5
0.4
0.3~
n ~ L
._
~ 0.1 _
O _
-0.1 _
-0.2 _
n .~ _
Body Force Cells and Body Force Vector (fx, fy)
| Pressure term only formula l
10
X
Figure 11: Body force cells and the body force vectors
that represent the foil according to the camberline
pressure model, equation (344; grid with 50x50 cells
(10 body force cells).
is plotted is the body force per unit mass and the
cell area where the force is to be distributed becomes
larger now.
In Figure 12, the solution of the Euler equations
is shown, which is obtained with the body force dis-
tribution shown in Figure 11. The formulation of the
two-dimensional Euler solver is similar to that of the
three-dimensional problem including the way the ar-
tificial dissipation is defined. The domain size and
the boundary conditions are the same to those ap-
plied in the previous case. A relatively coarse grid
with 50x50 cells is used because the number of body
force cells in this grid is similar to the number of
cells along the chord of a blade in the propeller flow
problems. The jagged pressure contours near the foil
NACA 0005, a=0.8, f/c=0.02, a=5°; p=0.08, a2=2.0, a4=1.0, At=0.5 CFE, 50x50 cells
.............. ,~: ~
A O''"'' ' .. _
<4~ If
Or....
~_~
.................... ~~ I\
. ~~w
I, I ' I, I I, ' i1 1 1 ~ 1
-1 -0.5 0 0.5 1
X
A
n :2~
0 25
0 15
0 05
-0 05
-0 15
-0 25
-0 35
-0 55
-0 65
-0 75
Figure 12: Stream lines and pressure contours
of the Euler solution (camberline pressure model);
NACA0005 thickness, NACA a=0.8 camberline, 2%
maximum camber, 5° angle of attack.
could be smoothed out by using increased artificial
dissipation. However, the present result with the
moderate artificial dissipation is preferred because
the present artificial dissipation seems to be adequate
in the region of interest, i.e. not too close to the cam-
berline. The pressure field and the streamlines ob-
tained from the Euler solver are similar to those of
the potential flow solution shown in Figure 10, even
though the foil thickness effect is completely lost in
the Euler solution.
The difference of the flow fields computed by the
potential solver and the Euler solver is shown in Fig-
ure 13. The difference of the two solutions is com-
puted by subtracting the potential flow field (Figure
10) from the Euler solution flow field (Figure 124.
For the most part of the field excluding very near the
foil surface and the wake surface, the difference be-
tween the two solutions is less than 1% of the inflow
velocity.
The calculation of the effective wake requires the
evaluation of velocities on a plane at about 0.3 pro-
peller radius upstream of the propeller plane. If a
section of a propeller blade at around 70% of the
propeller radius is considered, such plane would cor-
respond to a straight line as the dashed line shown
in Figure 13, which is sloped relative to the nose-tail
line of the foil section. Along this line, the velocity
magnitudes of the Euler solution and the potential
solution are compared in Figure 14. Along this line,
the maximum change of the velocity is about 7% of
the inflow velocity, and the maximum error of the Eu-
ler solution is about 1%. The computing time with
this 50x50 cell mesh is less than 3 minutes on a DEC
Alpha Station 600 5/266.
OCR for page 625
2
~ If. .
-no
-1 5~ q ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
·-2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3
X
Figure 13: Difference of the velocity magnitude be-
tween the Euler solution (camberline pressure model)
and the potential solution, i.e. (Euler solution - po-
tential solution); The dashed line from (-2,0) to
(1,1.5) is the cut along which the two solutions and
the difference are compared in Figure 14.
Compared to the run on a 200x200 mesh with
surface distribution model presented in the previous
section, the current run on a 50x50 mesh with the
camberline pressure model requires only 1/100 of the
computing time, while the predicted velocity field at
locations of interest has the same 1% maximum er-
ror. A three-dimensional version of the camberline
pressure model with relatively coarse mesh is applied
to compute the body force in the unsteady propeller
problem.
4.2 Bocly Force Representation of an
Actuator Disk
In actuator disk theory (Breslin & Andersen 1994),
the (dimensional) thrust, T. is assumed to be uniform
over the disk area, Ao The thrust loading coefficient,
CT, which is based on the ship speed, TOO, is defined
as follows.
T
-pAoU
If the body force is distributed over the dimensionless
distance /\~ in the axial direction, the dimensionless
body force per unit mass, fit, is expressed as follows.
C2~\~ (36)
When the thrust loading coefficient, CT, is spec-
ified, the body force, fit, can be readily computed by
equation (364. The axial distance /\~ is chosen to be
two cell thickness, after it has been confirmed that
NACA 0005, a=0.8, f/c=0.02, a=5°; p=0.08, a2=2.0, a4=1.0, At=0.5 CFL, 50x50 cells
1.15
Magnitude of velocity along a line y=1 +0.5x
| NACAOOOS + a=0.8 (f/c=0 02), 50x50 cells; pressure term only b.f. |
-: 1 ~
1 1
1 05
nq.~ _
_.__ _ ~ Potential solution
_ ~ ~ D ifference
- ~ ~ 1 ~ ~ . ~ ~
0.9 i i i i ~ i i i i ~ ~ i i i i
Distance from (-2,0)
n 1
o of
-n no
· · -0.1
Figure 14: Comparison of the velocity magnitudes
predicted by the Euler solver and the panel method
along the dashed line shown in Figure 13; 50x50 cells,
body force from the camberline pressure model, using
equation (344.
there is only local difference in the solution if this
distance is varied from one cell to four cell thickness
for the usual mesh of which the cell size near the disk
is 0.04.
4.3 Steady Bocly Force Model
In the steady three-dimensional Euler solver, the
body force which represents the propeller can be a
function of space, i.e. ffr,04. The variation in the
axial direction is neglected and the body force is dis-
tributed over the disk area of the propeller within
two cell thickness axially. The body force distributed
in this way can simulate the time averaged spatially
varying loading, such as in the case of a propeller
subject to inclined uniform inflow.
The body force distribution over the disk area is
calculated in the lifting line sense. First, the instan-
taneous blade circulation distributions F(r) along the
lifting lines at each blade angle are computed by the
cavitating potential flow solver (Kinnas et al 19984.
Then, the circulation distribution r'(r,0) over the
disk area is obtained by superimposing all of the in-
stantaneous blade circulation distributions F(r), con-
sidering the skew. Finally, it is assumed that the
radial and circumferential distribution of the thrust
and the torque follow the radial and circumferential
distribution of the circulation r'(r, B). In this work,
only the axial and tangential body forces, fit and fit
are applied. The details of the steady body force
model can be found in (Choi & Kinnas 20014.
OCR for page 626
o
. ~
cd
o
cd
inflow
rid, ~
Figure 15: Unsteady body force model for a propeller
blade; the cylindrical surface at a radius is expanded;
a view from top.
4.4 Unsteacly Bocly Force Mocle}
This is done for the 60 blade angles spaced by 6°,
so that the pressure force and the camber surface
geometry can be found by interpolation from these
60 sets of data for the blade angle corresponding to
the non-dimensional time clicking in the Euler solver.
After the time-interpolation of the pressure force has
been completed, the camber surface geometry and
the pressure force are interpolated in the radial direc-
tion first. For a certain radius, the camber surface ap-
pears as a camberline overlapped by the Euler solver
mesh as shown in Figure 15. Then, the body force for
the cells for which the camberline goes through can
be computed by collecting the pressure force over the
segment of the camberline contained by each cell.
5 NUMERICAL RESULTS
Several convergence studies were performed to find
the effect of the grid resolution and the effect of the
artificial dissipation coefficient, as reported in (Choi
2000, Choi & Kinnas 20014. A grid with 80 cells in
the axial direction and 50 cells in the radial direc-
tion, covering from four propeller radii upstream to
four propeller radii downstream and expanding five
propeller radii radially, was found to be appropri-
ate. Also, zero second order artificial dissipation and
fourth order artificial dissipation coefficients between
0.5 and 1.5 were found to be appropriate for the ex-
amined cases. Concerning the circumferential num-
ber of cells, it was found that at least 12 cells were re-
quired to resolve each harmonic of the effective wake;
i.e. 60 cells were required to resolve 5 harmonics.
In the unsteady solution, the body force must be a
function of time as well as space. Since the blade-
to-blade flow field in the vicinity of the propeller is
resolved in the unsteady solution, the blade should
be represented with more accuracy than that used in
the steady solution.
Three unsteady body force models were consid-
ered in this work; (a) the body force distributed
over the flat propeller-plane-projected blade area,
proportional to the blade bound vortex distribu-
tion F(r,0, t), (b) the body force model utilizing the
bound vortex distribution on blade which includes
the axial variation, F(~, r,0, t), and (c) the unsteady
three-dimensional extension of the camberline pres-
sure model explained in section 4.1.2. The model (c) 5.] Actuator Disk
was finally selected because numerical studies showed
that the models (a) and (b) gave unsatisfactory to-
tal velocity (about +5% error) at the plane where
the effective wake is evaluated. The selected model
has the advantage that it is based on the blade pres-
sure distribution which is directly related to the blade
force and, thus, more rigorous (also more accurate)
than the other two models based on the blade circula-
tion distribution. Note that a time averaged version
of this method is conceptually similar to the steady
body force model of (Toda et al 19904.
The unsteady body force for propeller blades is
obtained by an extension of the two-dimensional cam-
berline pressure model. The pressure, the unit nor-
mal vector, and the area of the BEM element are
known for each side of the blade from the propeller
potential flow solution. The product of the pressure,
the unit normal vector, and the area is calculated
first for each side, and then the two force vectors
from each side are added to obtain the force vector
corresponding to the element on the camber surface.
An actuator disk with thrust loading coefficient
CT=3.0 iS analyzed using the unsteady Euler solver.
In uniform flow, the uniform axial body force is ap-
plied suddenly on the cells at the disk at the initial
dimensionless time t=0 and the time evolution of the
solution is examined. The final solution reached at
the end of the transient flow field should agree with
the steady actuator disk result.
To solver this unsteady actuator disk problem,
the grid with 80x50x30 cells is used with the fol-
lowing parameters; 5000 time steps with /\t=0.002,
~2=O, and ~4=1. The run takes 26.4 CPU hours on a
DEC Alpha Station 600 5/266, and the convergence
is shown in Figure 16. The curves labeled as /\n etc.
stand for the maximum (over all cells) value of the
change between time steps of the axial velocity ~ etc.
at each node. The curve labeled as QM in this fig-
ure represents the maximum error in satisfying the
continuity equation, i.e. the residual of the continu-
ity equation. The dimensionless time, t, is scaled
OCR for page 628
y
WAKE FF-3 U (Pressure Correction Method) l
T = 1.0 h x
Actuator disk (CT=3),
In uniform flow,
SOX50X30 cells,
IpcMax - , IppIter 2,
At=0.002, 5000 Time steps,
a =0.5 a4=1.0,
26.4 CPU hours
WAKEFF-3U (Pressure Correction Method)
Actuator disk (CT=34,
In uniform flow,
SOX50X30 cells,
IpcMax - , IppIter 2,
At=0.002, 5000 Time steps,
aPP=0.5 64=1.0,
26.4 CPU hours
Z VOrt;C;tY
= 5
_ 4.2
_ 3.4
_ 2.6
............. 1 .8
............. 1
....... 0.2
-0.6
-1 .4
-2.2
3
-38
-46
. 54
-
y
T=2.0 Lx
Z VOrt;C;tY
= 5
_ 4.2
_ 3.4
~ 2.6
... 1 .8
...... 1
O 2
22
..~ -
3
-38
46
Figure 19: Time evolution of the z-vorticity field;
t=1.0 (top), t=2.0 (bottom).
be observed that the axial velocity curve reaches its
final shape at around t=4. It can be observed from
the figure that the converged pressure and velocity
along the center line agree well with the analytic val-
ues.
The time evolution of the solution field on the
center plane is shown in Figure 18. It is observed
that the entire flow field has reached almost a steady
state near t=4. In Figure 19, the z-vorticity (Liz)
fields are shown on the center plane corresponding
to the dimensionless times t=1.0 and t=2.0. It is
clearly shown in the figure that the start-up vortex
ring which originates from the rim of t,h~ ~,ct,~n,t,or
1- 1 ~ ~ ~ ~ 1 ~ -
disk at t=U.u, evolves forming a contracting vortex
tube. Since the axial velocity inside the slip stream
is faster than that outside, the z-vorticity, Liz, has
nosit,ive v~,l,~e in the upper half (y > 0) of the `1~-
1—
main, and negative in the lower half of the domain,
as can be seen in the figure. As the start-up vortex
convects downstream, it diffuses due to the coarser
cells downstream and due to the artificial dissipation.
There are two kind of iterations in the present
Table 1: CPU hours on a DEC Alpha Station 600
5/266 for combinations of IpcMax and IppIter.
| IpcMax | | IppIter |
1 11 2 1 4 1 6 1 8 T 12 1 18 1
2 11 1.0 1 1.0 1 11 1 x T X 1
6 1 2.7 1 2.9 1 3.2 1 x 1 x 1 x
10 1 4.5 I x I x I x I I
12 1 1591 1 T 1
18 18.018.91951 x T X 1
20 1901 1 1 T 1
(x means the run diverged.)
unsteady method; one is the pressure correction iter-
ation, of which the maximum number is IpcMax, and
the other is the inner iteration to solve the pressure
correction equation of which the maximum number
is IppIter. To find the optimum number of iterations,
several combinations of IpcMax and IppIter are tried
on the 50x30x30 grid, and the CPU times for each run
are compared in Table 1. It is observed from Table 1
that too many IppIter can cause the run to diverge.
The solution and the convergence history of the con-
verged runs are found to have little dependence on
IppIter. On the other hand, IpcMax is found to have
direct effects on the level of QM, though the gain in
accuracy becomes very small for IpcMax > 10. Com-
promising between the run time and the convergence,
the optimum values have been chosen as IppIter=2
and IpcMax=10.
5.2 Uniform Inflow
Due to the fact that the body force cells are placed
along the blade camber surface rather than on the
"exact" blade surface, the blade thickness effect on
the flow field is completely lost in the unsteady Euler
solver, as observed in the two-dimensional tests. This
means that the flow field obtained by the unsteady
solver is not consistent with the flow field predicted
by the propeller potential flow solver simply because
the propeller potential flow solver does include the
blade thickness effect. Placing the body force cells
on the "exact" blade surface, using the present dis-
cretization in the Euler solver, is practically impos-
sible due to the fact that the resulting size of a cell
is comparable to the blade thickness. Therefore, it is
necessary to quantify this inconsistency because the
difference between the two flow fields is directly re-
lated to the accuracy of the predicted effective wake.
The same inconsistency should affect the steady
Euler solver since it uses the disk-like distribution of
the body force based on the blade circulation distri-
OCR for page 629
Total
.~1 .05
lo
-
.~5
lo.
ACES=
_ A
A
AA
O C,Oo
oo0~0
One-bladed DTMB 4148, at r=0.7200 (j=8), At=0.004, a4=1.0, 80x50x60 cells
1.1 _
| Original Thickness |
~ WAKEFF-3U
A MPUF-3A
Effective Wake
Comparison of the Velocity Fields at X=-0.3R OriginalThickness1
One-bladed DTMB 4148, Js=0.8, fully wetted, uniform inflow, 80x50x60 cells, IpcMax=10, At=0.004, app=0.5, 64=1.0
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
-200 -150 -100 -50 0 50 100 150 200
Angular Position (deg)
Figure 20: Axial velocities predicted by the unsteady
Euler solver and by the potential flow solver for the
One-bladed DTMB 4148 propeller; ~eff=~0.323.
button. However, even though the blade thickness
effect cannot be properly treated by the steady body
force model, the validation test for the uniform in-
flow, presented in (Choi & Kinnas 2001), shows that
the difference between the two flow fields, predicted
by the steady Euler solver and the potential flow
solver, is as small as 1% of the inflow velocity. This
small error can be explained from the fact that the
time averaged effect of the blade thickness is small as
compared to the local instantaneous effect which can
be larger. In other words, the blade thickness effect
becomes more of an issue only in the unsteady Eu-
ler solver because the instantaneous blade-to-blade
flow field, which is affected by the thickness of blades
nearby, is determined in the unsteady Euler solver.
In this section, only the effect of the blade thick-
ness on the blade itself is examined by considering
an One-bladed propeller, so that the influence from
the thickness of the other blades of the propeller is
avoided. The One-bladed propeller is made by re-
moving two of the three blades from the DTMB 4148
propeller, while keeping the geometry of the remain-
ing blade (including the thickness) the same. The
geometry of DTMB 4148 propeller can be found in
(Choi 20004. The inflow condition is unit uniform
inflow so that the velocity field predicted by the po-
tential theory can be compared directly with that
predicted by the Euler solver. The propeller is work-
ing at the advance ratio, Js _ r'D/UOO=0.8, under
the fully wetted condition, where r' is the number of
propeller revolution per second and D is the propeller
diameter. The grid has 80x50x60 cells in the domain
—4 _ ~ < 4 and 0.2 < r < 5. The dimensionless_
time step size for the unsteady solution is /\t=0.004,
and the artificial dissipation coefficients are (72=0 and
Effective
~1 ~ T~1~L 11
no UK: 095 097 099101 103105
0.5
z
z
- 1
z
Figure 21: Comparison among the total, the pro-
peller induced, and the effective wake as predicted by
the unsteady Euler solver at 0.323 radius upstream
of the One-bladed DTMB 4148 propeller.
cT4=1.
The axial velocities along the circumference at
r=0.72 and axial location ~=-0.323 at the instance
the blade is at 0° position (top, or on y-axis) are
shown in Figure 20. The velocity computed by the
unsteady Euler solver (WAKEFF-3U) agrees very
well with the velocity computed by the potential
solver (MPUF-3A) at the angular positions away
from the blade. The agreement between the two ve-
locities becomes worse at the blade locations within
+50°, where flow deceleration and acceleration occur
when the leading edge and the trailing edge of the
blade, respectively, are aligned with the considered
point. This is attributed to the fact that the camber-
line pressure model can not simulate accurately the
effect of blade thickness on the flow near the blade,
as was also demonstrated in section 4.1.2. The dis-
agreement of the velocities results to an error in the
effective wake ranging from -1% to 2%.
In Figure 21, the three velocity fields at a loca-
tion 0.323 radius upstream of the propeller are shown
in order to observe the distribution of the total, the
propeller induced, and the effective wake velocities
on this plane. The velocities are plotted at the mo-
ment when the key blade is at the 0° (top) position.
While the total and the propeller induced velocities
agree very well globally, the two velocities near the
blade differ by 1 to 2%. This difference results in
the locally erroneous effective wake, as shown at the
right of Figure 21.
The results shown up to now are computed for
the propeller with the original thickness. To see the
effect of blade thickness, a similar computation is per-
formed with a zero-thickness propeller. The axial ve-
OCR for page 630
One-bladed DTMB 4148, at r=0.7200 (j=8), At=0.004, a4=1.0, 80x50x60 cells
1.1 ~
| Zero Thickness |
~ WAKEFF-3U
_ A MPUF-3A
_ ~ Effective Wake
.~1 .05
lo
-
.~5
Comparison of the Velocity Fields at X=-0.3R Zero Thickness I
One-bladed DTMB 4148, Js=0.8, fully wetted, uniform inflow, 80x50x60 cells, IpcMax=10, At=0.004, app=0.5, 64=1.0
A
_ ~
_ ~~ ~~ ~ _
AA E]
~3~9—_
~ ] ~
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
-200 -150 -100 -50 0 50 100 150 200
Angular Position (deg)
Figure 22: Axial velocities predicted by the unsteady
Euler solver and by the potential flow solver for zero-
thickness One-bladed DTMB 4148 propeller; ~eff=-
0.323.
locities along the same circumference are compared
in Figure 22. The agreement between the two veloci-
ties is very good. As a result, the effective wake has a
smaller error ranging from -0.7% to 0.2%. When the
MPUF-3A results in the two Figures 20 and 22 are
compared to each other, the blade thickness effect
is more prominent in the deceleration region (near
B=-25°) than in the acceleration region. It is very
interesting that the maximum error of the original
thickness result shown in Figure 20 occurs also near
the deceleration region. In the deceleration region,
the velocities predicted by the Euler solver are found
to be affected very little by the blade thickness, while
the potential solutions show appreciable effect of the
blade thickness. Therefore, it can be concluded that
the disagreement near the deceleration region shown
in Figure 20 is mostly due to the defect of the blade
thickness effect in the Euler solver.
In Figure 23, the velocity fields on the plane at
0.323 radius upstream of the propeller are shown for
the zero-thickness propeller case. It is observed that
the total and the propeller induced velocities have
very similar contour shapes, and that the effective
wake has maximum error near the hub ranging +1%.
The error near the hub is contributed by the differ-
ence of boundary condition at the hub. The zero
normal velocity boundary condition is applied on the
hub in the Euler solver, while the hub is not included
in the potential flow solver. Even though the poten-
tial flow solver has the ability to consider the hub
via a simplified image model, the hub is turned off in
the present potential solver computation because the
simplified image model has not been found to pro-
duce the correct flow field upstream of the propeller
Total
0.5
Effective
~1 ~ T~1~L 11
no UK: 095 097 099101 103105
Am'
z
z
- 1
z
Figure 23: Comparison among the total, the pro-
peller induced, and the effective wake as predicted by
the unsteady Euler solver at 0.323 radius upstream of
the zero-thickness One-bladed DTMB 4148 propeller.
plane.
To confirm the blade thickness is responsible for
the inconsistency shown in Figure 20 again, the fol-
lowing test is performed: When the propeller in-
duced velocity is computed, the velocity induced by
the dipoles (equivalent to vortex loops) and sources
on the blades and the associated trailing wakes are
summed. The dipoles are responsible for the lift of
the blade, while the sources represent the blade thick-
ness and the cavity thickness, if any. If the velocity
influence due to the sources is omitted in the sum-
mation, the blade thickness effect on the induced ve-
locity is excluded, and thus the potential flow model
becomes consistent to the camberline pressure model
in the Euler solver. The axial velocities predicted
in this way are shown in Figure 24. Note that the
agreement between the Euler solution and the po-
tential flow solution is as good as in the case of zero-
thickness propeller, shown in Figure 22.
In conclusion, the effect of the blade thickness is
found to be the major source of error in the calcula-
tion of the effective wake. Once the blade thickness is
eliminated either by using a zero-thickness propeller
or by omitting the velocity influence of the thickness
sources, the error in the effective wake is less than
1% of the inflow velocity.
5.3 Propeller DTMB 4842 in Uniform
Inflow
For the unsteady analysis of the five bladed DTMB
4842 propeller, working at Js=0.9 in uniform inflow,
a grid with 80x50x120 cells which covers the domain
—4 < ~ < 4, 0.323 < r < 5 is used. The propeller
OCR for page 631
One-bladed DTMB 4148, at r=0.7200 (j=8), At=0.004, a4=1.0, 80x50x60 cells
1.1 ~
| Original Thickness l ~ WAKEFF-3U
l No Source Effect in Induced Velocity | A MPUF-3A
~ Effective Wake
.~1 .05
lo
-
.~5
T
AN
_ AM
AK] ~
_ _ Ed
~ Al
A 1~ ~
Zi
"El
;9
_~ _ ~~ ~ ~ J~ ~ ~ ~
~~ 0~o
-
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
-200 -150 -100 -50 0 50 100 150 200
Angular Position (dog)
Figure 24: Axial velocities predicted by the unsteady
Euler solver and by the potential flow solver for the
one-bladed DTMB 4148 propeller at ~eff=~0.323.
When the propeller induced velocity is computed, the
influence of the blade thickness is excluded.
geometry is given in (Choi 20004. The number of
cells in the circumferential direction is chosen to be
120 because this propeller has five blades and more
variation is expected in the circumferential direction.
The blade thickness is reset to zero. One revolution
of the propeller is divided into 900 time steps using a
time step size /\t=0.002, the unsteady computation
parameters are IpcMax=10, IppIter=2, alpp = 0.5,
and the artificial dissipation coefficients are a2=0 and
c~4=5. The numerical simulation for the one revolu-
tion of the propeller took 39.5 CPU hours on a DEC
Personal Workstation 500au2.
In Figure 25, the axial body force distribution
that represents the propeller in the unsteady Euler
solver is shown at the instance the key blade is at
the top position (0°~. The body force cells are rep-
resented by the wire frame shaded according to the
strength of the axial component of the body force at
each cell.
The axial velocities along the circumference of
radius 0.7265 on this plane 0.323 radius upstream of
the propeller are shown in Figure 26. The velocities
are shown at the moment the key blade is at the top
position (0°~. Since there are five blades, the five
velocity peaks are shown in the figure. The mean
acceleration at this location is about 16% of the in-
flow velocity, while the variation between blades is
about 5.5% of the inflow velocity. The velocities pre-
dicted by the unsteady Euler solver (WAKEFF-3U)
and the potential flow solver (MPUF-3A) agree very
well leaving only -0.4% to +0.7% error in the effective
2DEC Personal Workstation 500au is approximately equiv-
alent to a 1 GHz Pentium PC.
WAKE F F-3 U Zero Thickness DTMB 4842 in uniform nominal wake;
SOx50x120 cells, IpcMax=lO, IppIter=2, otpp=0.5,
At=0.002, 900 steps, 62=O, 64=5.O, 39.5 WAKE-CPU hours ~ X
Figure 25: DTMB 4842 propeller represented with
body force cells; key blade at 0°.
wake velocity prediction.
5.4 Nominal Wake with Only the
Thirc! Harmonic
QX
14
13
12
1 1
10
9
8
7
6
5
4
3
2
1
o
In this section, a nominal wake which only has a third
cosine harmonic is considered with the zero-thickness
DTMB 4148 propeller. Since the propeller has three
blades and the inflow wake has only the third har-
monic, the resulting propeller loading should have
dominant third harmonic components. The nominal
wake has only the axial velocity which can be de-
scribed by the following equation.
{ 0 7 + 0.3 Cos(36) for r < l o (37)
The propeller working condition considered in
this section is Js=0.8 and fully wetted flow. The
fluid domain covering—4 < ~ < 4 and 0.2 < r < 5
is discretized with 80x50x60 cells. Since the dimen-
sionless time step size is determined to be /\t=0.004,
400 time steps correspond to one full revolution of
the propeller. The applied artificial dissipation coef-
ficients are (72=0 and a4=1.
The convergence with revolution number is
shown in Figure 27. The solution (n, v, w, p) at
(~,r,04=~-0.323,0.683,0.0) when the key blade is at
0° (top) position of the ninth revolution is subtracted
from that at the same location for the same key blade
position at the sixth, the seventh, and the eighth
revolution. Since the resulting differences of the so-
lutions from those of the ninth revolution converge,
the solution of the ninth revolution can be practically
regarded to fully developed; i.e. the solution repeats
itself for consecutive revolutions.
OCR for page 632
1.2 r
1.15
·<, 1 . 1
o
-
~ 1 .05
Zero Thickness DTMB 4842, at r=0.7265 (j=7), At=0.002, ad=5.0' 80x50x120 cells
~ Em, it,,—~
- Cad E2Y ED ~ ED
— EIA ~ MA ~ MA [a EIA
_
- Ad, ENS,,\ ~ ~~,~ Ens,,\
or ~ ~ ~ ~
~ _=
0 95, 1 1 1 1 1 1 1 1 1 1
-200 -150 -100
Residual of the Continuity Equation and Propeller Force during One Revolution
Zero-thickness DTMB 4148 in A3 Inflow Wolfe
1 o-2 0.7
3f:0 270 Blade ~ 8gle (deg.) SOx50x60 cells
M QM ~ 400, 10, 2, 20
T KT ~ 0.004, l.e-9, 1., 1., 0.5
o 10KQ ~ 0.0, 1.0
~ 12.2 CPU hours
1 n-3 L
.
10-'
10-0
10-6
10-7
~ WAKEFF-3U 1 o-8
A MPUF-3A
o Effective Wake
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
-50 0 50 100 150 200
Angular Position (deg.)
Figure 26: Axial velocities predicted by the unsteady
Euler solver and by the potential flow solver for the
zero-thickness DTMB 4842 propeller; ~eff=~0.323.
10-2
1n-3 _
Convergence with respect to Revolution
Zero-thickness DTMB 4148 propeller in A3 inflow wake
solution at (x,r,q)=(-0.323, 0.683, 0.0)
when the key blade is at Do position
CL
~ 10-5—
_
_
10-6—
1 n-7
:
_ ~ i
U AU
_ V AV
W AW _
P AP
~ \ \:
: : : : : : me\\ : : : :
you = 1 U - Ug 1
~V= V-Vg
_ Law= IW-W9 1
Ap=lP-P
_
7
5 6 7 8 9
Revolution
Figure 27: The convergence with respect to the rev-
olution for zero-thickness DTMB 4148 in the wake
with only the third harmonic.
~~m;~ ~
10-9 1 1 1 1 1 1 1 1 1 1 1 1 1 1 - Q
0 0.4 0.8 1.2 1.6
Time (=R/U)
0.6
0.5
0.4
o
0.3
0.2
0.1
Figure 28: The residual of the continuity equa-
tion and the propeller thrust and torque coefficients
during the last (ninth) revolution of zero-thickness
DTMB 4148 in the wake with only the third har-
mon~c.
Figure 28 shows the variation of the thrust and
torque coefficients, AT and KQ, of the propeller dur-
ing the ninth revolution. As expected, the three
propeller loading peaks, of which the amplitudes are
about 50% of the mean loading, are clearly observ-
able. The curve labeled QM in the figure is the resid-
ual of the continuity equation, which is less than 10-5
during the ninth revolution.
To observe the circumferential variation of the
velocities, the velocities along the circumference at
radius 0.7200 on the plane ~=-0.3 for the moment
corresponding to the blade angles, B=0°, are plot-
ted in Figure 29. Along the circumference, the to-
tal velocity is varying between 0.59 and 1.08, while
the propeller induced velocity varies from 0.8 to 1.4.
Notice that the maximum of the propeller induced
velocity occurs about 20° ~ 30° after the blade. The
effective wake varies from 0.48 to 0.96 and this ampli-
tude is smaller than that of the nominal wake which
varies between 0.4 and 1.0. In addition to these ve-
locities, the steady effective wake computed by the
steady Euler solver is plotted together. It is very in-
teresting that the unsteady effective wake at this mo-
ment (key blade at B=0°) is very close to the steady
effective wake.
In Figures 30 and 31, these velocities are com-
pared for different moments corresponding to blade
angles -24° and -72°. Notice that the nominal wake
and the steady effective wake (WAKEFF-3D) are,
by definition, always the same. It can be observed
from the comparison of velocities at different mo-
3A right-handed propeller is used in the present application.
That is, the propeller rotates in the direction of decreasing 0.
OCR for page 633
Zero Thickness DTMB 4148 in A3 Wake. 80x50x60 cells
Zero Thickness DTMB 4148 in A3 wake, at r=0.7200 (j=8), At=0.004, a4=1 .0, 80x50x60 cells
1.1 _
08
0.3 _
0.2 _
0 1
: : ~ ~
..
* * * *
: ** ~ : : ** ~
O. 1 1 1 1 ~ 1 1 1 1 ~ 1 1 1 1 ~ 1 1 1 1 ~ 1 1 1 1 ~ 1 1 1 1 ~ 1 1 1 1 ~ 1 1 1 1 ~
-200 -150 -100 -50 0 50 100 150 200
Angular Position (deg.)
Figure 29: Circumferential variation of the velocities
at r=0.72; Blade angle = 0°.
Zero Thickness DTMB 4148 in A3 wake, at r=0.7200 (j=8), At=0.004, a4=1 .0, 80x50x60 cells
1.1
1 L
0.9
n ~
~ 0.7
O 0.6 _
`~5 0.5
X
0.4 _
0.3 _
0.2 _
1 ,_ ~ , at_ ~ - _,
_ ^,~~Q~ ~^,e,~~Q~ ~^Q~~~ ~^^
0, _ I I I :l I I I l: I I I I I :l I I I l: I I I I I :l I I I l: I I l: I I :l I I I l: I I
-200 -150 -100 -50 0 50 100 150 200
Angular Position (deg.)
n 1
: : ~~ :
13 El
::
: : *a :
GAG
: * : ~ :
: ~ ~ : ~ :
: : ~~ : :
~ 13
. .
: : *0~* : :
: * ~ : :
: ~ ~ ~ C
Aid: : :
hi: hi: : :
::
*~*: : :
0 ~
*: ~ : :
8: : hi: :
_ ~~ ~ 8 ~~ 3 (3
- g* : : *a: g* :
- 88* ~ ~ ~ *&*
- * * *
=~
Figure 30: Circumferential variation of the velocities
at r=0.72; Blade angle = -24°.
Zero Thickness DTMB 4148 in A3 wake, at r=0.7200 (j=8), At=0.004, 64=1.0, 80x50x60 cells
1.1
n q
08
03
n ~
Figure 31: Circumferential variation of the velocities
at r=0.72; Blade angle = -72°.
0 8
0 4
n ~
0 0.25 0.5 0.75
A3, Vx
Figure 32: Comparison of the mean and the third
harmonic of the time averaged unsteady effective
wake.
meets that the unsteady effective wake does not vary
very much with time, even though the total and the
propeller induced velocity do change. Also, it is
found from this comparison that the unsteady effec-
tive wake is always very close to the steady effective
wake.
:e ~~ 6: : : ~ hi:
:*: g*: : : :*~:
= - ~ ~ ~ Add
_ *:;:* *
Total Velocity
Q Propeller Induced _
Unsteady Effective Wake
_ The time averaged unsteady effective wake can
be calculated easily from the time series data of the
unsteady effective wake. In Figure 32, the circum-
ferential mean and the third harmonic amplitude of
the time averaged unsteady effective wake are shown.
The circumferential mean is about 0.03 greater than
that of the nominal wake, while the third harmonic
amplitude is about 0.06 smaller than that of the nom-
inal wake. In Figure 32, the mean and the third
harmonic amplitude obtained from the steady Euler
solver (WAKEFF-3D) are shown together. The mean
of the time averaged unsteady effective wake is 0.721,
at radius 0.72, which is very close to the mean of the
steady effective wake with value 0.715, at the same
radius. Also, the third harmonic amplitude of the
time averaged unsteady effective wake is 0.241, while
that of the steady effective wake is 0.247. It can be
observed again that the time averaged unsteady ef-
fective wake and the steady effective wake are very
close to each other for this case.
In conclusion, the unsteady effective wake is
found to vary little with time for the case of DTMB
4148 propeller and the nominal wake with only the
third harmonic. As expected, the time averaged un-
steady effective wake is also very similar to the steady
effective wake obtained by the steady Euler solver.
OCR for page 634
Cavity Shapes Prop. 4148, H20 Wake, Js=0.9087, an=2.576, Fn=9.159
WAKEFF-3D (Tunnel), 60x44x60 cells, a4=1.5, Xeff=-0.3 1.2
// 30 deg. \\
1 1 1 1 1
0.5 1
Figure 33: Predicted cavity shape, 60 circumferential
cells.
Figure 34: Photographs of the cavity shape observed
in CAPREX I.
5.5 Cavitation on Propeller DTMB
4148
Since it is well known that propeller cavitation is very
sensitive to the non-uniformity of the inflow, the ac-
curacy of the predicted effective wake can be indi-
rectly judged by the cavitation shape predicted for
a given effective wake. The cavitation can be pre-
dicted by the propeller potential flow solver once the
effective wake is obtained through the iterations be-
tween the Euler solver and the potential flow solver.
It should be noted that in the evaluation of the ef-
fective wake the propeller induced flow field has been
determined by excluding the influence of the thick-
ness and cavity source.
In this work, the cavitation shapes observed in
CAPREX I (Mishima et al 1995) are considered. The
effective wake is predicted using the nominal wake
measured in CAPREX II (Mishima et al 1995) and
propeller DTMB 4148 working at Js=0.9087 and the
cavitation number based on the number of revolu-
tion, con—(Pshaf t—Pvapor)/(O 5prl2D24=2.576 which
is the condition at which the photographs of cavita-
tion were taken in CAPREX I. Note that the tun-
Total Velocity (prop.4148, Js=1.05, Xeff=-0.323)
~~ ~
Hi
~ 0.4 _
O.
- - - ~ - - R=0.689 (measured)
R=0.683 (WAKEFF-3D)
R=0.683 (WAKEFF-3D, tunnel flow)
R=0.720 (time average of WAKEFF-3U)
2 _
_
O. 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
-200 -100 0 100 200 300
Angular Position (deg.)
Figure 35: Comparison of predicted axial velocity at
r=0.68R with the measured one in CAPREX II.
nel effect was included by placing the circular tunnel
wall boundary in the prediction of the effective wake.
These effective wakes were predicted using grids with
60 cells circumferentially. In Figure 33, the predicted
cavity shapes are shown. The photographs taken in
the experiment CAPREX I are shown in Figure 34.
The predicted cavity planforms are very similar to
those shown in the photographs.
5.6 Total Velocity in front of Pro-
peller DTMB 4148
The velocity field predicted by the Euler solver can
be directly compared with the velocity measured in
CAPREX II (Mishima et al 19954. In the experiment,
the propeller DTMB 4148 was working at JS=1.05
in the non-axisymmetric wake. The velocity mea-
surement at 0.32 radius upstream of the propeller is
selected to be compared with the predicted velocity
because this is the place where the effective wake is
usually evaluated.
The velocity distribution along the circumfer-
ence at radius 0.68 is shown in Figure 35, and that at
radius 0.98 is shown in Figure 36. The axial veloci-
ties predicted by the steady Euler solver (WAKEFF-
3D), by the same steady solver with the tunnel flow
simulation (WAKEFF-3D, tunnel flow), and by the
unsteady Euler solver (WAKEFF-3U) are compared
together. The velocity marked as "time average of
WAKEFF-3U" is obtained by averaging the time his-
tory of the velocity at a fixed point predicted by the
unsteady Euler solver. The time averaging is ex-
actly equivalent to the method in which the mea-
sured data were processed during the experiment.
In the unsteady Euler solution, the propeller load-
ing is obtained from the effective wake predicted by
OCR for page 635
1.2 r
1
`0.8 _
O _
0.6 _
_
~ _
0.4 _
0.2 _
Total Velocity (prop.4148, Js=1.05, Xeff=-0.323)
. _~~~;
---O -- R=0.984 (measured)
R=0.980 (WAKEFF-3D)
R=0.980 (WAKEFF-3D, tunnel flow)
R=1.017 (time average of WAKEFF-3U)
O' AN
u2 30 -100 0 100 200 300
Angular Position (deg.)
Figure 36: Comparison of predicted axial velocity at
r=0.98R with the measured one in CAPREX II.
the steady solver. This propeller loading includes the
blade thickness effect because the purpose is to ob-
tain the total velocity and the effective wake is not
evaluated in this computation.
Even though the predicted profiles cannot follow
the sharp changes in the measured velocity, the global
agreement between the computed and the measured
velocity is very good. The effect of the tunnel in this
case is found to be very small because the propeller
loading is very light. Also, it can be observed from
the figures that the time average of the unsteady ve-
locity field is very close to the steady velocity field.
6 CONCLUSIONS
Robust steady and unsteady Euler equation solvers,
coupled with a previously developed propeller poten-
tial flow solver, are developed for the prediction of
the steady and unsteady effective wakes. To the au-
thor's knowledge, the present unsteady Euler solver
is the first method that ever computes the unsteady
eJective wake with sufficient space and time resolu-
tion. Through extensively two-dimensional studies, it
is shown that the surface distribution method, with
sufficient number of cells, can simulate the flow field
very accurately. In addition, the camberline pressure
model is suggested for coarser grid applications and
used in the present unsteady three-dimensional Eu-
ler solver. The blade thickness effect is identified as
the major source of error (in the order of 2%) in the
prediction of the unsteady effective wake. Omitting
the velocity influence of the thickness source has been
found to improve the accuracy of the evaluated ve-
locities, suppressing the error under 1% of the inflow
velocity.
The present method is applied to a nominal wake
inflow which contains only the third harmonic. It
is found that the unsteadiness in the predicted on-
steady eJective wake is small, and that the time av-
erage of the unsteady eJective wake is very close to
the steady elective wake predicted by the steady Enter
solver method. This finding leads to the conclusion
that the three-dimensional steady Euler solver can
be a computationally efficient, yet accurate, practical
tool for the prediction of effective wakes. The cavity
shape predicted using the effective wake obtained by
the present method showed good agreement with the
photographs taken during the experiment. The cir-
cumferential distribution of the axial velocity com-
puted by the present method agreed also very well
with the one measured upstream of the propeller.
As a future research, an improved body force
model which includes the mass source to represent
the blade thickness can be considered. The inclusion
of the mass source would also allow for a consistent
representation of unsteady cavity shapes within the
Euler equation solver. Also, a way of linking the
unsteady effective wake with the propeller potential
flow solver should be devised in the future. This is
necessary to investigate if a truly unsteady way of
feeding the unsteady effective wake into the propeller
potential solver could lead to the amplification of the
unsteadiness.
Acknowledgment
Support for this research was provided by the fol-
lowing companies and research centers: American
Bureau of Shipping, Daewoo Heavy Industries Ltd.,
David Taylor Model Basin, El Pardo Model Basin,
Hamburgische Schif~bau-Versuchsanstalt (HSVA),
Hyundai Maritime Research Institute, Kamewa AB,
Michigan Wheel Corporation, Rolla SP Propellers
SA, VA Tech Escher Wyss (former Sulzer-Hydro
GmbH), Ulstein Propellers AS, Volvo-Penta of the
Americas, and Wartsila Propulsion.
References
ANDERSON, J. 1995 Computational Fluid Dynam-
ics: The Basics with Applications. McGraw-Hill,
Inc., New York.
BRESLIN, J. AND ANDERSEN, P. 1994 Hydrody-
namics of Ship Propellers. Cambridge University
Press.
CHEN, B., STERN, F. AND KIM, W. 1994 Compu-
tation of unsteady viscous marine propulsor blade
and wake flow. Proceedings: Twentieth Symposium
on Naval Hydrodynamics. August.
OCR for page 636
CHo~, J.-K. 2000 Vortical inflow - propeller inter-
action using an unsteady three-dimensional Euler
solver. Doctoral dissertation, Department of Civil
Engineering, The University of Texas at Austin.
August.
CHo~, J.-K. AND K~NNAs, S. 1998a A 3-D Euler
solver and its application on the analysis of cavi-
tating propellers. Proceedings: Proceedings of the
25th ATTC. American Towing Tank Conference,
Iowa City, Iowa, September.
CHo~, J.-K. AND K~NNAs, S. 1998b Numerical
water tunnel in two and three dimensions. Journal
of Ship Research, 42, 2, June, pp. 86-98.
CHo~, J.-K. AND K~NNAs, S. 2001 Prediction of
Non-axisymmetric Effective Wake by a 3-D Euler
Solver. Journal of Ship Research (to appear).
CHoR~N, A. J. 1967 A numerical method for solving
incompressible viscous flow problems. Journal of
Computational Physics, 2, pp.12-26.
COURANT, R., FR~EDR~cHs, K. AND LEWY, H.
1967 On the partial difference equations of math-
ematical physics. IBM Journal, Mar, pp.215-234.
DAMLE, S., DANG, T. AND REDDY, D. 1997
Throughflow method for turbomachines applicable
for all flow regimes. Journal of Torbomachinery,
119, Apr. pp.256-262.
HUANG, T. AND Cox, B. 1977 Interaction of
afterbody boundary layer and propeller. Proceed-
ings: Symposium on Hydrodynamics of Ship and
Offshore Propulsion System. March.
HUANG, T. AND GROVES, N. 1980 Effective wake
: Theory and experiment. Proceedings: 13th Sym-
posium on Naval Hydrodynamics. October.
HUANG, T., WANG, H., SANTELLI, N. AND
GROVES, N. 1976 Propeller/stem boundarylayer
interaction on axisymmetric bodies: Theory and
experiment. Technical report. DTNSRDC 76-0113.
DTNSRDC. December.
KERw~N, J., KEENAN, D., BLACK, S. AND Dings,
J. 1994 A coupled viscous/potential flow de-
sign method for wake adapted multi-stage, ducted
propulsors using generalized geometry. Trans.
SNA ME, 102.
KERw~N, J., TAYLOR, T., BLACK, S. AND
McHuGH, G. 1997 A coupled lifting-surface
analysis technique for marine propulsors in steady
flow. Proceedings: Propellers/Shafting '97 Sympo-
sium. Soc. Naval Arch. & Marine Engnrs., Virginia
Beach, VA, September, 1-15 (Paper No. 204.
K~NNAs, S., GRIFFIN, P., CHOT, J.-K. AND KosA~,
E. 1998 Automated design of propulsor blades
for high-speed ocean vehicle applications. Trans.
SNA ME, 106.
K~NNAs, S., CHOT, J.-K., KosA~, E., YOUNG, J.
AND LEE, H. 1999 An integrated computational
technique for the design of propellers with specified
constraints on cavitation extent and hull pressure
fluctuations. Proceedings: CFD99 - The Interna-
tional CFD Conference. June 5-7.
LARSSON, L. 1997 The state-of-the-art of CFD for
ship hydrodynamics. Proceedings: The Interna-
tional CFD Conference. May 29-31, 1-17 (Paper
No. 04.
M~sH~MA, S., K~NNAs, S. AND EGNOR, D.
1995 The CAvitating PRopeller EXperiment
(CAPREX), Phases I & II. Technical report. De-
partment of Ocean Engineering, MIT. August.
No, R.-H. 1982 A multiple-grid scheme for solving
the Euler equations. AIAA Journal, 20, 11, Nov.
pp.1565-1571.
PATANKAR, S. 1980 Numerical Heat Transfer and
Fluid Flow. Hemisphere Publishing, New York.
RHIE, C. AND CHOW, W. 1983 Numerical study of
the turbulent flow past an airfoil with trailing edge
separation. AIAA Journal, 21, 11, Nov. pp.1525-
1532.
SHTH, W.-Z. 1988 A combined Euler equa-
tion/surface panel solution to the shear interac-
tion problem of an open or ducted propeller. Doc-
toral dissertation, Department of Ocean Engineer-
ing, M.I.T. October.
STANTER, M. J. 1998 The application of RANS
code to investigate propeller scale effects. Proceed-
ings: 22nd Symposium on Naval Hydrodynamics.
August 9-14, 199-211.
STERN, F., KIM, H., PATEL, V. AND CHEN,
H. 1988a Computation of viscous flow around
propeller-shaft configurations. Journal of Ship Re-
search, 32, 4, pp.263-284.
STERN, F., KIM, H., PATEL, V. AND CHEN, H.
1988b A viscous-flow approach to the computa-
tion of propeller- hull interaction. Journal of Ship
Research, 32, 4, pp. 246-262.
STERN, F., KIM, H., ZHANG, D., KERw~N, J. AND
JESSUP, S. 1994 Computation of viscous flow
around propeller-body configurations: Series 60 cb
=0.6 ship model. Journal of Ship Research, 38, 2,
pp. 137-157.
TODA, Y., STERN, F., TANAKA, I. AND PATEL, V.
1990 Mean-flow measurements in the boundary
layer and wake of a Series 60 CB=0.6 model ship
with and without propeller. Journal of Ship Re-
search, 34, 4, pp.225-252.
OCR for page 637
DISCUSSION
J. Kerwin
Massachusetts Institute of Technology, USA
This paper presents an extremely careful and
systematic approach to the problem of hull-propeller
interaction. I agree with the authors that getting the
effective wake right is particularly important in the
case of predicting unsteady propeller cavitation, since
it is known that small changes in the inflow can result
in large changes in cavity volume. The authors' work
is therefore particularly important since the weak link
in the overall process of predicting both unsteady
cavitation and propeller-induced hull vibratory
excitation is the uncertainty in the inflow.
As the authors point out in their introduction, we
have been working for some time on developing the
coupling concept for both steady and unsteady
propeller flows. In the case of steady flow, the global
flow solver that we employed was initially a steady,
axisymmetric RANS code. While this procedure
works well, we found that almost all of the
computing time was devoted to the RANS solution,
with only a small part being associated with the
vortex-lattice blade solution. In addition, RANS
"ridding issues proved to be the major source of time
delays in obtaining a solution.
We therefore decided to explore the feasibility of
using MTFLOW, an Euler code developed by
Professor Mark Drela at MIT, in place of a RANS
solver. This procedure is therefore closer to the
steady-flow method presented by the authors. Our
experience so far has been that both computing times
and user-preparation times have been substantially
reduced, while the accuracy of the solution is
comparable.
Figure 1 shows a preliminary result obtained for a
mixed-flow waterjet pump as part of an ongoing
Naval Engineer's thesis being carried out by C.J.
Hanson. The circumferential mean streamlines are
plotted, together with contours of swirl velocity.
Of principal interest is the comparison of rotor thrust
and torque predicted by the Euler code and the RANS
code. The results are quite close, in spite of the fact
that this is a particularly complex flow case. The
reduction in computing time from 12 hours to 30
minutes is very substantial.
Rotor Stator
M~BD RAN-S Method
IT = 1.97 KT = 1.87
KQ=0.415 KQ=0.396
Time = 30 min Time = 12 hours
Figure 1. Comparison of results for a mixed-flow
waterjet using a coupled EulerlLifting surface code
(MTFLOWIPBD) and a coupled RANSlLifting
surface code.
I am particularly interested in the authors'
comparison of their computed time-averaged
effective wake and their computed steady effective
wake as shown in Section 5.4. We make the
assumption, at the outset, that an axisymmetric
nominal wake remains axisymmetric. Thus, in the
steady flow case, one only needs to couple the lifting-
surface code with an axisymmetric RANS or Euler
solver. The authors' conclusion would seem to justify
this assumption. Could the authors elaborate on
whether this assumption is valid for all practical
cases, or whether there are situations where a full 3-D
unsteady Euler solution is required?
AUTHOR'S REPLY
The authors would like to thank Prof
Kerwin for taking time to review this paper and for
making valuable comments.
It is very encouraging that the resulting
thrust predicted by the coupled Euler/lifting surface
code is quite close to the result obtained by a coupled
RANS/lifting surface code. The fact that the propeller
thrust and torque predicted by the two methods are
quite close implies the closeness of the predicted
effective wake inflows. When we decided to solve
the Euler equations, we neglected the effect of
viscosity in the interaction between the vertical
inflow and the propeller potential flow, based on the
following two assumptions. (a) The effect of
viscosity is limited very near the wall boundary such
as a hub, and (b) the major player in the interaction is
the inviscid vorticity dynamics. These assumptions
seem to be valid even for a very complex flow inside
a mixed-flow water-jet according to Prof. Kerwin's
numerical experiment.
1
OCR for page 638
non axisymmetric crurlysis of The effective i flow
should be separated from the discussion of the steady
versus unsteady effective wake As c tlrlcl result of
the current unsteady non~xisymmetric Euler solver,
w obtain the unsteady nonoxi mmehic effective
wake i flow, q~(r,t),t), while the steady non-
axisymmehic effective wake, q~(r,tl), cm be
obtained from the iteration betw en the steady non-
axisymmetric Euler solver Ed the lifting su face
code Note that the body force used in this steady
analysis is The time averaged propeller force, which is
varying in pace here is yet mother prediction by
using en axisymmetric Euler solver, m which The
velocity flow field is c function of only the radius,
q~(~) Inthe axisymmetric crurlysis, the body force
is obtained by circumfe~enticlly averaging the
propeller force One of The conclusions m the paper is
that the time average of q~(r,t),t) is close to
q, (r, tl) for the case tested m the current work
In The course of our development of the
axisymmetric Ed non axisymmeh ic Euler codes, w
w re also intere ted in The question if the
cacumfere ticl average of q~(r,O) would be close
to q~(~) in The se- eras cases we have tested, we
have found that the circumferential average of
q~(r,O) is close to q~(~) herefme, w believe
that the axisymmeh ic crurlysis is also quite useful in
predicting the me m performance of propulsors
Finally, w carmot yet generalize our conclusion that
thetimeaverageof q,(~,l9',t) iscloseto q~(r,tl)
bee mse w have very limited mmmber of test cases
One obstacle to reach c final conclusion on this issue
is the fact that we carmot input The unsteady effective
i flow into The cunent Ifftmg surface crurlysis due to
the time invari mt i flow assumption built m the
medhod In my case, the unfreed. Euler solver is still
useful for Inherently unsteady flow cases such es
acceler t i on, dece lerct ion, Ed m cneuver i g of c ship
Representative terms from entire chapter:
euler solver
