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 882

lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
COMBINING ACCURACY AND EFFCIENCY WITH ROBUSTNESS IN SHIP STERN FLOW COMPUTATION 882
Combining Accuracy and Effciency with Robustness in Ship Stern
Flow Computation
A.van der Ploeg (Maritime Research Institute, The Netherlands)
L.Eça (Instituto Superior Técnico, Portugal)
M.Hoekstra (Maritime Research Institute, The Netherlands)
ABSTRACT
Using a mature RANS-solver for ship stern flows as a starting point, we consider some possibilities to enhance
efficiency and robustness without sacrificing accuracy. Among them are alternative domain decompositions and
vectorizable preconditioning techniques. It is shown that a good balance between flexibility, robustness and efficiency can
be obtained and that it is justified to neglect diffusion in main stream direction and some of the convective terms.
1 INTRODUCTION
For the use of viscous flow computations in the design process of ships, robustness of the numerical solution
procedure is a prerequisite. Practically this means that under reasonable demands on grid quality, convergence of the
iterative solution is warranted in flow simulations under typical ship operating conditions. However, robustness is in
several methods obtained by sacrificing accuracy (for example by using first-order discretisation). A dangerous route,
because the results of a robust code have no practical value if not a certain level of accuracy can be achieved. It is one
reason why considerable emphasis is put nowadays on verification and validation. Moreover, in view of usual time
constraints in ship design, it is desirable that a good quality solution can be obtained at low cost in a small amount of time
(efficiency). This paper aims to show how a good balance between robustness, accuracy and efficiency can be obtained in
the computation of the flow around a ships stern by solving the RANS-equations in a non-conventional way.
Most methods for the computation of the viscous flow around ships employ either the pressure-correction or the
artificial compressibility method. Recognizing that particularly for flows at very high Reynolds numbers these iterative
methods have difficulties in, respectively, restoring the coupling between momentum and mass conservation equations
and reaching a steady-state solution, an alternative approach has been chosen. This has grown out to the code
PARNASSOS, which is MARINs proprietary RANS-solver for ships (Hoekstra & Eça 1998, Hoekstra 1999). This code is
currently in use for quality assessment of hull designs on request of ship yards, navies and other customers. In our
approach the coupling between the equations is maintained in the iterative solution, an important factor for robustness.
Reduction of the size of the equation system is achieved by dividing the computation domain into sub-domains. Each sub-
domain is a grid plane roughly perpendicular to the main stream along the hull. In the iterative process these sub-domains
are visited in downstream direction, which is obviously the best choice. By a special update of the pressure field after
such a sweep through the domain, the convergence rate is enhanced. Between 40 and 100 sweeps are usually sufficient to
reduce the maximum norm of the changes of the non-dimensional variables between successive iterations to below 10−4.
Moreover elaborate verification studies have been done at model scale Rn by Hoekstra & Eça 1999 and at full scale Rn by
Eça & Hoekstra 1999.
After a brief outline of the original features of PARNASSOS in Section 2, this paper focuses on a number of new
elements which bring us closer to an optimum in combining robustness, accuracy and efficiency.
1. First, the possibility is included to choose larger sub-domains. The single-grid-plane sub-domain has
originally been chosen to minimize memory requirements, but on present-day machines memory is amply
available. So it is now possible to choose the size of the sub-domains in the range from one grid plane to even
the complete domain. In Section 3.1 we will describe this technique in more detail.
2. The systems of linear equations are solved with preconditioned GMRES. The preconditioning
the authoritative version for attribution.

OCR for page 882

lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
COMBINING ACCURACY AND EFFCIENCY WITH ROBUSTNESS IN SHIP STERN FLOW COMPUTATION 883
technique, which performs well on a vector processor and is suitable for implementation on parallel
machines, is partly responsible for efficiency and flexibility: it allows to choose the size of the sub-domains
between only one grid plane and the entire computational domain. It will be described in some detail in
Section 3.3.
3. Originally, we discarded the diffusion in the mainstream direction and also a part of the convection terms in
reversed-flow zones. With larger sub-domains as described under 1 these terms can readily be included
implicitly, while in the plane-by-plane version they can be incorporated explicitly. (Section 3.4). It will be
shown that these terms are insignificant indeed in representative circumstances, but that nevertheless the
behaviour of the iterative solution is improved by including them.
In Section 4 we illustrate the effects of the above additions to PARNASSOS by showing results of application to the
flow around the stern of the HSVA and Mystery tankers. Section 5 describes some current developments about how a free
surface can be incorporated, and in the final Section 6 the main achievements are summarized.
2 ORIGINAL SOLUTION STRATEGY
As mentioned before, PARNASSOS solves the fully coupled steady momentum and mass conservation equations in
primitive-variable form, without resorting to pressure correction methods or the artificial compressibility approach.
Additional transport equations associated with the turbulence model are treated as uncoupled from mass and momentum
equations. The grids are assumed to be body-fitted, but maybe generally non-orthogonal. They are stretched towards the
hull in order to resolve the gradients in the boundary layer. The governing equations are integrated down to the wall (no
wall functions used) even for full-scale Reynolds numbers. The inherently very high aspect ratio of the cells near the hull
puts high demands on the solver for the linear systems, which is one of the reasons to maintain the coupling between the
equations in the iterative solution.
For turbulence modelling in PARNASSOS, the concept of the isotropic eddy viscosity is used. It is possible to
choose from a wide range of turbulence models, varying from algebraic models to several two-equation turbulence
models (Eça & Hoekstra 2000).
A detailed description of the mathematical model, the computational grid and the PDEs in curvilinear coordinates is
given in Ref Hoekstra 1999. Here we restrict ourselves to a short description of the discretization. In PARNASSOS the
following finite difference approximations are used:
• For the numerical evaluation of the grid metric terms we use second-order, central difference schemes.
• In the continuity equation we use a second-order, three-point backward scheme for the mainstream derivative,
and a third-order four-point scheme with a fixed bias in normal and girth-wise direction.
• For the derivatives of the velocities that occur in the convective terms we use a second-order upwind scheme in
streamwise direction, and a third-order upwind scheme for the normal and girth-wise direction.
• For the gradients of the pressure in the momentum conservation equations we use a third-order scheme. For
stability reasons, the bias has to be opposite to that of the ‘corresponding' derivative in the continuity equation
(Hoekstra 1999).
• All second derivatives in the diffusive terms are discretized by second-order central-difference schemes.
Hence all terms in the momentum and mass conservation equations are discretized at least second-order accurately.
To avoid negative turbulence quantities, only in the transport equations in the turbulence models we use a first-order
upwind scheme for convection, while for diffusion the same discretization is applied as in the momentum conservation
equations.
Application of the afore-mentioned discretization leads to a huge set of non-linear algebraic equations, to which
quasi-Newton linearization is applied. In order to reduce the size of the discrete equation system, PARNASSOS uses a
marching solution scheme. The velocities and pressure of one grid plane across the mainstream direction (say, a
ξ=constant plane), are solved simultaneously. The grid planes are visited in downstream order, while the elliptic character
of the RANS-equations is numerically recovered by iteration. Each step of this iteration scheme includes not only the
downstream sweep through the computation domain, in which the eddy viscosity, the velocities and the pressure are
updated, but also an upwind sweep in which only the pressure is updated. In the sequel of this paper, this scheme will be
denoted as the ‘global iteration'.
The plane-by-plane marching yields of course enormous savings in memory requirements. But the
the authoritative version for attribution.

OCR for page 882

lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
COMBINING ACCURACY AND EFFCIENCY WITH ROBUSTNESS IN SHIP STERN FLOW COMPUTATION 884
price to be paid is that the upstream communication of downstream occurrences may proceed rather slowly via the global
iteration process. Nevertheless, a robust solution procedure has been achieved using a proper preconditioning technique.
3 NEW ELEMENTS IN SOLUTION PROCEDURE
3.1 Larger sub-domains
If memory is not a serious limitation, the size of the sub-domains can be increased. So, in the new approach, each
sub-domain is defined by a number (say, g-g≥1) of consecutive grid-planes. In a sub-domain the variables are solved
simultaneously. We construct a system of linear equations, from which we obtain a correction for the pressure and
velocities at all grid nodes in the current sub-domain. The resulting coefficient matrix is g times bigger than in the
traditional plane-by-plane approach in PARNASSOS, and contains terms representing the coupling between the g planes
covered by the sub-domain.
In the plane-by-plane solution strategy we are used to apply several steps of the quasi-Newton iteration in a sub-
domain in order to let the global iteration safely converge. It turned out that with g>1 it is possible to make only one step
of the quasi-Newton iteration before going to the next sub-domain. Hence the global iteration accounts completely for the
nonlinearity, while a significant saving in computation time is obtained.
Let us denote the number of grid nodes in streamwise, normal and girthwise direction by NX, NY and NZ,
respectively. In the extreme case of g NX (sub-domain=complete domain), the global iteration becomes the quasi-Newton
iteration for the whole problem. In that particular case, there is no need for the pressure-correcting upstream sweep, which
is normally part of a global iteration cycle.
We further note that, even if g NX, the turbulence model equation(s) is (are) still treated as uncoupled from the other
equations. Suppose that the three velocity components, the pressure and the eddy viscosity with grid indices i, j, and k are
denoted, respectively, by uijk, vijk, wijk, pijk and (vt)ijk. Then the complete non-linear problem that is obtained after
discretization can be written as
in which F is a vector-valued function with 5×NX× NY×NZ components. If n counts the number of global iteration
steps, one cycle of the global iteration process consists of the following steps:
1. compute the eddy viscosity as a function of and In other words: solve
approximately from
If a one- or two-equation turbulence model is used, the transport equation(s) are still solved using a plane-
by-plane strategy.
2. linearize the RANS-equations, using computed above, and compute and by solving the
resulting system of linear equations to a prescribed accuracy.
Hence we cannot expect a quadratic speed of convergence, not even when all convective terms have been linearized
second-order accurately.
3.2 The linear system solver
As has been mentioned above, the coefficient matrix is g times larger than in the plane-by-plane version of
PARNASSOS, and it contains coefficients which account for the coupling between the separate ξ=constant planes. This
matrix will be denoted by A. Let the entries of A be grouped in blocks so that all elements multiplying the variables in a
ξ=constant plane form a block. Let such a block, which itself is a square matrix of size 4×NY×NZ, be represented by one
entry Ai,l, in which i and l are row and column indices, with i for convenience taken equal to the grid-index of the relevant
ξ=constant plane within the current sub-domain. As a result of the chosen discretization, the block Ai,l can only be non-
zero if |l−i|≤2. Hence A has the following penta-diagonal structure:
(1)
Usually, the coefficients ei and fi only contain elements coming from the discretization of the pressure derivative in
main stream direction. Only in case of flow separation, these matrix elements can also contain small contributions from
the derivative of the velocities in main stream direction. We are faced now with the problem to solve the system of linear
equations Ax b. Direct methods are out of the question since they require too much storage and floating-point operations.
the authoritative version for attribution.
Therefore, we prefer to use an iterative method, suited for non-symmetric matrices, like

