|
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 669
Paper submitted to
TheTwenty-Fourth Symposium on Naval Hydrodynamics
July 8-13, 2002, Fukuoka, Japan
VALIDATION OFTHE FLOW AROUND ATURNING SUBMARINE
Chao-Ho Sung, Ming-Yee Jiang, Bong Rhee, Scott Percival,Paisan Atsavapranee, and In-Young Koh
David Taylor Model Basin
NSWCCD, MD 20817-5700, U.S.A.
ABSTRACT
A numerical procedure based on solving the incompressible
Reynolds-averaged Navier-Stokes equations in a steady rotating
coordinate system has been developed for the prediction of the
flow around a turning submarine. Computed results are com-
pared with the measured data obtained from an unclassified sub-
marine called ONR Body-l. The performance of the standard
k—w turbulence model and a realizable k—~ turbulence model
is compared and discussed. Several improvements in numerical
methods including higher-order spatially accurate schemes, lo-
cal refinements, and new wall function are suggested for more
accurate predictions.
INTRODUCTION
Before full scale sea trials are performed, experimental
data obtained by straight-line captive-model rotating arm experi-
ments(RAE) and radio-controlled models(RCM) provide the cru-
cial information for the design of maneuvering submarines. From
the point of view of the validation of Computational Fluid Dy-
namics(CFD) techniques, the validation of CFD with RAE data
is most fundamental. This is because RAE can be treated as a
time independent problem in a steady rotating coordinate system
thus avoiding complications by time dependency and other oscil-
latory flows. Yet RAE maintains most of the major flow physics
unique to a maneuvering submarine such as vertical flow interac-
tions with hull boundary layer and flow separations. Simplicity
of the flow problems in RAE makes it easier to identify CFD
difficulties and then to make necessary improvements. The pur-
pose of this paper is to present the current predictive capability
of a RAE model through validation with available experimental
data. Improvements to increase the prediction accuracy will also
be discussed.
The outline of the paper is the following. The govern-
ing equations of the incompressible Reynolds-averaged Navier-
Stokes(RANS) equations and the turbulence models in a steady
rotating coordinate system will be described. As expected, turbu-
lence models play an important role in the accurate prediction of
the flow around a maneuvering submarine. The numerical method
used will be discussed next. Only some special features of the nu-
merical schemes will be highlighted omitting most of the details.
The experimental performed on an unclassified submarine-like
body with a sail and four stern appendages (named ONR Body-1)
in the Rotating Arm Basin at the David Taylor Model Basin will
be discussed. The six components of hydrodynamic forces and
moments were measured for only one combination of roll, pitch
and yaw angle, rudder angle, sternplane angle, yawing angular
velocity and ahead speed. This is defined as an experimental grid
point. The six components of forces and moments are the most
crucial information for the design of maneuvering submarines.
The cross flows on the vertical planes at several locations down-
stream of the sail were also measured using the Particle Image
Velocimetry(PIV) technique. The measured data just described
will then be compared with the computed results and discussed.
In the conclusion, current status on the accuracy on the prediction
of the six components of forces and moments will be discussed
and future efforts to improve accuracy of prediction will also be
described.
GOVERNING EQUATIONS
The incompressible RANS equations in a steady rotating co-
ordinate system are
1
. .
OCR for page 670
= o
~ + ~[ninj +p*5ij - 2Q2(X2 +y2)5i
+2(Q x r)iui] = ~ (I,~—~ij )
(1)
(2)
~ * ~ * ~ *
As = ~coS[3 cos 1 (~ij i:)jk ;'ki )]
L
U* = j(S*)2 + (Q*)27 S* = ~S`j*SiJ*'
Q* = ~Qij*Qij*
where x; and Uj are the Cartesian coordinate and velocity S * _ 1 ( dui + Bu
componentrespectively,p*isthepressurepdividedbyaconstant t} —2 dxj
density p, r is a position vector, Q is steady rotation rate in the z-
direction, dij is the Kronecker delta, z' is the kinematic viscosity
and ~rij is the Reynolds stress tensor.
Qij* = 2(~ - ~) -
Both the k—w and the k—~ turbulence models are consid-
ered. The k—~ model by Wilcox[2] is
0dk + ~ (Uj k)—~ [(~ + ak ~t ) ~] =
_7ij ~ - p*wk (3)
bt + ~ (Uj~)—~ [(~ + ~ut) ~] =
—O' k rij ~ _ >~2 (4)
and the k—~ model by Shih[3] is
+ ~ (0jk)—~ [(U + aklJt) 7~] =
—rij ~—£ (5)
~ + ~ (Uj£)—7~ [(~ + ~t) ~] =
ClflS£—C2f2 k+~ + ~Vt(~)(~) (6)
tJt = C~f~k(k + ~) (7y
C —, 1
~ 4.0+As e
As = ~cos ¢, ~ = 3 cos-1 (~W*)
Cl=max(0.43,5+77), 9=Sk
W* Sij Sjk Ski U* = ;(S*) + (Q )
S* = ~Sij*St,*, Q* = ~Qij*Qij*
(9)
- 36ijk Qk ( 1 O)
/
Sl~f~ * ~ *
= ~ ~ij Qij
(8)
In the k—w model, k is the turbulent kinetic energy, w is the
specific dissipation rate and z/~ is the eddy viscosity given as k.
The values of model parameters are given by ~ * = 0.09, ~ = 4i,
0l = g and (Jk = ~W = 2. For the definitions of the additional
model parameters used in the k—£ model, readers are referred to
the paper by Shiht33. For practical applications, it is important
to point out that both models can be integrated right down to the
wall without using the normal distances to the wall. It is also
important to point out that the Shih's k—~ model is realizable in
the sense that negative turbulent kinetic energy components can
not be producedt4] as in the k—w model.
The forms of the two-equations in the rotating coordinate
system are the same as in the inertial coordinate system. This is
expected since it can be proved that turbulence models must be
form invariant under arbitrary translational acceleration of the ref-
erence frame and should only be affected by rotation through the
mean vorticity tensor[Si. The rotational effect is thus introduced
explicitly only in non-linear models through the mean vorticity
tensor defined in (10) in the expression for rij. For linear two-
equation models, rij is a function of mean strain tensor defined
in (9) only and rotation is not explicitly introduced. However,
rotational effect can still be introduced to the linear two-equation
models through the definitions of model parameters. Since all the
model parameters in the standard k—~ model given by equations
(3~-~4) are constant, the rotational effect is not introduced in this
model. However, the rotational effect is introduced to the linear
k—~ model, defined by equations (5~-~10), through C~ defined
by equation (8) which is a function of the mean vorticity tensor.
NUMERICAL METHOD
The incompressible RANS equations are solved by the arti-
ficial compressibility approach first proposed by Chorin t6] and
subsequently generalized and improved by Turkel [71. A finite
volume method is used. The mean flow (i.e., Eqns (1) and (2~) is
spatially discretized by a second-order accurate central difference
2
OCR for page 671
method with fourth-order accurate dissipation terms. The logic
for using the fourth-order accurate dissipation terms is twofold.
During the early stage of time stepping to reach convergence, the
fourth-order accurate dissipation terms act to suppress spurious
oscillations. But once the convergence is achieved, the contri-
bution from them to the solution is negligible because they are
fourth order accurate compared to the second-order accurate spa-
tial discretization scheme. The turbulent flow (i.e., Eqns (3~-~4)
and (5~-~10~) is discretized by one of the several upwind schemes
suggested by Yee [83. The reason for using an upwind scheme,
not a central difference scheme, to solve the turbulent flow equa-
tions is that the flux matrix is already diagonal; therefore there is
no additional cost in doing characteristic formulation. The time
stepping is based on an explicit one-step multi-stage Runge-Kutta
method to reach a steady-state solution. This approach is not just
applicable to the steady state solution but can also be extended
to solve the time dependent equations in a very simple manner.
Some discussion of this extension can be found in the papers by
Jameson [9] and Liu et al.Ll03. Several convergence acceleration
techniques including multigrid, local time step, implicit resid-
ual smoothing, preconditioning and bulk viscosity damping have
been implemented. To handle complex geometry, the multi-block
grid structure is adopted. These numerical techniques have been
implemented in a code named IFLOW, the acronym for Incom-
pressible FLOW. IFLOW is intended to be a general- purpose pro-
duction code for solving 2D, 3D, steady and unsteady problems.
The code is highly modular in structure so different turbulence
models and newly developed numerical schemes such as higher
order schemes can be easily implemented. Some special features
of the numerical schemes used will be highlighted. However,
detailed derivations will be mostly omitted.
PRECONDITIONED METHOD
The preconditioned method is developed based on a system p-1 =
of hyperbolic equations, but the idea goes back to the effort to
reduce the condition number of a matrix in linear algebra. For by- .
perbolic equations, the objective is to make the various speeds of
different wave modes more or less the same so that convergence
can be significantly accelerated. This is particularly important
when the artificial compressibility approach is adopted to solve
incompressible flows. The reason is that the sound speed, which
is one of the wave modes in the incompressible flow, propagates
much faster than the fluid speed. The result is a very slow conver-
gence as often encountered in the attempts to compute low Mach
number flows using any compressible flow code.
The preconditioned mean flow (i.e., Eqns (1) and (2~) can be
written in the conservative form as
To 1qt + FX + Gy + Hz = 0 (11)
where the preconditioned matrix PO and the three compo-
nents of fluxes F. G. and H are defined as
t1 + 7~3 2 yp-2u 73-2v 73-2W
p_l _ t1 + ~ + y'>-2U 1 + y3-2U2 )~-2UV yp-2uw
o - . (l+~+~-2V )~-2VU 1+~-2V2 y3-2vw
t1 + ~ + y)~-2W )~-2WU )~-2WV 1 + :~-2w2
(12)
q=
G=
IU F=
V
W
v
UV—Tyx
V2 +p* _ ~ ,
VW—Tyz
U
U2 +p* _
UV—[Xy
uw—[xz
?.,
UW—[ZX
VW - [Zy
w +p —[zz
where Ct,3-2 and ~ are preconditioning parameters,
rij,i, j = x,y,z, are Reynolds stresses. For mathematical anal-
ysis, it is easier to write Eqnfll) in a non-conservative form.
Neglecting the viscous terms, it can be derived as:
p lQt + Aqx + Bqy + Cqz = 0 (14)
The explicit forms of matrices A, B. C are omitted here. The
preconditioning matrix P-t is different from the previous one
Po ~ and is given by
(1 + ~y)3-2 ~y>-2n 3-2v ,JB-2W-
a>-2U
oL3-2V
oL>-2W
1 0 0
n 0 1
The condition or = ~y ensures the system of partial differential
equations is well posed. However, the implication of well posed-
ness in the case of numerical solution is not clear. In this paper,
itissetor=~=Oand
,3-2 = nlaX(|u|2, 6), ~ = 0.7. (16)
Through rather tedious algebraic manipulations, the eigen-
values and left and right eigenfunctions can be found. The max-
imum of eigenvalues is used to define a local time step. The
eigenfunctions are of no use to central difference schemes except
for establishing a nonrehecting boundary condition at far field.
The final system of equations to be solved in the conservative
3
OCR for page 672
form is
Po 1qt + FX + Gy + Hz
=
(po 1IPAIqxxx)x+
(po 1 IPBlQyyy)y + (pO 1 1PCIqzzz)z
Since the formulation is based on hyperbolic equations only,
viscous terms should be added to the fluxes F. G and H given by
Eqn (13~. The right-hand side of Eqn (17) is the fourth-orderma-
trix dissipation terms. The matrix dissipation gives the most accu-
rate solution but is less stable because a smaller amount of dissipa-
tion is added. As a compromise between accuracy and robustness,
vector dissipation is adopted in this paper. For vector dissipation,
matrices PA, PB and PC are replaced by their corresponding
radius spectra. In the curvilinear coordinates, (i,i = 1,2,3, the
maximum eigenvalue in the i-direction is given by
)\maX = - (|uiI + ~/ui2 +4,32|ai12)
ui=ue~, Pi
MULTIGRID METHOD
(18)
Multigrid is one of the most effective methods to accelerate
the rate of convergence and should be used routinely in every
production code. The approach in IFLOW follows mostly the
ideas of Brandt t11] and Jameson [12~. Several variations in-
cluding V-, W- and F-cycles have been implemented. In general,
the W- and F-cycles are more efficient but not significantly so.
More levels of multigrid are more efficient but cost a little more.
For simplicity, most computations performed with IFLOW use
three levels of multigrid in V-cycle. The multigrid method is
used routinely in IFLOW. For large scale computations on com-
plex geometries, the starting of computation is often jumpy. To
obtain a smoother start, a multigrid starting procedure is used.
Consider a 3-level multigrid computation. A 2-level multigrid
consisting of the medium and coarse grids is run for about 50
cycles and then the solution is extrapolated to the fine grid to start
the 3-level multigrid computations. In general, a solution ade-
quate for engineering applications can be achieved in 300-500
multigrid cycles. This level of efficiency is at least as good as the
best Computational Fluid Dynamics (CFD) codes available but
is far from the Textbook Multigrid Efficiency (TME) advocated
by Achie Brandt [133. The goal of TME is to obtain an engineer-
ing solution in about 10 multigrid cycles. Achieving TME is a
noble endeavor and will no doubt have a significant impact on
engineering applications of CI;D.
INITIAL CONDITIONS
The velocity relative to the fixed inertial coordinate system,
vi, is related to the velocity relative to the rotating coordinate
system, vr, by
(17) vi = vr + Q x r = (u + Q-y)ex + (v—Qx)ey + wez (19)
\ t:7, _ -l, , ~ _
where, as before, Q in the z-direction is the steady rotation
and r = (x, y, z) is the position vector with respect to the rotating
coordinate system, and ~ ex, ey, ez) are unit vectors . Since the
fluid is stationary in the fixed coordinate system, vi = 0. This
provides the following initial conditions for the mean flow in the
rotating coordinate system as
u = -Qy
v=Qx
w =0
(20)
Other initial conditions are defined as the following. The
pressure is assumed 1 everywhere. The turbulence kinetic energy
k is set somewhere between (0.01~2 and (0.0001~2. For the k—
model, w is set approximately 100 and then obtain the initial
eddy viscosity according to the relation zJ~ = ~ . For the k—
model, the initial eddy viscosity z'~ is set somewhere between 0.01
and 0.001 and the dissipation rate ~ is obtained from the relation
1~ = 0.09 ke .
BOUNDARY CONDITIONS
Only the solid wall and the farfield boundary conditions need
to be discussed. At the wall, the three components of velocity and
the turbulent kinetic energy k are set equal to zero, the pressure
p is derived by taking the the component normal to the wall of
the momentum equation (21. By neglecting the viscous term, this
component gives
~ vp* = Q2x(ex ri) + Q2y(ey ri) (21)
where n is the unit vector normal to the wall. It is noted that
in the absence of rotation, the usual boundary condition for the
pressure is recovered. The wall boundary condition of the specific
dissipation rate ~ originally given by Wilcox (p. 148 in [21) is
modified as
I, _ aOo
Ww— ~ ~W ,3 = 40
(22)
where Qw is the vorticity at the wall and aO is a constant
during calculation that can be varied from a value of 6 to 20
given by Wilcox. The choice of aO may vary the convergence rate
slightly but once the convergence is achieved the solution is about
the same. The motivation in deriving the modified wall bound-
4
OCR for page 673
ary condition (22) is to get rid of the requirement that the first
grid normal distance must be given. The presence of the normal
distance creates some difficulty in the coarser grids in multigrid
cycle. In general, two-equation turbulence models require that
the non-dimensional first grid normal distance y + be in the order
of 1. If the finest grid satisfies this requirement of a value of
1, then the y+ of the coarser grids will be greater than 4. The
consequence is that a proper turbulent boundary layer may fail
to develop. With Eqn (22), the normal distance does not appear
and the y+ of the first grid normal distance of the finest mesh
should be in the order of 1 to 3. The wall boundary condition for
the dissipation rate ~ in the k—~ model follows a suggestion by
Shiht4] as
,:_ 0.251u4
-
where u~ is the skin friction velocity.
(23)
In the absence of rotation, the gradients of the three compo-
nents of the Cartesian velocity are set to zero. In the presence of
rotation, there is a background flow described by equation (20~.
Thus the gradients of the perturbed flow field are set to zero. The
gradients of the turbulence quantities k, ~ and ~ are set equal
to zero. The pressure is obtained by a non-reflecting condition
discussed by Hedstrom [141, Rudy and Strikverda [151 and Sung
[163. This is one of the most important boundary conditions for
external flows and it can significantly affect both accuracy and
convergence. The idea is based on the characteristic formulation
of hyperbolic equations. The condition is formulated such that
the outgoing solution modes will not be reflected back into the
computational domain to corrupt the solution. For details readers
are referred to the mentioned references.
TWO-EQUATIONTURBULENCE MODELS
Several two-equation turbulence models have been imple-
mented in a general-purpose code IFLOW mentioned earlier. It
is well known that the convergence of the two-equation turbu-
lence models is rather temperamental. Two techniques have been
used in IFLOW and as a result the same rate of convergence as
using the Baldwin-Lomax turbulence model has been achieved.
The first technique is the point-implicit method for source terms.
Here, the positive part of the source term on the right-hand side is
treated explicitly and the negative part is treated implicitly. This
technique is, in fact, widely used. The second technique is to
establish a lower bound for the specific dissipation rate ~ and ~
by applying the Schwartz inequality. To illustrate the method, it
is sufficient to consider a linear Reynolds stress model
lid = Ui'Uj' = 3 kdij—C~f~*w (2Sij) (24)
where Sij is defined in (9).
By Schwartz inequality, it can be shown that
'Ui/uj/2 ~ 4k2
Taking square of the both sides of Eqn (24) gives
Fiji = 43 k2 + 2(C~fy Sky ) po' Po _ 25ijSij
(25)
(26)
A lower bound for ~ is then obtained by combining Eqns
(25) and (26) as
A>
(27)
The proportionality factor in Eqn (27) can be taken as a value
in the neighborhood of 2. Different values for this factor can affect
the convergence rate. But once the convergence is achieved, they
all give again about the same answer. The value used in this paper
is 2.1. The lower bound of ~ for the k—~ model is derived in a
similar manner and is given by
r
~ > C~f~ky/-Po
DESCRIPTION OF EXPERIMENT
(28)
A Rotating arm experiment was performed on an unclassified
generic submarine model called ONR BODY-1 in the Rotating
Arm Basin at NSWCCD as reported in Reference [11. Figure 1
shows a sting-mounted submarine for rotating arm testing. The
model was attached to the sting through two sets of block gages
which measured the longitudinal, lateral, and normal forces. The
center of gravity of the model is located at L, = 0.4646 for the case
investigated. The model is an axisymmetric body with a sail and
four stern appendages. The shapes and locations of the four ap-
pendages were deliberately made to differ from real submarines.
The sail has a NACA 0014 section with an aspect ratio of 0.27.
The four stern appendages are identical and have a NACA 0018
section with an aspect ratio of 1.2. The total length of the model
is 17.0ft and the diameter is l.55ft. The radius to the model
reference point at the center of gravity (CG) is 32.206ft. The
angular velocity of the rotating arm is 0.163 radios correspond-
ing to a linear speed of 5.26ft/s or 3.1 1 knots. Using the model
length of 17.0ft and the linear speed of 5.26ft/s as the units
for length and speed respectively, the non-dimensional angular
velocity is r' =0.53 and the Reynolds number is 8.2 x 106. The
model was oriented such that it was pitched at 2.0 °(bow up), rolled
2.1 °(to starboard) and yawed 9.5 °(to starboard). The rudder was
s
OCR for page 674
Figure 1. Sting-mounted submarine for rotating arm testing (Ref.[1])
maintained at a fixed angle of -20° for a starboard turn and the
sternplanes were set to -1.0° for bow up. Three components of
forces and three components of moments were measured in the
body-fixed coordinate system. A total of 20 laser sheets were also
taken as the model passed through a location fixed in the inertial
system using Particle Image Velocimetry(PIV). These laser sheets
provides instantaneous cross flows on 20 cross sections along the
model. Details of the experiment appear in [11.
DISCUSSION OF RESULTS
The size of the computational domain is severely restricted in
the simulation of the rotating arm experiment. The restriction is
particularly severe in the outer flow boundary. This is because the
computational domain should not cross over the center of rotation
implying that the outer flow boundary should not exceed the ro-
tating arm length, which is 1.88 ship lengths. The computational
domain used here is about one ship length in outer flow boundary
and about 3 ship lengths downstream from the tail of the model.
Three grids were used. The fine, medium and coarse grids having
a total of 2.89, 0.36 and 0.075 million grid cells, respectively. For
computer runs in this paper, the non-dimensional normal distance
to the wall from the first grid cell center, y +, is approximately 1
for the fine grid, 3 for the medium grid and 6 for the coarse grid.
A total of 12 blocks were used. The grid about the body is C-type
and the grid about the appendages are H-type. A typical grid is
shown in Figure 2. Notice that a small gap exists between body
and rudder to allow for a rudder deflection of 20 °.
IFLOW is run on Cray SV machines provided by The De-
partment of Defense High Performance Computing Moderniza-
tion Office(DOD-HPCMC) at NAVO. The code is parallelized by
Open-MP. Using 4 CPUs, the efficiency varies from 2.0 to near
perfection of 3.95 depending on how busy the machine is being
used. A typical production run takes about 500-600 multigrid
cycles. For an efficiency of 3.5 in 4 CPUS, the elapsed time is 3
hours and CPU time is 10 hours for 600 multigrid cycles using
the medium grid. For the fine grid, the run time is approximately
Figure 2. A typical grid used in the computation of flow around a ONR
BODY-1 fully appended model. (Notice a gap between body and rudder to
allow for a rudder deflection of 20°.)
Figure 3. Convergence history of root-mean-square of pressure
coarse grid
~ medium gric
—--- fire grid
100 200 300 400 500
multigrid cycles
multiplied by a factor of 8. Convergence histories of the root-
mean-square of pressure residues for three grids are shown in
Figure 3. The residues drop only about 3 orders of magnitude
and rather bumpy. The drop is smaller than the usual practice
with IFLOW.
6
OCR for page 675
Figure 4. Convergence history of force X'
o.o
-1 .o
-2.0
-3.0
do
X ~4 0
X
-5.0
-6.0
-7.0
~ =
-A,— '
;~ -
t;
_ ~ test data
A coarse grid
a medium grid
- fine grid
. _-
,
_
-80 F.,,, 1,,,, 1,,,, 1 1,,,, L . .
0 1 00 200 300 400 500
multigrid cycles
_~
Thus the quality of the grids used here need to be improved.
The convergence histories of the six components of forces(X ', Y'
and Z') and moments(K', M' and N') together with the mean
measured data are shown in Figures 4 through 9. Here, X',Y'
and Z' are the forces in the x, y and z directions of a coordinate
system fixed at the model and they are non-dimensionalized by
2 pUoo2 L2. K', M' and N' are the corresponding moments about
the center of gravity of the model at at I;, = 0.4646 and non-
dimensionalizedby 2pUoo2L3. Ingeneral,theconvergences are
quite satisfactory and the forces and moments converged to steady
values after about 200 multigrid cycles. The normal force Z '
and the pitching moment M' are very sensitive to the interaction
between the vortex shed by the sail and the cross flow in the hull
boundary layer caused by the rotation. The interaction is one of
the most important and fascinating flow physics in the flow about
a turning submarine and will be further discussed later.
Since there is only one experimental data grid point (an ex-
perimental grid point is defined in introduction), comparisons
between computations and the measured X ', Y', Z', K', M' and
N' are presented in tabular form in Table 1. At a first glance,
the comparison looks disappointing. But after further considera-
tion of experimental difficulties and uncertainties, the computed
results appear encouraging. For example, a small misalignment
in the sternplane pitch angle or rudder deflection by less than 1 °
can easily change the values of M' and 1\7' significantly. Should
this happen, the predicted M' and N' will be within experimental
error. To compensate for this uncertainty, multiple experimental
data grid point should be taken so that the trend of measured data
as angles of attack vary can be observed. Another uncertainty is
Figure 5. Convergence history of force Y'
30.0
20.0
15.0
lo
x
~ 10.0
0.0
2.0
1.0
0.0
-1 .0
lo
X -2.0
-
N
-3.0
-4.0
-5.0
1
25.0 ~ _
_
,
5.0 it_
t_-
- --- 1 .
~ 1
~ -,,,, .,,, I L
-5.u ~ ~ 100
· test data
coarse grid
0 medium grid
fine grid
~ ,
1 - 1 — ~ 1
1,,,, 1 1 1 1 1 1 1 1 1 ~ I I
200 300 400 500
multigrid cycles
Figure 6. Convergence history of force Z'
t---- 1
1
L
-
-6.0 ~ ~ 100
of
· test data
coarse grid
0 - - medium grid
fine grid
1
1 1 1
1 1 1 1 1 1 1 ~ I I I I , . , ,
200 300 400 500
multigrid cycles
the time dependency ofthe flow at such a high value of r ' = 0.53. A
small value of r' represents a mild turn with a large turning radius
and a large value of r' represents a sharp turn. For a sharp turn,
not only the possibility of instrument misalignments increases but
there is always a possibility that the flow may become unsteady.
This is a potential problem since both the experiment and com-
putation assume that the flow is steady in the rotating coordinate
7
OCR for page 676
Figure 7. Convergence history of force K'
o.s
n7
O.6
0.5
O.2
0.1
0.0
-0.1
10 n
6.0
an
on
o.o
-en
-4.n
-6.0
a,0
. ~ ~
~ _
· test dat _-
coarse grid
o medium grid
_— fine grid
._— _-
, . 1,,,,1 1,,,,1 ,,~l_-
0 100 200 300 400 500
multigr~d cycles
Figure 8. Convergence history of force M'
.
_ · test data
~ coarse grid
_ a medium end
..... ;~ fine grid
~~ 0 100 200 300
multigrid cycles
400 500
system. For these reasons, the historically measured data for r' in
the range between 0.5 and 0.6 contain very large scatters even for
a much simpler geometry like a body of revolution. More reliable
data have been obtained for r' in the range between 0.2 and 0.4.
In Table 1, the mean values of the measured
X',Y',Z',K',M' and N' are presented in the first row,
the corresponding standard deviations are presented in the
Figure 9. Convergence history of force N'
10.0
O.O
-on
lo
x -20 0
at
-30.0
-40.0
Ann
-60 0
0 100 200 300
multigrid cycles
.
5
_
_ '1
_ ~ test data _
-- I-- coarse grid
0 - - medium grid
_ ............ fine grid _
1 ~ 1 ~ 1 1 1 1 1 ~ 1 1 1 1 1_ ~ 1 1 1 1
400 500
second row and the standard deviations in percentage deviations
from the mean values are presented in the third row. The
predicted values with the coarse grid, medium grid and fine grid
using the realizable k—~ turbulence model are then presented.
The percentage deviations from the mean measured values are
listed in parentheses to give some idea of the agreement between
data and prediction. The use of percentage deviation could
create a false impression of excessive error when the mean
value is nearly zero as in the case of M' here. Finally, the
predicted results using the k—w turbulence model in medium
grid are presented. Looking at the measured data, the percentage
deviations in general fall within 10 to 20% with the exception
of Z' which is about 90% and M' nearly 700% . These large
percentage deviations are due to the small mean values of Z' and
M' at that specific experimental grid point. From historically
measured data of similar configurations, it is known that for
a slightly different experimental grid point, say yaw angle of
8° instead of 9.5°, the magnitudes of Z' and M' can increase
significantly to the orders of 10-3 and 10-4, respectively. This
is almost one order of magnitude increase in Z' and two orders
of magnitude increase in M'. This is the reason that more
than one experimental grid point should be taken as discussed
earlier. This is also the reason that the predicted values of Z' and
M' are rather sensitive to the change of grid size. With these
considerations in mind, the prediction is good. In particular, the
predictions of X', Y', Z' and K' are quite good.
It is expected that turbulence models will play an important
role in the flow around a turning submarine. From Table 1, it can
be seen that the prediction using the standard k—w model is con-
8
. . . .
OCR for page 677
Table 1. Comparison of calculated and measured Hydrodynamic
Forces and Moments of Fully Appended ONR Body on Rotating
Arm Captive Model Experiment (r' = 0.53, V = 3.11 knots, L = 17
ft. Re = 83 x106, o = -2.1°, ~ = 2°, ~ = 9.5°, sternplane rise 1°,
rudder deflection 20°
(x10-3) 1 (xlO-3) 1 (xlO-3) 1 (xlO-4) 1 (xlO-4) 1 (xlO-4)
:~:Mean~ ~ |~ ~ IO
Std. Dev. | 0.04
Std. Dev. in % | ~ 7
calc. k-e -1.53
coarse (44 DO )
calc. k-e -1.21
medium (10 % )
calc. k-e -1.29
fine and (17%) _
5.29
0.14
5
3.28
(-38 % )
5.06
(4%)
5.85
(11%)
_ _
-0.15
.-
0.07
.-
+91
.-
-0.39
(160%)
-0.33
(120%)
-0.11
_ (-27%) _
0.38
0.037
20
0.07
(-82 %
0.48
(26
0.53
(40 %
-0.05
0.176
690
2.32
(4700% )
0.85
(-1800% )
1.54
(-3180% ) .
- 1.50
0.103
14
-1.74
(16%)
-1.57
(5%)
-2.90
(93%)
calc. k-~ -1.57 3.73 0.13 0.29 2.44 2.61
medium (43 % ) (-29 % ) (-187 % ) (-24 % ) (4980 % ) (274 %
sistently not as good as the prediction using the realizable k—£
model. To understand the problem, a simple body of revolution
called SUBOFF at incidence is considered. It is clear from Figure
10 for the prediction of Z' that all models predict well at angles
of attack less than 4° and then deviate from the measured data
at higher angles. It can be seen that the realizable k—£ model
performs much better at higher angles of attack. The prediction
of the moment M' is shown in Figure 11. In general, all models
perform well but the realizable k—£ model still performs bet-
ter. Better prediction of M' by the realizable k—£ model is not
surprising since M' is essentially derived from the distribution
of Z'. A better prediction of Z' gives a better prediction of M'.
A conclusion can be drawn from this investigation. For a turbu-
lence model to perform well in the prediction of the flow about a
turning submarine it is essential that it can predict accurately the
force and moment of a bare hull at incidence.
One of the most important flow physics about a turning sub-
marine is the interaction between the vertical flow shed from the
sail and the cross flow on the hull boundary layer caused by ro-
tation. This interaction results in a bow pitch-up moment of the
submarine. For a barehull without a sail, there is no vertical flow
shed and the flow about the body is symmetric. As a result, there
is no normal force Z' nor pitch moment M'. In the presence of a
sail, a vertical flow is generated. The flow direction of the cross
Figure 10. Prediction of normal force Z' vs. incidence angles
(94x48x64, Re= 1 4x 1 o6 )
· SUBOFF TEST: +a
o SUBOFF TEST: -a
8 00 ~ — SUBOFF IFLOW(k-co) _
<~ . - SUBOFF IFLOW (k-£)
5.00
To 4-00
AX
N 3.00
2.00
1.00
0.00
~ no
=
.
~I,, ~,,,1
i/ ,)
/ /
(/ ~
id'
I/
/
- ..vvt ~ 1 2 16 20
a (deg.)
Figure 11. Prediction of pitch moment M' vs. incidence angles
(94x48x64, Re= 1 4x1 o6 )
.~ an
3.00
2.50
2.00'
o
_ ~
X ~
~ .
1.50t
1.00
O.50
0.006
· SUBOFF TEST: +a
a SUBOFF TEST: -a
SUBOFF IFLOW(k-~)
I ~ 1 SUBOFF~IFLOW (ken
~ 1 1 1
~ turf''
~ ~ i'
At''
1~, · , , , 1 , , , 1 1 1 1
4
~ 1
~~ 1
8 12 16 20
a (deg.)
flow of the shed vortex is opposite to that of the cross flow in
the hull boundary layer. Cancellation of these two cross flows
in opposite directions results in a lower velocity consquently a
higher pressure on the upper hull in comparison to the lower hull.
This downward normal force in the aft body causes the bow of
9
OCR for page 678
Figure 12. Longitudinal distribution of normal force Z',r',=0.4, ,3=9° for
both bare hull and bare hull with sail
0.3
0.1
_ n lo, __
X v.v~ _ .
-0.1 _
-0.2 _
I I I 1 1 1 1 ~ I ~ .
. ~
. ~
. .
. ;~-
,...
_ ~
or—
or. .-1 it., . rid ~
Cal wr ~ 1 Or Its ~
_ 1~ 1 1~
O barehull
- sail
~ 1
-on ,,,, I I,,,, ,,,, I,,, ,,,,,, I,,,, I I,,,, I,,,, I I,,,, I, ~ ~ .
t.0 0.1 . 0 3 0.4 0 5 0 6 0 7 0 ~ 0.9 ~ 0
the submarine to pitch up. This interesting flow physics is shown
in Figures 12 and 13 where the longitudinal distributions of Z'
and M' are shown respectively for a barehull and a barehull with
a sail in rotation with A' = 0.4 and ,B = 9°. For a barehull, both
Z' and M' are zero along the body. In the presence of a sail, a
downward force(+Z') is generated on the hull aft of the sail as
shown in Figure 12 resulting in a bow pitch-up moment as shown
in Figure 13.
For the ONR Body-1 fully appended submarine, the sail is
located between I;, = 0.2 and 0.3. Figure 14 compares the mea-
sured and computed crossbow velocities generated by the vortex
shed by the sail at A, = 0.34. The location of I, = 0.34 is approx-
imate since the laser sheets of the PIV techniques were taken in
the fixed inertial frame while computation is done in the rotating
frame. Thus the computed flow field needs to be projected onto
the laser sheet plane for the comparison. It can be seen that the
vortex shed by the sail is predicted well although the predicted
location of the center of vortex is slightly lower and further out-
board than the measured one and the strength is not as strong.
The interaction between the vortex shed by the sail and the hull
boundary layer resulting in a downward force producing a sub-
marine nose-up phenomenon is shown in Figure 1 Sa. At A, = 0.34
where the vortex has just been formed, the upper hull is domi-
nated by the sail-root high pressure field. As the vortex travels
downstream, the interaction becomes stronger and a downward
normal force is generated as shown at I;, = 0.48 in Fig 15b. This
observation is consistent with the longitudinal distribution of Z'
shown in Figure 12.
Figure 13. Longitudinal distribution of pitch moment M',r',=0.4, p=9°
for both bare hull and bare hull with sail
0.5
0.3
0.2 _
0.1
O Or ~ _ ,,.= ,
. ~ ~a,
-0.1 _^
-0.2 - _
-0.3
-0.4
-o.5o 0 0 1
. == _
ix
_~_
_
_ ~
Jam Illl~llr!
—~ barehull
~ sail
~ 1
the
1
t1111111111!111~1
I.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
CONCLUSIONS
A numerical procedure for the prediction of the flow around
a turning submarine has been developed. The procedure is based
on solving the incompressible Reynolds-averaged Navier-Stokes
equations coupled with a standard k—~ turbulence model and
a realizable k—~ turbulence model in a steady rotating coordi-
nate system. Computed results are compared with the measured
data obtained from an unclassified submarine called ONR Body-
1. ONR Body-1 consists of an unclassified submarine-like body
with a sail and four stern appendages. The non-dimensional an-
gular velocity is r' = 0.53 and the Reynolds number is 8.2 x 10 6.
The model was oriented such that it was pitched at 2.0 ° (bow up),
rolled 2.1°(to starboard) and yawed 9.5°(to starboard). The rud-
der was maintained at a fixed angle of -20 ° for a starboard turn and
the sternplanes were set to -1.0° for bow up. Only one data grid
point listed above was taken. Three components of forces(X ', Y'
and Z') and three components of moments (K', M' and N') were
measured. Laser sheets using Particle Image Velocimetry(PIV)
techniques were also taken to track the vortex shed by the sail.
Comparing with the measured data, the prediction of the forces
is in much better agreement with data than the prediction of the
moments. The reason for the poor prediction of moments can not
be ascertained since only one experimental data point was taken
as explained in the discussion. The predicted crossbow of the
vortex shed by the sail at L, = 0~34 compared well with the PIV
laser sheet. The interaction between the vortex shed by the sail
and the hull boundary layer results in a downward force in the
aft hull and consequently a pitch-up of the submarine nose. This
10
OCR for page 679
Figure 14. Crossflow Velocities at I, = 0.34
(right behind the trailing edge of sail)
-0.6 -
.2 0 0.2 0.4 0.6 0.8 1
Y(~)
Computed Crossflow Velocities for Plane #7 (x/L = .341 )
Figure 15. Calculated Pressure Coefficient
(a). at =~ = 0.34(right behind the trailing edge of sail)
(b). at L. = 0.48
-
-
-
-0.4 t
~ . _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ .
. . . _ _ _ _ _ _ _ _ _ _ _ _ _ _ .
a, . . . _ _ _ _ _ _ _ _ _ _ _ _ _ .
in, . . _ _ _ _ _ _ _ _ _ _ _ _ .
x. . _ _ _ _ _ _ _ _ _ _ _ _ .
5.26 fats Reference vector: W:: ~ - - - -:
~ \. :::::~::::::
~ . . . .
`... . .
~ · , ~ ~ , · , ~ , , , ~ , , , ~ , , , lo-,--,,, 1,-.-
-0.2 0 0.2 0.4 0.6 0.8 1
y (ft.)
important flow physics has been illustrated and explained with
the computed flow field.
As expected, turbulence models play an important role in the
accurate prediction of the flow around an maneuvering submarine.
A standard k—~ model and a realizable k—~ model were used in
this investigation. For a simple body of revolution at incidence,
the realizable k—~ model gives a much more accurate predictions
of lift and pitch moment than the standard k—~ model. It is not
surprising that the realizable k—~ model also gives a consistently
more accurate prediction of the flow around a turning submarine.
Since this investigation is the first time that the prediction of
the flow around a fully appended submarine in rotation has been
attempted based on the RANS solutions and the results compared
with the measured data, several improvements in numerical meth-
ods are needed for more accurate prediction. Since the accurate
tracking of the vortex shed by the sail is important for accurate
prediction of Z' and M', higher-order spatially accurate schemes
should be used. The conventional second-order spatially accurate
-
-
-
-
-2
y (ft.)
Cp
0.04
0.03
0.02
o.ol
Coo
-0.01
-0.02
-0.03
-0.04
-0.05
-0.06
-0.07
-0.08
.os
-0.10
schemes used in the present investigation produce too much nu-
merical dissipation which rapidly diffuses the shed vortex. Higher
orcer schemes have smaller numerical dissipation therefore vor-
tices will not be diffused as rapidly as second-order schemes.
Local refinement technique is also very effective in tracking the
vertical flows. It is also clear from this investigation that a grid
with 3x 1 o6 cells is far from sufficient for accurate prediction. The
majority of the grid cells are distributed near the wall because con-
ventional turbulence models require the first y + be approximately
1. There will be a great saving of grid cells if a new wall function
method can be developed so that a y+ in the range from 30 to
100 can be used. The new numerical methods just mentioned
including higher-order spatially accurate schemes, local refinem-
nent and new wall function method will be used for the prediction
of the flow around a turning submarine in the next attempt.
11
OCR for page 680
ACKNOWLEDGMENT
This work was partially funded by the Office of Naval Re-
search, Code 333, monitored by Dr. Patrick Purtel and Dr. Ki-
Han Kim and partially funded by the ILIR program at David
Taylor Model Basin monitored by Dr. John Barkyoumb. Com-
puter resources provided by The Department of Defense High
Performance Computing Modernization Office(DOD-HPCMC)
at NAVO is gratefully acknowledged. Useful discussions with
Mr. Tom Moran and Dr. Jerry Feldman are gratefully acknowl-
edged.
REFERENCES
1. Fu, Thomas, C., Paisan Atsavapranee and David E. Hess,
"PIV Measurements of a Turning Submarine Model (ONR BODY
1) Part 1: Experimental Setup", NSWCCD-50-TR-2002/019,
2002, Hydromechanics Directorate, Research and Development
Report.
2. Wilcox, D. C., Turbulence Modeling for CFD, DCW In-
dustries, Inc. CA, 1993
3. Shih, Tsan-Hsing, Jiang Zhu, William Liou, Kuo-Huey
Chen, Nan-Suey Liu and John Lumley, "Modeling of Turbu-
lent Swirling Flows", NASA Technical Memorandum 113112,
ICOMP-97-08; CMOTT-97-03, August 1997.
4. Shih, Tsan-Hsing, Jiang Zhu and John Lumley, "A New
Reynolds Stress Algebraic Equation model", NASA Technical
Memorandum 106644, ICOMP-94 -15; CMOTT-94-08, August
1994.
5. Speziale, Charles G.,"Turbulence Modeling in Noniner-
tial Frames of Reference", Theoretical and Computational Fluid
Dynamics 1, pp 3-19, 1989.
6. Chorin, A. J., "A Numerical Method for Solving In-
compressible Viscous Flow Problem", Journal of Computational
Physics, vol. 2, 275, 1967.
7. Turkel, E., "Preconditioned Methods for Solving the In-
compressible and Low Speed Compressible Equations", Journal
of Computational Physics, vol. 72, 277, 1987.
8. Yee, H. C., "A Class of High-Resolution Explicit and Im-
plicit Shock- Capturing Methods", NASA Technical Memoran-
dum 101088, February, 1989.
9. Jameson, A., "Time Dependent Calculations Using Multi-
grid with Applications to Unsteady Flows Past Airfoils and
Wings", AIAA 91-1596, June 1991.
10. Liu, C., X. Zheng and C. H. Sung, "Preconditioned Multi-
grid Methods for Unsteady Incompressible Flows", Journal of
Computational Physics, vol. 139, 35-57, 1998.
11. Brandt, A., "Multigrid Techniques: 1984 Guide, withAp-
plications to Fluid Dynamics", 1984,191 pages, ISBN-3-88457-
081 - 1; GMD-Studien Nr 85; Available from GMD-AIW, Postfach
1316, D-53731, St. Augustin 1, Germany, 1984.
12. Jameson, A., "Multigrid Algorithms for Compressible
Flow Calculations", Lecture Notes in Mathematics, No. 1228,
Proceedings of the Second European Conference on Multigrid
Methods, Cologne, pp. 166-201, October 1-4, 1985.
13. Brandt, A., "Barriers to Achieving Textbook Multigrid
Efficiency (TME) in CFD", NASA/CR-1998-207647, ICASE In-
terim Report No. 32, April 1998.
14. Hedstrom, G. W., "Nonreflecting Boundary Conditions
for Nonlinear Hyperbolic System", Journal of Computational
Physics, vol. 30, pp. 222-237, 1979.
15. Rudy, D. H., and J. C. Strikwerda, "Boundary Conditions
for Subsonic Compressible Navier-Stokes Equations", Comput-
ers and Fluids, vol. 9, pp.327-338, 1981.
16. Sung, C. H.,"An Explicit Runge-Kutta Method for 3D
Incompressible TurbulentFlows", DTNSRDC/SH -1244-01, July
1987.
17. Jameson, A. and L. Martinelli, "Mesh Refinement and
Modeling Errors in Fluid Simulation", AIAA Journal vol.36, No.
5, May 1998.
12
OCR for page 681
DISCUSSION
Steven Turnock
University of Southampton, United Kingdom
Would the authors like to comment on whether
the impressive calculations they have carried out
yet have the grid resolution necessary to achieve
reliable results that can be used for design?
Experience at Southampton suggests that
typically four or five million cells on a multi-
block structured mesh are required to accurately
capture lift and drag on a single low-aspect ratio
control surface. A simplistic analysis would
therefore suggest at least an order more cells are
necessary. Have the group considered moving to
an unstructured mesh (hybrid cells) which should
significantly reduce the computational effort.
AUTHORS' REPLY
The number of grid cells required depends on the
degree of solution resolution demanded. Due to
complexity of the experiment, we set the goal of
achieving prediction within 10 to 20 % of the
measured data in this first attempt. Apparently
the 3 million grid cells used in this investigation
is insufficient to achieve grid independent
solution. We agree with the discusser that a far
greater number of grid cells is needed. Our next
attempt will still be based on the structured grid
but with more sophisticated numerical
techniques including local refinement, higher-
order spatial discretization schemes and a new
wall function method for the turbulence models.
Use of these techniques should make it feasible
to achieve high solution resolution with much
less grid cells. Also the hybrid MPI/Open-MP
technique will be used to maintain optimum
parallelization efficiency even with a large
number of processors in shared memory
computers. After the code based on the structure
grid has become standard production code, we
will consider the approach based on unstructured
grid.
Representative terms from entire chapter:
turbulence models
