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.
Recent Developments in a Ship Stern Flow Prediction Code M. Hoekstra Maritime Research Institute Netherlands Wageningen, The Netherlands Abstract This paper is concerned with developments in a numerical method for the prediction of the steady incompressible flow around the stern of a ship. The method is based on a finite-difference discretisation of the Reynolds-averaged Navier-Stokes equations on a boundary-fitted curvilinear grid. The attention will be focused on three aspects, viz. the choice of the velocity variables, the boundary conditions and the global (outer) iteration process for the pressure field. Various illustrations will be given of the improved performance of the method. 1. Introduction It is in the first place owing to computer hard- ware developments that a computational simulation of the ship stern flow field has become feasible. In- deed, the availability of adequate computer facilities may be considered as a prerequisite for an accurate calculation to be carried out. It is not a sufficient condition, however. For some years already, the nu- merical simulation of the steady viscous flow around a ship's stern - usually under neglection of free sur- face effects - has been a prominent research topic. Yet, something like a commands opinio on how a re- lial~le, efficient and robust code is to be constructed has not been attained. It indicates that developing a suitable algorithm for the solution of the Reynolds- averaged Navier-Stokes equations for incompressible flows around complex body shapes is not plain sail- ing. At the author's institute, work in the field of stern flow computation was initiated around 1980 and has continued steadily since. It has resulted in a compu- tational procedure that has been applied with some success to a number of ship forms and axisymmetric bodies. The main features of the method have been outlined previously ~1,2,3] and a brief recollection of these features (Section 2) will suffice here as a frame- work for the subsequent discussion of some aspects of the method. Where the title of this paper promises the presen- tation of recent developments, 'recent' is to be inter- preted in relation to the time of appearance of the above references. The length limitations to which this paper is subjected do not allow a presentation of all changes that have taken place in the code since 1985. Therefore three major items have been se- lected for discussion here. They are: the choice of the velocity variables in relation to regular solution behaviour near grid singularities in the wake (Sec- tion 3~; boundary conditions and their numerical im- plementation (Section 4) and a new global pressure relaxation scheme (Section 5~. Moreover, results of a representative application will be shown as an il- lustration of the progress that has been made (Sec- tion 6). 2. Outline of computational procedure We are interested in a numerical simulation of the behaviour of the flow around the stern of a ship, be- ing towed steadily through still water. The free sur- face disturbances created by the ship are assumed to be of little significance so that it is allowed to replace the undisturbed free surface by a symmetry plane. Thus, taking a reference frame moving with the ship, we consider actually the double model of the submerged part of a ship, held fixed in a uniform flow directed from bow to stern. The flow domain is now assumed to be unbounded and the fluid to be incompressible. Three subdomains are distinguished as shown in Fig. 1, viz. the external flow zone where the flow behaves as being inviscid and irrotational; the bc)undary-layer zone, covering the forward part of the hull surface, where first-order boundary-layer theory is supposed to provide an adequate description of the flow; the stern-flow-and-wake zone which is of 87
primary concern in this paper. In the latter zone we solve eRectively the Reynolds-averaged Navier-Stokes equations but we shall assume that a predominant flow direction exists. Among other things, it allows us to drop all diffusion terms associated with that pri- mary flow direction. The equations are supplemented by an algebraic turbulence model. The boundary con- ditions will be described in Section 4. III /////HULL /~ Fig. 1. Division of flow field into three zones The flow-governing equations are written in terms of general curvilinear coordinates. The contravariant formulation of the momentum equations is adopted, which gives three relations expressing momentum con- servation along each of the three coordinate lines. The velocity components and the pressure are the dependent variables but a particular choice for the velocity variables is made as described in Section 3. The continuity equation is solved in its original form and is not replaced by a Poisson equation for the pressure as in Marker-and-Cell type methods. The equations are discretised on a boundary-fitted curvilinear grid. It is a single-block regular grid of N V * NY * NZ nodes. The grid is partly non-ortho- gonal but conceptually simple. All variables are de- fined on the grid nodes, in other words we do not apply grid staggering of one kind or another. That is convenient in many respects but it poses certain de- mands on the discretisation of the pressure gradients and the continuity equation which we are careful to fulfill to avoid 'checkerboarding' A. The discretisa- tion is second or third order accurate with the excep- tion of the main stream pressure gradient which is of first order accuracy. All derivatives with respect to the coordinate associated with the main stream di- rection are modelled by upstream difference formulae except the pressure gradient. Newton linearization is applied to the convective terms. The solution is obtained iteratively where the it- eration sequence may be characterized as a multiple- sweep space-marching process which takes maximum benefit of the existence of a predominant flow direc- tion. In fact, two iteration cycli can be distinguished: the local (or inner) and the global (or ousters iteration prc)cess. The local process refers to the solution of a sub- set of difference equations applying to all grid nodes having the same value of the main-stream coordinate. It is based on the Coupled Strongly Implicit Proce- dure (CSIP) [5], an incomplete factorization scheme, and yields a simultaneous solution of all variables in a cross-section of the computation domain. Iteration is required both by the non-linearity of the differential equations and by the incompleteness of the factoriza- tion. The CSIP has been selected because we think that a high degree of implicitness is desirable if not necessary for an efficient solution algorithm and be- cause it retains the coupling between the momentum and the continuity equations in the numerical solu- tion. The global iteration process involves the evalua- tion of the solution in repeated sweeps from the up- stream to the downstream boundary of the compu- tation domain. This process allows the influence of downstream occurrences to be felt by the upstream flow, an inherent property of the mathematical model used. To improve its convergence rate, each down- stream sweep is followed by an update of the solution of the pressure field in a reversed sweep as will be explained in Section 5. A further enhancement of the computational efficiency is obtained by grid se- quencing, with which we mean that the grid used is initially very coarse in the mainstream direction and is subsequently refined in two or three stages. 3. Grid singularities Fig. 2 gives a sketch of the computation domain as it looks after the symmetry properties of the flow have been taken into account. In addition to the Carte- sian reference frame x,y,z a curvilinear boundary- fitted coordinate system 6, a, ~ is introduced. A trans- formed space can be imagined in which (, 7~, ~ forms a rectangular coordinate system and in which the hull surface is plane. The appearance of the (,7~,( sys- tem in physical space depends on the transformation relations between (,71,( and ~,y,z. PHYSICAL DOMAIN z ,/ REV (a) TRANSFORMED DOMA I N TV ~577 REV Fig. 2. (computation domain 88 (b)
To simplify the grid generation procedure we have chosen ~ = x so that transverse sections in the phys- ica1 space remain so in the transformed space. An orthogonal mesh is then constructed in selected trans- verse sections by a procedure which relies on confor- mal mapping although the final mesh is not conformal due to subsequent stretching t64. By connecting these grids, the 3-D mesh is completed. It leads inevitably to non-orthogonality of the if, 71 and if, ~ lines in a part of the computation domain. However, before this task can be carried out it must be decided where the four corners of the trans- verse section in the transformed space are to appear in the physical space. A glance at Fig. 2 will help to recognize that the choice is obvious in the upstream part of the domain but less so behind the stern. Re- ferring to Fig. 3, the location of point A correspond- ing with A' can be taken somewhere between B and D or even coinciding with either of these points. In practice the final choice is primarily governed by the requirement that the (-lines must roughly be aligned with the main flow. The points A are therefore usu- ally located on a straight line extending from the keel (cf. Fig. 2~. z oriented physical velocity variables are used, so that v and w are the 71 and ~ component of the physical velocity. Then we must have v=O, w=finite along AD (Fig. 3), while w=O, v=finite along AB. It leads to a conflict at A because v and w cannot both be zero in A. Although this problem can be circumvented by e.g. grid staggering, the accuracy of finite-difference approximations is likely to be low in the vicinity of the singularity due to the sudden changes of some of the variables. We alleviate the difficulty by a special choice of the velocity variables, which may be explained as fol- lows. The coordinate transformation near the grid singularity in the wake tends to y + itz&) = -in(Ry + it (1) where i = x/=T, R = y/~ is the grid cell aspect ratio, c' is a proportionality constant while the singu- larity is located a+ the point y = 0,z = z, . Because of the relation (1) the singularity is called parabolic. Formula (1) is a well-known type of transformation in potential flow theory. It can be used to transform a parallel flow to a corner flow. If the parallel flow is supposed to occur in the physical space and is de- scribed by a velocity potential ~ = VT (` ZZ3, Y. (` Via is constant ) , D - ~ SAC D' C we find with the aid of (1) / ~ = ~Ta(~2 _ R202y AL ~ A L I B V2; 0~ -2aR: Van; Fig. 3. Cross-section of physical (left ~ and trans- formed (right) computation domain behind hull The coordinate transformation is singular in such points: the Jacobian of the transformation vanishes. Singularities may also appear elsewhere on the bound- aries of the domain - e.g. along the waterline when the frame line of the hull does not meet the symmetry plane (undisturbed free surface) at a right angle - but the mapping technique used excludes their occurrence in the interior of the domain. Grid singularities may deteriorate the accuracy of the calculations or destabilize the iterative solution process ~7~. This is particularly true when there is a finite velocity in the singularity as is the case be- hind the hull. Suppose for example that coordinate- We see that the velocity components V2,V3 be- have nicely in the transformed space. They vanish at the corner and increase linearly away from it so that accurate discretisation is possi ble. It is important now to recognize that V2,V3 are covariant quantities since ~ is a scalar. Apparently, covariant velocity components behave regularly near a parabolic singularity in contrast to physical or con- travariant components as can readily be verified. For example the contravariant velocity component v2 is in an orthogonal grid (as we have behind the hull) given by V2 = 922 V2V2 /g22 = - 2OR1/Tr1/ and tends to infinity for 71 ' O ~9
Formerly, we employed the contravariant veloc- ity components directly as dependent variables in the mathematical formulation. Denoting them by Vi (i = 1,2,3), we define now the three velocity variables as U = Vi; V = 922V2; W = ~V3 instead. Notice that v and w are not strictly covariant velocity components. But the multiplication factors for both v2 and V3 involve the Jacobian ~ (922 = R4), which is sufficient to obtain regular behaviour near the singularity. If the y, z velocity components in the physical space must be determined from these variables, dif- ficulties are encountered at the singularity. There is no trouble, however, in evaluating them in neighbour- ing grid points whereupon the required values at the singular point can be derived by interpolation. 4. Boundary conditions The stern-flow-and-wake zone has eight bound- aries (Fig. 2~: the hull surface, an inlet and an outlet boundary and the external boundary, the remaining boundary surfaces being located in symmetry planes. The terminology is intentional: the velocity compo- nent u is assumed to be positive on both the in and outlet boundaries. On the other hand the flow may enter or leave through the external boundary. Like the Navier-Stokes equations themselves, the Reynolds-averaged Navier-Stokes equations for incom- pressible flow require three boundary conditions on all boundaries [e.g. 84. However, if diffusion along the coordinate associated with the predominant flow di- rection is neglected, the equations exhibit Euler char- acter in that direction with a consequent change in the boundary condition requirements. Three condi- tions are still required at the inflow boundary, but only one condition must be imposed at the outflow boundary t94. It is obvious what conditions should be applied on the natural boundaries of the computation domain, viz. the hull surface and the symmetry planes. On the hull surface the no-slip and the impermeability condition give us three Dirichlet conditions for the velocity components. On a symmetry plane there is one Dirichlet condition for the normal velocity com- p<:,nent and a Neumann condition for each of the two other components. It is perhaps worth mentioning here that we im- pose the hull boundary conditions directly by calcu- lating the flow down to the hull and not indirectly via the use of wall functions as an approximate descrip- tion of the near-wall flow behaviour. The remaining boundaries are artificial bound- aries and it is not immediately clear what are suitable conditions. On the outlet boundary only one condi- tion must be prescribed, either for the pressure or for the normal velocity (i.e. u) t94. In view of our dis- cretisation it is natural to choose a condition for the pressure. A Neumann type condition is to be pre- ferred because it is less restrictive and allows some pressure variation over the outlet plane which may be important in view of the frequent occurrence of longitudinal vortices in a ship's wake; the existence of such vortices demand the pressure to be somewhat lower in the vortex core than in its surroundings. Mostly we use a vanishing longitudinal pressure gradient at the outlet plane, but sometimes a non- zero value, derived from potential flow calculations, is preferred. A remark should be made here about the numer- ical implementation of the pressure boundary condi- tion. In a space-marching scheme the discretisation of the streamwise pressure gradient involves a pressure value of the current as well as one of the preceding global iteration, in our discretisation: (` Pi + ~ pi ~ / /\ x i where rl counts the global iterations and i is the -station index. When the solution on the outlet plane at station any is to be determined, a value for pnN-~-+1 must be derived from the boundary condition, Pc = pgrad say. One might feel tempted to simply replace the pressure gradient as discretised above by the given boundary condition on the outlet plane: (\PN,\ +1 PNh )//\ XN,Y = pgrad It leads to unquiet solution behaviour near the out- let plane which eventually may deteriorate the global convergence. The correct procedure is to apply the Neumann condition for the pressure for each iteration level separately. Thus (\PN ~ +]PN V )//\X~,\- = pgrad for all n which implies that the pressure gradient term in the (-momentum equation must be 1nodelled on the out- let boundary as: pgrad -F (PN V PN V )//\~N V i.e. t he discretisatior, must include the pressure change wit l1 resl)( ct t`-) the prev;(:)us sweep. The three conditions to be specified on the inlet boundary must, again in view of the discretisation, apply to) the velocity components. Dirichlet condi- 90
tions may be chosen but since inlet conditions are obtained either from usually incomplete experimental data or from thin-boundary-layer calculations carried out in an approximate grid geometry, they are often imperfect, leading to a non-smooth solution near the inlet plane. This is the reason why we mostly replace the Dirichlet condition for v by a Neumann condition. Let us next turn to the boundary conditions on the external boundary. The old practice was to spec- ify the pressure and the two tangential velocity com- ponents. They were obtained from a potential flow calculation for the bare hull, neglecting the displace- ~ent effect of the boundary layer on the external flow. Such conditions are good enough if the external boundary is chosen far (several boundary-layer thick- nesses) from the hull. But that is an unattractive op- tion if one aims at computational economy because the mathematical model is then unnecessarily com- plicated in a considerable part of the domain. On the other hand, the conditions are evidently not ex- act when the extent of the computation domain is reduced. It may result in the formation of a weak non-physical boundary layer on the external bound- ary t94. After a close inspection we have indeed found ample evidence of its occurrence in our calculations in which we have always tried to use relatively small domains. A remedy is of course to correct the boundary con- ditions for viscous-inviscid interaction in the course of the solution process. Stern et al. t10] have com- pared a procedure of that type with a large domain approach and found it to be effective for two rela- tively simple test cases. Surprisingly, they applied the displacement body concept to update the exter- nal flow; surprisingly because of the ambiguity in the definition of the displacement thickness and the ap- proxima.te nature of the displacement body concept in viscous-inviscid interaction studies. A more accu- rate and unarnbiguc~us procedure is to use the normal velocity on the external boundary resulting from the viscous flow calculations as a boundary condition for repeated potential flow calculations for the exterior of hull plus stern-flow-and-wake zone. A practical implementation might be as follows. The potential flow around the bare hull may be obtained by, say, a. discrete source distribution on the hull surface; the source strengths are found by im- posing the impermeability condition and they are de- noted here by a,. Their induction yields as it were ,, ~ a first guess of the flow, being inviscid and irrota.- tional in all three zones of Fig. 1. Moreover, it gives a certain velocity fin the t'`:,unda.rics between zones I anal It and Ones I and Ill, respectively, Oldish inlay be decomposed into a normal (V,~,po~) and a tangen- tia.l (`V',pOt,`) component. Keep in mind now that the part of this basis flow covered by rep;:ion 1 of flog. 1 ~ ~ · T ~ I- ~ ~' V millet I,e reproduced by a source distribution a2 on the boundary of the union of hull, zone lI and zone III, if the normal. velocity (V,,,po`) would be used as a. bounty c`.'nditi`:,n. The viscous flow calculations in zones II and III, with (v`,,pot,) as a boundary condition, will yield a dif- ferent normal velocity (vnvi,C) on the boundaries be- tween the zones because of boundary-layer displace- ment effects. The related adjustment of the external flow might be calculated by using (vnvi`,C) as a new boundary condition to yield a new source distribution CJ3 defined on the same boundary as a2. However, as will become clear in a. moment, it is better to oper- ate with /\(T2 = CT3cr2, the correction of ~2 due to viscous-inviscid interaction. So we suggest to calcu- late /\CJ2' associated with a correction potential Icon' under the condition I31)cor on = Vn,visC Vn,pot, to evaluate the tangential velocity induced by them and to add it to (`v~,pO','). The pressure follows from Bernoulli's law which completes the set of new bound- ary conditions for the next viscous flow calculation. The process can be repeated until the solutions in the three zones match well enough. The advantage of operating with /\<J2 appears as soon as further approximations are introduced. In first-order boundary-layer theory the displacement ef- fects caused by boundary-layer formation in zone II hats a. negligible effect on the external flow. One may therefore simplify the procedure by assuming that i\~2 is non-zero on the external boundary of zone III only. The number of source panels required in a practical calculation is thereby drastically reduced. Although we think the above method to be practi- cable, we propose here an alternative that may lead to a. more elegant procedure because it does not involve new potential flow calculations. Instead of prescrib- ing fixed values for pressure and tangential velocity components, we set the two tangential vorticity com- ponents to zero and we let the pressure satisfy the Bernoulli equation on the external boundary in the viscous How computations. It implies that two of the three Dirichlet conditions are replaced by conditions of Neumann type. These conditions are exact on any external boundary, provided it is chosen in the effec- tively inviscid part (:,f the flow. Viscous-inviscid in- tera.ction should automatically be taken care of with these boundary conditions. As far as the author is 91
aware, the use of vorticity boundary conditions in primitive vartable formulations has nest been applied earlier although it has been suggested by Roe ills. One might wonder whether such conditions guar- antee a unique solution. After all, the conditions are valid for any external potential flow while no informa- tion on the velocity at infinity is conveyed by them. How is the correct solution selected from the many possible ones? Although a rigorous analysis is difficult and has not been completed yet, we think that the new bound- ary conditions give the solution appropriate for un- bounded flow conditions. It is clear that an out- side influence, disturbing the external flow, such as a nearby wall, would not be felt, except via the con- ditions specified on the inlet boundary. But in the absence of outside disturbances, the establishment of a certain external flow is fully governed by what hap- pens inside the stern-flow-and-wake zone (displace- ment effects of boundary layer and wake); there is no need for a communication of particulars of the exter- nal flow to the stern-flow-and-wake zone other than that it is an irrotational flow. The conditions on the inlet boundary make sure that we obtain the flow at the correct Reynolds number. At this moment, we have a limited experience with the application of vorticity boundary conditions. The test results seem to confirm that the correct so- lution is obtained but unfortunately they have all shown so far a reduced convergence rate. After a while the global convergence is fully governed by the changes on the external boundary which decrease very slowly. An example is given in Fig. 4, showing global convergence results for the Wigley hull. As a mea- sure for the rate of convergence we use the rate of decay of the maximum pressure change between suc- cessive iterations anywhere in the computation do- main. Besides two lines representing the results for the old and the new boundary conditions, respec- tively, a dashed line is shown indicating the maxi- mum pressure change on the external boundary. It is obvious from the figure that a grid refinement has been applied between iteration numbers 3 & 4 and 6 & 7; the new boundary conditions were imposed on the finest grid only. We may note that the slower convergence is not the result of a slow drift of the pressure level; the pressure adjustments continue to be partly positive and partly negative. I TH 01 R ~ CHLET BOUNDARY CONO ~ T I ENS W I TH VORT I C I TY BOUNDARY COND I T I ONS - xN ~ , Q ~ _ - C3 ~r~ O ,_ - - . , _ u'_ 1 0 4 8 12 ITERATION NO. \1\ y \N \\\ \~ it. 1 1 1 16 20 Fig. 4. Influence of boundary conditions on the global convergence for the flow around the Wigley hull 1 0.00 0.40 0.80 1.20 2X/L l l 1.60 2.00 Fig. 5. Calculated pressure change due to viscous- inviscid interaction 92
To illustrate that viscous-inviscid interaction is automatically taken into account with the new bound- ary conditions, we show a result of the calculated change of the pressure on the external boundary with respect to the data obtained from the potential flow calculation around the bare hull. In Fig. 5 the pres- sure changes on the intersection of the external bound- ary and the horizontal symmetry plane (undisturbed free surface) for the flow around the Wigley hull are presented. The figure shows the typical effect of vis- cous-inviscid interaction, viz. a pressure drop in a region in the vicinity of the stern bordered on both the upstream and the downstream side by a region of pressure rise. The slight increase of the pressure change towards the outlet station has been verified in an axisymmetric flow case to be a consequence of the inexactness of the boundary condition pa = 0; it disappears when a slightly negative pressure gradient is imposed as a boundary condition. We may conclude that the application of vortic- ity boundary conditions is worth to be investigated further. The only deficiency that we have observed is the reduced convergence rate caused by it. Attempts to remove this deficiency have not yet been made, however. As an aside we note that it is questionable whether with the alternative technique of repeated potential flow calculations the matching accuracy can ever be as good as /`Cp = 10-3 - achieved in Fig. 4 after 17 iterations - in view of the use of two totally different discretisation methods for, respectively, the viscous and inviscid part of the flow. 5. Global relaxation After discretisation and linearization, the differ- ence equations for all grid nodes can be put together in a matrix/vector equation Ad= b where A is a sparse square matrix, ~ is the vector of unknowns and b is a known vector (but containing previous iterates of the velocity components due to the non-linear convection terms). Let the entries of A be grouped in blocks so that all entries associated with the unknown variables in a certain transverse plane form a block. If such a block is represented by one element of a new matrix B and a corresponding grouping is carried through for the vectors ~ and b, we get BE = c (2) where B is a N,X' * N,X' square matrix with (in our discretisation) four non-zero diagonals and ~6 and c are vectors of size N ~ The it.er~.t.i~r'> ;:r~lilt.imn of thin system has been called in Section 2 the global or outer iteration process, to be distinguished from the local or inner iteration process needed to solve the block of equations associated with a certain transverse sec- tion. What will be called in this paper the 'basic method' to solve (2) is a Gauss-Seidel type iteration process: B is split into a sum of a lower (L) and an upper (U) triangular part where the elements of L are identical with the corresponding entries of B while U contains the upper diagonal of B only; the solution is updated via L~bn = cUt,bn-~ To enhance the convergence rate of this process we have used earlier a source term in the (-momentum equation, as originally proposed by Israeli and Lin t12] and described for our application in t13. It has in many cases proved to be a useful artifice, but our experience tells that it may have undesirable side- effects, particularly when used in combination with an algebraic turbulence model. Notably in the ini- tial phase of the calculations the source term may assume appreciable values, generating a significant overshoot in the velocity profiles or provoking flow separation. In the first case the turbulence model may show pathological behaviour in the determina- tion of the outer length scale, which deteriorates the global convergence. In the second case the required change of the discretisation formulae will have an in- fluence to the same effect. We have now abandoned the source term scheme in favour of an alternative convergence accelerator with a superior performance. It is described below. When an iterative process is used to solve a time- independent problem, it is often preferable to oper- ate with a transient form of the mathematical model and to try to find the steady state solution t134. Fol- lowing this suggestion, we add a quasi-time deriva- tive of the pressure to the (-momentum equation. If for the sake of compact notation we consider the equations valid in a Cartesian coordinate system, the momentum equation for the dominant flow direction becomes L(u, v, w) + PI = IT (3) where L is a differential operator incorporating the convection and diffusion terms but no 'time' deriva- tives. Ltu,v,wJ + pm = pt In the above Gauss-Seide] process the discretised -rno~nentu~n equation is ~^ w.~^v (Pin+ Pin-~)//~\Xi -- (pinpin-~)/~tRHSi (4) 93
where Lo is a. finite-difference analogue of L, RHSi is a. known right-hand side which does not contain variables at iteration level n and /\t is taken equal tc, i\x. Next, one may conceive a complementary discretisation of (3~: LA (`U, V, ?JJ~)n + (Pi +~ Pi + )//\xi(`pn+ipa)//\ t = RHSi (5) such that half of the sum of (4) and (5) yields a dis- cretisation of all terms of (3) at time (or iteration) level n, provided that pn = I tpn+~1 + pn~ ~ `6' is considered as a. new value for pa. From the combination of the formulae (4~5~6) above, a predictor-corrector scheme for the pressure is readily constructed. The predictor step is the basic method, the Gauss-Seidel scheme: L z~ (u, v, w in + ~ pin+ l1pi ~ / /\ hi = RH Si , giving new values u=,vn~wn for the velocity compo- nents and a first approximation p* to a new value of the pressure. The difference of (4) and (5) yields a simple algebraic relation to determine a fictitious pn+i The improved guess for pn is then obtained from the mean of this pn+i and the old pressure pun-. The latter operations may be combined in the correc- tor step pn = 1 tpn~ + path + I pan _ pa1' ~7y Because of the appearance of pa+ in the right-hand side, pa must be evaluated in an upstream marching sequence. Thus the basic process of a sweep from inlet to outlet station is followed by reversed sweep in which. an improved value for the pressure is obtained via the simple algebraic relation (7~. As a start -for this reversed sweep we have n _ * PN V ~~ PN,\- a.s in~media.lely follows from the above formulae when the mainstream pressure gradient is zero on the outlet boundary. It is easily verified by repeated substitution of (7) that a. pressure change at station xi can be expressed as where /\ pa = [of ( 2 ) ~f\Pk + ( 2 ) /\PNX~ pro pn _ pn--1 gyp* p* _ pa1 It implies that /\ pn is determined by a fraction (diminishing with distance) of downstream pressure changes resulting from the basic method. Thus we have incorporated the desirable feature of an infinite propagation speed of pressure influences in upstream direction (a property of the continuum equations) in the numerical scheme. In both the basic method and the Isra.eli/Lin source term scheme the propagation speed is only /`xi per iteration. Readers familiar with the recent literature on so- lution methods for the Reynolds-averaged Navier- St<'kes valuations in flows characterized by a dominant flow direction will have observed that the predictor- corrector scheme introduced here has some similarity with a. procedure suggested by Davis et al. in t144. The differences are essential, however. In t14] the two discretisations (4) and (5) appear also but the solu- tion princess is continued with An+. That would be all-right if the x-momentum equation stood on its own but when it is coupled with other equations in which the pressure appears as well it is a completely unsatisfactory approach. For /\t = /\x an essentially unstable scheme results. Also for smaller values of /\t we have not been able to produce acceptable results with their method. BRSIC SCHEME (GRUSS-SEIDEL) SOURCE TERM SCHEME (ISRRELI/LIN) + NEW PREDICTOR-CORRECTOR SCHEME ~2 Q (A N_ - C3 J I I I I I 1 1 1 1 1 4 8 12 16 20 24 I TERRT I ON NO . Fig. 6. Maximum pressure change between succes- sive sweeps in global iteration process for three methods 94
As an illustration of the excellent performance of our proposal we give global convergence results for the flow around an axisymmetric body as obtained by the basic scheme, the source term scheme and our new scheme in Fig. 6. In all three convergence curves two jumps appear which correspond each to a grid refinement, being executed when /`Cpma~ ~ 0.02 and < 0.01 respectively. Most striking is the gain in con- vergence rate with respect to the basic method which differs from the new method by the correction step (7) only. The more so because the correction step is extremely simple, requiring a completely insignif- icant amount of computation time. Similar results have been obtained for other cases. 6. Application Since the appearance of our earlier publications t1,2,3], the number of applications of our method has considerably increased. Some of these applications were repeated calculations for 'old' test cases (e.g. Wigley hull, SSPA 720), others concerned new explo- ratic~ns. Among them were computations including a representation of a propeller and/or a duct by a spec- ified external force field. Instead of giving results of some of those exercises or improved results for the test cases that appeared in previous publications, we have selected for presentation some data obtained for a rather demanding test case: the HSVA tanker. What makes this hull form to a difficult case for flow computations is in the first place the fullness of the hull (block coefficient = 0.85), implying a high viscous pressure resistance and presumably complex flow with a risk of boundary-layer separation. The body plan of the hull is given in Fig. 7. A 2.74 m model of it has been subjected to detailed measurements in a slotted-wall wind tunnel by Wieghardt and Kux t15] and a data tape with a part of their results has been kindly made available to us. On the basis of the geometry information provided on the data tape, a reconstruction of the hull shape was carried out with our hull Airing system. Intermediate frame lines were obtained by spline interpolation. Calculations were made at the Reynolds number of the measurements: Rn = 5 * 106. Neither the wind tunnel wall nor the rod extending from the back of the wind tunnel model were included in the simulation. The lengthwise extension of the stern-flow-ancl-wake zone was chosen as ().5 < 2;~/L < 1.65, where L is the length between perpendiculars and x _ O midships. The results to be presenter] were obtai~ecl on a 4~; * zl9 * 29 grid. Potential flow and boundary-layer calculations pro- vided us with the non-trivial part of the boundary conditions. On the external boundary, conventional conditions (prescription of tangential velocity com- ponents and pressure from the potential flow) were applied, while a zero pressure gradient condition was imposed on the outlet boundary. The initial guess for the pressure field was obtained by assuming pa = 0. One of the aims of the calculation of the stern flow field is to acquire information on the velocity distribution near the location where the propeller is to be mounted. Let us therefore first have a look at some contour plots of the axial velocity field near the stern. Fig. 8 shows the comparison of measurements and predictions in six transverse sections. In the mea- surements, the flow disturbance caused by a support wire in the wind tunnel is clearly visible. When that anomaly is taken into consideration there is gener- ally good agreement between both sets of data in the outer part of the boundary layer or wake. Near the hull, on the other hand, in the region just below the concave part of the frame lines, the measurements indicate a much stronger retardation of the flow. S shaped iso-velocity contours show up which do not clearly appear in the computed results. Usually the formation of that kind of wake distri- bution is connected with the development of a longi- tudinal vortex. How well that phenomenon is repro- duced by the calculations may be judged from Fig. 9 It shows vector plots of the transverse velocity com- ponents at various stations. From the measurements as well as the calculations only a part of the available information has been plotted but such that the vector lengths are directly comparable. ~ =~;I! At 1E it\ \ \I4\ql 119Y7 1 1 t o 95 l\` \ \ \ 1 / I I I~ \\\ ~ 11 /~J//JJ: Fig. 7. Body plan HSVA tanker ,? to 9
- - - 1 1 1 1 ~ -art K-~' ran I No L 4' --------~/l -----~-~---- r Calculated at 2~/l = 0.880 L Calculated at 2x/L = 0.90 - Measured at 2~/L = 0.885 _~ Measured at 2~/L = 0.90 7 L Calculated at 2x/L = 0.94 _____ Measured at 2~/L = 0.94 I____ it ! ://J~ J 1----- L Calculated at 2~/L = 0.980 ______ Measured at 2/L = 0.977 r ~ i ~ 1 " ' 1 Calculated at 2z/L = 1.02 ~ Calculated at 2/L = 1.10 ______ Measured at 2~/L = 1.02 ------ Measured at '22/L = l.O9 Fig. 8. Contour plots of axial velocity distributions 96
There is no clear evidence that the vortex for- mati<-'n is predicted to start at a wrong position but the vortex core seems to stay closer to the vertical symmetry plane than in the measurements, especially when is taken into account that the measurements are not truly symmetric with respect to the geomet- ric symmetry plane. This is in accordance with the stronger deceleration of the axial flow. The question remains what is the cause of the discrepancy in the axial velocity distribution. We may observe that there is little reason to suppose that a gross misprediction of the pressure is responsible. :'~ Calculated at 2~/L = 0.880 Ad_ Calculated at 2x/L = 0.90 Fig. 10 shows the girthwise pressure distribution on the hull at station 2~/l = 0.88 with an excellent agreement between measurements and calculations while the potential flow results deviate considerably. A preliminary exercise with vorticity boundary conditions on the external boundary indicated that the inclusion of viscous-inviscid interaction does not bring about either a change in the near-hull flow to the extent required for a reconciliation of measure- ments and predictions. Still to be investigated are the effects of grid refinement and of a variation in the turbulence model. tTtfTtftiltN t ~ ~ Measured at 2~/L = 0.885 ItittilltNtN ~ ~ Measured at 2x/L = 0.90 Fig. 9. Vector plots of transverse velocity components 97 \ \
f ~~ ill_ ~ :-- V~\ \ ~ \ ~ lilititt\Ntt ~ \ \ \ (3alc~llate:d at '2~/L -- 0.920 Measured at 2x/L = 0.926 rig. · ~ ~ ~ Calculated at 2x/L = 0.980 Measured at 2x/L = 0.977 Calculated at 2x/L = 1.020 Measured at 2~/L = 1.015 Fig. 9 (continued) 98
CRLCULRTED IN VISCOUS FLOW A ID ° CRLCULRTEO IN POTENTIAL FLOW A_ Al ~ MERSUREO ~ , / /// / / - o 00 0.20 0.40 0.60 O.80 1 .00 NORM . D I STANCE RLONG FRAME Fig. 10. Girthwise pressure distributions at station 2x/L = 0.88 It is not apparent from Fig. 8, but our results indi- cate boundary-layer separation at two locations: in a small region just above the level of the imaginary pro- peller axis, starting from about station 2x/L = 0.9, and near the waterline (horizontal symmetry plane) aft of station 2x/L = 0.97. In both cases the thick- ness of the separation bubble is very small. Only an oil flow experiment for the visualization of the lim- iting streamlines could shed some light on the faith- fulness of that prediction. However, the girthwise skin friction distribution at station 2x/L = 0.88, just ahead of one of the separation regions, is in fair agree- ment with the experimental data t16], as Fig. 11 shows. It may be concluded that many features of the experimental data are satisfactorily reproduced. In- deed, the present results are perhaps the most realis- tic canes among the few that have been reported so far for this test case. Yet, they can hardly be considered good enough, if the available set of measurements is considered to satisfy high quality standards. Where our method has been set up with the idea that it should give ultimately reliable results for a case like the HSVA tanker, further improvements are needed. They are to be found in a fine-tuning of the turbu- lence model en cl the use of finer grids rather than in the numerical scheme. 1 0 0 °- \ c) / CRLCULRTED MEASURED of ._ 1 1 1 1 1 1 1 1 1 No 00 0.20 0.40 0.60 0.~ 1.00 NORM . D I STANCE RLONG FRAME Fig. 11. Girthwise skin-friction distributions at station 2x/L = 0.88 7. Concluding remarks . In this paper we have given a detailed description of new developments in certain aspects of a computer code for the prediction of the stern flow field of a ship: · We have explained how coordinate-oriented ve- locity variables may be chosen so as to maintain a high numerical accuracy near the grid singu- larity in the wake; ~ We have presented some results of our first ap- plication of a new set of boundary conditions on the external boundary of the computation domain, which allow this boundary to be cho- sen close to the edge of the boundary layer and wake i We have shown how the convergence rate of the global relaxation of the solution may be improved by an easy algebraic update of the pressure field after every iteration cycle. Moreover results have been given of the application of the code to a difficult test case which indicate that the numerical scheme can cope with the complex flow around a full tanker, although the correlation with the available experimental data should be further im- proved. 99
References 1. Hoekstra, M. and Raven, H.C., "Ship boundary layer and wake calculation with a parabolised Navier-Stokes solution system", 4th Internation- al Conference on Numerical Ship Hydrodynam- ics, Washington D.C. (1985~. 2. Hoekstra, M. and Raven, H.C., "Application of a parabolised Navier-Stokes solution system to ship stern flow computation", Osaka Interna- tional Colloquium on Ship Viscous Flow, Osaka (1985~. 3. Hoekstra, M., "Computation of steady viscous flow near a ship's stern", in Notes on Numer- ical Fluid Mechanics, Vol. 17 (ea. Wesseling) 11 (Vieweg) (1986~. 4. . Bube, K.P. and StriLwerda, J.C., "Interior reg- ularity estimates for elliptic systems of differ- ence equations", SIAM Journal of Numerical Analysis, Vol. 20, No. 4 (1983~. Rubin, S.G., "Incompressible Navier-Stokes and parabolised Navier-Stokes formulations and com- putational techniques", in Computational Meth- ods in Viscous Flows, Vol. 3 in the series Recent Advances in Numerical Methods in Fluids (ea. Habashi) (Pineridge Press) (1984~. 6. Hoekstra, M., "Coordinate generation in sym- metrical interior, exterior and annular 2-D do- mains using a generalised Schwarz-Christoffel transformation", in Numerical Grid Generation in Computational Fluid Dynamics, (ea. Hauser 15 & Taylor) (Pineridge Press) (1986~. 7. Eriksson, L-E., "A study of mesh singularities and their effects on numerical errors", FFA TN 1984-10 (1984~. 8. Strikwerda, J.C., "Finite difference methods for the Stokes and Navier-Stokes equations", SIA M Journal of Scientific and Statistical Computa- tions, Vol. 5, No. 1, pp. 56-68 (1984~. 9. Oliger, J. and Sundstrom, A., "Theoretical and practical aspects of some initial boundary value problems in fluid dynamics", SIAM Journal of Applied Mechanics, Vol. 35, No. 3, pp. 419-446 (1978~. 10. Stern, F., Yoo, S.Y. and Patel, V.C., "Interac- tive and large-domain solutions of higher-order viscous-flow equations", AIA A Journal, Vol. 26, No. 9 (1988~. . Roe, P.L., "Remote boundary conditions for unsteady multidimensional aerodynamic com- putations", Computers & Fluids, Vol. 17, No. 1, pp. 221-231 (1989~. 12. Israeli, M. and Lin, A., "Numerical solution and boundary conditions for boundary-layer like flows", 8th International Conference on Numer- ical Methods in Fluid Dynamics, Aachen (1982~. 13. Roache, P.J ., "Computational Fluid Dynam- ics", Hermosa Publishers, Albuquerque, N.M. 87108 (1972~. 14. Davis, R.T., Barnett, M. and Rakich, J.V., "The calculation of supersonic viscous flows using the parabolised Na.vier-Stokes equations", Comput- ers & Fluids, Vol. 14, No. 3, pp. 197-224 (1986~. Wieghardt, K. and Kux, J., "Nomineller nach- strom auf Grund von Windkanalversuchen", Jahrbuch STG 74, pp. 303-318 (1980~. 16. Hoffmann, H.P., "Untersuchung der 3-dimen- sionalen, turbulenten Grenzschicht an einem Schiffsdoppelmodel im Windkanal", IFS Bericht No. 343 (1976~. 00
DI SCUSS ION by Y. Kodama You assumed a predominant f low direction and omitted a few terms. I'd like to hear your opinion about the ef feet of those neglected terms . Author's Reply The effect of those neglected terms on the final solution is of no practical significance in ship stern f lows. They do however inf luence the boundary condition requirements on the out l e t boundary. 101
