Below is the uncorrected machine-read text of this chapter, intended to provide our own search engines and external engines with highly rich, chapter-representative searchable text of each book. Because it is UNCORRECTED material, please consider the following text as a useful but insufficient proxy for the authoritative book pages.

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.

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.

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

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,

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.

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.

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.

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:

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.

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

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.

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.

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.

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.

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.