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.
New Viscous and Inviscid CFD Techniques for Ship Flows L. Larsson~ and L. Broberg FLOWTECH International AB4, Sweden K. J. Kim Daewoo Shipbuilding & Heavy Machinery LTD, Korea D. H. Zhang Wuhan University of Water Transportation, China ABSTRACT The research on CFD for ship flows carried out at SSPA and Chalmers Uni- versity of Technology in recent years is summarized. Methods for potential flow calculations including first and higher order theories with linear or non-linear free surface boundary con- ditions are presented. The importance of higher order effects and non- linearity is discussed. A viscous flow method based on the Reynolds-averaged Navier-Stokes equations is also intro- duced. The method may be run either in the partially parabolic mode or fully elliptically. A comparison is made between results of the two modes. Grid dependence studies and validation against experiments are presented for all methods. 1. INTRODUCTION Computational Fluid Dynamics (CFD) has been a major research area at SSPA for a number of years. During the seventies several computer programs for calculating laminar and turbulent boundary layers were developed. These methods turned out to be very useful for thin boundary layers but failed completely near the ship stern [1]. During the eighties a strong effort has therefore been made to develop more - Also Chalmers University of Techno- logy Presently at FLOWTECH International AB Presently at SSPA Maritime Consul- ting AB FLOWTECH International AB is a SSPA subsidiary specialized in CFD for external flows. 185 accurate methods, such as the higher order boundary layer integral method [2], [3], the ADI method in body-fitted coordinates [4] and the streamline curvature method [5]. All of these produced improved results near the stern, but none of them could accura- tely predict the near wake and the viscous/inviscid interaction. A more accurate approach to these problems was obviously required. The situation was anticipated around 1980, when the development of a Navier- Stokes method was started. After testing and analysing various routes to find the best possible coordinate sys- tem and solution procedure a method was developed in 1985-1988. Several diffe- rent programs were written based on somewhat different approximations of the Navier-Stokes equations [6], [7] [8] and [9]. Introduction of the pro- peller effect has been made recently [10]. In the early eighties interest was also directed towards the free surface potential flow, and the first methods, based on Dawson's theory, were deve- loped during 1983-1986 [11], [12], [13]. The call for more accurate solu- tions in certain cases prompted the development of methods with a more exact free surface boundary condition and better numerical accuracy. These methods are presented in detail in [14], [15], [16], [17], [18] and [19]. The purpose of the present paper is to introduce the new development, to give some examples of the application of the programs and to compare the results from different levels of approximation, in the potential flow as well as in the viscous flow.
2 . POTENT IAL FLOW METHODS 2.1 Governing equations In this class of methods, the flow is considered steady, inviscid, irrota- tional and incompressible. A right- handed coordinate system Oxyz is employed with the origin on the mean free surface, x positive in the direction of the uniform flow, and z positive upwards. A ship, piercing the free surface, is assumed to be in a uniform onset flow of velocity U.-. (See Fig 1.) The flow field around the ship may be described by a velocity potential ¢, which is generated by a distribution of sources on a surface S and by the uniform onset flow in the x- direction. z fx u=~ ~ p{fj,llj,Cj) Fig 1 Coordinate systems in poten- tial flow methods ¢(x,y,z) = ~ o(q)/r(p,q) dS + U-x S where o(q) is the source density on the surface element dS and r(p,q) is the distance from the point q to the field point p(x,y,z) where the potential is being evaluated. The potential ~ given in Eq (1) is governed by the Laplace equation. V2¢ = 0 in the fluid domain and satisfies the regularity condition at infinity (2) Vat = > (U-,O,O) as r -- (3) The source density ~ is to be determined from the boundary condi- tions on the hull and free surface. On the wetted hull surface the solid boundary condition is On = 0 where n denotes the outward normal. (4) At the free surface, two boundary conditions must be imposed, i e the flow must be tangent to the surface xhx + Achy _ ~ = 0 and the pressure constant gh + 2 ( Vat Vat - UP ) = 0 (6) Further, no upstream waves must be generated. The exact problem formulated above is nonlinear, since the equations (5) and (6) are nonlinear and are to be satis- fied on the wavy surface z=h(x,y), which is unknown. In the present work an iteration procedure is applied and the free surface boundary condition in each iteration is linearized about the previous solution. When the process converges the exact solution is approached. Unknown sources ~ on the hull and wavy surface z=h(x,y) will induce a potential ~ and a wave elevation h, which satisfy the boundary conditions (5) and (6). Dl(~'h) = ~xhx + Why - Liz - O D2(o,h) = h ~ 2g [ UP - (¢x+¢y +¢z ) ] = 0 Small perturbations 5o and ah of the previous solution are introduced, and the equations expanded to first order in a Taylor series. Dl(a,h) ~ Dl(~°,h°) + ADl(~'h°) + ~Dl(~°,h) ~ Dl(~°,h° ) + aa Dl(Cl,h° )~o ah D1(~ ,h)6h ~ 0 D2(a,h) ~ D2(o ,h ) + ~D2(~'h ) + AD2(~°,h) ~ D2(o ,h ) + at D2(~,h° ) + a D2(~°,h)ah - O 186 (8)
Here the superscript, °, corresponds to the previous solution, which is assumed to be known. Thus, D1(~° h°) = ~xhx + Why - fez D2(~° h°) = h° - 2g [Up - (~x + By + Liz)] where ~ is the potential induced by a° on the free surface z = h°(x y). Introducing the perturbation 5¢, the partial increments of D1 and D2 may be written AD1(o h°) = 6¢xhx + 5¢yhy Fez (10) AD1(~° h) = ~X5hx + Why + (~xzhx + Achy - ~zz)5h AD2(o h ) = g(~x5¢x + byway Z z (11) (~° h) = ah + g(~x~xz + ~Y.YZ + ~z~zz)6h Therefore the linearized free surface boundary conditions are xhx + Why - fez + ~X5hx + Achy + (~Xzhx+~yzhy ~zz)5h = 0 (12) { 1 + g(.x~xz+.y~yz+~z~zz) }ah = 1 { U 2 ~ 2 ~ 2 ~ 2 (~x6¢x+~y5¢y+~z6¢Z) } -h° (12) and (13) are to be applied on z=h°(x,y) (13) To start the iteration process the double model case is first computed. This yields ~ and ho, and a linear solution ¢, h may be obtained from (12) and (13). However, for this particular case Liz, as well as z- derivatives of fix and by are equal to zero, for symmetry reasons, so (12) and (13) are simplified to h ° + ~ h ° - ~ + ~ ah + Why = 0 (14) Sh - {u2+~x+~y - 2(~x~x+~y~y)} / 2g These are the equations solved in the linear methods. They may be applied either on z=0 [11], [19] or on z=h°(x,y) [15], the so called Bernoulli wave. To obtain a solution satisfying the exact boundary conditions (non-linear in ~ and h) the iterations continue from the linear solution using (12) and (13). After 5-10 iterations ah and 5¢ are usually sufficiently small, and the process may be considered converged. 2.2 Numerical Solution 2.2.1 Free Surface Grid In the numerical implementation the first step is to discretize the integral domain S in (1) into a large number of small panels. The domain consists of the hull surface and part of the free surface. In principle the free surface panel distribution should be extended over the entire region where there are significant waves, but for practical reasons the region of free surface panels has to be limited to the area near the hull. The influence of the truncation will be discussed later. To generate the free surface grid, two different principles have been used. Following Dawson [26], Xia Ill] generated a grid based on the stream- lines of the double model flow. This approach has the advantage that a somewhat simpler form of the boundary conditions (14) and (15) may be used in the linear case. Dawson showed that the equations can be rewritten as ordinary differential equations in l, when l is the arc length along a streamline. Using a streamline grid thus enables simple numerical differentiation (Dawsons derivation was not quite correct, but the error is small).There is, however, one serious disadvantage of the streamline grid, as noted by Xia & Larsson [12]: the grid gets very coarse at the ends of the hull, where a fine grid is required. Therefore another approach has been chosen in the more recent methods t12]-[19]. The longitudinal lines are now generated numerically by interpolation between the edge of the panelled region and the hull. The transverse lines can be 187
either hyperbolas [11]-[17] or straight lines [18], [19]. A typical grid is shown in Fig.2. For non-linear calcu- lations the free surface panels are moved to the wavy surface in each iteration. tax Fig 2 Panel distribution. SSPA Ro-Ro ship, medium bulb To obtain derivatives of a function in the x and y directions, derivatives along the L (longitudinal) and T (transverse) directions are obtained from gx = aXlgT - aX2gL gy = gT where all = fL r ~ ax2 = - i1 + fL2 and y = fL(x) is the equation for a longitudinal line. (16) Derivatives with respect to L and T are obtained either from three point or four point finite difference operators gLi = CAigi + CBigi_NL + CCigi-2NL gTi = GAigi + GBigi+l + GCigi+2 or gLi = CAigi + CBigi_NL + CCigi-2NL + CDigi_3NL gTi = GAigi + GBigi+1 + GCigi+2 + GDigi+3 (17) where NL is the number of longitudinal strips on the free surface. The coefficients of the upstream operator (CA, CB, CC and CD) and the downstream operator (GA, GB, GC, GD) are calculated from the distance bet- ween successive control points on the free surface panels, see [11]. Dawson conjectured that the use of upstream differences along the stream- lines should prevent upstream waves. This has been true of all cases com- puted by the present methods even though the streamlines have been re- placed by the body-fitted longitudinal lines. During the iteration process in the nonlinear procedure the calculation domain is changed in each iteration. The free surface sources are moved to the previous wavy surface in every step and the boundary condition is applied there. This means that the hull surface has to change as well, so a special panelization procedure has been intro- duced in the program. A typical panel distribution is shown in Fig. 3. Fig 3 Hull panel distribution in the non-linear method. Series 60, CB = 0.60. Fn = 0.32 2.2.2. Higher Order Effects Xia's methods tll]-[14] were of first order i.e. the panels on the hull and free surface were flat and with a con- stant source strength. To obtain a more accurate solution without increasing the number of panels Ni tl5]-[18] introduced a higher order technique in which the panels are parabolic and the source distribution linear. The tech- nique, which is similar to the one pro- posed by Hess t20], has been further improved by Kim [19]. A parabolic panel is defined in following form (' = ZO + Ao: + BOn + po:~2 + 2Qo:'nt + ROn'2 (18) where a panel element co-ordinate system (:', n ' , ~ ~ ) is constructed using the four corner points. Thus, the 188
origin is defined as the average of the corner points, while the ;'-direction is normal to the two diagonal vectors. The six parameters (Zo,Ao,Bo,Po,Qo and Ro) are determined by: (a) requiring the panel to pass through the corner points of the panel (four conditions), and (b) requiring the panel to pass as closely as possible, in a least squares sense, to the eight additional input points for the adjacent panels (two additional conditions). In the present method the control point is the point on the curved panel closest to the average point of the four input points. This control point is characterized by the condition that the vector from the origin is parallel to the local normal vector [: ei+n ej+( (:',p')ek] x t-(' t -A 'n .ej+ek] = 0 Eq (19) is equivalent to the two scalar equations G(:' ,n ' ) = ~ '+; '(a ' ,n ~ )~: ~ (;',n ' ) = 0 (20) H(~t,n ~ ) = n t+` ~ (a ~ ,n ~ )¢n ~ (a ~ ,n ' ) = 0 These nonlinear equations are solved by Newton-Raphson iteration. Once the control point is determined, the panel co-ordinate system (I', n ' , ~ ' is transformed to a new projected flat panel which is tangent to the parabolic panel and the tangent point is both the control point and the origin of new panel co-ordinate system (gin, ~ ) see Fig 1. The equation of the panel may be written F(~'n'() = (_ [p:2 + 2Q:n + Rn2] = 0 The source density distribution is assumed to be linear on each panel. () = ~° + ~~: + Ann (22) This is fitted in the least squares sense to the values of source density at the control points of the M adjacent panels. Thus the source derivatives on the panel in question may be expressed in terms of the unknown values of the source density on the adjacent panels in the form M (I) =ok k M (a ) = ~ Ck Ok n k=o (23) The desired source density coefficients Cal) and Cin) are obtained by minimizing the error function Min ~ M (l9) ct?)Ckn ) k-o[ k (~°+~:k+~n~k)]2 (24) According to Eq (1) the perturbation potential Bid at the i-th field point ( ~ i' n i';i) induced by the j-th panel on which the source density ~ is distributed is ij = |Aj(~/r) dS (25) Since in the higher order method the integration is to be carried out over a curved surface, expressions for 1/r and dS are different from the corresponding first order ones. The expression for Aid is further complicated by the variation of ~ on each panel. In [18] the following expression for Ail is derived hit = .(j°) j + { P.ijP) + 2Q~.jQ) + Ouija) }aj + ~(jl6)~ + ~(1~)~ = A* ~ + ~ { C( ~ )~( 14 ) (21) J J k=1 k lo + Ckn)~ij1n)}~k (26) Here M is the number of adjacent panels. 189
All the ~ terms in (26) can be calculated numerically and they may be interpreted as follows: () corresponds to a flat panel with a constant source density, .(~), ~(Q) and () are caused by the parabolic panel shape, () and ~(\n] come from the linear variation of the source density. All these terms depend only on the geometry of the j-th panel, but the curvature terms P. Q. R and the source derivative coefficients Cog) and con) depend on the surrounding panels. as NE · = ~ X1 j=1 The individual It's in Eq (26) are ~yi Fiji) = iAj r dads Fiji ) = jAj (i;~/rf dads Fiji) = iAj (i~j~j/rf dada (27) (R) = iJ |Aj (ing/r3f dads (14)= |Ajfj/rf duds ~ij(ln) = lAj nj/rf dads where rf is the distance between (ti' ni' (i) and the j-th control point ( Id, nj, O) on the flat panel (see Fig 1). Then the velocity induced by the panel is obtained by taking the gradient of the corresponding potential. At= ~~aj + E {ck(~)~(l4)+ck(n)~(ln)}~k n nisi +k[-l{ck(~)~(l()+ck(n)~(ln)}~ (28) E1{Ck(~)~(l;)+Ck(n)~(lM) The subscript ij is omitted in the equation for simplicity, *and new velocity terms (~*, .*, ~ ) are intro- duced. They are velocity Components induced at the i-th control point by a unit source density on the j-th panel. All the induced velocities in (28) can be obtained analytically. The velocities the reference coordinate can be calculated from the coordinate transformation :~x) all a21 all ~ 4] Bye = al2 a22 a32 ~ ~~ t (29) zJ al3 a23 a33 l rJ These velocities may also be written am a21 a31 a19 all aR2 Xi j (Sj ~ Uco NE ~ Yi-~- j=1 J J NE ~ Zing. j=1 J J (30) Here the velocity influence coeffi- cients Xij, Yij and Zij are the velocity components in the reference coordinate system (x, y, z) at the i-th control point, induced by a source density which is unity at the control point of the j-th panel. In the first order method the velocity induced by a panel depends only on the panel itself. The essentially new feature of the higher order method is that the velocity induced by a panel depends on the values of source density also at the control points of adjacent elements. Thus the influence coefficients for a panel depend not only on the geometry of that panel but also on the geometry of adjacent panels. The second order derivatives of Fig required in (12) and (13) are more complicated than the first order ones. In the present work, the magnitude of the second derivative terms related to the curvature and linear variation of source density are assumed to be small and vanish rapidly during the iteration. Therefore the second order derivative terms are calculated from a source velocity which corresponds to the flat panel with constant source density. xz NE also ~yz j-1 a23 J fizz i a33 j | all a21 a3l ~13 a23 an 190 fall a22 a32 ~ nt Ann Ant (31) ; din ~~ ij
NE is here again the total number of panels. 2.2.3 Solution Procedure The boundary conditions (4) and (12), (13) or (14), (15) may be transformed into a set of linear equations in ~ using the relations above for the first and second derivations of A. [ A ] { ~ } = { B } (32) where tA] is a NE x NE matrix and {B} is the right-hand side array. The upper part of the A-matrix and B- vector, corresponding to the hull boundary condition, is obtained in a straightforward way, while the lower part, corresponding to the free surface, is considerably more complex. The equations are given in [18] and [19]. Unfortunately, the lower part of the matrix is not diagonally dominant, so iterative methods cannot be used in the solution. This is therefore found by a Gaussian elimination procedure. 2.2.4 Wave Resistance Having obtained the source strengths, velocities may be computed from (30) and the pressures from the Bernoulli equation. To obtain the resistance, two principally different methods have been employed. Either the force is found by pressure integration over the hull or by integration of the sources over the free surface. In the former case, different formulas are used depending on whether or not the higher order technique is employed, i.e. NB NB W i CPi NXi Nisi / ~ Nisi (33) where NB is the number of panels and Nx and ASi are the x-component of the normal and the area respectively of a panel, or NB JC i E ( P° PC ~ Purl) (a11N(; + a21Nrl + a31N; ) (34) |NB ( 1+2~2 )d~dq| ~ |Ei ( 1+2¢ 2 )d~dn The latter formula is derived assuming that the pressure varies linearly over the panel in a similar way as the source strength. Cp: and Cp~ are the slopes of the pressure variation, all, a21 and a31 are the first column of the transformation matrix between the panel and reference coordinate systems, N; E, N. E and N: E are the components of the Formal vector in the panel coordinate system and ~ 2 is a term which takes into account the change in panel area due to curvature. The formula for obtaining the resistance from the free surface sources is derived from the momentum theorem applied to a control volume outside the hull, under the undistrubed free surface and extending to infinity sidewards, downwards and longitudi- nally. The detailed derivation is given in [18]. It yields NF NB Cw = (4n ~ UBici~Si) / (USA ASi) (35) 1 i It should be noted that this formula is valid only for the linear case. 2.3 Validation The methods presented have been validated in a number of ways. A grid dependence study was first carried out for the linear first order method to investigate the sensitivity of the solution to the panelization on the free surface. Most likely the results from this study are valid also for the other methods. The improvement in accuracy when introducing the higher order method was then tested against two analytical results without a free surface. Comparison with measurements have been made for a number of cases such as the Wigley hull, the Series 60, CB=0.60 hull, the HSVA-tanker, the high speed ship ATHENA, a Ro-Ro ship designed at SSPA and two sailing yachts. All these comparisons are described in detail in [11]-[19]. As an example only the calculations for the Ro-Ro ship will be presented here. This case is parti- cularly interesting since a study of several bulb alternatives was made. It is also the only case for which calcu- lations have been carried out using the first and higher order linear, as well as the higher order non-linear method. Interesting comparisons between all three methods may thus be made. 2.3.1 Grid Dependence Study Since one of the most important questions in connection with the present methods is the dependence of the solution on the panel distribution, particularly on the free surface, a systematic variation was carried out. 191
The test case was the Wigley hull. Employing the coordinate system defined in Fig 1, the equation for the hull surface reads y = 10 ( 1 - x2) ( 1 ~ Z64 ) (36) In table 1 (last page) the number of panels on the hull and on the surface, the extent of the panelled part of the free surface in the x- and y- directions, the Froude number and the wave resistance coefficient are given. In run No 1, only approximately 1/4 of the half ship length, 1, was used for the free surface extent in the y- direction. The calculation gave an unresonably high wave resistance at Fn = 0.266 and broke down at Fn = 0~45. In runs of No 2, 3 and 4, the free surface region was enlarged and divided into 5, 10, and 15 strips respectively. At Fn = 0.266, the resistance changed signi- ficantly from run to run, but at Fn = 0.45 there was no change between the last two runs. The free surface extent in the x- direction was tested in runs No 3, 5, 6, 7 and 8. From No 3 to No 5 the free surface region was extended by 1/4 1. The wave resistance did not change very much, but as the free surface region was enlarged by 1/2 from run No 6 to No 7, the wave resistance changed, particularly for the highest Froude numbers. From run No 7 to run No 8, part of the free surface region was moved from the bow to the stern. The wave resistance did not change at Fn = 0.266, whereas it was increased by 10% viva at Fn = 0 45 There were no systematic changes in the panel distribution on the hull, but two slightly different hull panel distributions were used in the successive tests. This should not cause large changes in the final results. The considerable difference in the wave resistance between runs No 3 and No 6 is likely to have been caused by the change of the panels on the free surface. More even panels were placed on the free surface for run No 6 than those for case No 3. The results obtained imply that the choice of the panels depends on the Froude number to be evaluated. Smaller panels should be used for low Froude numbers, while a larger free surface portion is required for high Froude numbers. It is wise to vary the panel arrangement for different Froude numbers. A panel distribution similar to case No 9 of Table 1 has been used in most work afterwards. 2.3.2 Analytical Test Cases To investigate the improvement in accuracy when introducing the higher order technique the flow around a sphere was calculated. In Fig 4 the velocity distribution along a generator through the stagnation point is shown. Two grids were tested, one with a distribution similar to a typical ship case, i.e. with 22 longitudinal and 40 transverse strips on one quarter of the sphere and one with 5xlO strips. It is seen that both higher order solutions are extremely accurate, while there is some error also for the 880 panel first order case. Small as this error may seem, it may have a significant in- fluence on the drag. Integrating the pressure on one half of the sphere yields the following results for the drag coefficient CX First order Higher order Exact Cx 0.0550 0.0622 0.0625 Error 12% 0.5% - ~ Anolytic solution + First order solution 0 Higher order solution First order solution O Higher order solution .0 0.5 r' B80 effective panels 50 effective panels ,. V:,': ~ v< , , ~ 30 60 (3 90 Fig 4 Velocity distribution on a sphere 192
The higher order calculation is thus considerably more accurate. It should be noted that the computer time for a given number of panels increases by about 10% for a case of this size. z.3.4 The SSPA Ro-Ro Ship The body plan of the SSPA Ro-Ro ship is shown in Fig 5 with a small, a medium and a large bulb. A comparison between the measured residuary resistance and the computed wave resistance (pressure integration) using two linear (first and higher order) versions and one non-linear (higher order) version of the same program [19] is made for the medium bulb in Fig 6. The panel distribution is shown in Fig 2. Fig 5 Body plan, SSPA Ro-Ro ship with different bulbs C W X 1 0 3 3.0, 2.0 1.0 O _ Fig 6 0 0.1 - Measurements Linear solution (first order) · Linear solution (higher order) · Non - I inear solution _. .1 ,~^ 0.2 En i l ,, 0.3 Predicted wave resistance compared with measured residual resistance. SSPA Ro-Ro ship, medium bulb A common problem in most evaluations of wave resistance calculations is that the wave pattern resistance is seldom measured. In the present case the form factor has been determined using a Prohaska plot, and the viscous resistance based on the ITTC-57 corre- lation and the form factor has been subtracted from the total resistance to get the residuary resistance. Cwx10-3 0.6 r 0.5 ~ J 0.4 03 02 01 F., 0.15 021 025 1 ~ Medi um bulb Otto 0~0 O 1 2 3 4 5 6 7 I terati on Fig 7 Convergence history. Wave resistance of SSPA Ro-Ro . ship, medium bulb As appears from Fig 6 the difference between the linear and non-linear higher order methods is quite small for this case. This is also apparent from Fig 7 where the convergence history of the wave resistance is shown. The final converged value after 7 iterations is not very much different from the first value, corresponding to the linear solution. Larger differences have been noted for bluffer ships [16] and non- linearity should be extremely important for hulls with large overhangs such as sailing yachts. 193
2h Calculated 0.02 .. Bulb /; />\\ ~ ~ _, -0.01 No bulb L Measured 0.02 Hi o.olL I I I ,0 2X/L 0;0 2 X/L Fig 8 Wave profiles at En = 0.21. SSPA Ro-Ro ship CWX 10-3 an, . 2.0 1.0 Type Measurement Calculation No bulb Small bulb ___ Medium bulb Large bu l b 1. 11' hi: /i,1 :/'1 ~/~/// 0.1 an Fig 9 Wave resistance. SSPA Ro-Ro ship The first order solution of Fig 6 is not as accurate as the higher order solutions, but even the first order results are not grossly in error. Interesting comparisons between results for the three bulbs and the case without a bulb can be made in Fig 8 and 9 based oh higher order linear solutions. The measured wave profiles of the four cases are very much different as appears from Fig 8, and the differences are quite well predicted. The absolute values of the wave resistance, Fig 9, are not as accurate in all cases, but the method is able to rank the cases in the right order. 3 VISCOUS FLOW METHODS Although r as indicated in the Introduction' a number of different methods for computing the viscous flow (boundary layer/wake) around ships have been developed at SSPA/CTH, only the Navier-Stokes methods will be presented in the present review. 3.1 Coordinate systems and grid In the SSPA Navier-Stokes methods the coordinate systems of Fig 10 are employed. These systems, which differ from the ones of the potential flow methods, are more convenient for tensor notation which is required for the fully transformed Navier-Stokes equations in curvilinear nonortho- gonal coordinates. ;~ are the global Cartesian reference coordinates, while xi are the grid oriented ones, as appears from Fig 10. The relation between the two systems is of the form V2xi = fi (i = 1, 2, 3) (37) where v2 is the Laplacian.operator in orthogonal coordinates (ha) and fir are control functions which may be assigned appropriate values to yield the desired stretching of coordinate surfaces. The equations are inverted by exchanging the dependent and independent variables siiaiaj:1 + fkak~l = V251 (38) The Laplacian operator may be written 2 1 3 hlh2h3 ~ V hlh2h3 j z- 1 a a; j ( 2 - ~ 194 (39)
: ~ ~ . 1 1 1 1 1 1 1 1 Fig 10 Coordinate systems in Navier- Stokes methods here hi, h2, h3 are metric coefficients in the frame id. After selecting 41(xl) and assigning proper values for the control functions fin, equations are solved numerically. A more detailed description of the grid generation is given in [9]. Fig 11 shows the coordinate system for the SSPA 720 Model. 2 Governing equations The governing equations solved are the Reynolds equations obtained from the Navier-Stokes equations by averaging over the largest, turbulent time scale. In this way the unsteadiness due to turbulence is removed at the expense of having to introduce a turbulence model for the Reynolds stesses. In the coordinate system above the incompressible momentum equations may be written (see [6] for a derivation). Otui + h(k)ukuik = - h(i)gikakP + (gjkulj + h(i)h(k)g u,j)SkvE VE{gjkakajui + gikOk(ulQjil) + gjkQ gj kr1 ui } J , (40) q ~ F l 1 l 1 2 1 0 l n l 0 ~ 0 4 0.2 View of the water plane Fig 11 Longitudinal and transverse views of the grid. SSPA 720 Model The incompressible continuity equation may be written J 1ak(Jh(k)uk) = 0 (41) where J is the Jacobian us are the physical components of the velocity vector, i.e. the components based on unit vectors parallel to the covariant base vectors go, whose magnitude igi~ is denoted h(i), gij is tlhe contravariant metric tensor and r i j is a Christoffel symbol of the second kind. a i is to be interpreted as a partial derivative with respect to xi and ul,j is defined from the relation 195
ui,j = a jui + QjilU1 (42) where Qjl in(j) Qjl and (43) Qji1 = h(J) in(l) in(i) Fji1 _51i h(i)ajh(i) where bli is the Kronecher delta and 1 k = in(k) Ok An index written within a bracket means that it is not considered in the summation rule. The quantity vE, finally, is the effective viscosity 1 VE = R + it (46) where Rn is the Reynolds number and it is the eddy viscosity which is obtained from k2 it = Cp The following standard values are used for the constants: c = 0.09, Cal = 1.44, cat = 1.92, o~ = 1.0 and a~ = 1.3 Two turbulence models are being used. If no wall functions are employed the turbulence model is of the two-layer type. Thus, in the inner region an equation for the turbulent kinetic energy k is solved together with a prescribed length scale. In the outer part the length scale is replaced by an equation for the rate of dissipation of kinetic energy,- . If wall functions are used as the inner boundary condition for the flow, the k-E model is used throughout. The transport equations for k and e are as follows auk + h(i)uiaik = [g jaiajk ~ giiFiiDlk] (48) ~ · . ok g a jvE a ik + G - ~ at6+h(i)U aid = HE [siiaia ~ giJr1 a ~ gijajVEai~ + (49) Cal k G - ce2 k The generation term can be written G = vt(h(m)h(i'U iUlm + g gmjh(m)h(~)u iU]n) If no wall functions are used the rate of dissipation and the eddy viscosity are obtained from k3/2 ~ = Vt Cp ik 1p respectively. The following length (47) scales are used 1 ~ = Cly ( 1 - e~RY/Ac ) 1p = Cly ( 1 _ e-Ry/Ap ) (50) (51) (52) where Ry = Rn ik y , y is the distance in the normal direction to the surface. The following constants are used: C1 = ~Cp3/4 , Ad= 2C1 , Ap = 70 A= 0.418 , as reported in Chen and Patel [23] The velocity components ui, the turbulent quantities k and 6, the time t and the pressure p are made dimen- sionless in the usual way by the free- stream velocity UO, the reference length L and the density P. 3.3 Numerical solution The Reynolds equations and the transport equations for k and ~ can be rearranged into the following general form 196
gllalal~ + g22a2a2~+ g33a3a3¢ = 2A a 3¢ + 2B~ a 2¢ + D¢~¢ + Ed a to + SO where ( ) 2A = Reff(a~u3/h(3) - g3jDjvE) 2B~ = R ff(a~u2/h(2) - g2jDjvE) DO = Reff(a~u1/h(l) - g jDjvE) EN a¢Reff so = sit, - 2(gl2ala2¢ + gl3ala3¢ + g23a2a3~) al = 1 for ~ = u1, u2, u3 and ak = ok aE = ~E Reff is the effective Reynolds number and is defined as Reff = 1/VE Finally, the source terms so are sui = Refft(h (m) um - g i a jVE)Qmlu +h(i)giiDj (p + 2 k) (l (i)h(m)g u jamvE] - gems Qi 1 2 imQ1 ulj + gimQJlQlmun+ gj FJmu, Sk = gi FJmalk - Reff(G - E) 54) so = gimrJmOl~ ~ Reff(c~1 k G-CE2 k ) In most ship flows there is a predominant flow direction usually approximately aligned with the x1 direction of the grid. Under these circumstances the term a 1a 1¢ may be neglected and the elliptic equations (53)become parabolic. The advantage of this approximation is that a marching technique may be used in the solution. On the other hand, the fully elliptic equations are more accurate for complex flows, particularly if recirculation of the flow occurs. At SSPA two different programs have been developed: NASELL for the fully elliptic equations and NASPAR for the parabolized equations. The discretization scheme by Chen & Patel [21] has been adopted with some minor modifications due to the presence of the term gllDlal¢. Within each cell the coefficients are treated as con- stant. The finite-analytical scheme is used in the cross-plane and a second- order central and first-order upwind difference approximation is used for the second and first-order derivative in the xl-direction, respectively. In NASELL, this results in an eleven-point discretization formula that takes the following form using compass notations: UP NBANB¢NB + At¢P - Rob NB = D,U,N,S,E,W,NE,NW,SE,SW (55) NB CNB/p, At = CpE¢/(pdt), Rib = CpS¢/p ~ = 1 + Cp(jD¢~ + 2gl + E¢/dt) CD = (~D¢~°)max + g = (D¢,O)max + g Cs = A/g(Bk)' CN = PA - Cs CW = PB/g(Ah)' CE PB cw SW (1 - PA - PB)/g(Ah)/g(Bk) SE e Csw, CN = e~2BkC C -2Ah Cp = h (1 - PA)/f(Ah)/~ - k2( 1 - PB)/f(Bk)/2 PA = E2e g(Ah)g(Bk)f(Ah) tow PB = 1 + (h/k) (PA - l)f(Bk)/f(Ah) 197
~ - ( -1 )m( )`mh ) 2 m-1 [(Ah)2 + (~mh)2] 2cos~(A2+B2+> k h = ~ h = (m ~ 2)n g(x) = l + e~2X, f(x) = xg(x)/(l - e~2X) h = 1/igp , k = 1/iOp , A = h(A¢)p B = k(B¢)p In NASPAR,downstream points, index D, are excluded, so the discretization formula contains ten points. The compass notation and the location of points are described further in Fig 12. The superscript n-1 denotes the previous time level, while the superscript n for the present time level has been left out. The time step is denoted aft. Asymptotic expressions for PA and PB are used for large cell Reynolds numbers, following Chen & Patel t21]. ~ ~ PN ~ ~ V n ~ ups NW ~ N . ~ W ' P ~ ~ i. SW 5 I NE 1 ~ E 1 ~ SE Fig 12 Staggered grid and compass notations The scheme is not symmetric with respect to the coordinate directions, but treats the x1-direction in a simpler way. Therefore special care has to be taken with solutions produced in regions where the velocity vector is at large angles to the x1-coordinate line. In the present application to ship stern flows this problem is rarely en- countered. Longitudinal vortices occur, but without any larger deviation between the velocity and coordinate- line directions. In very few cases small local separations with reverse flow may occur in the stern-most part. (56) The flow in the transverse plane is very essential for stern flows, since it controls phenomena such as the transverse distribution of fluid and the curvature of the streamlines, which determines the pressure. The more accurate finite-analytical scheme is therefore applied in this plane. The pressure-velocity coupling is treated with the SIMPLER-algorithm, Patankar [22]. Details from a derivation of the equations for the particular coordinate system and velocity components are given in t6]. The resulting equations are Ud = Ud ~ Bd(PD ~ PP) Un = Un ~ Bn(PN ~ pp) (57) u3 = ue - Be(pE - pp) where the ui is the solution of the momentum equations using the pressure from the previous time level and B' is the coefficient for the orthogonal part of the pressure gradient in the i-th momentum equation. . . . effh(i)g Cp/' (no summation) The compass notation for the location of the velocity components in the staggered grid are shown in Fig 12. The pressure correction p is determined from the equation 198
* appp = ~ aNBPNB (58) where aNB = (JB)nb' ap = ~ aNB NB = D,U,N,S,E,W D = [u J/h(l)]u + [u J/h(2)]s + [us J/h(3)]e The pressure, finally, is calculated from the pressure equation: appp = ~ aNgPNB (60) (59) where D = [u J/h(1)]u + [u J/h(2)]s + [U3J/h(3)]e (61) The pseudo-velocities u equal the sum of the terms on the right-hand side of the corresponding momentum equation (55) , except for the source term, which is modified in such a way that the orthogonal part of the pressure gradient is taken out: ui = ~ AnbUnb + At ( Ui ) p - Ri nb = u,d,n,s,e,w,ne,nw,se,sw and Ri = Ri + B a iP (no summation) (62) In NASPAR the pressure correction equation has to be parabolic so the downstream values are taken from a previous time step. This may be expected to slow down convergence, but does not influence the final result. The methods can handle unsteady flows, but only steady flows have been considered so far. The time history is therefore merely a route to the steady state. Time accuracy is irrelevant and the largest possible time step should be used. The different steps in NASPAR, involved in the calculation of veloci- ties, turbulence quantities and pres- sure at a time level can be described as follows: 199 - The momentum, pressure correction and turbulence equations are parabolic and can therefore be solved in one x1- constant plane at a time from inlet to outlet plane. This is done as follows: o The momentum equations(55)are solved using pressure and coefficients based on velocities from previous time level. The algebraic equation systems are solved using a fixed number of line by line sweeps with a tri-diagonal matrix algorithm. Lagging and zero pressure field are used in the first time level. o The coefficients for the parabolized pressure-correction equation(58) are calculated and the algebraic equations are solved using a fixed number of line by line sweeps with a tri-diagonal matrix algorithm. The velocities are corrected using (57). O The divergence of the pseudoveloci- ties u3(61) and the coefficients of the pressure equation are stored. o The transport equations for the tur- bulent kinetic energy and the rate of turbulent energy dissipation (54)are solved using a fixed number of line by line sweeps with a tri-diagonal matrix algorithm. - When the marching procedure is com- pleted, i e the last x1 = constant plane is reached, the elliptic pressure equation (60) is solved. Two techniques have been used to solve the resulting algebraic equations. One in which the pressure has been updated in one plane at a time, sweeping upstream from the outlet to the inlet plane. A line by line sweep technique is then used in each plane. The other solution tech- nique is complete line by line sweep with a tri-diagonal matrix algorithm. - The calculation continues at the next time level and all the above steps are repeated. In the elliptic program NASELL the marching technique cannot be employed. The steps are then as follows - Velocities, pressures and turbu- lence quantities in the solution domain are estimated, usually from a NASPAR solution. - The momentum equations (55) are solved based on the known pressure using a line by line iteration technique. - The velocity field is corrected to become divergence free using (57) and (58) Stones strongly implicit algorithm is used.
- The pressure is computed from (60) using the same algorithm. - The transport equations for k and are solved using the line by line tech- nique. - The same steps are repeated at the next time level. In both methods the solution is assumed converged when the global divergence is below a critical value. 3.4 Boundary Conditions The calculation domain is shown in Fig 11. (1) Upstream and initial condition The profiles of the velocity components us and the turbulent quantities k and ~ are required at the inflow section. Only standard boundary layer profiles for ul, k and ~ (U2 = us = 0) have been used in the calcualtions so far. No condition for the pressure is required. Lagging has been used during the first time step for the velocity and the turbulent quantities and the pressure field is assumed to be con- stant and equal to zero. (2) Outflow condition In the present calculations the outflow section is located in the far wake and a zero pressure gradient condition alp = 0 is used. (3) Outer edge The outer plane may be located far from the body, where the flow is nearly undisturbed and thus p = 0, ul = 1,u3 = 0 , a2k = a2 ~ = 0 and the normal velocity component u2 is calculated from the local continuity. Alternatively, p, u1 and us may be taken from a potential flow solution. (4) Symmetry planes The undisturbed water surface and the centerline plane are treated as symme- try planes, i e us = 0, 33u1 = 33u2 = ask = a3E = 0 (5) Wake center plane The following conditions have been applied at wake center planes that are degenerated to a single line 200 u2 = us = 0, a2u1 = a2k = a26 = 0 otherwise U2 = 0, 32U1 = a2u3 = a2k = a2c = 0 For the body boundary condition there are two alternatives: either the computational domain may be extended into the viscous sublayer and the no- slip condition applied (this requires a two-layer turbulence model) or the wall law may be used. In the former case the conditions are simply ul = u2 = us = k = 0 If wall functions are used a procedure similar to the one of Chen & Patel [21] has been employed. The procedure requires that the two points closest to the wall are located in the logarithmic layer. It is assumed that the velocity at the innermost point is parallel to the wall which yields u2 = 0 at x2 = 2 and that the velocity vector does not rotate between the first and second point in the plane parallel to the wall. The last assump- tion yields ( q ~ x2=2 =(q ~ X2=3 (q X2=2 =(q- x2=3 (63) where q is the magnitude of the velo- city. The procedure is iterative and starts by calculating the magnitude of the velocity q2 at x2 = 2 with an assumed skin friction velocity us from the wall law q2 = ~ 1 lny+2 + B where K= 0.42,B = 5.45 and y+ = Re use. The velocity components at x2 = 2 can then be calculated from equation(63). The components serve as boundary condi- tions for the momentum equations. q3 can be calculated from the solution and an updated u can be calculated from the wall law applied at point x2 = 3. The procedure is repeated until conver- gence. The boundary conditions at x2 = 2 for the k- and E-equations are uT2 k2 = ~ up 62 = ~Y2
3.5 Validation 3.5.1 Grid dependence study A study was carried out to investi- gate how the solutions depend on the grid and what the proper number of clusters is in each coordinate direction for ship stern flow calcu- lations. Due to computer limitations the study had to be rather restricted, however. The systematic variations were made using four different grids to investi- gate the grid dependence in the three coordinate directions. For all four grids the calculation domain was the same, that is the inlet plane was placed at the midship, 2x/L = 0.0; the outlet plane at 2x/L = 10 in the far wake; the outer edge boundary on a cir- cular cylinder located one ship length from the centre line; and the two innermost grid points were located within the logarithmic layer. In case 1 the grid has a number of clusters equal to LLxMMxNN = 60 x 19 x 9, where LL, MM and NN2are the number of clusters in the x1, x and x3 directions respectively. The number of clusters in the transverse direction x3 in case 1 might be too small to capture the variation of the cross flow and the rapid pres- sure change near the region of the keel. Thus in case 2, NN was increased to 15. The grid in case 2 has the number of clusters equal to LLxMMxNN = 60 x 21 x 15. In case 3, the grid was designed to test the grid dependence in the pre- dominant flow direction x1, LL was decreased to 30. The grid had the number of clusters equal to LLxMMxNN = 30 x 21 x 15. The velocity profiles and the pro- files of other physical quantities in the normal direction might be quite influenced by the number of clusters in the normal direction x2. Based on the grid of case 3, we increased the number of clusters in the normal direction to MM = 41, getting the grid case 4: LLx MMxNN = 30 x 41 x 15. A summary of the four cases is given in Table 2. Case LL 1 60 2 60 3 30 4 30 MM NN 19 9 21 15 21 15 41 15 Total number of clusters 10 260 18 900 9 450 18 450 Table 2. The four grid variation cases The calculations for the four cases were carried out using standard boun- da2ry later profiles for ul, k and ~ (u = u = 0) as the inlet profiles. \ Of 2X/L = 0.8 o N l n __ i a 80X19X9 - - 80X21X15 30XZiXi5 ...... 30X41X15 , l l l .00 0 . 05 0. 10 0. 15 0 . 20 GIRTH LENGTH 2X/L= 0.8 60XI9X9 BOX21X15 · 30X21X15 ..... - poX41X15 o . Go 0.05 0. 10 0. 15 0.20 GIRTH LENGTH Fig 13 Grid dependence. Distribution of Cp and up at 2x/L = 0.8 The girthwise pressure and friction velocity variations at 2x/L = 0.8 for the four grids are shown in Fig 13. The influence of the cluster numbers in the transverse direction can clearly be seen from the figure. Comparing the grid case 1, 60 x 19 x 9, with other grids it is obvious that nine clusters are too few to resolve the problem in the transverse direction. The Cp variation is almost the same for the other three grids, since the same number of clusters was used in the transverse direction. Owing to the 201
limitation of the computer size we did not increase the number of clusters in the transverse direction to more than 15. So it is difficult to say whether or not the solution would be improved if more than 15 clusters were used. This number seems to be the minimum required, however. O _ . ~X4iX15 - - 30~1Xi5 o 0° Waterl i ne N _ ~ O DO ~ 1.00 U 0 - - X41X15 - - ~~lX15 45° / 0 00 0.~ t.00 U - ~X41Xi5 _ ~~iXi5 90° Keel o.oo 0.50 t.oo u it, ~ Grid 60 x l9 x 9 I' ///// / - / Grid 60 x 21 x 15 ////// / Grid 31 x 41 x 15 /~/~ ///// / 11// A// Grid 31 x 21 x 15 Fig 14 Grid dependence. Velocity profiles at 2x/L = 0.8 Fig 15 Grid dependence. Iso-wakes at Waterline, mid-girth and keel 2x/L = 0.9 202
The influence of the number of clusters in the normal direction can be tested by comparing the results of grid 3 and grid 4. The total velocity pro- files at 2x/L = 0.8 along the water line, keel and at the bilge are shown in Fig 14. There are small discrepan- cies between the results for the two grids. Also in Fig 13 the girthwise up variation of the grid 30 x 41 x 15 is only slightly different than the others. This means that there is a very small grid sensitivity in the normal direction. If we look at the isowakes from the results of the four grids at 2x/L ~ 0.91 shown in Fig 15, they correspond quite well with each other. The pressure and friction velocity ° variations along the water line and the keel from the results of grid 2 and 3 are shown in Fig 16. Grids 1 and 4 are not shown in this figure, since their (xl, x2)-planes are slightly different. 60X2iXt5 30X2iXt5 C1 CO o n CO O _ o O _ tu O _ ,,\ Watertine 1 1 1 1 1 1 1 1 1 0.0 O.Z 0.4 0.6 O.B 1.0 1.2 I.4 i.B 1.8 2.0 2X/L 60X2iX15 · 30X21X15 Keel 0 0 0.2 0.4 0.6 O.B 1.0 t.2 1.4 1.6 1.B 2.0 2X/L ~ 1 1 0.0 O.Z 0.4 0.6 O.B ZX/L o o * ,_, o ~ A o CU o 80XZlXIS 30XZlX15 Waterline 1 1 1 1.0 60XZiX15 · 30XZlX15 Keel 1 1 1 0.4 0.6 ZX/L O.B 1.0 Fig 16 Grid dependence. Cp and up along the waterline and keel The influence of the number of clusters in the xl-direction is tested from this comparison. The results are almost the same in the whole region. Increasing the number of clusters from 30 to 60 gives only very small changes. It is therefore concluded that 30 clusters is enough in this direction. These limited systematic studies of grid-dependence indicate that there is some grid sensitivity in the cross section, but not in the longitudinal direction. Special attention must be paid to regions of abrupt change in geometry, such as the stern region and near the keel. Concentration of clusters to these regions may be an adequate alternative to fine grids. Considering this possibility a grid of 30 x 20 x 15 would be proper for practical use. 203
3.5.2 Location of the outer edge of the computational domain As indicated above, the outer edge of the computational domain can be located so far from the body that free stream conditions may be specified. Another alternative is to obtain velocities (except for the normal component) and pressure from a potential flow solu- tion. However, outside the boundary- layer there is a region in which the viscous and inviscid flow interact. To apply the boundary condition here an iterative procedure between the poten- tial and viscous solutions must be in- troduced. To investigate the location of the region of viscous-inviscid interaction for the SSPA Model a series of calculations was carried out using different locations of the outer edge. The largest distance to the edge mea- sured from the keel is 1.08 (dimension- less by L/2) for a grid with 60 x 19 x 15 clusters. The distance was varied by dropping clusters in the normal direction so that the distance was re- duced to 0.51, 0.24, 0.11 and 0.05 respectively. The last grid consists of 60 x 11 x 15 clusters. For each calcu- lation the boundary conditions at the outer edge and at the inflow (outside the boundary layer) were calculated from the potential flow solution. The resulting pressure-fields are normalized at a point by subtracting a constant in the following comparisons. Fig 17 shows pressure profiles at four different x-stations, distinguished by symbols. The symbols are placed at the corresponding outer edge. The line types distinguish different calcula- tions, i e different locations of the outer edge. The profiles shown in Fig 17 (N = 14) are from the region close to the keel where the boundary layer is very thin, and the profiles shown in Fig 17 (N = 7) are from bilge region, where the boundary layer is thick. Both figures show that the profiles from the calculation with the outer edge located at 0.05 are different from the others, indicating that the boundary conditions have been applied in the interaction region. Thus, the interaction region extends out to a surface which is located between the outer edges of the calculation domain 0.11 and 0.05. A calculation of the potential flow and the laminar and thin turbulent boundary-layer normally precedes the stern-flow calculation. The cost for using the potentialflow condition is therefore low. The gain would be a reduction of clusters in the normal direction and a convenient way to handle blockage-effects in tunnel simulations. Variation of outer surf ace (N-14) 2X/L-0.7 ZX/L-O .8 2)(/L-0 .8 2X/L-0.95 · 2Are/L-t . 03 2Are/L-0.51 2Are/L-0.24 - 2Are/L-O.~t ~ 2Are/L-0.05 0.0 O. 1 O.Z Cp N Var I at ion of outer sun f ace tN-7) - x. Or ID 0 _ ~ _ x 2X/L-0 .7 ~ 2X/L-D . ~ + 2x/L-o . g 0 2X/L-0 .85 · 2Are/L-I. OB 2Are/L-0 . 51 2Are/L-0.24 2Are/L-0 . ~ 1 - 2Are/L-0 . 05 Be.! 0.0 0.t 0.2 Cp - Fig 17 Variation in outer edge location. Cp distribution in the x2 direction 3.5.3 Comparison between oarabolized and elliptic solution The comparion between the parabolized and elliptic methods was carried out for the SSPA 720 Model with a grid consisting of 60 x 19 x 15 clusters. The inlet profiles were standard boundary-layer profiles, normal and transverse velocity 204
components were zero and a constant velocity u1 = was outside the boundary- layer. The solution after one sweep from the partially-parabolic method was used as input for the fully-elliptic algorithm. The pressure and friction velocity variation in the girthwise direction at 2x/L = 0.8 are shown in Fig 18. The partially-parabolic and fully-elliptic algorithms give almost the same results. No noticeable difference between the solutions can be seen. This is not surprising since the only simplification made in the parabolized algorithm, as mentioned before, is that the second derivative ( a 1 a 1 ¢) in the predominant direction is neglected. This higher order term is of signi- ficant importance only for bluff stern flow with separation, but not for slender sterns, like the one of the SSPA model. Thus the two algorithms should give the same results. The CPU time for the partially- parabolic algorithm is about 70% of that of the fully elliptic one. O SSPA720 Cp AT 2X/L - O.B 0 0 0 . 1 0 . ~ 0 . 2 0 . 2 N GIRTH LENGTH T 0.0 0.2 0.4 0.6 O.B 1.0 1.2 t.4 1.6 I.8 2.0 2X/L cat o l 0 lo lo 0 lo lo - O 0 _ ELLIPTIC Pi"B~IC SSPA720 Ut AT ZX/L - O.B ~ ELLIPTIC PARABOLIC cam L, ' l l 1 lo 0 0 0.! 0.1 0.2 0.2 ,o GIRTH LENGTH r Fig 18 Comparison between results of elliptical parabolized methods 3.5.4 Comparison with experimental data The grid 60 x 21 x 15 was used to make the final calculations. The inlet profile at midship (2x/L = 0.0) for the longitudinal velocity component U within the boundary layer was deter- mined from the wall law using the measured boundary layer thickness and friction velocity from Larsson. The normal component V and the transverse component W were assumed to be zero. Outside the boundary layer the distribution of U was calculated using Hess & Smith's potential flow calcula- tion method. The k and ~ within the boundary layer were calculated according to Klebanoff & Bradshaw as presented in [24] and set to zero out- side. The calculated results for the Reynolds number Re = UoL/ = 5.0 x 106, were compared with the experimental data of Larsson [25]. Owing to the different coordinate systems used in the calculations and the experiments it is difficult to make direct comparisons of all quantities of interest. Therefore only limited comparisons are presented here. fir o SSPA720 Cp AT O DE6. mu lo o 0 _ ~ I~ tic \ - 60X21Xt5 0 SSPA720 Cp AT 90 DE6. ru 0 _ o 0 _ 60X2IXI5 l 0.0 0.2 0.4 0.6 O.B I.0 1.2 I.4 1.6 t.8 2.0 2XIL Fig 19 Comparison with experiments. Cp distribution along water- 1lne and keel 205
The pressure distribution on the hull surface along the water line and the wake centerline as well as along the keel are prsented in Fig 19 along with the measurements. The original experimental data of Larsson[25] are subject to windtunnel blockage effects. Corrections have been introduced in the figures using the formula given by Larsson[25], which was formulated after carrying out potential-flow calcu- lations with and without tunnel constraints. Fig 20 shows the girthwise pressure distributions at three stations. It is clear from these figures that the agreement between calculations and measurements is quite good. o ° r SSPA720 Cp AT 2XJL - O .70 l ._ ~/? c MEASUPEH£NT 60X2IXt5 O 00 0.05 0.10 0.15 O. 20 GIRTH LEN6TH SSPA720 Cp AT 2X/L - O.BO - ... 1° _ to N O' I _ 0.00 0.05 ~ MEASUREMENT - ~nx2lxts _ . _ 0.10 0.15 0.20 GIRTH LENGTH SSPA720 Cp AT ZX/L - O.S o to o can to o 0 _, o to up SSPA720 Ut AT O DEG. BOX21X15 0.2 0.4 0.6 0.8 1.0 Z)(/L SSPA720 Ut AT 90 DEB. · Rnx~~Xls 010 O'Z 0'4 0~6 ·11 110 2X/L Fig 21 up distribution along water- line and keel Fig 21 shows the calculated wall shear velocity distributions along water line and keel and in Fig 22 the girthwise distributions at three transverse cross sections are given along with the measurements. It can be seen from the figures that in the stern region the calculations overpredict the wall shear velocity along the water line, but underpredict it along the keel. This indicates that the calculated velocity profiles in these regions might differ somewhat from the measured profiles. It is difficult to make comparisons of velocity profiles, owing to the different coordinate systems used in calculations and experiments, as mentioned above. SSPA720 Ut AT 2X/L - O . 70 g o to _~ \ lo/ ~ htE AS UREME NT \_/ 60X2IXt5 n, 00 O .05 0. 10 0 . 15 O. ZO BIRTH LEN6TH Fig 20 C distribution at three stations . 00 0 .05 0 . 10 0. 15 GIRTH LEN6TH 2~ 0.20
SSPA7ZO Ut AT ZX/L - O.BO to o o \ O MEASUREMENT 60XZiXt5 :~ ~ ~ O MEASUREMENT 0 60XZiXt5 to 0 I , I 0.00 0.05 0.10 0. t5 0.20 6IR1H LENGTH 172~ ~ AT 2Xh _ n P' Fig 23 s~ BOX2IX15 0.00 0.05 0.10 0.15 O. GIRTH LENGTH Fig 22 up distribution at three stations Wake contours are given for two stations in Fig 23. A reasonable correspondence with measurements i. noted at 2x/L = 0.8, while at 0.9 the outermost iso-wake is displaced outwards. Close to the keel the boundary layer thickness is over- predicted at both stations. The results were obtained with the wall law as the inner boundary condition. Preliminary results from the latest version of NASPAR, where the wall law is removed, indicate that the prediciton of the flow near the keel can be much im- proved in this way. Measurements Colculotion Iso-wakes measurements x=O.R u/U=0.6, 0.7, 0.8, 0.9, 0.95 calculation x=0.81 q/U=0.6, 0.7, 0.8, 0.9, 0.95 Measurements Colculation :~ / Iso-wakes measurements x= 0.9, u/U=0.6, 0.7, 0.8, 0.9, 0,95, 1.0 calculation x=0.&9, q/U=0.6, 0.7, 0.8, 0.9, 0.95 Iso-wakes at two stations 4 FINAL COMMENTS Methods for computing free surface potential flows as well as complex viscous flows have been presented. The development has taken place over a long period of time and includes five PhD projects. Several computer programs have been written, but they are now combined into one program for the potential flow and one for the Navier- Stokes flow. These flow programs, together with a boundary layer code (for laminar and turbulent flows as well as transition) represent a powerful tool, named SHIPFLOW, for investigating the flow and resistance properties of ships. The SHIPFLOW system is also equipped with state-of- the-art pre- and postprocessors, which makes it easy to use. For preliminary project investigations the system represents an alternative to tank testing, i.e. it could be used as a numerical towing tank. It should be noted that there are several features of the programs which have not been presented in this survey. Such features are the transom stern option and the sinkage and trim calcu- lation in the potential flow as well as the propeller effect in the Navier- Stokes code. Quite interesting is also the most recent development, in which a mathematical optimization scheme is linked to the flow programs, enabling optimization of hull forms from a resistance point of view to be carried out automatically. Owing to space limi- tations these features have had to be left out, but they are explained in the references [10], Ill] and Elf]. 207
run no panel arrangement (half range) hull free surface Table 1. Grid dependence study. Wigley hull. extent of the free surface F (in half ship length) n x-direction y-direction 1 24 x 6 34 x 5 -1.50 to 1.50 0.28 0.266 0.450 2 24 x 6 3 24 x 6 34 x 5 -1.50 to 1.50 0.775 34 x 10 -1.50 to 1.50 0.775 4 24 x 6 34 x 15 -1.50 to 1.50 0.775 5 24 x 6 36 x 10 -1.75 to 1.75 0.775 6 22 x 6 36 x 10 -1.50 to 1.50 0.775 7 22 x 6 32 x 10 -2.00 to 2.00 0.775 8 22 x 6 40 x 10 -1.50 to 2.50 0.775 0.266 0.313 0.266 0.313 0.350 0.400 0.450 0.266 0.450 0.220 0.310 0.350 0.450 0.266 0.450 0.266 0.450 0.266 0.450 9 22 x 6 40 x 10 -2.00 to 2.00 0.775 0.266 References 1. Larsson, L. (editor), "The SSPA- 4. ITTC Workshop on Ship Boundary Layers 1980", SSPA Publication No. 90 (1981). 2. Larsson, L., & Chang, M.-S., "Numerical Viscous and Wave Resistance Calculations Including 5. Interaction", 13th Symposium on Naval Hydrodynamics, Tokyo (1980). Kim, K.-J. & Larsson, L., "Comparison Between First and Higher Order Methods for Computing the Boundary Layer and Viscous Resistance of Arbitrary Ship Hulls", International Symposium on Resistance and Powering Perform- ance, Shanghai (1989). v w x 103 1.204 0.598 1.437 0.785 1.495 1.212 1.849 3.054 0.847 3.054 0.603 1.551 1.199 3.046 0.859 2.788 0.931 3.317 0.931 3.662 0.913 (0.966) 3.168 (3.235) Broberg, L. & Larsson, L., "A Calculation Method for Ship Stern Flows Using and Analytic Body- Fitted Coordinate System", 15th Symposium on Naval Hydrodynamics, Hamburg (1984). Larsson, L. & Johans son, L.-E., "A Streamline Curvature Method for Computing the Flow near Ship Sterns", 14th Symposium on Naval Hydrodynamics, Ann Arbor (1982). Broberg, L., "Calculations of Partially Parabolic Stern Flows", SSPA Report 2801-2 (1988). 7. 208 Broberg, L., "Numerical Solution of the Time Averaged Navier- Stokes Equations for Ship Stern Flow", SSPA Report No. 2803-1 (1988).
8. Broberg, L., "Numerical Calcula- 18. tion of Ship Stern Flow", Ph.D. Thesis, Chalmers University of Technology (1988). 9. Zhang, D.H., Broberg, L. & Larsson, L., "A Numerical Method for Generating Smooth Body-Fitted Coordinate Systems Suitable for Ship Stern Flow Calculations", International Symposium on Resistance and Powering Perform- ance, Shanghai (1989). 10. Zhang, D.H., "Numerical Prediction of Propeller-Hull Interaction in Viscous Flow", SSPA Report No. 2963-1 (1989) . 11. Xia, F., "Calculation of Potential Flow with a Free Surface", SSPA Report No. 2912-1 (1984) . See also Report No. 65, Dept of Marine Hydrodynamics, Chalmers University of Technology, Gothenburg, Sweden (1984) . 12. Xia, F. & Larsson, L., "A Calculation Method for the Lifting Potential Flow around Yawed, Surface-Piercing 3-D Bodies", 16th Symposium on Naval Hydrodynamics, Berkeley ( 1986) . 13. Xia, F., "Numerical Calculation of Ship Flows with Special Emphasis on the Free Surface Potential Flow", Ph.D. Thesis, Chalmers University of Technology ( 1986) . 14. Xia, F., "A Study on the Numerical Solution of Fully Non-Linear Ship Wave Problems", SSPA Report No. 2912-3 (1986) . 15. Ni, S.Y., "A Higher Order Panel Method for Double Model Linearized Free Surface Potential Flows", SSPA Report No. 2912-5 (1987) . 16. Ni, S.Y., "A Method for Calculat- ing Non-Linear Free Surface Potential Flows Using Higher Order Panels", SSPA Report No. 2912-6 (1987). 17. Ni, S.Y., "Higher Order Panel Methods for Potential Flows with Linear or Non-Linear Free Surface Boundary Conditions", Ph.D. Thesis, Chalmers University of Technology (19 87) . 2 Ni, S.Y., Kim, K.-J., Xia, F. & Larsson, L., "A Higher Order Panel Method for Calculating Free Surface Potential Flows with Linear Surface Boundary Condi- tions", International Symposium on Resistance and Powering Perform- ance, Shanghai (1989). 19. Kim, K.-J., "Ship Flow Calcula- tions and Resistance Minimiza- tion", Ph.D. Thesis, Chalmers University of Technology (1989). 20. Hess, J.L., "A Higher Order Panel Method for Three-Dimensional Potential Flow", Douglas Report N 62269-77-C-0437 (1979). 21. Chen, H.C. & Patel, V.C., "Calculation of Trailing Edge, Stern and Wake Flows by a Time- Marching Solution of the Partially Parabolic Equations", Iowa Inst. of Hydraulic Research, Report No. 285 (1985). 22. Patankar, S.V., "Numerical Heat Transfer and Fluid Flows", McGraw Hill (1980). 23. Chen, H.C. & Patel, V.C., "Practical Turbulence Models for Complex Flows Including Separa- tion", AIAA 19th Fluid Dynamics, Plasma Dynamics and Lasers Conference, Honolulu (1987). 24. Markatos, N.C., "The computation of Thick Axisymmetric Boundary Layers and Wakes around Bodies of Revolution", Proc. Inst. Mech. Engineers, Vol 198C (1984). 25. Larsson, L., "Boundary Layer of Ships. Part III: An Experimental Investigation of the Turbulent Boundary Layer on a Ship Model", SSPA, Goteborg, Report No. 46 (1974). 26. Dawson, C.W.: "A Practical Computer Method for Solving Ship- Wave Problems", Proceedings of the Second International Conference on Numerical Ship Hydrodynamics, Berkeley (1977).
DISCUSSION by C.M. Lee This is one of the rare papers which cover both inviscid and viscous flows about ships. I believe that a rational approach to ship flow problems should be based on an intelligent combination of potential and viscous flow solutions. This paper apparently shows the efforts in that direction. The authors deserve congratulation for their efforts. Bulbous bows certainly can be cited as one of the useful inventions achieved by the ship hydrodynamic research community. However, we are still uncertain of its effects on flow around ship bows, particularly on the viscous flows. Does a bulbous bow really contributes to reducing the total resistance of a fat ship like a super tanker at low speeds? I think it is about the right time to examine the bow flow for ships with bulbous bow in detail to find the stagnation point on the bow as a function of ship speed in order to develop a correct model for computing ship flows more scientifically. Author's Reply We agree with Prof. Lee that a thorough investigation of the flow around the bulb of a fat hull would be a most interesting exercise. Ideally, such an investigation should include an experimental and a numerical part. The experiments should include detailed velocity measurement, such as the ones by Fry and Kim, presented at the 15th Symposium on Naval Hydrodynamics using an LDV, and flow visualization to find the stagnation point. The numerical part should include free surface as well as viscous flow calculations. For this detailed investigation the best approach would be to use a free surface Navier-Stokes method, like the one developed by Prof. Miyata and his co-workers. In such a case the wave breaking, which is likely to be important in this case, could also be considered. DISCUSSION by S. Ogiwara There are several ways in prediction of wave resistance by the panel method, namely, the pressure integral over a hull surface, computation of momentum change and so on. If you have experiences to compare the results between them, especially for the case of the series models with protruding bulb, would you comment about that? Author's Reply We have compared the two ways of computing the resistance for several hulls. Results are presented in reference [19] of the paper. For the thin Wigley hull the differences between the two methods are very small, while for the Series 60, CB=0.60 hull, differences of about 10% were noted. The momentum approach yielded smaller values than the higher order pressure integrations. The Series 60 results are in our experience typical for ships of moderate thickness. However, for the very full HSVA tanker we experienced the opposite trend with the momentum results higher than the ones from the pressure integration. The differences were again around 10% . DISCUSSION by F. Stern Fig.17 of your paper is very surprising and contradictory to our work on viscous- inviscid interaction, ie, our work indicates a much greater extent of the interaction. Is the inviscid-flow calculation for the bare or displacement body? If for the former, please comment on the lack of dependency in your solution to the location of the outer boundary even for outer-boundary locations near the boundary-layer edge. Author's Reply It should be noted that the distance to the outer edge of the grid is measured from the keel, where the boundary layer is very thin. Fig.A1 may be helpful. We claim in the paper that there is an interaction at 2x/L = 0.05 but not at 2x/L = 0.11. As appears from Fig.A1, the latter edge is about three times the maximum boundary layer thickness away from the hull. When reading the paper one might get the impression that the 2x/L = 0.11 grid extends only to the boundary layer edge, in which case the results would have been most questionable. \ Fig.A1 210 11\ N: \ | Boundory \// \\\\ ~: \ No os - \~1 ----~_ \, ~