OCR for page 882

lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
COMBINING ACCURACY AND EFFCIENCY WITH ROBUSTNESS IN SHIP STERN FLOW COMPUTATION 885
Bi-CGSTAB(l) (van der Vorst 1992) or GMRES(M) (Saad & Schultz 1986). Usually, the most expensive operations in
such methods are the computation of innerproducts, vector updates, and the matrix-vector multiplication. In our particular
situation, the matrix-vector multiplication is much more expensive than an innerpoduct or vector update. Hence we
should minimize the number of matrix-vector multiplications, and therefore we choose GMRES(M) as the linear system
solver (where (M) indicates that GMRES is restarted every Mth step).
3.3 Preconditioning
The speed of convergence of GMRES strongly depends on the eigenvalue spectrum of A. Ideally, the matrix should
have all eigenvalues clustered; then A would be close to the identity matrix and GMRES will rapidly converge. In order to
bring a given spectrum closer to that ideal situation, one can use a preconditioning technique. Instead of solving Ax b, one
solves M−1Ax M−1b, in which M is some approximation of A. The matrix M should be easy to construct and for a given
vector x the computation of M−1x should be cheap. There are several possibilities: one can try to approximate the inverse
directly, for example, by using a polynomial of the original matrix A. However, the total number of matrix-vector
multiplications will then not be less than with full GMRES. Therefore, we have chosen a different kind of preconditioning
technique: we exploit the fact that the equations are almost parabolic in the streamwise direction, by neglecting the terms
ei and fi in (1) during the construction of the preconditioner. We note that GMRES still has to solve the system of linear
equations with the coefficient matrix (1) that includes these blocks. This is the essential difference with the approach in
the plane-by-plane version of PARNASSOS, in which the systems of linear equations do not have any entries to account
for the coupling between the separate ξ=constant planes, since they have all moved to the right-hand side vector b.
The matrix M is constructed in such a way that it has the block structure
(2)
and it is an approximation of
(3)
Hence mi should be an approximation of the block di. To that end, we construct a coupled incomplete LU-
decomposition of di. Then mi LiUi, in which the factors Li and Ui are sparse and lower- and uppertriangular, respectively.
The incomplete decomposition is constructed on a 4×4 block-level. Each ‘entry' in Li and Ui consists of a 4×4 block.
Both the construction and the triangular solves using the factors Li and Ui should not cost too many floating-point
operations. Another requirement is that both the construction of Li and Ui and the triangular solves can be vectorized. The
last requirement implies that the sparsity pattern of Li and Ui should be regular, so that indirect addressing can be avoided.
Let the entries of di be grouped in square blocks of size 4 so that all elements multiplying the velocity components
and pressure in a grid point form a block. Let such a block be represented by one entry. The matrix di then has a block
sparsity pattern that corresponds to the following 9-point discretization stencil
(4)
The coefficients in CNW, CNE, CSW and CSE are relatively small, and therefore, those blocks are neglected during the
incomplete LU-decomposition of di. Furthermore, all fill-in blocks are neglected. Hence the matrix Li|Ui has a block
sparsity pattern that corresponds to the following 5-point stencil
(5)
Both in the construction of Li and Ui and in the triangular solves there are no recurrence relations when looping
along a diagonal of the grid. This means that by diagonal ordering, we can vectorize the loops. This can be regarded as a
kind of hyperplane ordering.
the authoritative version for attribution.
When g 1, we use a special preconditioning technique, in which only the coefficients in CNE and CSW,

OCR for page 882

lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
COMBINING ACCURACY AND EFFCIENCY WITH ROBUSTNESS IN SHIP STERN FLOW COMPUTATION 886
together with all fill-in blocks are neglected. This preconditioning technique has not been vectorized, since computations
in which g 1 are more likely to run on scalar machines, instead of on a vector processor or a parallel computer.
As mentioned above, it is not necessary to solve the system of non-linear equations for each domain: it suffices to
apply only one step of the quasi-Newton iteration if g≥2. As a result, the required CPU-time has been reduced by
approximately a factor of 10, compared to the plane-by-plane implementation of PARNASSOS. If an algebraic turbulence
model is used the reduction in CPU-time can be even more since the number of global iterations decreases by increasing
g. If a one-equation or two-equation turbulence model is used, this favourable effect has not yet occurred, which is
probably due to the strongly nonlinear source terms.
Increasing the size of the sub-domains (i.e. increasing the value of g) leads to higher memory requirements. Below
we estimate the required memory by listing the most memory-consuming fields of variables. Herein N denotes the total
number of grid points.
• The velocities and pressure, coordinates of the grid points, pressure of the previous sweep, quantities for the
turbulence model: 11N variables.
• The matrix:
• The preconditioner:
• Right-hand side and solution vector of linear system:
• For GMRES, if restarted every 25th step:
• The transformation from Cartesian to contravariant velocity components:
Together with several fields of variables that are not so large, we arrive at the following estimate of the number of
variables to be stored:
(6)
In order to increase the Mflop-rate on one processor of the Cray C90, memory bank conflicts should be avoided.
Therefore, if the leading dimension of large arrays is a power of two, on the Cray C90 it is increased by one. This gives an
additional increase of the memory usage by approximately 15%.
A very robust strategy is to avoid the use of a marching scheme by choosing g NX. However, from numerical
experiments it has appeared that the number of global iteration steps does not increase markedly if g equals NX/2 instead
of NX, while the memory requirements can almost be halved.
3.4 Including ‘negligible' terms
The flow around more or less slender bodies like ships and aircrafts at high Reynolds number is characterised by a
dominant flow direction. Diffusion in that direction is practically negligible, and we have discarded the associated terms
in the momentum conservation equations in setting up PARNASSOS. Also the convection terms for the dominant flow
direction have always been treated somewhat differently from the convection terms in the lateral directions. Although
PARNASSOS has been dealing with flows with boundary layer separation in the sense that there is a region where the
flow is against the dominant flow direction, we have assumed that the reversed-flow region will be of relatively small
extent and therefore with low convective transport against the main flow. As a consequence we have neglected the
corresponding convection term.
While in a plane-by-plane marching scheme the ‘negligible' terms have to be included explicitly, i.e. in the right-
hand side vector, in a sub-domain of greater extension the neglected terms can readily be included implicitly.
Accordingly, we have investigated the role of these terms to confirm that they are negligible indeed and to verify their
effect on the convergence behaviour.
4 NUMERICAL RESULTS
The present comparisons are performed for the calculation of the flow around the HSVA-tanker at model scale
Reynolds number, Rn 5×106, and around the Mystery Tanker both at model scale Reynolds number, Rn 5×106, and at full
scale Reynolds number, Rn 2×109. Except for the HSVA-tanker, where Cebeci-Smith's algebraic turbulence model is used
(Cebeci & Smith 1984), the results shown in this section were generated using Menter's one-equation model (Menter
1997). The grid around the HSVA-tanker consists of 137 x 81 x 37 nodes in longitudinal, wall-normal and girthwise
direction, respectively. For the calculation of the flow around the Mystery tanker at model scale the grid consists of 161 x
81 x 41 nodes. For the full-scale computation an increased number of 121 grid nodes is used in wall-normal direction,
bringing the total number of nodes for that case to almost 1 million. Although the grids used to generate the results shown
in this section cover only
the authoritative version for attribution.

OCR for page 882

lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
COMBINING ACCURACY AND EFFCIENCY WITH ROBUSTNESS IN SHIP STERN FLOW COMPUTATION 887
the region around the afterbody of the ship, the current implementation of PARNASSOS allows to extend the domain to
include the bow part also.
The inlet and outlet plane are x=constant planes. The inlet plane is located at x 0.5L (midship) and the outlet plane at
x 1.25L. The external boundary is an elliptical cylinder, given by:
The remaining boundaries are the still water surface, plane z 0.056L, the symmetry plane of the ship, y 0, and the hull
surface.
On the hull surface we specify Dirichlet conditions for the three velocity components, as dictated by the no-slip and
impermeability conditions. On the symmetry plane of the ship, a Dirichlet condition for the velocity component normal to
the boundary is used, and Neumann condition for both tangential velocity components. In this section, free surface effects
are not taken into account, hence on the still water surface boundary we use similar symmetry conditions (double-body
computation). Section 5 presents some current developments about how to incorporate free surface effects. Both at the
outlet plane and the external boundary, an exact condition is not available for a computation domain of finite extent.
Hence the location of these boundaries determines for a great deal the quality of the boundary conditions. At the external
boundary we impose Dirichlet boundary conditions for the pressure and for the velocity components tangential to the
boundary, derived from an inviscid-flow computation. At the outlet plane, we take as the boundary condition that the
gradient of the pressure in longitudinal direction vanishes. At the inlet plane, we use Dirichlet boundary conditions for the
velocity components. If diffusion in main stream direction is not taken into account, this completes the set of boundary
conditions. If all diffusive terms are taken into account, we use as extra boundary condition at the outlet plane that the
normal derivatives of the velocities vanish.
If a one-equation turbulence model is used, suitable values for the eddy viscosity (derived with the aid of an
algebraic turbulence model) are imposed on the inlet boundary, while vt is set to zero on a no-slip boundary as well as on
the external boundary. At the other boundaries, we use similar boundary conditions in the turbulence models as for mass
and momentum equations.
An illustration of a possible effect of choosing large subdomains is shown in Fig. 1 and 2, which show the
convergence behavior for the computation of the flow around the HSVA tanker on model-scale Reynolds number. If g
NX/2 the whole separation region is contained in one sub-domain, instead of being distributed over more, and a
considerable improvement in the convergence behavior can be obtained. We like to point out that on the vertical axis the
maximum (L∞)-norm and not the more common L2-norm (which would have suggested a much smoother convergence
behavior) is measured. The CPU-time on one processor of a Cray C90 for the computation corresponding with Fig. 2 is
approximately 10 minutes.
Figure 1: Convergence behavior for HSVA-tanker, Cebeci-Smith's turbulence model. g 4. Rn 5×106.
Next, we study the effect of taking into account all convective terms on the convergence behavior of the global
iteration. Fig. 3 shows the convergence behaviour for the computation of the flow around the Mystery tanker on model-
scale Reynolds number that is obtained when in reversed-flow zones convection in main stream direction is dropped,
whereas Fig. 4 shows the convergence behavior obtained when all convective terms are taken into account. It appears that
taking all convective terms into account has a stabilising influence on the global convergence behavior. This is caused by
the second-order upwind scheme for the gradients in the streamwise convection, which gives a positive contribution to the
main diagonal of the coefficient matrix. For comparison, Fig. 5 shows the convergence behavior in case g 1, also for
model-scale Reynolds number. It appears that the number of global iteration steps slightly increases if g is increased from
1 to 4. However, the required CPU-time is reduced considerably since at each do
the authoritative version for attribution.

