|
Items in cart [0] |
Questions? Call 888-624-8373 |
| Copyright © 2009. National Academy of Sciences. All rights reserved. Terms of Use and Privacy Statement |
Below are the first 10 and last 10 pages of uncorrected machine-read text (when available) of this chapter, followed by the top 30 algorithmically extracted key phrases from the chapter as a whole.
Intended to provide our own search engines and external engines with highly rich, chapter-representative searchable text on the opening pages of each chapter.
Because it is UNCORRECTED material, please consider the following text as a useful but insufficient proxy for the authoritative book pages.
Do not use for reproduction, copying, pasting, or reading; exclusively for search engines.
OCR for page 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 finite-difference discretisation of
the Reynolds-averaged Navier-Stokes equations on a
boundary-fitted 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 Navier-Stokes 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)undary-layer zone, covering the forward part
of the hull surface, where first-order boundary-layer
theory is supposed to provide an adequate description
of the flow; the stern-flow-and-wake zone which is of
87
OCR for page 88
primary concern in this paper. In the latter zone we
solve eRectively the Reynolds-averaged Navier-Stokes
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 flow-governing 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 Marker-and-Cell type methods.
The equations are discretised on a boundary-fitted
curvilinear grid. It is a single-block regular grid of
N V * NY * NZ nodes. The grid is partly non-ortho-
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 space-marching 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 main-stream 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
cross-section of the computation domain. Iteration is
required both by the non-linearity 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 89
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 3-D mesh is completed. It leads inevitably
to non-orthogonality 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 finite-difference
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 well-known 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. Cross-section 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 90
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 stern-flow-and-wake 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 Navier-Stokes equations themselves, the
Reynolds-averaged Navier-Stokes 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 no-slip 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 near-wall 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 space-marching 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 91
OCR for page 93
OCR for page 94
OCR for page 95
OCR for page 96
OCR for page 97
OCR for page 98
OCR for page 99
OCR for page 100
OCR for page 101
OCR for page 102
Representative terms from entire chapter:
boundary conditions
tions may be chosen but since inlet conditions are
obtained either from usually incomplete experimental
data or from thin-boundary-layer calculations carried
out in an approximate grid geometry, they are often
imperfect, leading to a non-smooth 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 boundary-layer 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
non-physical 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 viscous-inviscid 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 viscous-inviscid 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 stern-flow-and-wake 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 boundary-layer 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
viscous-inviscid 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 stern-flow-and-wake 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 stern-flow-and-wake 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
To illustrate that viscous-inviscid 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-
cous-inviscid 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 = 10-3 - 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 non-linear 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 non-zero 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 Gauss-Seidel 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 quasi-time 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 Gauss-Seide] process the discretised
-rno~nentu~n equation is
~^ w.~^v (Pin+ —Pin-~)//~\Xi -- (pin—pin-~)/~t—RHSi (4)
93
where Lo is a. finite-difference analogue of L, RHSi
is a. known right-hand 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 predictor-corrector scheme for the pressure
is readily constructed. The predictor step is the basic
method, the Gauss-Seidel 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 right-hand
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 _ pn--1 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 Reynolds-averaged 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
all-right if the x-momentum 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 (GRUSS-SEIDEL)
SOURCE TERM SCHEME (ISRRELI/LIN)
+ NEW PREDICTOR-CORRECTOR 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
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 boundary-layer 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 slotted-wall 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 stern-flow-ancl-wake
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 boundary-layer calculations pro-
vided us with the non-trivial 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 iso-velocity 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
- - -
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
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 viscous-inviscid interaction does not
bring about either a change in the near-hull 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
\ \
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
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 boundary-layer 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 fine-tuning 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 skin-friction 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 coordinate-oriented 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
References
1. Hoekstra, M. and Raven, H.C., "Ship boundary
layer and wake calculation with a parabolised
Navier-Stokes 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 Navier-Stokes 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 Navier-Stokes and
parabolised Navier-Stokes 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 2-D do-
mains using a generalised Schwarz-Christoffel
transformation", in Numerical Grid Generation
in Computational Fluid Dynamics, (ea. Hauser 15
& Taylor) (Pineridge Press) (1986~.
7. Eriksson, L-E., "A study of mesh singularities
and their effects on numerical errors", FFA TN
1984-10 (1984~.
8. Strikwerda, J.C., "Finite difference methods for
the Stokes and Navier-Stokes equations", SIA M
Journal of Scientific and Statistical Computa-
tions, Vol. 5, No. 1, pp. 56-68 (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. 419-446
(1978~.
10. Stern, F., Yoo, S.Y. and Patel, V.C., "Interac-
tive and large-domain solutions of higher-order
viscous-flow 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. 221-231 (1989~.
12. Israeli, M. and Lin, A., "Numerical solution and
boundary conditions for boundary-layer 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.vier-Stokes equations", Comput-
ers & Fluids, Vol. 14, No. 3, pp. 197-224
(1986~.
Wieghardt, K. and Kux, J., "Nomineller nach-
strom auf Grund von Windkanalversuchen",
Jahrbuch STG 74, pp. 303-318 (1980~.
16. Hoffmann, H.P., "Untersuchung der 3-dimen-
sionalen, turbulenten Grenzschicht an einem
Schiffsdoppelmodel im Windkanal", IFS Bericht
No. 343 (1976~.
00
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
