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 295
Computation of the Flow past Shiplike Hulls
J. Piquet and M. Visonneau
ENSM
Nantes, France
Abstract
The present work first details the implementation
problems of a numerical procedure for the solution of the
Reynolds-averaged Navier-Stokes equations in boundary-
fitted coordinates. The procedure is validated using well
documented experiments on the HSVA tanker and the
comparisons demonstrate the ability of the method to
predict ship stern flows.
I. Introduction
For a long time, wind tunnel or towing tank
experiments have been the own practical way to get
information on ship performances because of the
restrictions of inviscid or boundary layer models. With the
development of computers, it becomes possible to obtain
numerical information with an answer turnover and a level
of details which are hardly possible with experiments. The
prediction of shiplike afterbody flows, which is considered
hereafter, is necessary for the knowledge of viscous
resistance and, more importantly, for propeller and control
surface design as both are immersed in regions where
viscous effects are significant. Classically, such problems
were handled solving first a potential problem with panel
methods. The inviscid flow field on the body was used as
a boundary condition for thin boundary layer equations
the solution of which allowed the determination of viscous
corrections. Unfortunately, it was demonstrated during the
SSPA-ITTC workshop on ship boundary layers [1] that
such a procedure failed in the stern region because of the
rapid thickening of the boundary layer in this zone. The
limited success of generalizations of thin boundary layer
equations involving high order corrections, see e.g.
[21~33~4], was subsequently demonstrated so that the
tendency towards computing the full solution of the
Navier-Stokes equations became stronger and stronger for
Increased computer ressources became more and more
available at continuously decreasing costs. However,
significant obstacles to progress in the treatment of such
threedimensional flows remain, the most important being
the difficulty in describing with enough resolution the
geometry of the physical domain in which the body is
immersed, even if, as here, no free surface effect is
considered. Another important difficulty is the need of
reliable data to validate computations insofar as the
comparisons between the experiments and the
computations involve systematic interpolation procedures.
Table 1 summarizes the main characteristics of the
presently available methods which share several common
295
features as well as differences among which the following
appear important. (i), used coordinates and grid generation
techniques, adaptivity of the mesh; (ii), type of master
equations and retained approximations -all of which depart
from thin boundary layer assumptions-; (iii), turbulent
models, how are they used and with what specific
boundary conditions; (iv), discretization schemes
pressure velocity coupling methods, iterative modes
solution methods used for the linear systems; (v),
acceleration means; (vi), initial and boundary conditions.
All the methods (including this one) that have been
published, but [51~61~7l, rest on the use of Reynolds
averaged Navier-Stokes equations (hereafter called
RANSE) presented in §2 and written in a body-fitted
curvilinear coordinate system. Except works
[83~91~101~1 1], only the independent variables are
transformed from the physical coordinates x (horizontal
along the axis of the body), y (horizontal, orthogonal to the
axis of the body), z (vertical) to the curvilinear coordinates
~ = (\, A, () while the dependent variables are retained in
the physical domain: the velocity components are the
components in a Cartesian frame with the noticeable
exception of t121~131~14l, where cylindrical velocity
components are preferred. In order to allow more easily
comparison with experiments, the coordinate planes
\=const. are identified with x=const. planes. In the
following (§2.2), ~ will refer to the coordinate away from
the wall such that the equation of the wall is ~ = 1; ~ will
refer to the girthwise coordinate. The other alternative
[81~91~10~11] uses contravariant coordinates. This
simplifies the continuity equation but introduces Christoffel
symbols in the momentum equations, the full treatment of
which is cumbersome (usually a lot of terms are omitted)
and either time or storage consuming.
Now, considering the master equations, ie. those
which are discretized numerically, they usually involve a k-
£ type of closure with a wall-function approach which
avoids numerical troubles often associated to the sink-
source terms in the £ equation while saving storage for the
nodal unknowns. Some works depart from this trend by
the use of either a classical mixing length model [81~91 or
of a subgrid scale model [ 151. Considering more
sophisticated Reynolds stress models does not seem to
improve significantly the results [161. In the following, a
k-£ model is retained but the wall function approach is
avoided.~§2.4) Only the £ equation is bypassed close to the
wall. For a significant increase in n~,mencal troubles and in
OCR for page 296
~ ~ o ~ ~ 2 · ~ c
o o ° o ° ° °—C
- o 4J ~ - > o ~ c
-
_ ~ ~ ~ ~ _ ~ s ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ E
-
-
-
c
O ~ ~ ~ 0 0 C
~ `~ C~ ° cn c c °, ° c c E ~ E ° E
c 0 ~ r 2. ~ E 0 o E . w c E ° C c E c ~ C ° c
a, ~ ~ ~n ~ ~ ~ ~ d, O > ~ ~ O > ~ ~ L} ~ ~
~ .~ c ° ° 0 · O (' c c° c
G 0 1~ ~ _
O e/ ~ ,,. O
~ ~ ~ O E E
~ O
_ . . _ 2 _ _
O ~ O E ~, 0 0 °
~ n e c · ~ ~ ~ O
>. a~ i~ ~ ~ ~o ~ Q w
I ~q t~ ~ .. _ u~ Ci) ~
a 0 J v ~, 0
`, tO Q > C
~ . b E E . E.
O) a, ~ ~ ~ ~ ~
b ~ Y Y Y ~ I
a ~n _ ~~,
~D
n
296
w
D
-
a,
o
0
C)
-
-
-
C'
S
as
w
D
a,
OCR for page 297
computing time, the delicate problem of the
threedimensional orientation of the log law specification,
dealt with only in t14], is avoided, as well as significant
overestimations of the velocity profiles, close to the wall
in the propeller disk region.~171~183.
From the closed RANSE, the master equations are
almost invariably obtained -except [143~193- using the so-
called partially parabolic approximation in which
longitudinal \-transport due to viscosity and turbulence is
neglected (de, if ~ is the advected quantity, only the ¢)
terms are dropped). This approximation is usually
developped for viscous flow fields in which pressure is the
dominant transmitter of influence in the upstream direction.
In certain conditions, it does not forbid the presence of -
seldomly present- longitudinal flow reversal (along \)
which may occur close to transom sterns or when a
propeller is in operation. In this work, the master equations
do not involve such an approximation for reasons which
will be made clear; therefore the fully elliptic equations are
used.
Turning now to the numerical aspects which are
studied in §3, a general prerequisite is a method for
generating the body-fitted gnds. The most common choice,
Including ours already detailed in §2.3, rests on elliptic
solvers [2~)] because they are easy to use and because they
allow a priori controls of the quality of the grid close to the
body. A possible alternative has been proposed in
[81~9] [101~19] [21]where twodimensional ends, generated
via conformal transformation designed, for instance, from
the Schwarz-Christoffel theorem, are stacked in the ~
direction. Such grids build-in transverse orthogonality to
the price of distorsion problems in the ~ direction, close to
the aft sections. Other algebraic techniques have also been
used [22] [23] [24] [25] [26] but they appear often
abandonned.
The most important aspect of the numerics lies in the
way the incompressibility constraint is enforced. This may
be done in three ways. A first category of methods
includes those where the pressure-velocity coupling is
simulated with a suitable modification of the continuity
equation. Compressible flow methods used for low Mach
numbers belong to this category, although severe
difficulties are usually found in the limit of vanishing Mach
numbers. The artificial compressibility method [27],
systematically used in some institutes (see e.g. [283), and
applied to ship flows in [293~30], gives also rise to
convergence problems and needs a careful selection of
numerical parameters such as the pseudocompressibility
factor.The second category includes methods in which the
continuity equation is satisfied identically at each iteration.
Apart from unsegregated methods [31], such methods
follow the boundary layer practice in which the flow
domain is swept from upstream to downstream, implicit
differencing of momentum being used for marching
stability. Upstream influence through the pressure field has
to be accounted for, and this is done introducing some
form of forward differencing for the streamwise pressure
gradient which allows departure-free behavior [32~. For
ship applications, such methods have been mainly
developped in the framework of the partially parabolic
approximation in [81~91~101. In this work, we have
followed (§3.2) the most common so-called pressure
correction methods in which the pressure velocity coupling
is solved iteratively, updating successively the velocity
variables In the momentum equations and the pressure in a
pressure equation, in such a way that solenoidality is
satisfied at convergence. In the ship applications, both the
advection and the diffusion of momentum are implicit in
time; the updating procedure being handled on each
equation by the well known "SIMPLE" t33l, "SIMPLER"
[34], or "PISO" [35] methods. One of the dominant
aspects of such techniques lies in the way the pressure
equation Is handled as the pressure solver controls the
global convergence of the method (and therefore, its
robustness on fine grids), and as it generates the greatest
part of the computational cost. The present work strongly
differs from similar computations [121~131~14] by the
characteristics of the pressure and pressure correction
solvers which become critical as soon as grid clustering is
strongly increased by the removal of the the wall function
approach.
Such modifications are necessary for the following
reason. Because only the independent variables are
transformed, the directions of the velocity components
differ from the direction of the curvilinear coordinates. As
a first consequence, the projection phase of the pressure-
velocity coupling does not connect always the driving
component of the pressure gradient to the corrected
velocity component . Because of the strong grid clustering
the stiffness of the "threedimensional" pressure matrix is
so much increased that standard relaxation techniques
usually retained do not converge anymore. This argument
will be clarified in (§3.3) As a second consequence, the
possibility of pressure odd-even oscillations is recovered
even if the grid is staggered (the most common choice, see
Table 1~. They are classically avoided, as when, as here,
colocated grids are used, with weighted averages of
coefficients defined during the projection phase. For the
methods belonging to the third category, the discretization
of the advective terms rests often on strongly diffusive
classical hybrid schemes. In this work, we use (§2.1) the
so-called finite-analytic method of [37] applied, for
instance, in [133~143~151~181.
This quick review of methods would stay incomplete
without consideration of initial and boundary conditions.
They will be specified in §44, when the results are studied.
While the iterative procedure is usually initialized from
rest, from uniform flow or from potential theory, the
considered domain boundaries usually extend from a
midship station to several lengths behind the body.
[51~61~71~153~381~39] depart from this choice on model
ships and start from a uniform flow field far upstream. In
this work, we follow the common practice which considers
that the inlet flow at midship stations can be characterized
by some simple integral parameters and that the crossflow
at midship can be neglected. These integral parameters are
obtained either from experiments or from thin boundary
layer assumptions used upstream. This allows to save
computer ressources and to avoid transition problems in
the fore part of ship models. The flow domain extends
usually to some ship lengths downstream and to one ship
length away from the axis of the ship. Symmetry
conditions are retained on the vertical plane -the so-called
keel plane- (only one side of the ship is handled) as well as
on the horizontal plane -the water plane-. Only the so-
called double-model problem is therefore treated, except in
[51~61~73~151~38 1 where free surface effects are
considered. Related to the flow domain choice, is the wake
topology which involves, the most often, a singularity line
either of the polar type or of the parabolic type. In the
former case, transverse grids are topologically equivalent
to concentric circles Art = const.) and their radii if=
297
OCR for page 298
coasts. In the latter case, the singular line starts from the
bottom of the ship, in its vertical plane of symmetry into
the wake, defining a line of separation between the upper
part of this vertical plane of symmetry (the ''fictive rudder")
which is a (symmetry) surface ~ =1 that extends the wall
surface (rl=l) in the wake, and its bottom part which is a
(symmetry) surface (=2.
Among the available experimental information which
has been extensively reviewed in t40][41], we shall mainly
focus (§4) on the HSVA tanker [44]~451. For this case,
very detailed mean flow measurements are provided in
several x = const. planes; this minimizes interpolation
processing. Moreover, the data are free from blockage
effects because a wind tunnel with slotted walls was used.
Uncertainties on the inlet conditions appear therefore to be
an admissible penalty. The paper is outlined as follows. §2
presents the equations and the main features of the used
method; §3 focusses on the numerical aspects; §4
discusses the results and presents significant comparisons
with the HSVA experiments.
2. The Master Equations
2.1. The Primitive Form.
The full RANSE (2.1 ) for the mean velocity field U
and for the mean pressure deviation P
V.U=O, DU/at + (U.V)U + VP = Re~1V2U - V.uu (2.1)
need a closure assumption for the Reynolds stresses -uu
which is taken under the isotropic form (2.2)
uu = 2k 1/3 - vT [V U + VTU] (2.2)
where VT = C~k2/e is the turbulent viscosity, k is the
turbulent kinetic energy and ~ the so-called dissipation of
turbulent kinetic energy. k and ~ are provided by the
standard following transport equations (2.3) and (2.4)
Ok/3t + U.Vk = G - ~ + V.[ Rk-1Vk ~ (2.3)
8e /Ot + U.Ve = Cal 1 ~ G/k - Cal 2/k +V IRE - 1 Ve ~ (2.4)
G is the rate of production of turbulent kinetic energy
uu: VU, Reff is the effective Reynolds number defined
by Reff -1 = Re~1 + vT while Rk-1 = Re~1 + vT / ok and
Rat -1 = Re~1 + vT/c;~ . Unless specified, model
"constants" are given by their standard values:
Cp =.09,C~ 1 = 1.44,C~ 2 = 1.92,c~k = 1,c,~ = 1.3 (2.5)
2.2 Physical space and computational space
2.2.1 The problem . The physical domain Q. where the
flow is studied, is bounded by the ship hull B of surface
DB, by the outer boundary £, by the inlet station U. by the
outlet station D and by the vertical V and horizontal H
planes of symmetry. Equations (2.1 to 4) which classically
involve the Cartesian velocity components
{U°C}-{U,V,W} are partially transformed from the
Cartesian rectilinear coordinates {xa}-ix, y, z} in the
physical space, to the curvilinear coordinates hi - {k,ll,`}
in the so-called computational space. The coordinate
transformation is designed in order that the boundaries OB
and £ become ~ = const. surfaces, that the boundaries U
and V become ~ = const. surfaces and that the planes of
symmetry become ~ = const. surfaces. The domain Q
298
becomes a parallelepiped in the computational space in
which the discretization consists in stacked unit cubes of
sides AX = A~ = A: = 1. Each unit cube of the
computational space is a curvilinear hexahedron in the
physical space, the "sides" of which are measured by the
modulii of the covariant vectors al = ar/aXi.
2.2.2.Useful Relationships. The transformation
necessarily involves byproducts from the basis aj. Of
particular interest are [201:
(i) The area vectors bJ = ajxak (i,j,k in cyclic order) which
measure the oriented area of a small surface of unit sides
. .
along ~J and ~ ~ on a \~ = const. surface in the
computational space. bl appears as built with two small
triangle-likes surfaces in the physical space.
(ii) The Jacobian J of the transformation from the
computational space of the coordinates {\i} to the physical
space of the coordinates {xa). J measures the "physical"
volume of a parallelepiped of unit sides in the
computational space; this parallelepiped appears as an
hexaedron-like volume in the physical space.
(iii) The covariant and contravariant metric tensors
.. . .
gij = ai.aj; glJ = bl.bJ/detgij where the determinant of g
is the square J2 of the foredefined Jacobian.
An important relationship to be used is the following
simple restatement (2.6) of the chain-rule derivative
formula:
a~k/8xa=J-1 teak (2.6)
which allows the computation of the standard following
operators:
. . . . .
V.U = J~ D(JUl)/~1 = JO a[balu(x]l3~l (2.7)
where UP and [V°]a are the physical components in the
Cartesian orthonormalized physical space.
[V¢~]a = J~ 1 teak DO fink (2.8)
V o = div [g akk ~ = J akl [ Jg 3\k ~
(2.9)
2.2.3.The Grid generation method. Because the values of
the curvilinear coordinates hi are specified on am, the
values of ki inside Q result from a boundary value problem
in which the dependent variables are the {bit and the
independent (unknown) variables are the {xa3, the
correspondence between {x(X} and (\i} being one-to-one.
If (2.9) is written for V2x = V2y = V2z = 0, the resulting
partial differential equations (2.10) may be used to generate
the coordinates:
2
V2`q, glk a ~ + 1 [ g ] ~ with Ax —(x,y,z) (2 10)
aX aX aX aX
Because the coefficient of the first order derivative,
in (2.10), is the laplacian of the curvilinear coordinate:
k 2 k 1 a [Jg ] ik (2. 1 1 )
aX
The common practice is to specify it as a grid clustering
control fk = v2\k. In this case, (2.10) becomes (2.12)
OCR for page 299
vie, lk a ~ + fj at (2.12)
afar ax
Now, if v2xa is computed with (2.11) rather than with
(2.10), fk must be equal to the rhs. of (2.12~. Therefore,
when (2.12) is used to compute the diffusion terms (for
¢=U,V,W,k,~ ), the iterative method for solving the mesh
must be solved to convergence. In all that follows, only
the special case ~ = ~ (x) is considered so that
fl= gllX~/x~__2agll t121~46] results from the
assumed longitudinal grid distribution x(~. Thus, only the
equations for y and z need to be solved, once control
functions f2 and f3 have been specified.
~ i, ..
f~andf~ are connected to giJ in such a way that
2b=-f2/g22 - Sr~/Sr~ and 2c=-f3/g33 - 0~/~; (as well
as 2a) are frozen during the iterative use of (2.12) for
¢=y,z. They are prescribed according to (2.13) whore S is
the curvilinear abscissa along the ~ = cte lines and ~ is the
polare angle in the x = const. correction.
S
f =FA(11,`) for \<\A; FA= ~g S I OVA (2.13a)
S
f2=FB(q,`) for \2^B; FB= g22 S~'1 1 \=\B (2.13b)
f =FC(\,11,`) for ~A<~<\B ;
(\—~ B FAR-—~A)FB
CF=
~ -\
A B
a,b,c are fixed from an initial guess taylored to the
problem (see §3~. The numerical solution of the grid
equations uses exponential schemes similar to (3.6) -see
infra- in both directions; it considers crossed second-order
derivatives as source terms which are discretized, as well
as metric components, with centered schemes. The
solution of the resulting linear systems endly results from
SLOR alternate ~ and ~ sweeps. The grid generation
procedure thus defined has two aims: it first allows a
correct clustering of grid points close to the wall, and, as
important, it smoothes the grid. So, zero machine
convergence of the procedure is not necessary. The f2 and
f3 functions that result from (2.13) after SLOR sweeps are
not used in the Navier-Stokes solver. Rather, the discrete
equations corresponding to the Poisson operators are
considered as a 2x2 linear system whose solution gives the
needed f2 and f3.
2 3 The Flow Equations
Apart from the continuity equation -see (2.7)-, the
transport equations for the mean momentum, for k and ~
can be written under the following Master Equation (2.14)
form:
gl 1~:,~` + g22~r~ + g33q>~` = 2A¢~:
+ 2Bq>~, + 2Co'DX + R¢Ot + SO (2.14)
where
2A`p, = Jeff l b3 as ~ - f3
2B~ = Jeff tb2 At ~ _ f2
2C<, = J—Ebl a¢] -
Red = a¢Reff; S =s—2 lgl241, + gi3~' + gads
and
a =a<, U-J Ebb VT ;~+b~ VT r,+b~ VT `] (2. 1 6a)
tab—alp, V-J l:b2 VT ;~+b2 VT T~+b2 VT `] (2. 1 6b)
a =a~ W-J [b3 VT A+b3 VT rl+b3 vT,`] (2.1 6c)
The index ~ refers to any of the advected quantities,
namely U. V, W. k, e; the indices X, it, ~ refer to the
partial derivatives; coefficient ad, is equal to 1 except
ak = Cork, an = o~ . Source terms so for U. V, W. k, ~ are
obtained from (2.17) written with Cartesian derivatives
indexed x, y, z.
SU=Reff [PX- gradVT. a
au
SV=Reff [py- gradvT · ]
By
(2.13c) au
SW=Reff [pz- gradvT · I]
(2.17a)
(2. 17b)
(2.17c)
Sk=~C7kReff (Gee) (2.17d)
is ~ e2)
So ~C,~Reff Decal k G - Cur 2 k ) (2.17e)
299
2.4 The Turbulent Model
Two reasons for the eviction of the wall function
approach have been already mentionned in the
introduction, namely the problem of specifying the
directional dependence of Up and the practical
overestimation of velocity profiles. A third reason is the
noticed strong dependence of the pressure level on the
location of the grid points close to the surface. As, at least
the ~ equation is avoided close to the wall, it is necessary
to specify the closure assumptions to supply. For the
HSVA case (§3), a one equation is used in the inner zone
in order to make the closure problem free from the
estimation of Up (whose sign may vary in the domain).
The closure assumptions (2.18, 2.19) follow the choices
of [47]
~ = k3/2/1~ : vT = Cook 1~.~ (2.1 8a,b)
where the scales of the two equations differ [481:
1p = C1 n ~ 1 - exp (-Ret/ k n / April;
1~ = C1 n ~ 1 - exp (-Reek n / An )] (2.19)
The constant C1 results from log law assumptions, n is the
normal distance to the wall, which give C1 = KC~,-3/4; the
OCR for page 300
OCR for page 301
OCR for page 302
OCR for page 303
OCR for page 304
OCR for page 305
OCR for page 306
OCR for page 307
OCR for page 308
OCR for page 309
OCR for page 310
OCR for page 311
OCR for page 312
Representative terms from entire chapter:
velocity components
asymptotic behavior of vT when n - O fixes An = 2C1; Ap
= 70 allows to recover the 5.45 constant of the flat plate
log law. The (straightforward) patching of the k inner
model with the outer k-e model is located at
Reek n _ 200-250.
3. The Numerics
3.1 The Transport Equations
The numerical method for the treatment of (2.14)
closely follows [123~14~. Coefficients of ~ derivatives are
evaluated at the center P of the control volume and
coordinates are normalized in the following way.
>* = >/4gll A* = 1l/
1 +CC(CU+CD)
{E CnbVnb+CC [CUvu +CDvd +—(Vn -Vn) Sv ] }
We- 1
1+CC(CU+CD)
TZ Cnbwnb+cc [CuWu +CDwd +—(We -We)-Sw ] }
CC is the finite analytic coefficient computed at the
colocated velocity mode (in the standard MAC method [49]
Cd, for Ud, Cn, for Vn, Ce, for We are computed at the
staggered velocity nodes). The matrix E is built with
coefficients like l+Cc (Cu+CD), etc...The matrix A
gathers terms like Cnb, CcOu, CCCD, ...Terms like R / ~
are in f. The G matrix for pressure contributions is inside
source terms (Su), (Sv)...
-Step 2- Equation (3.15b) can be written
U U J D(blP'`+ UPS+ by)
~ _
V=V-TlD(b2PA+ b2P~+ b <;)
V=W J-lD(biP + ~+ b3P )
w. 1ere
(3.18)
CcR
1+ CC(CU+ CD)
The underlined terms in (3.18) are generated by the
m~salignment between the coordinate lines and the
directions of the velocity components. With respect to the
staggered grid approach, it is seen that the D coefficient for
the pressure gradient does not depend now on the
correspond~ng velocity component.
-Step 3 - Substituting (3.18) into the continuity equation
(3.19):
[b,U + b2V + b3W] + [b,2U + b22V + b32Wl~+
[b,3U + b2V + b33W] = 0 (3.19)
gives the pressure f~eld (3.16) which can be viewed also as
the solution of the continuous equation (3.20)
—[aij_ ] = div V; aij = DJgij = DT~, bk bk (3.20)
ax ax k=1
Equation (3.20) is now discretized with a centered
scheme. Usually, this gives a 27 node pressure molecule
(9 in the upstream plane, 9 in the central plane, 9 in the
downstream plane). Seven points: U. D, N. S. E, W. P
come from second order derivatives 32 / axia~i, i=;; the
other twenty points come from the crossed second - order
derivatives 82 / 8kiD)J, i ~ j. Because of the retained
discretization, the pressure molecule involves only 19
points (five points in U and D planes). Among the whole
set of contributions to the pressure matrix, it is necessary
to keep all the fluxes in order to allow a complete
construction of the aij coeff~cients, i = j, because of the
misalignment problem. Extradiagonal coefficients i ~ j
arising from crossed second order derivatives are
connected to the level of non orthogonality of the grid.
They are not retained implicity but rather treated as source
terms. As a result diagonal dominance and symmetry of the
discretized pressure matrix is easily obtained. 18 velocity
components hav~e to be specified in (3.20) for the
evaluation of div V in each control volume
div V = (b~U + b2V + b3W),d - (b~U + b2V + b3W)U
+ (b2U + b2V + b32W) - (bi2U + b2V + b3W)S
+ (blU + b2V + b3W)e - (b~U + b2V + b3W)w (3.21)
Because of the definition of the control volume, the
whole set of velocity components has to be reconstructed
from the known centered velocity components. A simple
interpolation is used for that purpose. Because of the
structure of the retained discrete molecule, the discrete
pressure matrix is symmetric.
-Step 4- The computation of the non solenoidal V*
f~eld uses (3.17), the pressure field resulting from step 4
being present in source terms. The correction procedure
starts only when the field V* has been computed
th~ughout the complete domain.
3.4 Solution of the linear Systems
3.4.1 The Problem
Before discussing the technical solutions which have
allowed correct solutions to be obtained, let us point out
the origin of the difficulties. Pressure equations and
pressure correction equations are almost invariably solved
with SLOR type methods applied to the seven point
molecule induced by the treatment of under lined terms in
(3.18,19) as source terms. The main interest of this
approach is, as already seen, that the P matrices are now
symmetric and, usually, diagonally dominant. SLOR
methods can therefore be expected to converge. Moreover
the coding of the needed tridiagonal matrix algorithm is
straightforward. Unfortunately, due to the problem of
misalignment, dominant terms are omitted in the system
and treated as source terms. As a consequence, global
convergence of the metrod is made more diff~cult. Unless
strong departure from orthogonality is admitted for the
grid, this is not dramatic when wall functions are used
because grid clustering in the direction ~ remains rather
weak: the first points away from the wall being located so
that 50 < y+ < 200. When wall fuctions are avoided,
strong grid clustering implies a dramatic increase of the
pressure matrix stiffness so that, for threedimensional
problems, divergence of the iterative procedure occurs,
mainly through a too important pressure correction in the
wake. If the pressure matrix is enriched, [18], (using e. g.
nine points at the current station ~ = ki) it becomes more
difficult to invert. Thus, improving at least the pressure
solver is necessary.
3.4.2 The pressure solver
It has been found [50] that convergence was
systematically achieved in practice with an incomplete LU
preconditionned biconjugate gradient (PBCG method).
Conjugate gradient techniques are used as acceleration
methods in which the elements of the pressure matrix
(hereafter called A) are generated at each iteration, without
any need of storing the full matrix. The CG method
satisfies a minimisation property so that a monotone
decrease of the error is possible as well as a convergence
disruption to which is associated a monotonic decrease of
the error. This is not verified in practice because it is not
the full pressure matrix that is solved. Because the
spectrum of the A matrix is continuously distributed
preconditionning is necessary as a means of improving the
convergence of the BCG method: the system to solve Ax
= b is replaced by the equivalent system By = c where
K(B) << K(A), where K iS the condition number of the
301
matrix, so that eigenvalues of B are more clustered thaw
those of A. Wha^; is needed is an approximate inverse A
for A such that AAx = Ab or By = c where B = AA; c =
fib. B being "not to far" from the identity, the CG method
for B converges quickly.
The A matrix is selected so as to preserve the
sparsity of A. Here, an incomplete LU decomposition [51]
[52] -no fill in- is used - the so-called LU (1, 2)
decomposition performed worse than the retained classical
LU (1, 1) decomposition - . The solution of the LU
systems involves recursions which cannot be easily
vectorized so that the solutions of these triangular systems
is time consuming. Fortunately, a Neumann series method
applied to the L and U systems has been found to converge
nicely. The argument is presented for L. If L is written
under its Jacobi splitting L = DL + EL, L = DL [1 - 1 ]
where 1 = DL-1EL. It appears in practice that 1 is a
convergent matrix. A classical theorem of linear algebra
then indicates that 1-1 is regular and (1-1~-1 = 1+1+12+...
The Neuman serie is truncated to m terms and the resulting
approximate inverse of 1-1 can be computed with the
Homer scheme. The whole set of computations can be
vectorized with a speedup ratio higher than 30. Further
information on the PBCG algorithm are given in [511.
4. The Results. HSVA Tanker.
The flow domain covers 0 < X < 6 where X is
adimensionalized with L, the ship being from X = -1 to X
= 1.05; rs < r < L where r is the radial distance from the
axis of the ship y = z = 0 and rs is the radial location of the
hull. The flow domain is discretized with 80, 40, 31 points
in the axial, radial and girthwise directions. Starting from
an a-priori specified surface grid distribution, a volumic
grid is firstly generated using a transfinite interpolation
procedure [203~55] which is modified in order to allow
clustering of points in areas of strong concavity. The
resulting grid is irregular as the discontinuities on the body
surface are propagated throughout the domain. The elliptic
grid generation method is used to smooth out the grid
irregularities. Fig.3 presents a perspective view of the grid
which appears to behave correctly.
Inlet conditions, written at X = 0, are estimated in
such a way that results at X = .291 roughly match with
experiments [11. The main Indetermination lies in the fact
that specified data [13~44] involve 81 1, 61, Cf and H12,
where:
tan
o
^00
r--
Ue-U U d U-U
[ Qe ~ Qe 1 | [ Q
J
o
Hl2=8llI6l ;Cf=tW/ 2pU2
from which needed values for 6, Up and Qe are estimated
(,Bw - O). These estimations allow the generation of inlet
velocity profiles by the method of Coles & Thompson
t531. It is therefore necessary to briefly mention the results
of the sensitivity analysis of 6, Us and Qe data. It is found
that, if the influence Up is very weak, the result of the
computations depend on the level of the overshoot Of Qe
with respect to 1, and of the location ~ of this overshoot.
We consider that this is the most serious weakness of this
approach and this justifies recents attempts [6] [29] t30]
[39] to compute the flow past a complete full shape,
starting from far upstream. Apart from difficulties
connected to the transition specification -a problem on ship
models, not on full scale ships-, the flow becomes then
somewhat sensitive to the level of the turbulent viscosity
outside the viscous zone so that consideration of these
aspects is left for a future study. Other boundary
conditions are either wall-type conditions or symmetry
conditions. Neumann boundary conditions are used for the
pressure, except for r = L where P is taken equal to zero.
120 global iteration are usedwith a local time
stepping procedure. Local time step is divised to ensure a
fixed amount of diagonal dominance will respect to
discretized momentum equations.
Fig.4 presents longitudinal pressure distributions. A
sensible grid effect in the longitudinal direction is present
in this case. Girthwise distributions agree reasonably well
with experiments although some more resolution should be
necessary (fight. Comparison of secondary velocities with
experiments [44] is given in fig.6. The convergence of
streamlines, which starts to develop at midship, increases
more downstream and gives rise to a counterclockwise
primary vortex in the hollow of the stern. The primary
vortex is close to the keel plane both in the computation
and in the experiments which indicate a "blind" zone for
AL > .95, where the longitudinal velocity component is
too small to allow significant measurements. The
longitudinal vortex where the hull terminates is also found
to agree with experiments. If isowakes are now considered
(fig. 7) -they are defined as isovalues of the component U
in crossplanes-, it is found that the experimental trends are
correctly captured (except, of course, for string effects
which are not present in the computation, see e.g. x = 38
mm, O mm, 100 mm, in [44] data-~. It must be mentionned
that a bilinear interpolation is used to specify isowakes so
that it is misleading to compare isowakes .95, 1, or 1.05
as corresponding velocity profiles are usually very "flat"
close to U = 1. The apparent overestimation of the viscous
thickness close to the waterline is probably due to the fact
that the inlet viscous thickness does not depend on girth,
while it should thicken from the keelplane to the
waterplane.
Velocity components U. V, W as well as pressure
data can be compared more extensively (fig. 8a to c). U and
V are presented on the upper left and right sides,
respectively, W and P are presented on the lower left and
right sides, respectively. For each series of plots, the
evolution is considered with respect to y at a given station
X = const., for several depths z = const.. In general, it is
seen that the calculations are in correct agreement with the
data, especially on V and W. except, for U. in small
regions close to the hull or in the neighborhood of the
wake centerplane. This discrepancy is believed to be due to
the hull shape modifications. Due to the large amount of
grid points which would have been necessary to account
for the hub support boss located for -19 mm. < x < 0 mm.
(in [401 data), this detail of the geometry has been
eliminated (fig. 9~. Also, coordinate lines (L): ~ = 1, (=
const. evolve continuously from the ship hull into the
vertical plane V of symmetry of the wake. They define a
family of curves that describe the hull surface for ~ <
kw (X) where kw (X) defines parametrically the
intersection of the hull-surface DB with V; this intersection
is located for .95 < X < 1.05. For ~ < kw(X), lines (L)
are almost rectilinear in V, ie. zeal = 1, ~ given) depends
only weakly on \. In V, the line ~ < kw(X) usually
302
consists of a vertical line located at X = .95 for Zb < Z <
Zr. For Zr < z < 0, the line ~ ~ kw(X) (which is a
piecewise linear function of X) is rather close to a straight
line, the equation of which can be written:
Z/Zr = 1 +10 [.95 - X ~
so that the hull has disappeared in the cross plane X
= 1.05. Now, the approximate hull profile in V is not
described exactly by ~ < kw(X) but rather by a "stair"
function (fig.9) consisting of a succession of lines
X=const., each of them separated by aline close to
z = const.. The no-slip condition is evidently written on the
resulting "stair surface" DB s which is a crude
approximation of the "longitudinal geometry" while
crossplanes are correctly body-fitted. As a result, the
description of the viscous layer is rather poor.
Nevertheless, a comparison with [14] -where a similar
treatment is present- indicates that the slight improvement
obtained on the secondary flow can be attributed to the fact
that wall factions are not used here. The discrepancy
between the computed and measured U components is
important, especially in the wake. Apart from the fact that
this discrepancy is already present upstream, another
possible explanation could be the incorrect account of the
loss of no-slip which results from the turbulent model.
Fig. 10 presents the comparisons between the
computed and measured longitudinal vorticity component
isolevels . In spite of the uncertainties involved in
processing both the measured and computed data, such a
comparison allows examination of differences of velocity
gradients rather than velocity themselves. This test is
therefore more stringent. As a result, it is found that the
computations reproduce both the general features of the
contours and the magnitudes. It is also evident that the
diffusivity which is present in the viscous region is more
correct, with respect to experiments, than that which would
result from the use of wall functions [141. The differences
in the shape of the isovalues, in the wake, could be
foreseen from the differences in the velocity distributions.
Fig. 1 1 presents a color picture of lines of wall friction, the
hull is coloured with the pressure field, high Cp and low
Cp being associated to red and blue colours respectively.
The flow which looks clearly twodimensional at midship
enters an oblique low pressure region which deviates the
wall streamlines. The resulting convergence of these
friction lines implies a thickening of the viscous layer,
close to the aft, that makes thin boundary layer methods
breakdown [11. A vertical convergence separation line
seems also present, close to ~ = kw (=.95~. Due to the
approximate characteristics of the hull geometry and to the
fact that the U velocity component is very low here, the
materiality of this result remains to be confirmed from
experiments. The maximum pressure levels are found on
the hull X ~ 1-1.025. More downstream, the vertical
structure of the flow is confirmed: streamlines associated
to the secondary velocity components V and W are drawn
for two crossections X = const.
S. Con cl usi on
A fully elliptic numerical method for the solution of
the full Reynolds-averaged Navier-Stokes equations has
been applied to the computation of the flow past a shiplike
hull. The method uses a system of numerically generated
curvilinear coordinates and retains the Cartesian velocity
components as independent variables. The turbulent
closure of the equations is handled with the well-known
303
k-e model. The use of wall functions is avoided; instead,
a one equation model is used close to the wall. Bypassing
the wall functions implies a very high grid resolution in the
direction normal to the wall as well as high aspect ratios for
the discrete molecules. Due to the resulting increased
stiffness of the pressure matrix, it has been found
necessary to improve the pressure solver in order to get rid
of convergence problems. The numerical method has been
described and evaluated on one ship hull for which
experimental data are available. From the presented results
the following conclusions can be drawn.
With respect to the modelling of physical flow
characteristics, the main features of ship stern flows appear
correctly captured. In particular, the viscous inviscid
interaction gives rise to a pressure field in good agreement
with experiments. The method allows also a correct
description of the mean velocity flowfield in the thin
boundary layer, provided suitable initial conditions are
available [51~. In the thick region, close to the aft, where
thin boundary layer methods systematically breakdown,
the overall features of the flow are also captured:
secondary velocities and even longitudinal vorticity
components are correctly predicted. Where discrepancies
are found, they may be attributed to the technique which
has been used to approximate geometrical details in the
framework of a monoblock structured grid. Improvements
in the description of the flow close to the propeller disk
region therefore require a more detailed resolution of the
geometry.
With respect to the methodology, ship stern flows do
not appear to involve too strong a viscous-pressure
interaction. Therefore, the segregated approach seems well
adapted. Apart from geometrical grid aspects
improvements could result from one of the following three
points.
(i) The convection diffusion scheme. The finite analytic
scheme only roughly accounts for the correct physical
imbalances: pressure gradients vs. diffusion in the near
field, pressure gradient vs. convection in the far field, for
all the possible convection velocity directions. More
particularly, the hypothesis of local uniformity of the
influence coefficients and of the source terms, which
allows locally (finite) analytic solutions to be computed, is
probably too restrictive, especially in the direction away
from the wall, where the convection velocity and the
source terms vary rapidly, owing to strong gradients of
turbulent quantities.
(ii) The Poisson solver. This is the most time-consuming
part of the method, especially on threedimensional grids
where the aspect ratios can be very high. The conjugate
gradient technique appears very well adapted to the
problem in spite of its induced storage increase. It can be
noticed that it may be somewhat unreasonable to compute
the pressure -which is regular and does not vary
significantly close to the walls- on a grid which appears
rather well suited for the velocity distributions and for the
turbulent quantities in that they especially vary close to the
walls. Because less points are needed for the pressure, the
use of a multigrid procedure should efficiently cut down
cpu times. Because relaxation is a local procedure, a
multigrid technique might be also a way of improving the
solution in the boundary-unfitted region close to the aft.
(iiiJ The Turbulent Model is also a potential source of
uncertainty for the description of near-wall flows where
strong curvature effects are present. Although the flows
considered are mainly pressure-controlled, improvements
in the turbulence models might lead at least to improved
length scale predictions (the turbulent kinetic energy is
known to be overpredicted in thick boundary layers [40]).
Nevertheless, it is considered that improved results will be
rather an outcome of a better treatment of the geometry and
of a better convection-diffusion algorithm.
Acla~owledgrnents: Authors gratefully acknowledge
financial support of DRET through contract 86 / 104. Cpu
on Cray2 was attributed by the Scientific Committee of
CCVR. Cpu on VP200 was attributed by DS / SPI. The
presented method has been developped from an earlier
code version of the method [121.
preferences
[1] L.Larsson. (Ed.) Proceedings of SSPA-IITC on
Ship Boundary Layers. SSPA Report N°90 (1980)
[2] L. Larsson and M.S. Chang, Numerical Viscous
and Wave Resistance Calculations Including Interaction.
Proc. 13th ONR Symp. Naval Hydrodynamics, Tokyo,
(1980)
t3] T. Nagamatsu, Calculation of Ship Viscous
Resistance and its Application. Journ. Soc. Naval Arch.
Japan 157, 47 (1985) - See also Proc. 2nd lnt Symp.
Viscous Flow. SSPA-Goteborg, (1985)
t4] S. Soejima, Calculation of Three-dimensional
Boundary Layers Around Ship Hull Forms. Proc. 2nd.
Symp. Num. Phys. Aspects of Aerod. Flows, Long
Beach, CA. (1985)
[5] H. Miyata and S. Nishimura, Finite Difference
Simulation of Non linear Ship Waves. J. Fluid Mech.
157, 327-357 (1985)
[6] H. Miyata, S. Nishimura and A. Masuko, Finite
Difference Simulation for Non linear Waves Generated by
ships of Arbitrary Threedimensional Configuration.
J. Comps Phys. 60, N°3, 391-396 (1985)
[7] H. Miyata, H. Kajitani, S. Akifuji, M. Baba,
M.Kanal, Wave Formation about a Ship Bow advancing in
Head Seas. Spring Meeting Soc. Nav. Arch. Japan (1987)
[8] M. Hoekstra, H.C. Raven, Application of a
Parabolized Navier-Stokes Solution System to Ship Stern
Flow Computation. Proc. Osaka Int. Colloq. on Ships
Viscous Flow, 125-142 (1985)
[9] M. Hoekstra and H.C. Raven, Ship Boundary
Layer and Wake Calculation with a Parabolized Navier-
Stokes Solution System. Proc. 4th Int. Conf. Num. Ship
Hydrodynamics, Washington D.C. 470-491 (1985)
[10] H.C. Raven and M. Hoekstra, A Parabolized
Navier-Stokes Solution Method for Ship Stern FLow
Calculations. Proc. 2nd Int. Symp. Ship Viscous
Resistance, SSPA Goteborg (1985)
t11] L. Broberg, Numencal Calculation of Ship
Stern Flow. PhD Thesis, Univ. Of Chalmers (1988~; also
SSPA Rept. 2801-2 & 2803-1
[12] H.C. Chen and V.C. Patel, Calculation of
Trailing-Edge, Stern and Wake Flows by a time-Marching
Solution of the Partially Parabolic Equations. IIHR Rept
N°285 Iowa City, la (1985)
[ 13] H.C. Chen and V.C. Patel, Numerical
Solutions of the Flow over the Stern and in the Wake of
Ship Hulls. Proc. 4th. Num. Conf. on Numerical Ship
Hydrodynamics, Washington D.C. 492-511 (1985)
[14] V.C. Patel, H.C. Chen and S. Ju, Ship Stern
and Wake Flows: Solutions of the Fully-Elliptic
Reynolds-Averaged Navier-Stokes Equations and
Compansons with Expenments. IIHR Rept. 323 Iowa
City,la(1988)
[15] H. Miyata, T. Sato and N. Baba, Finite
Difference Solution of a viscous Flow `, i.h Free-surface
Wave about an Advancing Ship. J. Comp. Phys. 72,
393-421 (1987)
[16] G. Tzabiras, On the Calculation of the 3D
Reynods Stress Tensor by Two Algorithms. Proc. 2nd.
Int. Symp. Ship Viscous Flows, Goteborg, Pap. 15
(1985)
[17] J. Piquet, P. Queutey and M. Visonneau,
Computation of Viscous Flows past Axisymmetric Bodies
with and without Propeller in Operation. Numerical
Methods in Laminar and Turbulent Flows (C. Taylor,
W.G. Habashi, M.M. Hafez, Eds.) 5, Part 1, Pineridge
Press, 644-655 (1987)
[ 18] J. Piquet and M. Visonneau, Steady
Threedimensional Viscous Flow Past a Shiplike Hull.
Proc. 7th. GAMM Conf. Num. Methods in Fluid
Dynamics. Notes on Numerical Fluid Dynamics,
(M. Deville, Ed.) 20, Vieweg, 294-301 (1988)
[19] G.D. Tzabiras, Numerical and Experimental
investigation of the Flow Field at the stern of Double Ship
Hulls. PhD Thesis, Nat. Tech. Univ. Athens (1984)
[20] J.F. Thompson, Z.U.A. Warsl and C.W.
Mastin, Numerical Grid Generation, Foundations and
Applications. North Holland Publ. New York (1985)
[21] G. Tzabiras and T.A. Loukakis, A Method for
Predicting the Flow Around the Stern of Double Ship
Hulls. Int. Ship. Progress N°345 (1983)
[22 ~ A.M. Abd El Meguid, N.C. Markatos,
D.B. Spalding, K. Muraoka, A Method of Predicting
three-dimensional Turbulent Flows around Ships'Hulls.
Proc. Ist. Int. Symp. Ship Viscous Resistance, SSPA
Goteborg (1978)
[23] A.M. Abd El Meguid, N.C. Markatos,
K. Muraoka, D.B. Spalding, A comparison between the
Parabolic and Partially-parabolic Solution Procedures for
Three-dimensional Turbulent Flows around Ship's Hulls.
App. Math. Modelling 3,249 (1980)
[24] C.E. Jansson and L. Larsson, Ship Flow
Calculations using the Phoenics Computer Code. Proc.
2nd. Int. Symp. on Ship Viscous Resistance SSPA,
Goteborg (1985)
[25] K. Muraoka, Examination of a two-Equation
Model of Turbulence for Calculating the Viscous Flow
around Ships. J. Society Naval Arch. Japan, 147, 35
(1980)
[26] K. Muraoka, Calculation of thick Boundary
Layer and Wake of Ships by a Partially Parabolic Method.
Proc. 13th. ONR. Symp. Naval Hydrodynamics Tokyo
601-616 (1982)
[27] R. Peyret and T.D. Taylor, Computational
Methods for Fluid FLow. Springer-Veriag Ed., New-
York. (Senes in Comp. Phys.)
[28] D.C. Kwak, J.L. Chang, S.P. Shanks,
S. Chakravarthy. An incompressible Navier-Stokes Flow
Solver in Threedimensional Curvilinear Coordiante System
Using Primitive Variables. AIA,4 J. 24, 390-396 (1986)
[29 1 Y. Kodama, Computation of 3-D
Incompressible Navier-Stokes Equations for Flow Around
a Ship Hull Using an Implicit Factored Method. Proc.
Osaka Int. Symp. on Ship Viscous Flow, 109-124
[30] Y. Kodama, Computation of High Reynolds
Number Flow past a Ship Hull using the IAF Scheme.
Proc. Soc. Naval Arch. Japan, Spring Meeting (1987)
[31 ] S.P.Vanka, Computations of Turbulent
Recirculation Flows with Fully Coupled Solution of
Momentum and Continuity Equations. Argonne National
Lab. Rept ANL-83-74 (1983)
[32] S.G. Rubin, Incompressible Navier-Stokes and
Parabolized Navier-Stokes Formulations and
Computational Techniques. Computational Methods in
Viscous Flows, 3; Series: Recent Adv. in Num. Meth.
Fluids, Habashi, Ed., Pinendge Press (1985)
304
[33] S.V. Patankar and D.B. Spalding, A
Calculation Procedure for heat, Mass and Momentum
Transfer in Three-Dimensional Parabolic Flows. Int. J.
Heat Mass Transfert 15, 1787 (1972)
t34] S.V. Patankar, Numerical Heat Transfer and
Fluid Flow, Hemisphere Publishing Corp. New-York.
[35] R.I.Issa, Solution of the Implicity Discretized
Fluid Flows Equations by Operator-Splitting, J. Comp.
Phys., 62, N°1, 40-65 (1986)
[36] M. Israeli and A. Lin, Iterative numerical
solutions and Boundary Conditions for the Parabolized
Navier-Stokes. Comp. & Fluids 13, N°4, 397-410 (1985)
[37] C.J. Chen and H.C.Chen, Finite Analytic
Numerical Method for Unsteady Two-Dimensional Navier-
Stokes Equations. J. Comp. Phys. 53, N°2, 209-226
(1984)
t383 T. Hino, Numerical Simulation of a Viscous
Flow around a Ship Model, Proc. Soc. Naval Arch. Japan,
Spring Meeting (1987)
[39] A. Masuko, Y. Shirose and S. Ishida
Numerical Simulations of the Viscous Flow around Ships
including Bilge Vortices, Proc. 1 7th. Symp. Naval
Hydrodynamics, Den Haag (1988)
[40] V.C. Patel, Some Aspect of Thick Three-
Dimensional Boundary Layers, Proc. 14th. ONR Symp.
Naval Hydro. Ann Arbor, 999-1040 (1982)
[41] V.C. Patel, Ship Stern and Wake Flows:
Status of Experiment and Theory, Proc. 17th. Symp.
Naval Hydrodynamics, Denn Haag (1988)
[42] L. Larsson, Boundary Layers on Ships, PhD
Thesis, Univ. Chalmers (1975)
[43] L. Lofdahl and L. Larsson, Turbulence
Measurements Near the Stern of a Ship Model, J. Ship
Research 28, N°3 (1984)
f44] K. Wieghardt and J. Kux, Nomineller
Nachstrom auf Grund von Windkanalversuchen, Jahrb.
Der Schiffbautechnischen Gesellschaft (STGJ, 303-318,
(1980)
[45] K. Wieghardt, Zur Kinematik einer
Nachlaufstromung, Z. Flugwiss. Weltraumforsch. 7,
Heft 3 (1983)
/ -
Fig.lc
Fig. la - Physical Domain
Fig.lb - Computational Box
Fig.lc - Transverse plane in the Computational Domain.
[46] H.C. Chen and V.C Patel, Calculation of Stern
Flows by a Time Marching Solution of the Partially
Patabolic Equations, Proc. 15 th. ONR Symp. Naval
Hydrodynamics, Hamburg (1984)
[47] H.C. Chen and V.C. Patel, Near-Wall
Turbulence Models for Complex Flows including
Separation, AIAA Paper 87-1300
[48] M. Wolfshtein, The Velocity and Temperature
Distribution in One-Dimensional Flow with Turbulent
Augmentation and Pressure Gradient, Int. J. Heat Mass
Transfert 12, 301-318 (1969)
[49] F.H. Harlow and J.R. Welsh, Numer-ical
Calculation of Time Dependent Viscous Incompresible
FLows of Fluid with a Free Surface, Phys. Fluids 8
2182-2189 (1965~. See also Los Alamos Rept. LA-3425
(1966)
[50] J. Piquet, P. Queutey and M. Visonneau,
COmputation of the Threedimensional Wake of a Shiplike
Body, Proc. 11th. Int. Conf. Num. Meth. Fl. Dyn., Lect.
Notes Phys. Springer-Verlag (1989)
[51] J. A Meijerink and H.A. Van Der Vorst. An
Iterative Solution Method for Linear Systems of which
the coefficient Matrix is a Symmetric M-Matrix, Maths.
Comp. 31, N°137, 148-162 (1977)
[52] J.A. Meijerink and H.A. Van Der Vorst.,
Guidelines for the usage of Incomplete Decompositions in
Solving Sets of Linear Equations as they occur in Practical
Problems, J. Comp. Phys. 44 (1981)
[53] T. Cebeci, K.C. Chang and K. Kaups, A Three
Dimensional General Method for Threedimensional
Laminar and Turbulent Boundary Layers on Ship Hulls,
Proc. 12 th. ONR Symp. Naval Hydrodynamics,
Washington (1978)
t54] K.H. Mori & N. Ito, Near-Wake Computations
by Solving the Vorticity Transport Equation in a Body-
fitted Coordinate System, Proc 4th. Int. Conf. Num. Ship
Hydro., Washington D.C., 512-528 (1985,)
[55] L.E. Erikkson, Generation of Boundary
Conforming grids around Wing-Body Configurations
using Transfinite Interpolation, AIAA Journ., 20, N°10
1982.
305
NW tNC NE
_ WC = r E C 5
_ _ , ~ . ,
SWI ISC h 15£
1-' -1' - - -1
N t
n
., (i)
he'd 1._ __~
~ ' ' : 1 _/ '
1 1 1 ~ ,l,~
_~. _: ~ D
In it T .
~ ~ _ ~ __l__,L__ _ it
1 {/ I.~/1 ~
Fig.2 - Colocated grid. a, location of velocity components
and of pressure.
cat
\
~1
n
.5 .S .7 .3 .9 1.~ I.t 1.: 1.3 1.' 1.; 1.i 1.7 1.] t.9 2.0 2.1 2.2 2.3 2.' 2.5
!~' Fig.3 - Perspective view of the grid surface j = 1; 11 = 0.
~ \\
p
F
:~
Fig.4 - Longitudinal Pressure Distributions.c,, experiments
—, computations.
,,.
.
1
I,.
.
111
,,.
Fig.5 - Girthwise Pressure Distributions. , experiments;
, computations.
306
Al
. -
lI,
·:
i:
- ~
al~
~J
11 . 1
l l ~ l ~ l
0 .04 .o~ .12 .16 .20 .20
x
or
(o
no.
A.
co
N
(O
LD
CO
tO
O o
C'
.n
a;
_
a....
(
t
!~
~i
O . I
Fig.S - Girthwise Pressure Distributions. Q. experiments;
, computations.
U
W
=2
o
N
~0
~1
O-
o
~1
~ :
_ _
_
·
0 .04 .08 .12 .16 .20 .24
x
tO I
(~ 1
~r
p
O- . . , · .
0 .04 .08 .12 .16 .20 .24
Fig.8a - Velocity components and Pressure at a given
crosssection X/L = .944, for several depths Z = const., as
a function of Y.a, experiments; , computations.
307
Fig.8b - W
1 1
. 1 1 1 1 1
0 .04 .08 .12 .16 .20 .24
0~
_-
~ .
._
~o
a, -
Fig.8b - U
~ ~ I ~ ~ ~ ~ ~ ~ ~ ~
0 . Oq . 08 . 1 2 . 1 6 . 20 . 24 0 . 04 . 08 . 1 2 . 1 6 20 . 24
Fig.8b - P
~ , , , , , , 1
0 .04 .08 .12 .16 .20 .24
0 .04 . 08 . 1 2 . 1 ~ .20 . 24
o '
a)
o
~D
(o
rn
r~
_~.
Fig.8c - P
.
o- ~
o
,
o
-
~:
(.D- ~
~.
O
o.
N
~,
(D '
_
-
(~1:
_ :: ~
t0 _ ~ O _ A A, e
, mA
(~0 — ~ ~ A A A A o
O
-
~ _ _
.04 .08 .12 .16 .20 .24
Fig.8c - W
w~ ..
ll ~ ~ ~ ~ ~ ~
O ~0~ A~
1 , ~ I · I
0 .04 .00 .12 .16 .20 .24 0 .04 .08 . i2 . i6 .20 .24
Fig.8b, c - Velocity components and Pressure as a function of Y. at given crosssections: XIL =
.962 and XJL = 1.007, for several depths Z = const. , experiments; , computations.
308
ins
-
i~''2_22,'2
Fig.6 - Secondary velocities at several crosssections X =
const.; left, experiments; right, computations.
wit by G
~~ ~~ O
X/L = 0.942
o,/~7
X/L = 0.955
~n,.l .~y
o-
lo
(D
lo
.
At_
tool
.~
Fig.7 - Isowake lines; left, experiments; right,
computations.
309
~ x = -157mm
-
i /, /,,
~ /, /, ~
/ 32^~"
_
~ ,
~ ,
1 _', '
~ ,
i
;
x = -53mm
/'
i.,
it,,
/'',''
/, A,
/ , - .,
/,
~ ,
,,
, ~ ;. s .
/,` ~ ,'
/' ',2 ~
/';
,, i,
~4 ,,"-
~; =_` ,,
~ _ ~,~ _ _ _ _ _
/~/ "
/,' ~
/'
1
!
oc .
I C :
r .
0c ~
_
1
G~
. ~
~Jt
Ot
1 ~
_~_
ot
.(
CO
C~`
r~ ~
.~ CO~
'Dt
I L
t
O . Q2 . 0~ 06 . 08 . 1 0 . 1 2 . 1 4
X
X1L=.949 ,~~ - ~-
= , , ~ -r—- r—
O .02 .O] .06 .08 .10 . 1;? . 11
1
',. X/L = .962
~ ''— , I r I I I
G . 02 .01 . 06 . 013 . 1 Ct . 12 . 1'1
Fig. 10 - Contours of isolevels of the longitudinal component of the vorticitv fi~la in tr~ncv`~r~-
crosssections. Left, experiments [443; right, computations.
3~0
, _~_ .. ~ ~—v · _^ v~
o-
x=l~ /
.
_~~ ~ ~ e
~5
=
1
1
1
. ... _ _
:'°~ /
~ '0 ~
~ I
WS
~-
_~ \\
/
/
0 .02 .~] .06 .~S .10 .12 .11
Fig. 10 - Contours of isolevels of the longitudinal component of the vorticity field in transverse
crosssections. Left, experiments [441; right, computations.
Fig.9 - Detailed view of the aft part of the HSVA tanker. ---, \(x) exact intersection between V
and ~B; , present approximate "stair" representation of this line (note: the sting is omitted
in the computations).
~Groupe de 1~1ode1 isat ion Nu~erique RSCETE - . ~ ~ . -- ~
., , =, ; ' ', , -, ,_ ' - ~ .+ - - ^ - '; . _i __ ' - ' ' ' ' ·
~ A;~ .~— ~ ' j ~ ~ , ; , A ~ y* ii ~ 't;; ' ~ ~ . ! '7n ~ - i i T. i
· r 7 ~ ~ ~ ;; ~ , ., 7 7 ; ~ ~ ; ~ . ~
.._. .- ~ ' _ '~.; 7.~. 2' ~ i; ~ 7 ~ ~ ; . ;~i _ - ~ ~ —+~ ~ ~ ~ ! j7~ 7~A j1 ~ _~_;7 ~ t-; - ~ ~ ~ ·— ~+ j ~ !~ .~ ~ 7
_ ^;, ~ .: 2 ~ , ,. - h,, _,, 7 ~, ~ ~ ~ ,~ ~ _ ! ~ ~ ~ ~ ~ 7 ; ~ T i 7 i i . _
;~ - ~~ ~ i '' ~ ~; -t ~ ; , ~~~ , ~, ~ i - ~ -' ' ~ , ~ , $~7 ; j r
~' .;~J ~' 2 ~ 't ; ' ; 1 ~i !
~'.'.~;'~~~ 72:.,,~,, .~' -=r-'_,
j ,_ + ~ —'';,, ~_
~ f. + -
,: _. ~. " '~, ~ , +, . :', ., T.': -
~ ;; ~+~ _3,~
+
- - - . - . . . - i i . . . : . _ . ? .: _ _ 7: :: ~ . L _ . 7 7 - ~ 7 : A
i~ _ =- -~ +; ' -i '~ - - ; ~ ;'i ;- - '- - ' . ' i- ii7 .; - ~ - ,.
Fig. 1 1 - Perspective view of the aft part of the HSVA tanker; skin-friction lines.
311
DISCUSSION
by Y. Kodama
You used a non-aligned grid system near
the stern, which greatly reduces load for grid
generation and improves grid orthogonality.
However, I think a special treatment should be
needed in order to satisfy conservation there.
Did you do those treatments there? If not,
what is your opinion about the effect of
possible violation of conservation there on
the computed results?
Author's Reply
Because of the use of an H-grid associate
to the constraint x=x(~),the grid which has
been built here is no more body-fitted in the
intersection between the overhang and the
plane of symmetry(y=0). This implies two
consequences;
(i) the exact geometry is not well
represented, so that the no slip
condition is not written at the right
place in the plane of symmetry.
(ii) the development of the boundary layer
below the overhang is not correctly
captured because of the lack of
discretization points and a distribution
law which is not in agreement with the
strong variation of the velocity profiles in
the plane y=0. Due to the used topology,
this weakness is not easily surmountable
because the grid lines which are supposed to
describe the overhang symmetry line boundary
layer come from the surface of the hull.
It is believed that the correct
representation of the body surface is more
important than conserving momentum on a
geometry which is incorrectly reproduced in
the symmetry plane. The best remedy to this
fundamental weakness is to relax the
constraint x=x(~), working on more general
(and less orthogonal) grids.
312