OCR for page 882

lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
COMBINING ACCURACY AND EFFCIENCY WITH ROBUSTNESS IN SHIP STERN FLOW COMPUTATION 888
Figure 4: Convergence behavior Mystery tanker,
Figure 2: Convergence behavior for HSVA-tanker,
Cebeci-Smith's turbulence model. g NX/2. Rn 5×106. Menter's turbulence model. No terms neglected. g 4. Rn
5×106.
Figure 3: Convergence behavior Mystery tanker,
Menter's turbulence model. Streamwise convection Figure 5: Convergence behavior Mystery tanker,
neglected in reversed-flow zone. g 4. Rn 5×106. Menter's turbulence model. No terms neglected. g 1. Rn
5×106.
the authoritative version for attribution.

OCR for page 882

lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
COMBINING ACCURACY AND EFFCIENCY WITH ROBUSTNESS IN SHIP STERN FLOW COMPUTATION 889
main only one step of the quasi-Newton method is applied per global iteration step. Fig. 6 shows the convergence
behavior obtained by increasing the size of the subdomains (g 16). It appears that, in contrast with the results shown in
Fig. 2, the speed of convergence is not improved by choosing the value of g larger than 4. This can be explained by the
fact that in Menter's model the convergence behavior can be deteriorated by strongly non-linear source terms, whereas in
Cebeci-Smith's turbulence model the eddy viscosity is completely determined by the velocity field.
The CPU-times on one processor of a Cray C90 for the computations corresponding with Fig. 4 and 6 are
approximately one hour. We note that in practice it is sufficient to reduce the maximum of the change in the pressure to
below 10−5. Hence in that case the amount of required CPU-time is significantly smaller.
Figure 6: Convergence behavior Mystery tanker, Menter's turbulence model. No terms neglected. g 16. Rn 5×106.
Next, we study the influence of the Reynolds number on the global convergence behavior. Comparing Fig. 7 with
Fig. 4 shows that the speed of convergence decreases some what, but the increase in the number of iterations is not very
dramatic. This is in accordance with previous experience that the convergence speed of PARNASSOS is essentially
independent of Rn. In this case, the required CPU-time is approximately hour, the increase being mainly due to the
increased number of grid nodes.
In order to study the effect of the terms that have been neglected so far in PARNASSOS, some integral quantities
and maximum, minimum and mean values of the unknowns are given for the Mystery tanker in Tables 1 and 2. The
different versions of PARNASSOS have been indicated as follows
Figure 7: Convergence behavior Mystery tanker, Menter's turbulence model. No terms neglected. g 4. Rn 2×109.
• PNS: streamwise diffusion and convection in streamwise flow separation regions, (U1<0), are neglected.
• PNS: Only neglects streamwise diffusion.
• PNS: No terms in the RANS-equations are neglected.
The quantities given in the tables are:
• Friction resistance coefficient,
• Pressure resistance coefficient,
the authoritative version for attribution.
• Wake fraction, Wf:

OCR for page 882

lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
COMBINING ACCURACY AND EFFCIENCY WITH ROBUSTNESS IN SHIP STERN FLOW COMPUTATION 890
The integrals in the definitions above are calculated with Gaussian quadrature rules assuming a bi-linear variation of
the unknowns between the grid nodes. The area Ω for the calculation of the wake fraction is the propeller disc which has
been defined by:
• x=0.989L
• Axis of the propeller at z=0.0166L
• Hub radius equal to zero.
• Propeller radius R=0.015L.
The maximum, (max), minimum, (min), and mean, (med), values are obtained for the following flow variables:
• U1, axial velocity.
• Vw, cross-stream velocity,
• axial velocity at the outlet plane.
• Cp, pressure coefficient,
• vt, eddy-viscosity scaled by
Table 1: Model Scale, Rn=5×106 Table 2: Full Scale, Rn=2×109
PNS+ PNS++ PNS+ PNS++
Version PNS Version PNS
1.605 1.605 1.606 0.791 0.791 0.791
0.747 0.747 0.741 0.537 0.537 0.535
0.653 0.653 0.652 0.657 0.657 0.657
Vwmed×10 0.631 0.634 0.632 Vwmed×10 0.633 0.633 0.632
Cpmed×102 −0.361 −0.364 −0.361 Cpmed×102 −0.304 −0.304 −0.303
(vt)med×104 0.262 0.262 0.262 (vt)med×104 0.109 0.109 0.109
0.624 0.621 0.619 0.353 0.353 0.353
Wf Wf
−0.121 −0.071 −0.068 0.000 0.000 0.000
1.064 1.064 1.064 1.073 1.073 1.073
0.687 0.687 0.686 0.804 0.804 0.803
0.310 0.311 0.311 0.410 0.410 0.411
Vwmax Vwmax
−0.111 −0.111 −0.111 −0.155 −0.155 −0.155
Cpmin Cpmin
0.149 0.149 0.148 0.203 0.203 0.203
Cpmax Cpmax
(vt)max×104 1.267 1.265 1.261 (vt)max×104 0.626 0.626 0.626
The results in these tables show that the effect of the “missing” terms is very small, especially for full scale Reynolds
number. However, in case of flow reversal, the missing convective term has a larger effect on the maximum reversed flow
velocity
Fig. 8 shows the limiting streamlines computed by PNS++ both for model scale and full scale. A ○ in the figure
indicates where flow separation was detected. The results of PNS and PNS+ are not shown in these
Figure 8: Limiting streamlines Mystery Tanker, above: Rn=5×106, below: Rn=2×109
the authoritative version for attribution.

OCR for page 882

About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
the authoritative version for attribution.
dotted: PNS·, dashed dotted: PNS· ·. Rn 2×109.
dotted: PNS·, dashed dotted: PNS· ·. Rn 5×106.
Figure 10: Isolines U1 for x 0.989L, solid line: PNS,
Figure 9: Isolines U1 for x 0.989L, solid line: PNS,
COMBINING ACCURACY AND EFFCIENCY WITH ROBUSTNESS IN SHIP STERN FLOW COMPUTATION
dotted: PNS·, dashed dotted: PNS· ·. Rn 2×109.
dotted: PNS·, dashed dotted: PNS· ·. Rn 5×106.
Figure 12: Isolines Cp for x 0.989L, solid line: PNS,
Figure 11: Isolines Cp for x 0.989L, solid line: PNS,
891

OCR for page 882

lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
COMBINING ACCURACY AND EFFCIENCY WITH ROBUSTNESS IN SHIP STERN FLOW COMPUTATION 892
Figure 15: Log10 of the differences between PNS and
Figure 13: Isolines Log10(vt) for x 0.989L, solid line:
PNS· · in the magnitude of the transverse velocity vector
PNS, dotted: PNS·, dashed dotted: PNS· ·. Rn 5×106.
for plane x 0.989L. Rn 5×106.
Figure 14: Isolines Log10(vt) for x 0.989L, solid line:
PNS, dotted: PNS·, dashed dotted: PNS· ·. Rn 2×109. Figure 16: Log10 of the differences between PNS and
PNS· · in the magnitude of the transverse velocity vector
for plane x 0.989L. Rn 2×109.
the authoritative version for attribution.

OCR for page 882

lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
COMBINING ACCURACY AND EFFCIENCY WITH ROBUSTNESS IN SHIP STERN FLOW COMPUTATION 893
figures, because these limiting streamlines are virtually identical. Only in the neighbourhood of the propeller-plane,
the computed velocities and pressure differ slightly. Therefore, Fig. 9–14 show the computed results in the propeller plane
x 0.989L. These results show that the effect of the neglected terms is not large, and for model scale Reynolds number, the
differences are located mainly in the region that contains flow separation. From these figures it appears that for full-scale
Reynolds number the influence of the tradionally neglected terms in PARNASSOS is even smaller than for model-scale
Reynolds number.
Despite these small differences, the velocity fields computed by PNS and PNS· · show hardly any difference, as is
clearly demonstrated by Fig. 15 and 16. They show the isolines of the Log10 of the differences in the magnitude of the
transverse velocity vector. Fig. 15 shows the results for model-scale Reynolds number in the propeller plane, whereas
Fig. 16 shows similar results for full-scale Reynolds number.
5 CURRENT DEVELOPMENTS: INCORPORATING THE FREE SURFACE
By basing the modelling and the solution approach on the physics of the problem considered, an accurate, efficient
and robust method has thus been obtained. A similar philosophy underlies our current research on incorporation of the
free surface.
Viscous effects on the wave pattern of a ship are relatively limited, and are essentially confined to the stern area. If
we leave the stern waves and wake out of account for a moment, experience from validations shows that the nonlinear
free-surface potential flow model can produce accurate predictions of the wave pattern. A RANS-solution subject to the
kinematic and dynamic free-surface boundary conditions (FSBC) should reproduce the same wave pattern everywhere
except in the stern area. Conversely, if we first use a potential flow code to compute the wave pattern, and then solve for
the viscous flow under that given wave surface (imposing just the kinematic FSBC), we should get the same solution, i.e.
the dynamic condition should be automatically satisfied as well. Fig. 17 illustrates the result of this composite approach
(Windt & Raven 2000) based on the codes RAPID (Raven 1996) and PARNASSOS, for the Series 60 model at Fn 0.316.
Rn 3.4×106. The bottom half shows isolines for the wave elevation around the afterbody, the top half shows the
corresponding isolines for the error in the dynamic FSBC, expressed as a wave elevation difference. The error is
significant only in the stern area and is negligible everywhere else. Therefore, for most of the domain we have already
found a solution equivalent to a full RANS/FS solution, with an effort just marginally larger than that for a double-body
RANS solution. This easy and most efficient procedure will soon replace the double-body approach for most of our
practical computations.
Remains, however, the stern area where viscous effects on the wave pattern are present, as shown by the larger
pressure differences at the wave surface. Locally, a full RANS/FS method is still needed. While most such methods use a
time-dependent solution approach, unsteady waves may delay the approach to steady state considerably. We prefer to
solve the steady form of the free-surface problem, by iteration rather than by time-stepping. Refs Raven & van
Brummelen 1999 and van Brummelen & Raven 2000 derive a new form of the FSBC suitable for iterative solution, and
show the superior performance for a 2D test case. The iterative process converges very fast, and the total effort is just 2 or
3 times that for corresponding problem without free surface. The same approach is now being implemented in
PARNASSOS to be tested in 3D, and will hopefully supplement the composite procedure just described.
6 CONCLUSIONS
The results for the Mystery tanker have shown that the effect of the terms that have been neglected so far in
PARNASSOS is indeed very small. Especially at full scale Reynolds number, it is justified to omit diffusion in
mainstream direction altogether and to discard part of the convection terms in regions with flow separation. However,
especially the latter terms can have a positive influence on the convergence behaviour. Therefore, it is better to take all
convective terms into account. In the current production version of PARNASSOS, we also take the streamwise diffusion
into account. Although we have demonstrated that the effect of this term on both the convergence behavior and the
computed results is very small for the Mystery tanker, taking all terms into account removes any doubt in this respect.
The new feature in PARNASSOS, a space-marching solution with sub-domains of larger size than a single grid
plane has improved the efficiency and the flexibility: it is now possible to choose the size of the sub-domains in the range
from one-grid-plane to even the complete domain. The plane-by-plane solution strategy has the advantage that it is very
efficient with respect to memory usage. A full-scale computation using 1 million grid points can be done on a PC with
only 64 Mbyte of main memory.
the authoritative version for attribution.

OCR for page 882

lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
COMBINING ACCURACY AND EFFCIENCY WITH ROBUSTNESS IN SHIP STERN FLOW COMPUTATION 894
Figure 17: Composite solution for viscous flow under wave surface. Bottom half: wave elevation around stern (inviscid). Top
half: error in dynamic FSBC (viscous), expressed as wave elevation difference. Contour interval ∆z 0.0005L.
The multi-plane strategy, on the other hand, has the advantage that it is faster, mainly because the quasi-Newton
linearization process can be shifted to the global iteration. In practice, the best strategy can be chosen dependent on the
machine available for the computations. Treating 4 planes simultaneously seems a good choice in many practical
situations. In that case, the memory requirements increase only by approximately a factor of 4 compared to the plane-by-
plane solution strategy. Even a flow computation at full-scale Reynolds number can then still be done on a workstation or
a PC that contains at least 256 Mbytes of memory.
If enough computer memory is available and if an algebraic turbulence model like the Cebeci-Smith model is used,
increasing the size of the sub-domains is a good option since it reduces the number of global iteration steps. If Menter's
turbulence model is used, the number of global iteration steps does not decrease significantly by increasing g. In the near
future, we hope to remedy this phenomenon by treating the nonlinear source terms differently.
More than 80% of the CPU-time is spent on the solution of the systems of linear equations, using preconditioned
GMRES. The Mflop-rate that is obtained on a vector computer like the Cray C90 is quite high: on one processor of this
machine, approximately 300 Mflop/sec is attained. This is achieved by choosing a suitable ordering in the computation
which circumvents recurrences in the construction and application of the preconditioner.
We are confident that the enhancements incorporated in PARNASSOS as described in this paper have brought us
closer to the optimum combination of flexibility, accuracy, efficiency and robustness in ship stern flow computations.
ACKNOWLEDGEMENTS
The authors wish to thank H.Raven for many stimulating discussions and numerous suggestions for improving the
presentation of the paper. They also thank the university of Utrecht for the fortran routine that implements GMRES(M).
REFERENCES
[1] T.Cebeci and A.M.O.Smith, Analysis of Turbulent Boundary Layers. Academic Press, November 1984.
[2] L.Eça and M.Hoekstra, On the numerical verification of ship stern flow calculations. In paper presented at the MARNET workshop, Barcelona,
Spain, 1999.
[3] L.Eça and M.Hoekstra, Numerical Prediction of Scale Effects in Ship Stern Flows with Eddy-Viscosity Turbulence Models. In Twenty-Third
Symposium on Naval Hydrodynamics Applications to Ship Flow and Hull Form Design, Osaka, Japan, September 2000.
[4] M.Hoekstra and L.Eça, PARNASSOS: An efficient method for ship stern flow calculation. In
the authoritative version for attribution.

