Below are the first 10 and last 10 pages of uncorrected machineread 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, chapterrepresentative 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 87
Recent Developments in a Ship Stern Flow Prediction Code
M. Hoekstra
Maritime Research Institute Netherlands
Wageningen, The Netherlands
Abstract
This paper is concerned with developments in a
numerical method for the prediction of the steady
incompressible flow around the stern of a ship. The
method is based on a finitedifference discretisation of
the Reynoldsaveraged NavierStokes equations on a
boundaryfitted curvilinear grid. The attention will
be focused on three aspects, viz. the choice of the
velocity variables, the boundary conditions and the
global (outer) iteration process for the pressure field.
Various illustrations will be given of the improved
performance of the method.
1. Introduction
It is in the first place owing to computer hard
ware developments that a computational simulation
of the ship stern flow field has become feasible. In
deed, the availability of adequate computer facilities
may be considered as a prerequisite for an accurate
calculation to be carried out. It is not a sufficient
condition, however. For some years already, the nu
merical simulation of the steady viscous flow around
a ship's stern  usually under neglection of free sur
face effects  has been a prominent research topic.
Yet, something like a commands opinio on how a re
lial~le, efficient and robust code is to be constructed
has not been attained. It indicates that developing a
suitable algorithm for the solution of the Reynolds
averaged NavierStokes equations for incompressible
flows around complex body shapes is not plain sail
ing.
At the author's institute, work in the field of stern
flow computation was initiated around 1980 and has
continued steadily since. It has resulted in a compu
tational procedure that has been applied with some
success to a number of ship forms and axisymmetric
bodies. The main features of the method have been
outlined previously ~1,2,3] and a brief recollection of
these features (Section 2) will suffice here as a frame
work for the subsequent discussion of some aspects of
the method.
Where the title of this paper promises the presen
tation of recent developments, 'recent' is to be inter
preted in relation to the time of appearance of the
above references. The length limitations to which
this paper is subjected do not allow a presentation of
all changes that have taken place in the code since
1985. Therefore three major items have been se
lected for discussion here. They are: the choice of
the velocity variables in relation to regular solution
behaviour near grid singularities in the wake (Sec
tion 3~; boundary conditions and their numerical im
plementation (Section 4) and a new global pressure
relaxation scheme (Section 5~. Moreover, results of
a representative application will be shown as an il
lustration of the progress that has been made (Sec
tion 6).
2. Outline of computational procedure
We are interested in a numerical simulation of the
behaviour of the flow around the stern of a ship, be
ing towed steadily through still water. The free sur
face disturbances created by the ship are assumed
to be of little significance so that it is allowed to
replace the undisturbed free surface by a symmetry
plane. Thus, taking a reference frame moving with
the ship, we consider actually the double model of
the submerged part of a ship, held fixed in a uniform
flow directed from bow to stern. The flow domain is
now assumed to be unbounded and the fluid to be
incompressible. Three subdomains are distinguished
as shown in Fig. 1, viz. the external flow zone where
the flow behaves as being inviscid and irrotational;
the bc)undarylayer zone, covering the forward part
of the hull surface, where firstorder boundarylayer
theory is supposed to provide an adequate description
of the flow; the sternflowandwake zone which is of
87
OCR for page 87
primary concern in this paper. In the latter zone we
solve eRectively the Reynoldsaveraged NavierStokes
equations but we shall assume that a predominant
flow direction exists. Among other things, it allows
us to drop all diffusion terms associated with that pri
mary flow direction. The equations are supplemented
by an algebraic turbulence model. The boundary con
ditions will be described in Section 4.
III /////HULL /~
Fig. 1. Division of flow field into three zones
The flowgoverning equations are written in terms
of general curvilinear coordinates. The contravariant
formulation of the momentum equations is adopted,
which gives three relations expressing momentum con
servation along each of the three coordinate lines.
The velocity components and the pressure are the
dependent variables but a particular choice for the
velocity variables is made as described in Section 3.
The continuity equation is solved in its original form
and is not replaced by a Poisson equation for the
pressure as in MarkerandCell type methods.
The equations are discretised on a boundaryfitted
curvilinear grid. It is a singleblock regular grid of
N V * NY * NZ nodes. The grid is partly nonortho
gonal but conceptually simple. All variables are de
fined on the grid nodes, in other words we do not
apply grid staggering of one kind or another. That is
convenient in many respects but it poses certain de
mands on the discretisation of the pressure gradients
and the continuity equation which we are careful to
fulfill to avoid 'checkerboarding' A. The discretisa
tion is second or third order accurate with the excep
tion of the main stream pressure gradient which is of
first order accuracy. All derivatives with respect to
the coordinate associated with the main stream di
rection are modelled by upstream difference formulae
except the pressure gradient. Newton linearization is
applied to the convective terms.
The solution is obtained iteratively where the it
eration sequence may be characterized as a multiple
sweep spacemarching process which takes maximum
benefit of the existence of a predominant flow direc
tion. In fact, two iteration cycli can be distinguished:
the local (or inner) and the global (or ousters iteration
prc)cess.
The local process refers to the solution of a sub
set of difference equations applying to all grid nodes
having the same value of the mainstream coordinate.
It is based on the Coupled Strongly Implicit Proce
dure (CSIP) [5], an incomplete factorization scheme,
and yields a simultaneous solution of all variables in a
crosssection of the computation domain. Iteration is
required both by the nonlinearity of the differential
equations and by the incompleteness of the factoriza
tion. The CSIP has been selected because we think
that a high degree of implicitness is desirable if not
necessary for an efficient solution algorithm and be
cause it retains the coupling between the momentum
and the continuity equations in the numerical solu
tion.
The global iteration process involves the evalua
tion of the solution in repeated sweeps from the up
stream to the downstream boundary of the compu
tation domain. This process allows the influence of
downstream occurrences to be felt by the upstream
flow, an inherent property of the mathematical model
used. To improve its convergence rate, each down
stream sweep is followed by an update of the solution
of the pressure field in a reversed sweep as will be
explained in Section 5. A further enhancement of
the computational efficiency is obtained by grid se
quencing, with which we mean that the grid used is
initially very coarse in the mainstream direction and
is subsequently refined in two or three stages.
3. Grid singularities
Fig. 2 gives a sketch of the computation domain as
it looks after the symmetry properties of the flow have
been taken into account. In addition to the Carte
sian reference frame x,y,z a curvilinear boundary
fitted coordinate system 6, a, ~ is introduced. A trans
formed space can be imagined in which (, 7~, ~ forms a
rectangular coordinate system and in which the hull
surface is plane. The appearance of the (,7~,( sys
tem in physical space depends on the transformation
relations between (,71,( and ~,y,z.
PHYSICAL DOMAIN
z
,/
REV
(a)
TRANSFORMED DOMA I N
TV
~577
REV
Fig. 2. (computation domain
88
(b)
OCR for page 87
To simplify the grid generation procedure we have
chosen ~ = x so that transverse sections in the phys
ica1 space remain so in the transformed space. An
orthogonal mesh is then constructed in selected trans
verse sections by a procedure which relies on confor
mal mapping although the final mesh is not conformal
due to subsequent stretching t64. By connecting these
grids, the 3D mesh is completed. It leads inevitably
to nonorthogonality of the if, 71 and if, ~ lines in a part
of the computation domain.
However, before this task can be carried out it
must be decided where the four corners of the trans
verse section in the transformed space are to appear
in the physical space. A glance at Fig. 2 will help to
recognize that the choice is obvious in the upstream
part of the domain but less so behind the stern. Re
ferring to Fig. 3, the location of point A correspond
ing with A' can be taken somewhere between B and
D or even coinciding with either of these points. In
practice the final choice is primarily governed by the
requirement that the (lines must roughly be aligned
with the main flow. The points A are therefore usu
ally located on a straight line extending from the keel
(cf. Fig. 2~.
z
oriented physical velocity variables are used, so that
v and w are the 71 and ~ component of the physical
velocity. Then we must have v=O, w=finite along AD
(Fig. 3), while w=O, v=finite along AB. It leads to
a conflict at A because v and w cannot both be zero
in A. Although this problem can be circumvented by
e.g. grid staggering, the accuracy of finitedifference
approximations is likely to be low in the vicinity of
the singularity due to the sudden changes of some of
the variables.
We alleviate the difficulty by a special choice of
the velocity variables, which may be explained as fol
lows. The coordinate transformation near the grid
singularity in the wake tends to
y + it—z&) = in(Ry + it (1)
where i = x/=T, R = y/~ is the grid cell aspect
ratio, c' is a proportionality constant while the singu
larity is located a+ the point y = 0,z = z, . Because
of the relation (1) the singularity is called parabolic.
Formula (1) is a wellknown type of transformation
in potential flow theory. It can be used to transform
a parallel flow to a corner flow. If the parallel flow
is supposed to occur in the physical space and is de
scribed by a velocity potential
~ = VT (` Z—Z3, Y. (` Via is constant ) ,
D  ~ SAC D' C we find with the aid of (1)
/ ~ = ~Ta(~2 _ R202y
AL ~ A L I B V2; 0~ 2aR: Van;
Fig. 3. Crosssection of physical (left ~ and trans
formed (right) computation domain behind
hull
The coordinate transformation is singular in such
points: the Jacobian of the transformation vanishes.
Singularities may also appear elsewhere on the bound
aries of the domain  e.g. along the waterline when
the frame line of the hull does not meet the symmetry
plane (undisturbed free surface) at a right angle  but
the mapping technique used excludes their occurrence
in the interior of the domain.
Grid singularities may deteriorate the accuracy of
the calculations or destabilize the iterative solution
process ~7~. This is particularly true when there is
a finite velocity in the singularity as is the case be
hind the hull. Suppose for example that coordinate
We see that the velocity components V2,V3 be
have nicely in the transformed space. They vanish at
the corner and increase linearly away from it so that
accurate discretisation is possi ble.
It is important now to recognize that V2,V3 are
covariant quantities since ~ is a scalar. Apparently,
covariant velocity components behave regularly near
a parabolic singularity in contrast to physical or con
travariant components as can readily be verified. For
example the contravariant velocity component v2 is
in an orthogonal grid (as we have behind the hull)
given by
V2 = 922 V2—V2 /g22 =  2OR1/Tr1/
and tends to infinity for 71 ' O
~9
OCR for page 87
Formerly, we employed the contravariant veloc
ity components directly as dependent variables in the
mathematical formulation. Denoting them by Vi (i =
1,2,3), we define now the three velocity variables as
U = Vi; V = 922V2; W = ~V3
instead. Notice that v and w are not strictly covariant
velocity components. But the multiplication factors
for both v2 and V3 involve the Jacobian ~ (922 =
R4), which is sufficient to obtain regular behaviour
near the singularity.
If the y, z velocity components in the physical
space must be determined from these variables, dif
ficulties are encountered at the singularity. There is
no trouble, however, in evaluating them in neighbour
ing grid points whereupon the required values at the
singular point can be derived by interpolation.
4. Boundary conditions
The sternflowandwake zone has eight bound
aries (Fig. 2~: the hull surface, an inlet and an outlet
boundary and the external boundary, the remaining
boundary surfaces being located in symmetry planes.
The terminology is intentional: the velocity compo
nent u is assumed to be positive on both the in and
outlet boundaries. On the other hand the flow may
enter or leave through the external boundary.
Like the NavierStokes equations themselves, the
Reynoldsaveraged NavierStokes equations for incom
pressible flow require three boundary conditions on
all boundaries [e.g. 84. However, if diffusion along the
coordinate associated with the predominant flow di
rection is neglected, the equations exhibit Euler char
acter in that direction with a consequent change in
the boundary condition requirements. Three condi
tions are still required at the inflow boundary, but
only one condition must be imposed at the outflow
boundary t94.
It is obvious what conditions should be applied on
the natural boundaries of the computation domain,
viz. the hull surface and the symmetry planes. On
the hull surface the noslip and the impermeability
condition give us three Dirichlet conditions for the
velocity components. On a symmetry plane there is
one Dirichlet condition for the normal velocity com
p<:,nent and a Neumann condition for each of the two
other components.
It is perhaps worth mentioning here that we im
pose the hull boundary conditions directly by calcu
lating the flow down to the hull and not indirectly via
the use of wall functions as an approximate descrip
tion of the nearwall flow behaviour.
The remaining boundaries are artificial bound
aries and it is not immediately clear what are suitable
conditions. On the outlet boundary only one condi
tion must be prescribed, either for the pressure or for
the normal velocity (i.e. u) t94. In view of our dis
cretisation it is natural to choose a condition for the
pressure. A Neumann type condition is to be pre
ferred because it is less restrictive and allows some
pressure variation over the outlet plane which may
be important in view of the frequent occurrence of
longitudinal vortices in a ship's wake; the existence
of such vortices demand the pressure to be somewhat
lower in the vortex core than in its surroundings.
Mostly we use a vanishing longitudinal pressure
gradient at the outlet plane, but sometimes a non
zero value, derived from potential flow calculations,
is preferred.
A remark should be made here about the numer
ical implementation of the pressure boundary condi
tion. In a spacemarching scheme the discretisation of
the streamwise pressure gradient involves a pressure
value of the current as well as one of the preceding
global iteration, in our discretisation:
(` Pi + ~ pi ~ / /\ x i
where rl counts the global iterations and i is the
station index. When the solution on the outlet
plane at station any is to be determined, a value for
pnN~+1 must be derived from the boundary condition,
Pc = pgrad say. One might feel tempted to simply
replace the pressure gradient as discretised above by
the given boundary condition on the outlet plane:
(\PN,\ +1 —PNh )//\ XN,Y = pgrad
It leads to unquiet solution behaviour near the out
let plane which eventually may deteriorate the global
convergence. The correct procedure is to apply the
Neumann condition for the pressure for each iteration
level separately. Thus
(\PN ~ +]—PN V )//\X~,\ = pgrad for all n
which implies that the pressure gradient term in the
(momentum equation must be 1nodelled on the out
let boundary as:
pgrad F (PN V — PN V )//\~N V
i.e. t he discretisatior, must include the pressure change
wit l1 resl)( ct t`) the prev;(:)us sweep.
The three conditions to be specified on the inlet
boundary must, again in view of the discretisation,
apply to) the velocity components. Dirichlet condi
90
OCR for page 87
tions may be chosen but since inlet conditions are
obtained either from usually incomplete experimental
data or from thinboundarylayer calculations carried
out in an approximate grid geometry, they are often
imperfect, leading to a nonsmooth solution near the
inlet plane. This is the reason why we mostly replace
the Dirichlet condition for v by a Neumann condition.
Let us next turn to the boundary conditions on
the external boundary. The old practice was to spec
ify the pressure and the two tangential velocity com
ponents. They were obtained from a potential flow
calculation for the bare hull, neglecting the displace
~ent effect of the boundary layer on the external
flow. Such conditions are good enough if the external
boundary is chosen far (several boundarylayer thick
nesses) from the hull. But that is an unattractive op
tion if one aims at computational economy because
the mathematical model is then unnecessarily com
plicated in a considerable part of the domain. On
the other hand, the conditions are evidently not ex
act when the extent of the computation domain is
reduced. It may result in the formation of a weak
nonphysical boundary layer on the external bound
ary t94. After a close inspection we have indeed found
ample evidence of its occurrence in our calculations
in which we have always tried to use relatively small
domains.
A remedy is of course to correct the boundary con
ditions for viscousinviscid interaction in the course
of the solution process. Stern et al. t10] have com
pared a procedure of that type with a large domain
approach and found it to be effective for two rela
tively simple test cases. Surprisingly, they applied
the displacement body concept to update the exter
nal flow; surprisingly because of the ambiguity in the
definition of the displacement thickness and the ap
proxima.te nature of the displacement body concept
in viscousinviscid interaction studies. A more accu
rate and unarnbiguc~us procedure is to use the normal
velocity on the external boundary resulting from the
viscous flow calculations as a boundary condition for
repeated potential flow calculations for the exterior
of hull plus sternflowandwake zone. A practical
implementation might be as follows.
The potential flow around the bare hull may be
obtained by, say, a. discrete source distribution on the
hull surface; the source strengths are found by im
posing the impermeability condition and they are de
noted here by a,. Their induction yields as it were
,, ~
a first guess of the flow, being inviscid and irrota.
tional in all three zones of Fig. 1. Moreover, it gives
a certain velocity fin the t'`:,unda.rics between zones I
anal It and Ones I and Ill, respectively, Oldish inlay
be decomposed into a normal (V,~,po~) and a tangen
tia.l (`V',pOt,`) component. Keep in mind now that the
part of this basis flow covered by rep;:ion 1 of flog. 1
~ ~ · T ~ I ~
~' V
millet I,e reproduced by a source distribution a2 on
the boundary of the union of hull, zone lI and zone
III, if the normal. velocity (V,,,po`) would be used as a.
bounty c`.'nditi`:,n.
The viscous flow calculations in zones II and III,
with (v`,,pot,) as a boundary condition, will yield a dif
ferent normal velocity (vnvi,C) on the boundaries be
tween the zones because of boundarylayer displace
ment effects. The related adjustment of the external
flow might be calculated by using (vnvi`,C) as a new
boundary condition to yield a new source distribution
CJ3 defined on the same boundary as a2. However, as
will become clear in a. moment, it is better to oper
ate with /\(T2 = CT3—cr2, the correction of ~2 due to
viscousinviscid interaction. So we suggest to calcu
late /\CJ2' associated with a correction potential Icon'
under the condition
I31)cor
on = Vn,visC —Vn,pot,
to evaluate the tangential velocity induced by them
and to add it to (`v~,pO','). The pressure follows from
Bernoulli's law which completes the set of new bound
ary conditions for the next viscous flow calculation.
The process can be repeated until the solutions in the
three zones match well enough.
The advantage of operating with /\
aware, the use of vorticity boundary conditions in
primitive vartable formulations has nest been applied
earlier although it has been suggested by Roe ills.
One might wonder whether such conditions guar
antee a unique solution. After all, the conditions are
valid for any external potential flow while no informa
tion on the velocity at infinity is conveyed by them.
How is the correct solution selected from the many
possible ones?
Although a rigorous analysis is difficult and has
not been completed yet, we think that the new bound
ary conditions give the solution appropriate for un
bounded flow conditions. It is clear that an out
side influence, disturbing the external flow, such as
a nearby wall, would not be felt, except via the con
ditions specified on the inlet boundary. But in the
absence of outside disturbances, the establishment of
a certain external flow is fully governed by what hap
pens inside the sternflowandwake zone (displace
ment effects of boundary layer and wake); there is no
need for a communication of particulars of the exter
nal flow to the sternflowandwake zone other than
that it is an irrotational flow. The conditions on the
inlet boundary make sure that we obtain the flow at
the correct Reynolds number.
At this moment, we have a limited experience
with the application of vorticity boundary conditions.
The test results seem to confirm that the correct so
lution is obtained but unfortunately they have all
shown so far a reduced convergence rate. After a
while the global convergence is fully governed by the
changes on the external boundary which decrease very
slowly. An example is given in Fig. 4, showing global
convergence results for the Wigley hull. As a mea
sure for the rate of convergence we use the rate of
decay of the maximum pressure change between suc
cessive iterations anywhere in the computation do
main. Besides two lines representing the results for
the old and the new boundary conditions, respec
tively, a dashed line is shown indicating the maxi
mum pressure change on the external boundary. It
is obvious from the figure that a grid refinement has
been applied between iteration numbers 3 & 4 and 6
& 7; the new boundary conditions were imposed on
the finest grid only.
We may note that the slower convergence is not
the result of a slow drift of the pressure level; the
pressure adjustments continue to be partly positive
and partly negative.
I TH 01 R ~ CHLET BOUNDARY CONO ~ T I ENS
W I TH VORT I C I TY BOUNDARY COND I T I ONS

xN
~ ,
Q
~ _

C3 ~r~
O ,_


.
, _
u'_ 1
0 4 8 12
ITERATION NO.
\1\
y \N
\\\
\~
it.
1 1 1
16 20
Fig. 4. Influence of boundary conditions on the
global convergence for the flow around the
Wigley hull
1
0.00 0.40 0.80 1.20
2X/L
l
l
1.60 2.00
Fig. 5. Calculated pressure change due to viscous
inviscid interaction
92
OCR for page 87
To illustrate that viscousinviscid interaction is
automatically taken into account with the new bound
ary conditions, we show a result of the calculated
change of the pressure on the external boundary with
respect to the data obtained from the potential flow
calculation around the bare hull. In Fig. 5 the pres
sure changes on the intersection of the external bound
ary and the horizontal symmetry plane (undisturbed
free surface) for the flow around the Wigley hull are
presented. The figure shows the typical effect of vis
cousinviscid interaction, viz. a pressure drop in a
region in the vicinity of the stern bordered on both
the upstream and the downstream side by a region
of pressure rise. The slight increase of the pressure
change towards the outlet station has been verified
in an axisymmetric flow case to be a consequence of
the inexactness of the boundary condition pa = 0; it
disappears when a slightly negative pressure gradient
is imposed as a boundary condition.
We may conclude that the application of vortic
ity boundary conditions is worth to be investigated
further. The only deficiency that we have observed is
the reduced convergence rate caused by it. Attempts
to remove this deficiency have not yet been made,
however. As an aside we note that it is questionable
whether with the alternative technique of repeated
potential flow calculations the matching accuracy can
ever be as good as /`Cp = 103  achieved in Fig. 4
after 17 iterations  in view of the use of two totally
different discretisation methods for, respectively, the
viscous and inviscid part of the flow.
5. Global relaxation
After discretisation and linearization, the differ
ence equations for all grid nodes can be put together
in a matrix/vector equation
Ad= b
where A is a sparse square matrix, ~ is the vector of
unknowns and b is a known vector (but containing
previous iterates of the velocity components due to
the nonlinear convection terms). Let the entries of
A be grouped in blocks so that all entries associated
with the unknown variables in a certain transverse
plane form a block. If such a block is represented by
one element of a new matrix B and a corresponding
grouping is carried through for the vectors ~ and b,
we get
BE = c
(2)
where B is a N,X' * N,X' square matrix with (in our
discretisation) four nonzero diagonals and ~6 and c
are vectors of size N ~ The it.er~.t.i~r'> ;:r~lilt.imn of thin
system has been called in Section 2 the global or outer
iteration process, to be distinguished from the local
or inner iteration process needed to solve the block
of equations associated with a certain transverse sec
tion.
What will be called in this paper the 'basic method'
to solve (2) is a GaussSeidel type iteration process:
B is split into a sum of a lower (L) and an upper (U)
triangular part where the elements of L are identical
with the corresponding entries of B while U contains
the upper diagonal of B only; the solution is updated
via
L~bn = c—Ut,bn~
To enhance the convergence rate of this process we
have used earlier a source term in the (momentum
equation, as originally proposed by Israeli and Lin
t12] and described for our application in t13. It has
in many cases proved to be a useful artifice, but our
experience tells that it may have undesirable side
effects, particularly when used in combination with
an algebraic turbulence model. Notably in the ini
tial phase of the calculations the source term may
assume appreciable values, generating a significant
overshoot in the velocity profiles or provoking flow
separation. In the first case the turbulence model
may show pathological behaviour in the determina
tion of the outer length scale, which deteriorates the
global convergence. In the second case the required
change of the discretisation formulae will have an in
fluence to the same effect.
We have now abandoned the source term scheme
in favour of an alternative convergence accelerator
with a superior performance. It is described below.
When an iterative process is used to solve a time
independent problem, it is often preferable to oper
ate with a transient form of the mathematical model
and to try to find the steady state solution t134. Fol
lowing this suggestion, we add a quasitime deriva
tive of the pressure to the (momentum equation.
If for the sake of compact notation we consider the
equations valid in a Cartesian coordinate system, the
momentum equation for the dominant flow direction
becomes
L(u, v, w) + PI = IT (3)
where L is a differential operator incorporating the
convection and diffusion terms but no 'time' deriva
tives.
Ltu,v,wJ + pm = pt
In the above GaussSeide] process the discretised
rno~nentu~n equation is
~^ w.~^v (Pin+ —Pin~)//~\Xi  (pin—pin~)/~t—RHSi (4)
93
OCR for page 87
where Lo is a. finitedifference analogue of L, RHSi
is a. known righthand side which does not contain
variables at iteration level n and /\t is taken equal
tc, i\x. Next, one may conceive a complementary
discretisation of (3~:
LA (`U, V, ?JJ~)n +
(Pi +~ —Pi + )//\xi—(`pn+i—pa)//\ t = RHSi (5)
such that half of the sum of (4) and (5) yields a dis
cretisation of all terms of (3) at time (or iteration)
level n, provided that
pn = I tpn+~1 + pn—~ ~ `6'
is considered as a. new value for pa.
From the combination of the formulae (4~5~6)
above, a predictorcorrector scheme for the pressure
is readily constructed. The predictor step is the basic
method, the GaussSeidel scheme:
L z~ (u, v, w in + ~ pin+ l1—pi ~ / /\ hi = RH Si ,
giving new values u=,vn~wn for the velocity compo
nents and a first approximation p* to a new value
of the pressure. The difference of (4) and (5) yields
a simple algebraic relation to determine a fictitious
pn+i The improved guess for pn is then obtained
from the mean of this pn+i and the old pressure pun.
The latter operations may be combined in the correc
tor step
pn = 1 tpn—~ + path + I pan _ pa—1' ~7y
Because of the appearance of pa+ in the righthand
side, pa must be evaluated in an upstream marching
sequence. Thus the basic process of a sweep from
inlet to outlet station is followed by reversed sweep in
which. an improved value for the pressure is obtained
via the simple algebraic relation (7~. As a start for
this reversed sweep we have
n _ *
PN V ~~ PN,\
a.s in~media.lely follows from the above formulae when
the mainstream pressure gradient is zero on the outlet
boundary.
It is easily verified by repeated substitution of (7)
that a. pressure change at station xi can be expressed
as
where
/\ pa = — [of ( 2 ) ~f\Pk + ( 2 ) /\PNX~
pro pn _ pn1 gyp* p* _ pa—1
It implies that /\ pn is determined by a fraction
(diminishing with distance) of downstream pressure
changes resulting from the basic method. Thus we
have incorporated the desirable feature of an infinite
propagation speed of pressure influences in upstream
direction (a property of the continuum equations) in
the numerical scheme. In both the basic method and
the Isra.eli/Lin source term scheme the propagation
speed is only /`xi per iteration.
Readers familiar with the recent literature on so
lution methods for the Reynoldsaveraged Navier
St<'kes valuations in flows characterized by a dominant
flow direction will have observed that the predictor
corrector scheme introduced here has some similarity
with a. procedure suggested by Davis et al. in t144.
The differences are essential, however. In t14] the two
discretisations (4) and (5) appear also but the solu
tion princess is continued with An+. That would be
allright if the xmomentum equation stood on its
own but when it is coupled with other equations in
which the pressure appears as well it is a completely
unsatisfactory approach. For /\t = /\x an essentially
unstable scheme results. Also for smaller values of /\t
we have not been able to produce acceptable results
with their method.
BRSIC SCHEME (GRUSSSEIDEL)
SOURCE TERM SCHEME (ISRRELI/LIN)
+ NEW PREDICTORCORRECTOR SCHEME
~2
Q
(A N_

C3
J
I I I I I 1 1 1 1 1
4 8 12 16 20 24
I TERRT I ON NO .
Fig. 6. Maximum pressure change between succes
sive sweeps in global iteration process for
three methods
94
OCR for page 87
As an illustration of the excellent performance of
our proposal we give global convergence results for
the flow around an axisymmetric body as obtained
by the basic scheme, the source term scheme and our
new scheme in Fig. 6. In all three convergence curves
two jumps appear which correspond each to a grid
refinement, being executed when /`Cpma~ ~ 0.02 and
< 0.01 respectively. Most striking is the gain in con
vergence rate with respect to the basic method which
differs from the new method by the correction step
(7) only. The more so because the correction step
is extremely simple, requiring a completely insignif
icant amount of computation time. Similar results
have been obtained for other cases.
6. Application
Since the appearance of our earlier publications
t1,2,3], the number of applications of our method has
considerably increased. Some of these applications
were repeated calculations for 'old' test cases (e.g.
Wigley hull, SSPA 720), others concerned new explo
ratic~ns. Among them were computations including a
representation of a propeller and/or a duct by a spec
ified external force field. Instead of giving results of
some of those exercises or improved results for the
test cases that appeared in previous publications, we
have selected for presentation some data obtained for
a rather demanding test case: the HSVA tanker.
What makes this hull form to a difficult case for
flow computations is in the first place the fullness of
the hull (block coefficient = 0.85), implying a high
viscous pressure resistance and presumably complex
flow with a risk of boundarylayer separation.
The body plan of the hull is given in Fig. 7. A
2.74 m model of it has been subjected to detailed
measurements in a slottedwall wind tunnel by
Wieghardt and Kux t15] and a data tape with a part
of their results has been kindly made available to us.
On the basis of the geometry information provided on
the data tape, a reconstruction of the hull shape was
carried out with our hull Airing system. Intermediate
frame lines were obtained by spline interpolation.
Calculations were made at the Reynolds number
of the measurements: Rn = 5 * 106. Neither the wind
tunnel wall nor the rod extending from the back of the
wind tunnel model were included in the simulation.
The lengthwise extension of the sternflowanclwake
zone was chosen as ().5 < 2;~/L < 1.65, where L is the
length between perpendiculars and x _ O midships.
The results to be presenter] were obtai~ecl on a 4~; *
zl9 * 29 grid.
Potential flow and boundarylayer calculations pro
vided us with the nontrivial part of the boundary
conditions. On the external boundary, conventional
conditions (prescription of tangential velocity com
ponents and pressure from the potential flow) were
applied, while a zero pressure gradient condition was
imposed on the outlet boundary. The initial guess for
the pressure field was obtained by assuming pa = 0.
One of the aims of the calculation of the stern
flow field is to acquire information on the velocity
distribution near the location where the propeller is
to be mounted. Let us therefore first have a look at
some contour plots of the axial velocity field near the
stern. Fig. 8 shows the comparison of measurements
and predictions in six transverse sections. In the mea
surements, the flow disturbance caused by a support
wire in the wind tunnel is clearly visible. When that
anomaly is taken into consideration there is gener
ally good agreement between both sets of data in the
outer part of the boundary layer or wake. Near the
hull, on the other hand, in the region just below the
concave part of the frame lines, the measurements
indicate a much stronger retardation of the flow. S
shaped isovelocity contours show up which do not
clearly appear in the computed results.
Usually the formation of that kind of wake distri
bution is connected with the development of a longi
tudinal vortex. How well that phenomenon is repro
duced by the calculations may be judged from Fig. 9
It shows vector plots of the transverse velocity com
ponents at various stations. From the measurements
as well as the calculations only a part of the available
information has been plotted but such that the vector
lengths are directly comparable.
~ =~;I! At
1E
it\ \ \I4\ql 119Y7 1 1
t
o
95
l\` \ \ \ 1 / I I I~
\\\ ~ 11 /~J//JJ:
Fig. 7. Body plan HSVA tanker
,?
to
9
OCR for page 87
  
1
1
1
1
~ art
K~'
ran
I No
L 4'
~/l
~~
r
Calculated at 2~/l = 0.880 L Calculated at 2x/L = 0.90
 Measured at 2~/L = 0.885 _~ Measured at 2~/L = 0.90
7
L Calculated at 2x/L = 0.94
_____ Measured at 2~/L = 0.94
I____
it !
://J~ J 1
L Calculated at 2~/L = 0.980
______ Measured at 2/L = 0.977
r ~ i ~
1 " ' 1
Calculated at 2z/L = 1.02 ~ Calculated at 2/L = 1.10
______ Measured at 2~/L = 1.02  Measured at '22/L = l.O9
Fig. 8. Contour plots of axial velocity distributions
96
OCR for page 87
There is no clear evidence that the vortex for
mati<'n is predicted to start at a wrong position but
the vortex core seems to stay closer to the vertical
symmetry plane than in the measurements, especially
when is taken into account that the measurements
are not truly symmetric with respect to the geomet
ric symmetry plane. This is in accordance with the
stronger deceleration of the axial flow.
The question remains what is the cause of the
discrepancy in the axial velocity distribution. We
may observe that there is little reason to suppose that
a gross misprediction of the pressure is responsible.
:'~
Calculated at 2~/L = 0.880
Ad_
Calculated at 2x/L = 0.90
Fig. 10 shows the girthwise pressure distribution on
the hull at station 2~/l = 0.88 with an excellent
agreement between measurements and calculations
while the potential flow results deviate considerably.
A preliminary exercise with vorticity boundary
conditions on the external boundary indicated that
the inclusion of viscousinviscid interaction does not
bring about either a change in the nearhull flow to
the extent required for a reconciliation of measure
ments and predictions. Still to be investigated are
the effects of grid refinement and of a variation in the
turbulence model.
tTtfTtftiltN t ~ ~
Measured at 2~/L = 0.885
ItittilltNtN ~ ~
Measured at 2x/L = 0.90
Fig. 9. Vector plots of transverse velocity components
97
\ \
OCR for page 87
f ~~ ill_ ~ :
V~\ \ ~ \ ~ lilititt\Ntt ~ \ \ \
(3alc~llate:d at '2~/L  0.920 Measured at 2x/L = 0.926
rig. · ~ ~ ~
Calculated at 2x/L = 0.980 Measured at 2x/L = 0.977
Calculated at 2x/L = 1.020 Measured at 2~/L = 1.015
Fig. 9 (continued)
98
OCR for page 87
CRLCULRTED IN VISCOUS FLOW
A
ID
° CRLCULRTEO IN POTENTIAL FLOW A_
Al ~
MERSUREO
~ ,
/
///
/
/
 o 00 0.20 0.40 0.60 O.80 1 .00
NORM . D I STANCE RLONG FRAME
Fig. 10. Girthwise pressure distributions at station
2x/L = 0.88
It is not apparent from Fig. 8, but our results indi
cate boundarylayer separation at two locations: in a
small region just above the level of the imaginary pro
peller axis, starting from about station 2x/L = 0.9,
and near the waterline (horizontal symmetry plane)
aft of station 2x/L = 0.97. In both cases the thick
ness of the separation bubble is very small. Only an
oil flow experiment for the visualization of the lim
iting streamlines could shed some light on the faith
fulness of that prediction. However, the girthwise
skin friction distribution at station 2x/L = 0.88, just
ahead of one of the separation regions, is in fair agree
ment with the experimental data t16], as Fig. 11
shows.
It may be concluded that many features of the
experimental data are satisfactorily reproduced. In
deed, the present results are perhaps the most realis
tic canes among the few that have been reported so far
for this test case. Yet, they can hardly be considered
good enough, if the available set of measurements is
considered to satisfy high quality standards. Where
our method has been set up with the idea that it
should give ultimately reliable results for a case like
the HSVA tanker, further improvements are needed.
They are to be found in a finetuning of the turbu
lence model en cl the use of finer grids rather than in
the numerical scheme.
1 0
0 °
\ c)
/
CRLCULRTED
MEASURED
of
._ 1 1 1 1 1 1 1 1 1
No 00 0.20 0.40 0.60 0.~ 1.00
NORM . D I STANCE RLONG FRAME
Fig. 11. Girthwise skinfriction distributions at
station 2x/L = 0.88
7. Concluding remarks
.
In this paper we have given a detailed description
of new developments in certain aspects of a computer
code for the prediction of the stern flow field of a ship:
· We have explained how coordinateoriented ve
locity variables may be chosen so as to maintain
a high numerical accuracy near the grid singu
larity in the wake;
~ We have presented some results of our first ap
plication of a new set of boundary conditions
on the external boundary of the computation
domain, which allow this boundary to be cho
sen close to the edge of the boundary layer and
wake i
We have shown how the convergence rate of
the global relaxation of the solution may be
improved by an easy algebraic update of the
pressure field after every iteration cycle.
Moreover results have been given of the application
of the code to a difficult test case which indicate that
the numerical scheme can cope with the complex flow
around a full tanker, although the correlation with
the available experimental data should be further im
proved.
99
OCR for page 87
References
1. Hoekstra, M. and Raven, H.C., "Ship boundary
layer and wake calculation with a parabolised
NavierStokes solution system", 4th Internation
al Conference on Numerical Ship Hydrodynam
ics, Washington D.C. (1985~.
2. Hoekstra, M. and Raven, H.C., "Application of
a parabolised NavierStokes solution system to
ship stern flow computation", Osaka Interna
tional Colloquium on Ship Viscous Flow, Osaka
(1985~.
3. Hoekstra, M., "Computation of steady viscous
flow near a ship's stern", in Notes on Numer
ical Fluid Mechanics, Vol. 17 (ea. Wesseling) 11
(Vieweg) (1986~.
4.
. Bube, K.P. and StriLwerda, J.C., "Interior reg
ularity estimates for elliptic systems of differ
ence equations", SIAM Journal of Numerical
Analysis, Vol. 20, No. 4 (1983~.
Rubin, S.G., "Incompressible NavierStokes and
parabolised NavierStokes formulations and com
putational techniques", in Computational Meth
ods in Viscous Flows, Vol. 3 in the series Recent
Advances in Numerical Methods in Fluids (ea.
Habashi) (Pineridge Press) (1984~.
6. Hoekstra, M., "Coordinate generation in sym
metrical interior, exterior and annular 2D do
mains using a generalised SchwarzChristoffel
transformation", in Numerical Grid Generation
in Computational Fluid Dynamics, (ea. Hauser 15
& Taylor) (Pineridge Press) (1986~.
7. Eriksson, LE., "A study of mesh singularities
and their effects on numerical errors", FFA TN
198410 (1984~.
8. Strikwerda, J.C., "Finite difference methods for
the Stokes and NavierStokes equations", SIA M
Journal of Scientific and Statistical Computa
tions, Vol. 5, No. 1, pp. 5668 (1984~.
9. Oliger, J. and Sundstrom, A., "Theoretical and
practical aspects of some initial boundary value
problems in fluid dynamics", SIAM Journal of
Applied Mechanics, Vol. 35, No. 3, pp. 419446
(1978~.
10. Stern, F., Yoo, S.Y. and Patel, V.C., "Interac
tive and largedomain solutions of higherorder
viscousflow equations", AIA A Journal, Vol. 26,
No. 9 (1988~.
. Roe, P.L., "Remote boundary conditions for
unsteady multidimensional aerodynamic com
putations", Computers & Fluids, Vol. 17, No.
1, pp. 221231 (1989~.
12. Israeli, M. and Lin, A., "Numerical solution and
boundary conditions for boundarylayer like
flows", 8th International Conference on Numer
ical Methods in Fluid Dynamics, Aachen (1982~.
13. Roache, P.J ., "Computational Fluid Dynam
ics", Hermosa Publishers, Albuquerque, N.M.
87108 (1972~.
14. Davis, R.T., Barnett, M. and Rakich, J.V., "The
calculation of supersonic viscous flows using the
parabolised Na.vierStokes equations", Comput
ers & Fluids, Vol. 14, No. 3, pp. 197224
(1986~.
Wieghardt, K. and Kux, J., "Nomineller nach
strom auf Grund von Windkanalversuchen",
Jahrbuch STG 74, pp. 303318 (1980~.
16. Hoffmann, H.P., "Untersuchung der 3dimen
sionalen, turbulenten Grenzschicht an einem
Schiffsdoppelmodel im Windkanal", IFS Bericht
No. 343 (1976~.
00
OCR for page 87
DI SCUSS ION
by Y. Kodama
You assumed a predominant f low direction
and omitted a few terms. I'd like to hear your
opinion about the ef feet of those neglected
terms .
Author's Reply
The effect of those neglected terms on the
final solution is of no practical significance
in ship stern f lows. They do however inf luence
the boundary condition requirements on the
out l e t boundary.
101
OCR for page 87