OCR for page 882

lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
COMBINING ACCURACY AND EFFCIENCY WITH ROBUSTNESS IN SHIP STERN FLOW COMPUTATION 895
Third Osaka Colloquium on Advanced CFD Applications to Ship Flow and Hull Form Design, Osaka, Japan, 1998.
[5] M.Hoekstra and L.Eça, An example of error quantification of ship-related CFD results. In Seventh Conference on Numerical Ship Hydrodynamics,
Nantes, France, 1999.
[6] M.Hoekstra, Numerical Simulation of Ship Stern Flows with a Space-Marching Navier-Stokes Method. PhD thesis, MARIN, Wageningen, 1999.
[7] F.R.Menter, Eddy Viscosity Transport Equations and Their Relation to the k-ε Model. Journal of Fluids Engineering, 119:876–884, December 1997.
[8] H.C.Raven, A Solution Method for the Nonlinear Ship Wave Resistance Problem. PhD thesis, MARIN, Wageningen, 1996.
[9] H.C.Raven and H.van Brummelen, A new approach to computing steady free-surface viscous flow problems. In 1st MARNET-CFD workshop,
Barcelona. To be downloaded from http://www.marin.nl/publications/pg _resistance.html, 1999.
[10] Y.Saad and M.H.Schultz, A generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM J. Sci. Statist. Comput.,
7:856–869, 1986.
[11] H.van Brummelen and H.C.Raven, Numerical solution of steady free-surface Navier-Stokes flow. In 15th Int. Workshop on Water Waves and
Floating Bodies, Caesarea, Israel. To be downloaded from http://www.marin.nl/publications/pg_resistance.html, 2000.
[12] H.A.van der Vorst, Bi-CGSTAB: A fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems. SIAM J. Sci.
Statist. Comput., 13(2):631–644, 1992.
[13] J.Windt and H.C.Raven, A composite procedure for ship viscous flow with free surface. 2000.
the authoritative version for attribution.

OCR for page 882

lengths, word breaks, heading styles, and other typesetting-specific formatting, however, cannot be retained, and some typographic errors may have been accidentally inserted. Please use the print version of this publication as
About this PDF file: This new digital representation of the original work has been recomposed from XML files created from the original paper book, not from the original typesetting files. Page breaks are true to the original; line
COMBINING ACCURACY AND EFFCIENCY WITH ROBUSTNESS IN SHIP STERN FLOW COMPUTATION 896
DISCUSSION
A.Sánchez-Caja
VTT Manufacturing Technology, Finland
In your very interesting paper you mention that the outlet plane for your computations is located at x=1.25L, midship
being located at x=0.5L. This means that the position of the outlet plane relative to the aft-perpendicular is 0.25L, which
seems to be too small. Have you checked the influence of the outlet plane location on the results, for example on the value
of the drag coefficient?
AUTHOR'S REPLY
The choice of the location of the outlet plane is based on our experience and is a compromise between maximizing
grid density and minimizing the influence of approximate boundary conditions. The inexactness of the boundary
conditions on the outlet plane is in the first place related to the condition imposed for the pressure. Where in most
methods a Dirichlet condition (undisturbed pressure) is applied, we use a less restrictive Neumann condition that allows
the selection of a smallest domain size, without significant loss of accuracy.
As an illustration, we give some results for the Mystery tanker at model scale Reynolds number. The table below
summarizes a few meaningful results for three positions of the outlet plane. It appears that the computed data hardly
changes when the position of the outlet plane relative to the aft-perpendicular is moved from 0.25L to 0.365L.
Furthermore, the plots of the isolines of the pressure and velocity distributions obtained with the outlet plane at x/L=1.365
are virtually identical to the plots shown in Section 3.4, where the position of the outlet plane is at x/L=1.25. On the other
hand, by moving the outlet boundary to, say, x/L=2.0 and keeping the same number of grid nodes we would observe a
much more drastic effect on the solution. So positioning the outlet plane at x/L=1.25 is a deliberate choice.
Table 3: Computed values for several locations of the outlet plane. Rn=5x106.
1.135 1.250 1.365
x/L
1.605 1.606 1.606
0.688 0.741 0.744
0.624 0.619 0.619
Wf
U1min −0.069 −0.068 −0.068
U1max 1.064 1.064 1.064
0.310 0.311 0.311
Vwmax
−0.110 −0.111 −0.111
Cpmin
0.148 0.148 0.148
Cpmax
the authoritative version for attribution.