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.
Computation of the Flow past Shiplike Hulls J. Piquet and M. Visonneau ENSM Nantes, France Abstract The present work first details the implementation problems of a numerical procedure for the solution of the Reynolds-averaged Navier-Stokes equations in boundary- fitted coordinates. The procedure is validated using well documented experiments on the HSVA tanker and the comparisons demonstrate the ability of the method to predict ship stern flows. I. Introduction For a long time, wind tunnel or towing tank experiments have been the own practical way to get information on ship performances because of the restrictions of inviscid or boundary layer models. With the development of computers, it becomes possible to obtain numerical information with an answer turnover and a level of details which are hardly possible with experiments. The prediction of shiplike afterbody flows, which is considered hereafter, is necessary for the knowledge of viscous resistance and, more importantly, for propeller and control surface design as both are immersed in regions where viscous effects are significant. Classically, such problems were handled solving first a potential problem with panel methods. The inviscid flow field on the body was used as a boundary condition for thin boundary layer equations the solution of which allowed the determination of viscous corrections. Unfortunately, it was demonstrated during the SSPA-ITTC workshop on ship boundary layers [1] that such a procedure failed in the stern region because of the rapid thickening of the boundary layer in this zone. The limited success of generalizations of thin boundary layer equations involving high order corrections, see e.g. [21~33~4], was subsequently demonstrated so that the tendency towards computing the full solution of the Navier-Stokes equations became stronger and stronger for Increased computer ressources became more and more available at continuously decreasing costs. However, significant obstacles to progress in the treatment of such threedimensional flows remain, the most important being the difficulty in describing with enough resolution the geometry of the physical domain in which the body is immersed, even if, as here, no free surface effect is considered. Another important difficulty is the need of reliable data to validate computations insofar as the comparisons between the experiments and the computations involve systematic interpolation procedures. Table 1 summarizes the main characteristics of the presently available methods which share several common 295 features as well as differences among which the following appear important. (i), used coordinates and grid generation techniques, adaptivity of the mesh; (ii), type of master equations and retained approximations -all of which depart from thin boundary layer assumptions-; (iii), turbulent models, how are they used and with what specific boundary conditions; (iv), discretization schemes pressure velocity coupling methods, iterative modes solution methods used for the linear systems; (v), acceleration means; (vi), initial and boundary conditions. All the methods (including this one) that have been published, but [51~61~7l, rest on the use of Reynolds averaged Navier-Stokes equations (hereafter called RANSE) presented in §2 and written in a body-fitted curvilinear coordinate system. Except works [83~91~101~1 1], only the independent variables are transformed from the physical coordinates x (horizontal along the axis of the body), y (horizontal, orthogonal to the axis of the body), z (vertical) to the curvilinear coordinates ~ = (\, A, () while the dependent variables are retained in the physical domain: the velocity components are the components in a Cartesian frame with the noticeable exception of t121~131~14l, where cylindrical velocity components are preferred. In order to allow more easily comparison with experiments, the coordinate planes \=const. are identified with x=const. planes. In the following (§2.2), ~ will refer to the coordinate away from the wall such that the equation of the wall is ~ = 1; ~ will refer to the girthwise coordinate. The other alternative [81~91~10~11] uses contravariant coordinates. This simplifies the continuity equation but introduces Christoffel symbols in the momentum equations, the full treatment of which is cumbersome (usually a lot of terms are omitted) and either time or storage consuming. Now, considering the master equations, ie. those which are discretized numerically, they usually involve a k- £ type of closure with a wall-function approach which avoids numerical troubles often associated to the sink- source terms in the £ equation while saving storage for the nodal unknowns. Some works depart from this trend by the use of either a classical mixing length model [81~91 or of a subgrid scale model [ 151. Considering more sophisticated Reynolds stress models does not seem to improve significantly the results [161. In the following, a k-£ model is retained but the wall function approach is avoided.~§2.4) Only the £ equation is bypassed close to the wall. For a significant increase in n~,mencal troubles and in
~ ~ o ~ ~ 2 · ~ c o o ° o ° ° °C - o 4J ~ - > o ~ c - _ ~ ~ ~ ~ _ ~ s ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ E - - - c O ~ ~ ~ 0 0 C ~ `~ C~ ° cn c c °, ° c c E ~ E ° E c 0 ~ r 2. ~ E 0 o E . w c E ° C c E c ~ C ° c a, ~ ~ ~n ~ ~ ~ ~ d, O > ~ ~ O > ~ ~ L} ~ ~ ~ .~ c ° ° 0 · O (' c c° c G 0 1~ ~ _ O e/ ~ ,,. O ~ ~ ~ O E E ~ O _ . . _ 2 _ _ O ~ O E ~, 0 0 ° ~ n e c · ~ ~ ~ O >. a~ i~ ~ ~ ~o ~ Q w I ~q t~ ~ .. _ u~ Ci) ~ a 0 J v ~, 0 `, tO Q > C ~ . b E E . E. O) a, ~ ~ ~ ~ ~ b ~ Y Y Y ~ I a ~n _ ~~, ~D n 296 w D - a, o 0 C) - - - C' S as w D a,
computing time, the delicate problem of the threedimensional orientation of the log law specification, dealt with only in t14], is avoided, as well as significant overestimations of the velocity profiles, close to the wall in the propeller disk region.~171~183. From the closed RANSE, the master equations are almost invariably obtained -except [143~193- using the so- called partially parabolic approximation in which longitudinal \-transport due to viscosity and turbulence is neglected (de, if ~ is the advected quantity, only the ¢) terms are dropped). This approximation is usually developped for viscous flow fields in which pressure is the dominant transmitter of influence in the upstream direction. In certain conditions, it does not forbid the presence of - seldomly present- longitudinal flow reversal (along \) which may occur close to transom sterns or when a propeller is in operation. In this work, the master equations do not involve such an approximation for reasons which will be made clear; therefore the fully elliptic equations are used. Turning now to the numerical aspects which are studied in §3, a general prerequisite is a method for generating the body-fitted gnds. The most common choice, Including ours already detailed in §2.3, rests on elliptic solvers [2~)] because they are easy to use and because they allow a priori controls of the quality of the grid close to the body. A possible alternative has been proposed in [81~9] [101~19] [21]where twodimensional ends, generated via conformal transformation designed, for instance, from the Schwarz-Christoffel theorem, are stacked in the ~ direction. Such grids build-in transverse orthogonality to the price of distorsion problems in the ~ direction, close to the aft sections. Other algebraic techniques have also been used [22] [23] [24] [25] [26] but they appear often abandonned. The most important aspect of the numerics lies in the way the incompressibility constraint is enforced. This may be done in three ways. A first category of methods includes those where the pressure-velocity coupling is simulated with a suitable modification of the continuity equation. Compressible flow methods used for low Mach numbers belong to this category, although severe difficulties are usually found in the limit of vanishing Mach numbers. The artificial compressibility method [27], systematically used in some institutes (see e.g. [283), and applied to ship flows in [293~30], gives also rise to convergence problems and needs a careful selection of numerical parameters such as the pseudocompressibility factor.The second category includes methods in which the continuity equation is satisfied identically at each iteration. Apart from unsegregated methods [31], such methods follow the boundary layer practice in which the flow domain is swept from upstream to downstream, implicit differencing of momentum being used for marching stability. Upstream influence through the pressure field has to be accounted for, and this is done introducing some form of forward differencing for the streamwise pressure gradient which allows departure-free behavior [32~. For ship applications, such methods have been mainly developped in the framework of the partially parabolic approximation in [81~91~101. In this work, we have followed (§3.2) the most common so-called pressure correction methods in which the pressure velocity coupling is solved iteratively, updating successively the velocity variables In the momentum equations and the pressure in a pressure equation, in such a way that solenoidality is satisfied at convergence. In the ship applications, both the advection and the diffusion of momentum are implicit in time; the updating procedure being handled on each equation by the well known "SIMPLE" t33l, "SIMPLER" [34], or "PISO" [35] methods. One of the dominant aspects of such techniques lies in the way the pressure equation Is handled as the pressure solver controls the global convergence of the method (and therefore, its robustness on fine grids), and as it generates the greatest part of the computational cost. The present work strongly differs from similar computations [121~131~14] by the characteristics of the pressure and pressure correction solvers which become critical as soon as grid clustering is strongly increased by the removal of the the wall function approach. Such modifications are necessary for the following reason. Because only the independent variables are transformed, the directions of the velocity components differ from the direction of the curvilinear coordinates. As a first consequence, the projection phase of the pressure- velocity coupling does not connect always the driving component of the pressure gradient to the corrected velocity component . Because of the strong grid clustering the stiffness of the "threedimensional" pressure matrix is so much increased that standard relaxation techniques usually retained do not converge anymore. This argument will be clarified in (§3.3) As a second consequence, the possibility of pressure odd-even oscillations is recovered even if the grid is staggered (the most common choice, see Table 1~. They are classically avoided, as when, as here, colocated grids are used, with weighted averages of coefficients defined during the projection phase. For the methods belonging to the third category, the discretization of the advective terms rests often on strongly diffusive classical hybrid schemes. In this work, we use (§2.1) the so-called finite-analytic method of [37] applied, for instance, in [133~143~151~181. This quick review of methods would stay incomplete without consideration of initial and boundary conditions. They will be specified in §44, when the results are studied. While the iterative procedure is usually initialized from rest, from uniform flow or from potential theory, the considered domain boundaries usually extend from a midship station to several lengths behind the body. [51~61~71~153~381~39] depart from this choice on model ships and start from a uniform flow field far upstream. In this work, we follow the common practice which considers that the inlet flow at midship stations can be characterized by some simple integral parameters and that the crossflow at midship can be neglected. These integral parameters are obtained either from experiments or from thin boundary layer assumptions used upstream. This allows to save computer ressources and to avoid transition problems in the fore part of ship models. The flow domain extends usually to some ship lengths downstream and to one ship length away from the axis of the ship. Symmetry conditions are retained on the vertical plane -the so-called keel plane- (only one side of the ship is handled) as well as on the horizontal plane -the water plane-. Only the so- called double-model problem is therefore treated, except in [51~61~73~151~38 1 where free surface effects are considered. Related to the flow domain choice, is the wake topology which involves, the most often, a singularity line either of the polar type or of the parabolic type. In the former case, transverse grids are topologically equivalent to concentric circles Art = const.) and their radii if= 297
coasts. In the latter case, the singular line starts from the bottom of the ship, in its vertical plane of symmetry into the wake, defining a line of separation between the upper part of this vertical plane of symmetry (the ''fictive rudder") which is a (symmetry) surface ~ =1 that extends the wall surface (rl=l) in the wake, and its bottom part which is a (symmetry) surface (=2. Among the available experimental information which has been extensively reviewed in t40][41], we shall mainly focus (§4) on the HSVA tanker [44]~451. For this case, very detailed mean flow measurements are provided in several x = const. planes; this minimizes interpolation processing. Moreover, the data are free from blockage effects because a wind tunnel with slotted walls was used. Uncertainties on the inlet conditions appear therefore to be an admissible penalty. The paper is outlined as follows. §2 presents the equations and the main features of the used method; §3 focusses on the numerical aspects; §4 discusses the results and presents significant comparisons with the HSVA experiments. 2. The Master Equations 2.1. The Primitive Form. The full RANSE (2.1 ) for the mean velocity field U and for the mean pressure deviation P V.U=O, DU/at + (U.V)U + VP = Re~1V2U - V.uu (2.1) need a closure assumption for the Reynolds stresses -uu which is taken under the isotropic form (2.2) uu = 2k 1/3 - vT [V U + VTU] (2.2) where VT = C~k2/e is the turbulent viscosity, k is the turbulent kinetic energy and ~ the so-called dissipation of turbulent kinetic energy. k and ~ are provided by the standard following transport equations (2.3) and (2.4) Ok/3t + U.Vk = G - ~ + V.[ Rk-1Vk ~ (2.3) 8e /Ot + U.Ve = Cal 1 ~ G/k - Cal 2/k +V IRE - 1 Ve ~ (2.4) G is the rate of production of turbulent kinetic energy uu: VU, Reff is the effective Reynolds number defined by Reff -1 = Re~1 + vT while Rk-1 = Re~1 + vT / ok and Rat -1 = Re~1 + vT/c;~ . Unless specified, model "constants" are given by their standard values: Cp =.09,C~ 1 = 1.44,C~ 2 = 1.92,c~k = 1,c,~ = 1.3 (2.5) 2.2 Physical space and computational space 2.2.1 The problem . The physical domain Q. where the flow is studied, is bounded by the ship hull B of surface DB, by the outer boundary £, by the inlet station U. by the outlet station D and by the vertical V and horizontal H planes of symmetry. Equations (2.1 to 4) which classically involve the Cartesian velocity components {U°C}-{U,V,W} are partially transformed from the Cartesian rectilinear coordinates {xa}-ix, y, z} in the physical space, to the curvilinear coordinates hi - {k,ll,`} in the so-called computational space. The coordinate transformation is designed in order that the boundaries OB and £ become ~ = const. surfaces, that the boundaries U and V become ~ = const. surfaces and that the planes of symmetry become ~ = const. surfaces. The domain Q 298 becomes a parallelepiped in the computational space in which the discretization consists in stacked unit cubes of sides AX = A~ = A: = 1. Each unit cube of the computational space is a curvilinear hexahedron in the physical space, the "sides" of which are measured by the modulii of the covariant vectors al = ar/aXi. 2.2.2.Useful Relationships. The transformation necessarily involves byproducts from the basis aj. Of particular interest are [201: (i) The area vectors bJ = ajxak (i,j,k in cyclic order) which measure the oriented area of a small surface of unit sides . . along ~J and ~ ~ on a \~ = const. surface in the computational space. bl appears as built with two small triangle-likes surfaces in the physical space. (ii) The Jacobian J of the transformation from the computational space of the coordinates {\i} to the physical space of the coordinates {xa). J measures the "physical" volume of a parallelepiped of unit sides in the computational space; this parallelepiped appears as an hexaedron-like volume in the physical space. (iii) The covariant and contravariant metric tensors .. . . gij = ai.aj; glJ = bl.bJ/detgij where the determinant of g is the square J2 of the foredefined Jacobian. An important relationship to be used is the following simple restatement (2.6) of the chain-rule derivative formula: a~k/8xa=J-1 teak (2.6) which allows the computation of the standard following operators: . . . . . V.U = J~ D(JUl)/~1 = JO a[balu(x]l3~l (2.7) where UP and [V°]a are the physical components in the Cartesian orthonormalized physical space. [V¢~]a = J~ 1 teak DO fink (2.8) V o = div [g akk ~ = J akl [ Jg 3\k ~ (2.9) 2.2.3.The Grid generation method. Because the values of the curvilinear coordinates hi are specified on am, the values of ki inside Q result from a boundary value problem in which the dependent variables are the {bit and the independent (unknown) variables are the {xa3, the correspondence between {x(X} and (\i} being one-to-one. If (2.9) is written for V2x = V2y = V2z = 0, the resulting partial differential equations (2.10) may be used to generate the coordinates: 2 V2`q, glk a ~ + 1 [ g ] ~ with Ax (x,y,z) (2 10) aX aX aX aX Because the coefficient of the first order derivative, in (2.10), is the laplacian of the curvilinear coordinate: k 2 k 1 a [Jg ] ik (2. 1 1 ) aX The common practice is to specify it as a grid clustering control fk = v2\k. In this case, (2.10) becomes (2.12)
vie, lk a ~ + fj at (2.12) afar ax Now, if v2xa is computed with (2.11) rather than with (2.10), fk must be equal to the rhs. of (2.12~. Therefore, when (2.12) is used to compute the diffusion terms (for ¢=U,V,W,k,~ ), the iterative method for solving the mesh must be solved to convergence. In all that follows, only the special case ~ = ~ (x) is considered so that fl= gllX~/x~__2agll t121~46] results from the assumed longitudinal grid distribution x(~. Thus, only the equations for y and z need to be solved, once control functions f2 and f3 have been specified. ~ i, .. f~andf~ are connected to giJ in such a way that 2b=-f2/g22 - Sr~/Sr~ and 2c=-f3/g33 - 0~/~; (as well as 2a) are frozen during the iterative use of (2.12) for ¢=y,z. They are prescribed according to (2.13) whore S is the curvilinear abscissa along the ~ = cte lines and ~ is the polare angle in the x = const. correction. S f =FA(11,`) for \<\A; FA= ~g S I OVA (2.13a) S f2=FB(q,`) for \2^B; FB= g22 S~'1 1 \=\B (2.13b) f =FC(\,11,`) for ~A<~<\B ; (\~ B FAR-~A)FB CF= ~ -\ A B a,b,c are fixed from an initial guess taylored to the problem (see §3~. The numerical solution of the grid equations uses exponential schemes similar to (3.6) -see infra- in both directions; it considers crossed second-order derivatives as source terms which are discretized, as well as metric components, with centered schemes. The solution of the resulting linear systems endly results from SLOR alternate ~ and ~ sweeps. The grid generation procedure thus defined has two aims: it first allows a correct clustering of grid points close to the wall, and, as important, it smoothes the grid. So, zero machine convergence of the procedure is not necessary. The f2 and f3 functions that result from (2.13) after SLOR sweeps are not used in the Navier-Stokes solver. Rather, the discrete equations corresponding to the Poisson operators are considered as a 2x2 linear system whose solution gives the needed f2 and f3. 2 3 The Flow Equations Apart from the continuity equation -see (2.7)-, the transport equations for the mean momentum, for k and ~ can be written under the following Master Equation (2.14) form: gl 1~:,~` + g22~r~ + g33q>~` = 2A¢~: + 2Bq>~, + 2Co'DX + R¢Ot + SO (2.14) where 2A`p, = Jeff l b3 as ~ - f3 2B~ = Jeff tb2 At ~ _ f2 2C<, = JEbl a¢] - Red = a¢Reff; S =s2 lgl241, + gi3~' + gads and a =a<, U-J Ebb VT ;~+b~ VT r,+b~ VT `] (2. 1 6a) tabalp, V-J l:b2 VT ;~+b2 VT T~+b2 VT `] (2. 1 6b) a =a~ W-J [b3 VT A+b3 VT rl+b3 vT,`] (2.1 6c) The index ~ refers to any of the advected quantities, namely U. V, W. k, e; the indices X, it, ~ refer to the partial derivatives; coefficient ad, is equal to 1 except ak = Cork, an = o~ . Source terms so for U. V, W. k, ~ are obtained from (2.17) written with Cartesian derivatives indexed x, y, z. SU=Reff [PX- gradVT. a au SV=Reff [py- gradvT · ] By (2.13c) au SW=Reff [pz- gradvT · I] (2.17a) (2. 17b) (2.17c) Sk=~C7kReff (Gee) (2.17d) is ~ e2) So ~C,~Reff Decal k G - Cur 2 k ) (2.17e) 299 2.4 The Turbulent Model Two reasons for the eviction of the wall function approach have been already mentionned in the introduction, namely the problem of specifying the directional dependence of Up and the practical overestimation of velocity profiles. A third reason is the noticed strong dependence of the pressure level on the location of the grid points close to the surface. As, at least the ~ equation is avoided close to the wall, it is necessary to specify the closure assumptions to supply. For the HSVA case (§3), a one equation is used in the inner zone in order to make the closure problem free from the estimation of Up (whose sign may vary in the domain). The closure assumptions (2.18, 2.19) follow the choices of [47] ~ = k3/2/1~ : vT = Cook 1~.~ (2.1 8a,b) where the scales of the two equations differ [481: 1p = C1 n ~ 1 - exp (-Ret/ k n / April; 1~ = C1 n ~ 1 - exp (-Reek n / An )] (2.19) The constant C1 results from log law assumptions, n is the normal distance to the wall, which give C1 = KC~,-3/4; the
asymptotic behavior of vT when n - O fixes An = 2C1; Ap = 70 allows to recover the 5.45 constant of the flat plate log law. The (straightforward) patching of the k inner model with the outer k-e model is located at Reek n _ 200-250. 3. The Numerics 3.1 The Transport Equations The numerical method for the treatment of (2.14) closely follows [123~14~. Coefficients of ~ derivatives are evaluated at the center P of the control volume and coordinates are normalized in the following way. >* = >/4gll A* = 1l/<g22; `* = `/<g33 ~3.1' so that the normalized form (3.2) is obtained: O~ *>,* + O11*~* + ** = 2A By* + 2B Oil* + 2C Go* + R Ot + (Skip (3.2) with: A = (A¢)p/4g33p; B = (g¢)p/4g22p; C = (C¢)p/<gl lp R = (Rasp for a parallelepipedic volume (in the computational space) whose sides are given by (3.3~: < * = 1 = 1/4g l l; ^~, * = k = 1/4g22; A`* = h = 1/4g33 Now, the equation (3.1) is splitted as follows [141: 2C Go* - O~*~* + R At + Sq, = Gig*, at*, (*, t) (3.4) O~*~* + O`*~* - 2A O`* - 2B Oil* = G (3.5) Assuming that G and SO are constant on each element, a backward time derivative reduces (3.4) and (3.5) to a one dimensional and a twodimensional advection-diffusion equation, respectively. The solution of (2.4) is written as: a [exp(2ck*3 - 1] + by* + c (3.6) Substituting (3.6) into (3.4), one is left with: g = [ 2CO~* - ~~*~ * + R At + SO ] p = (CU+CD)~P - Cud>U - CD~ + -lip ~pn-ly+ +(S<:,Yp with C exp(Cl) C exp(-Cl) Cu = l sinh C1 ; CD = 1 sinh C1 ~ is the time step, the index n-1 refers to the known state, indices U and D refer to upstream and downstream nodal values respectively. The resolution of (3.5) is classical and the result can be written under the form (3.8) (3.7) OP = ~ Cnb 4)nb CP g (3.8) The index nb refers to neighbouring nodes NE, NW SE, SW, NC, SC, EC, WC (fig.l). If the relationship (3.8) characterizes compact schemes (the molecule of which involves at most three points per direction), the values of influence coefficients Cnb, Cu. Cp depend on the type of discretization in the plane ~,(. Following [121~371~40], we use the finite analytic discretization. The resulting values of Cnb are omitted here. Combining (3.7) and (3.8) gives (3.9~. [ 1+ Cp (CU+CD)+ Cp ~ Dip] = (3 9) ~ R n-1 CnbOnb+CP CU~U+ CDOD+ _ OP -(S~)P 3.2. The Pressure Velocity Coupling. 3.2.1. The algebraic problem. The discrete solution of RANSE is written (E-A)V+GP=f; DV=g (3.10) where V and P are the vectors of nodal unknowns for the velocity and the pressure. Boundary conditions and known iterates from time tn-1 are gathered in f and g. D and G refer to the discrete form of the divergence and of the gradient operator which may involve interpolation procedures so that, in general, D ~ GT . E-A gathers the implicit part of the advection and diffusion terms in (3.9), E being the diagonal matrix of the coefficients C at node P. The concrete form of E and A is not important for the argument; let us only mention that, with the finite analytic method, diagonal dominance implies that the spectral radius of E-1A is less than one so that an iterative method which gives V from P through (3. 11) V= (E-A)-1 (f GP) (3.11) converges quickly. If (3.11) is substituted into (3.10), the pressure is obtained with (3.12) [D(E-A)-1G] P = D(E-A)-1f g (3.12) Pressure correction techniques rest on the use of the (3~3) approximate inverse E-1 for (E-A)1. 3 . 2. 2. The "P150 " Pro ced ure [35 ~ . It con si sts in the following steps, starting from a guessed solution (Vn~l, pn-l) from which (fin, En) is computed in the following way. (i) prediction step. The velocity predictor V* is first computed solving -step 4 - (3.13~: (E-A)V* = f Gpn-1 (3.13) V* does not satisfy the continuity equation (DV* ~ g) and must be corrected. Two correction steps are used to generate An, pny (id f rst correction step. The first corrected field (V**, P*) to be computed is such that: EV** - AV* + GP* = f; DV** = g: step I (3.14a,b)1 The momentum equation (3.14) specifies V** from V defined by (2. 15a) V = E l (AV* + f (3.15a) according to the following projection form - step 2 -: V** = V -lop* (3.l5b)1 while the pressure equation for P* is obtained from a substitution of (3.15b) into (3.14b) - step 3 -: DE-1G P* = DV - g (3.16~1 (iii) second (mandatory) correction step. The second corrected field (V***, P**) - (fin, fin) satisfies EV*** - AV** + GP** = f; DV*** = g (3.14a,b)2 Therefore, equations (3.15a,b)2 and (3.16~2 result from (3.15a,b)1 and (3.16~1if the following substitution is performed: V* - V**, V** - V***, P* - P**. In our implementation, where only the steady state is looked for, the second correction step is not performed so that vn=V*, pn= p* 300 3.3 Implementation problems. -Step l- In the cell-centered colocated grid approach, where velocity components pressure and turbulent quantities are defined at the center of the control volume the momentum equations (3.9) can be written under the form (3.17~. U- 1 1+CC(CU+CD) {Z CnbUnb+CC ~ CUUU +CDUd +(Ud -Ud)-SU ~ ~
1 +CC(CU+CD) {E CnbVnb+CC [CUvu +CDvd +(Vn -Vn) Sv ] } We- 1 1+CC(CU+CD) TZ Cnbwnb+cc [CuWu +CDwd +(We -We)-Sw ] } CC is the finite analytic coefficient computed at the colocated velocity mode (in the standard MAC method [49] Cd, for Ud, Cn, for Vn, Ce, for We are computed at the staggered velocity nodes). The matrix E is built with coefficients like l+Cc (Cu+CD), etc...The matrix A gathers terms like Cnb, CcOu, CCCD, ...Terms like R / ~ are in f. The G matrix for pressure contributions is inside source terms (Su), (Sv)... -Step 2- Equation (3.15b) can be written U U J D(blP'`+ UPS+ by) ~ _ V=V-TlD(b2PA+ b2P~+ b <;) V=W J-lD(biP + ~+ b3P ) w. 1ere (3.18) CcR 1+ CC(CU+ CD) The underlined terms in (3.18) are generated by the m~salignment between the coordinate lines and the directions of the velocity components. With respect to the staggered grid approach, it is seen that the D coefficient for the pressure gradient does not depend now on the correspond~ng velocity component. -Step 3 - Substituting (3.18) into the continuity equation (3.19): [b,U + b2V + b3W] + [b,2U + b22V + b32Wl~+ [b,3U + b2V + b33W] = 0 (3.19) gives the pressure f~eld (3.16) which can be viewed also as the solution of the continuous equation (3.20) [aij_ ] = div V; aij = DJgij = DT~, bk bk (3.20) ax ax k=1 Equation (3.20) is now discretized with a centered scheme. Usually, this gives a 27 node pressure molecule (9 in the upstream plane, 9 in the central plane, 9 in the downstream plane). Seven points: U. D, N. S. E, W. P come from second order derivatives 32 / axia~i, i=;; the other twenty points come from the crossed second - order derivatives 82 / 8kiD)J, i ~ j. Because of the retained discretization, the pressure molecule involves only 19 points (five points in U and D planes). Among the whole set of contributions to the pressure matrix, it is necessary to keep all the fluxes in order to allow a complete construction of the aij coeff~cients, i = j, because of the misalignment problem. Extradiagonal coefficients i ~ j arising from crossed second order derivatives are connected to the level of non orthogonality of the grid. They are not retained implicity but rather treated as source terms. As a result diagonal dominance and symmetry of the discretized pressure matrix is easily obtained. 18 velocity components hav~e to be specified in (3.20) for the evaluation of div V in each control volume div V = (b~U + b2V + b3W),d - (b~U + b2V + b3W)U + (b2U + b2V + b32W) - (bi2U + b2V + b3W)S + (blU + b2V + b3W)e - (b~U + b2V + b3W)w (3.21) Because of the definition of the control volume, the whole set of velocity components has to be reconstructed from the known centered velocity components. A simple interpolation is used for that purpose. Because of the structure of the retained discrete molecule, the discrete pressure matrix is symmetric. -Step 4- The computation of the non solenoidal V* f~eld uses (3.17), the pressure field resulting from step 4 being present in source terms. The correction procedure starts only when the field V* has been computed th~ughout the complete domain. 3.4 Solution of the linear Systems 3.4.1 The Problem Before discussing the technical solutions which have allowed correct solutions to be obtained, let us point out the origin of the difficulties. Pressure equations and pressure correction equations are almost invariably solved with SLOR type methods applied to the seven point molecule induced by the treatment of under lined terms in (3.18,19) as source terms. The main interest of this approach is, as already seen, that the P matrices are now symmetric and, usually, diagonally dominant. SLOR methods can therefore be expected to converge. Moreover the coding of the needed tridiagonal matrix algorithm is straightforward. Unfortunately, due to the problem of misalignment, dominant terms are omitted in the system and treated as source terms. As a consequence, global convergence of the metrod is made more diff~cult. Unless strong departure from orthogonality is admitted for the grid, this is not dramatic when wall functions are used because grid clustering in the direction ~ remains rather weak: the first points away from the wall being located so that 50 < y+ < 200. When wall fuctions are avoided, strong grid clustering implies a dramatic increase of the pressure matrix stiffness so that, for threedimensional problems, divergence of the iterative procedure occurs, mainly through a too important pressure correction in the wake. If the pressure matrix is enriched, [18], (using e. g. nine points at the current station ~ = ki) it becomes more difficult to invert. Thus, improving at least the pressure solver is necessary. 3.4.2 The pressure solver It has been found [50] that convergence was systematically achieved in practice with an incomplete LU preconditionned biconjugate gradient (PBCG method). Conjugate gradient techniques are used as acceleration methods in which the elements of the pressure matrix (hereafter called A) are generated at each iteration, without any need of storing the full matrix. The CG method satisfies a minimisation property so that a monotone decrease of the error is possible as well as a convergence disruption to which is associated a monotonic decrease of the error. This is not verified in practice because it is not the full pressure matrix that is solved. Because the spectrum of the A matrix is continuously distributed preconditionning is necessary as a means of improving the convergence of the BCG method: the system to solve Ax = b is replaced by the equivalent system By = c where K(B) << K(A), where K iS the condition number of the 301
matrix, so that eigenvalues of B are more clustered thaw those of A. Wha^; is needed is an approximate inverse A for A such that AAx = Ab or By = c where B = AA; c = fib. B being "not to far" from the identity, the CG method for B converges quickly. The A matrix is selected so as to preserve the sparsity of A. Here, an incomplete LU decomposition [51] [52] -no fill in- is used - the so-called LU (1, 2) decomposition performed worse than the retained classical LU (1, 1) decomposition - . The solution of the LU systems involves recursions which cannot be easily vectorized so that the solutions of these triangular systems is time consuming. Fortunately, a Neumann series method applied to the L and U systems has been found to converge nicely. The argument is presented for L. If L is written under its Jacobi splitting L = DL + EL, L = DL [1 - 1 ] where 1 = DL-1EL. It appears in practice that 1 is a convergent matrix. A classical theorem of linear algebra then indicates that 1-1 is regular and (1-1~-1 = 1+1+12+... The Neuman serie is truncated to m terms and the resulting approximate inverse of 1-1 can be computed with the Homer scheme. The whole set of computations can be vectorized with a speedup ratio higher than 30. Further information on the PBCG algorithm are given in [511. 4. The Results. HSVA Tanker. The flow domain covers 0 < X < 6 where X is adimensionalized with L, the ship being from X = -1 to X = 1.05; rs < r < L where r is the radial distance from the axis of the ship y = z = 0 and rs is the radial location of the hull. The flow domain is discretized with 80, 40, 31 points in the axial, radial and girthwise directions. Starting from an a-priori specified surface grid distribution, a volumic grid is firstly generated using a transfinite interpolation procedure [203~55] which is modified in order to allow clustering of points in areas of strong concavity. The resulting grid is irregular as the discontinuities on the body surface are propagated throughout the domain. The elliptic grid generation method is used to smooth out the grid irregularities. Fig.3 presents a perspective view of the grid which appears to behave correctly. Inlet conditions, written at X = 0, are estimated in such a way that results at X = .291 roughly match with experiments [11. The main Indetermination lies in the fact that specified data [13~44] involve 81 1, 61, Cf and H12, where: tan o ^00 r-- Ue-U U d U-U [ Qe ~ Qe 1 | [ Q J o Hl2=8llI6l ;Cf=tW/ 2pU2 from which needed values for 6, Up and Qe are estimated (,Bw - O). These estimations allow the generation of inlet velocity profiles by the method of Coles & Thompson t531. It is therefore necessary to briefly mention the results of the sensitivity analysis of 6, Us and Qe data. It is found that, if the influence Up is very weak, the result of the computations depend on the level of the overshoot Of Qe with respect to 1, and of the location ~ of this overshoot. We consider that this is the most serious weakness of this approach and this justifies recents attempts [6] [29] t30] [39] to compute the flow past a complete full shape, starting from far upstream. Apart from difficulties connected to the transition specification -a problem on ship models, not on full scale ships-, the flow becomes then somewhat sensitive to the level of the turbulent viscosity outside the viscous zone so that consideration of these aspects is left for a future study. Other boundary conditions are either wall-type conditions or symmetry conditions. Neumann boundary conditions are used for the pressure, except for r = L where P is taken equal to zero. 120 global iteration are usedwith a local time stepping procedure. Local time step is divised to ensure a fixed amount of diagonal dominance will respect to discretized momentum equations. Fig.4 presents longitudinal pressure distributions. A sensible grid effect in the longitudinal direction is present in this case. Girthwise distributions agree reasonably well with experiments although some more resolution should be necessary (fight. Comparison of secondary velocities with experiments [44] is given in fig.6. The convergence of streamlines, which starts to develop at midship, increases more downstream and gives rise to a counterclockwise primary vortex in the hollow of the stern. The primary vortex is close to the keel plane both in the computation and in the experiments which indicate a "blind" zone for AL > .95, where the longitudinal velocity component is too small to allow significant measurements. The longitudinal vortex where the hull terminates is also found to agree with experiments. If isowakes are now considered (fig. 7) -they are defined as isovalues of the component U in crossplanes-, it is found that the experimental trends are correctly captured (except, of course, for string effects which are not present in the computation, see e.g. x = 38 mm, O mm, 100 mm, in [44] data-~. It must be mentionned that a bilinear interpolation is used to specify isowakes so that it is misleading to compare isowakes .95, 1, or 1.05 as corresponding velocity profiles are usually very "flat" close to U = 1. The apparent overestimation of the viscous thickness close to the waterline is probably due to the fact that the inlet viscous thickness does not depend on girth, while it should thicken from the keelplane to the waterplane. Velocity components U. V, W as well as pressure data can be compared more extensively (fig. 8a to c). U and V are presented on the upper left and right sides, respectively, W and P are presented on the lower left and right sides, respectively. For each series of plots, the evolution is considered with respect to y at a given station X = const., for several depths z = const.. In general, it is seen that the calculations are in correct agreement with the data, especially on V and W. except, for U. in small regions close to the hull or in the neighborhood of the wake centerplane. This discrepancy is believed to be due to the hull shape modifications. Due to the large amount of grid points which would have been necessary to account for the hub support boss located for -19 mm. < x < 0 mm. (in [401 data), this detail of the geometry has been eliminated (fig. 9~. Also, coordinate lines (L): ~ = 1, (= const. evolve continuously from the ship hull into the vertical plane V of symmetry of the wake. They define a family of curves that describe the hull surface for ~ < kw (X) where kw (X) defines parametrically the intersection of the hull-surface DB with V; this intersection is located for .95 < X < 1.05. For ~ < kw(X), lines (L) are almost rectilinear in V, ie. zeal = 1, ~ given) depends only weakly on \. In V, the line ~ < kw(X) usually 302
consists of a vertical line located at X = .95 for Zb < Z < Zr. For Zr < z < 0, the line ~ ~ kw(X) (which is a piecewise linear function of X) is rather close to a straight line, the equation of which can be written: Z/Zr = 1 +10 [.95 - X ~ so that the hull has disappeared in the cross plane X = 1.05. Now, the approximate hull profile in V is not described exactly by ~ < kw(X) but rather by a "stair" function (fig.9) consisting of a succession of lines X=const., each of them separated by aline close to z = const.. The no-slip condition is evidently written on the resulting "stair surface" DB s which is a crude approximation of the "longitudinal geometry" while crossplanes are correctly body-fitted. As a result, the description of the viscous layer is rather poor. Nevertheless, a comparison with [14] -where a similar treatment is present- indicates that the slight improvement obtained on the secondary flow can be attributed to the fact that wall factions are not used here. The discrepancy between the computed and measured U components is important, especially in the wake. Apart from the fact that this discrepancy is already present upstream, another possible explanation could be the incorrect account of the loss of no-slip which results from the turbulent model. Fig. 10 presents the comparisons between the computed and measured longitudinal vorticity component isolevels . In spite of the uncertainties involved in processing both the measured and computed data, such a comparison allows examination of differences of velocity gradients rather than velocity themselves. This test is therefore more stringent. As a result, it is found that the computations reproduce both the general features of the contours and the magnitudes. It is also evident that the diffusivity which is present in the viscous region is more correct, with respect to experiments, than that which would result from the use of wall functions [141. The differences in the shape of the isovalues, in the wake, could be foreseen from the differences in the velocity distributions. Fig. 1 1 presents a color picture of lines of wall friction, the hull is coloured with the pressure field, high Cp and low Cp being associated to red and blue colours respectively. The flow which looks clearly twodimensional at midship enters an oblique low pressure region which deviates the wall streamlines. The resulting convergence of these friction lines implies a thickening of the viscous layer, close to the aft, that makes thin boundary layer methods breakdown [11. A vertical convergence separation line seems also present, close to ~ = kw (=.95~. Due to the approximate characteristics of the hull geometry and to the fact that the U velocity component is very low here, the materiality of this result remains to be confirmed from experiments. The maximum pressure levels are found on the hull X ~ 1-1.025. More downstream, the vertical structure of the flow is confirmed: streamlines associated to the secondary velocity components V and W are drawn for two crossections X = const. S. Con cl usi on A fully elliptic numerical method for the solution of the full Reynolds-averaged Navier-Stokes equations has been applied to the computation of the flow past a shiplike hull. The method uses a system of numerically generated curvilinear coordinates and retains the Cartesian velocity components as independent variables. The turbulent closure of the equations is handled with the well-known 303 k-e model. The use of wall functions is avoided; instead, a one equation model is used close to the wall. Bypassing the wall functions implies a very high grid resolution in the direction normal to the wall as well as high aspect ratios for the discrete molecules. Due to the resulting increased stiffness of the pressure matrix, it has been found necessary to improve the pressure solver in order to get rid of convergence problems. The numerical method has been described and evaluated on one ship hull for which experimental data are available. From the presented results the following conclusions can be drawn. With respect to the modelling of physical flow characteristics, the main features of ship stern flows appear correctly captured. In particular, the viscous inviscid interaction gives rise to a pressure field in good agreement with experiments. The method allows also a correct description of the mean velocity flowfield in the thin boundary layer, provided suitable initial conditions are available [51~. In the thick region, close to the aft, where thin boundary layer methods systematically breakdown, the overall features of the flow are also captured: secondary velocities and even longitudinal vorticity components are correctly predicted. Where discrepancies are found, they may be attributed to the technique which has been used to approximate geometrical details in the framework of a monoblock structured grid. Improvements in the description of the flow close to the propeller disk region therefore require a more detailed resolution of the geometry. With respect to the methodology, ship stern flows do not appear to involve too strong a viscous-pressure interaction. Therefore, the segregated approach seems well adapted. Apart from geometrical grid aspects improvements could result from one of the following three points. (i) The convection diffusion scheme. The finite analytic scheme only roughly accounts for the correct physical imbalances: pressure gradients vs. diffusion in the near field, pressure gradient vs. convection in the far field, for all the possible convection velocity directions. More particularly, the hypothesis of local uniformity of the influence coefficients and of the source terms, which allows locally (finite) analytic solutions to be computed, is probably too restrictive, especially in the direction away from the wall, where the convection velocity and the source terms vary rapidly, owing to strong gradients of turbulent quantities. (ii) The Poisson solver. This is the most time-consuming part of the method, especially on threedimensional grids where the aspect ratios can be very high. The conjugate gradient technique appears very well adapted to the problem in spite of its induced storage increase. It can be noticed that it may be somewhat unreasonable to compute the pressure -which is regular and does not vary significantly close to the walls- on a grid which appears rather well suited for the velocity distributions and for the turbulent quantities in that they especially vary close to the walls. Because less points are needed for the pressure, the use of a multigrid procedure should efficiently cut down cpu times. Because relaxation is a local procedure, a multigrid technique might be also a way of improving the solution in the boundary-unfitted region close to the aft. (iiiJ The Turbulent Model is also a potential source of uncertainty for the description of near-wall flows where strong curvature effects are present. Although the flows considered are mainly pressure-controlled, improvements in the turbulence models might lead at least to improved length scale predictions (the turbulent kinetic energy is
known to be overpredicted in thick boundary layers [40]). Nevertheless, it is considered that improved results will be rather an outcome of a better treatment of the geometry and of a better convection-diffusion algorithm. Acla~owledgrnents: Authors gratefully acknowledge financial support of DRET through contract 86 / 104. Cpu on Cray2 was attributed by the Scientific Committee of CCVR. Cpu on VP200 was attributed by DS / SPI. The presented method has been developped from an earlier code version of the method [121. preferences [1] L.Larsson. (Ed.) Proceedings of SSPA-IITC on Ship Boundary Layers. SSPA Report N°90 (1980) [2] L. Larsson and M.S. Chang, Numerical Viscous and Wave Resistance Calculations Including Interaction. Proc. 13th ONR Symp. Naval Hydrodynamics, Tokyo, (1980) t3] T. Nagamatsu, Calculation of Ship Viscous Resistance and its Application. Journ. Soc. Naval Arch. Japan 157, 47 (1985) - See also Proc. 2nd lnt Symp. Viscous Flow. SSPA-Goteborg, (1985) t4] S. Soejima, Calculation of Three-dimensional Boundary Layers Around Ship Hull Forms. Proc. 2nd. Symp. Num. Phys. Aspects of Aerod. Flows, Long Beach, CA. (1985) [5] H. Miyata and S. Nishimura, Finite Difference Simulation of Non linear Ship Waves. J. Fluid Mech. 157, 327-357 (1985) [6] H. Miyata, S. Nishimura and A. Masuko, Finite Difference Simulation for Non linear Waves Generated by ships of Arbitrary Threedimensional Configuration. J. Comps Phys. 60, N°3, 391-396 (1985) [7] H. Miyata, H. Kajitani, S. Akifuji, M. Baba, M.Kanal, Wave Formation about a Ship Bow advancing in Head Seas. Spring Meeting Soc. Nav. Arch. Japan (1987) [8] M. Hoekstra, H.C. Raven, Application of a Parabolized Navier-Stokes Solution System to Ship Stern Flow Computation. Proc. Osaka Int. Colloq. on Ships Viscous Flow, 125-142 (1985) [9] M. Hoekstra and H.C. Raven, Ship Boundary Layer and Wake Calculation with a Parabolized Navier- Stokes Solution System. Proc. 4th Int. Conf. Num. Ship Hydrodynamics, Washington D.C. 470-491 (1985) [10] H.C. Raven and M. Hoekstra, A Parabolized Navier-Stokes Solution Method for Ship Stern FLow Calculations. Proc. 2nd Int. Symp. Ship Viscous Resistance, SSPA Goteborg (1985) t11] L. Broberg, Numencal Calculation of Ship Stern Flow. PhD Thesis, Univ. Of Chalmers (1988~; also SSPA Rept. 2801-2 & 2803-1 [12] H.C. Chen and V.C. Patel, Calculation of Trailing-Edge, Stern and Wake Flows by a time-Marching Solution of the Partially Parabolic Equations. IIHR Rept N°285 Iowa City, la (1985) [ 13] H.C. Chen and V.C. Patel, Numerical Solutions of the Flow over the Stern and in the Wake of Ship Hulls. Proc. 4th. Num. Conf. on Numerical Ship Hydrodynamics, Washington D.C. 492-511 (1985) [14] V.C. Patel, H.C. Chen and S. Ju, Ship Stern and Wake Flows: Solutions of the Fully-Elliptic Reynolds-Averaged Navier-Stokes Equations and Compansons with Expenments. IIHR Rept. 323 Iowa City,la(1988) [15] H. Miyata, T. Sato and N. Baba, Finite Difference Solution of a viscous Flow `, i.h Free-surface Wave about an Advancing Ship. J. Comp. Phys. 72, 393-421 (1987) [16] G. Tzabiras, On the Calculation of the 3D Reynods Stress Tensor by Two Algorithms. Proc. 2nd. Int. Symp. Ship Viscous Flows, Goteborg, Pap. 15 (1985) [17] J. Piquet, P. Queutey and M. Visonneau, Computation of Viscous Flows past Axisymmetric Bodies with and without Propeller in Operation. Numerical Methods in Laminar and Turbulent Flows (C. Taylor, W.G. Habashi, M.M. Hafez, Eds.) 5, Part 1, Pineridge Press, 644-655 (1987) [ 18] J. Piquet and M. Visonneau, Steady Threedimensional Viscous Flow Past a Shiplike Hull. Proc. 7th. GAMM Conf. Num. Methods in Fluid Dynamics. Notes on Numerical Fluid Dynamics, (M. Deville, Ed.) 20, Vieweg, 294-301 (1988) [19] G.D. Tzabiras, Numerical and Experimental investigation of the Flow Field at the stern of Double Ship Hulls. PhD Thesis, Nat. Tech. Univ. Athens (1984) [20] J.F. Thompson, Z.U.A. Warsl and C.W. Mastin, Numerical Grid Generation, Foundations and Applications. North Holland Publ. New York (1985) [21] G. Tzabiras and T.A. Loukakis, A Method for Predicting the Flow Around the Stern of Double Ship Hulls. Int. Ship. Progress N°345 (1983) [22 ~ A.M. Abd El Meguid, N.C. Markatos, D.B. Spalding, K. Muraoka, A Method of Predicting three-dimensional Turbulent Flows around Ships'Hulls. Proc. Ist. Int. Symp. Ship Viscous Resistance, SSPA Goteborg (1978) [23] A.M. Abd El Meguid, N.C. Markatos, K. Muraoka, D.B. Spalding, A comparison between the Parabolic and Partially-parabolic Solution Procedures for Three-dimensional Turbulent Flows around Ship's Hulls. App. Math. Modelling 3,249 (1980) [24] C.E. Jansson and L. Larsson, Ship Flow Calculations using the Phoenics Computer Code. Proc. 2nd. Int. Symp. on Ship Viscous Resistance SSPA, Goteborg (1985) [25] K. Muraoka, Examination of a two-Equation Model of Turbulence for Calculating the Viscous Flow around Ships. J. Society Naval Arch. Japan, 147, 35 (1980) [26] K. Muraoka, Calculation of thick Boundary Layer and Wake of Ships by a Partially Parabolic Method. Proc. 13th. ONR. Symp. Naval Hydrodynamics Tokyo 601-616 (1982) [27] R. Peyret and T.D. Taylor, Computational Methods for Fluid FLow. Springer-Veriag Ed., New- York. (Senes in Comp. Phys.) [28] D.C. Kwak, J.L. Chang, S.P. Shanks, S. Chakravarthy. An incompressible Navier-Stokes Flow Solver in Threedimensional Curvilinear Coordiante System Using Primitive Variables. AIA,4 J. 24, 390-396 (1986) [29 1 Y. Kodama, Computation of 3-D Incompressible Navier-Stokes Equations for Flow Around a Ship Hull Using an Implicit Factored Method. Proc. Osaka Int. Symp. on Ship Viscous Flow, 109-124 [30] Y. Kodama, Computation of High Reynolds Number Flow past a Ship Hull using the IAF Scheme. Proc. Soc. Naval Arch. Japan, Spring Meeting (1987) [31 ] S.P.Vanka, Computations of Turbulent Recirculation Flows with Fully Coupled Solution of Momentum and Continuity Equations. Argonne National Lab. Rept ANL-83-74 (1983) [32] S.G. Rubin, Incompressible Navier-Stokes and Parabolized Navier-Stokes Formulations and Computational Techniques. Computational Methods in Viscous Flows, 3; Series: Recent Adv. in Num. Meth. Fluids, Habashi, Ed., Pinendge Press (1985) 304
[33] S.V. Patankar and D.B. Spalding, A Calculation Procedure for heat, Mass and Momentum Transfer in Three-Dimensional Parabolic Flows. Int. J. Heat Mass Transfert 15, 1787 (1972) t34] S.V. Patankar, Numerical Heat Transfer and Fluid Flow, Hemisphere Publishing Corp. New-York. [35] R.I.Issa, Solution of the Implicity Discretized Fluid Flows Equations by Operator-Splitting, J. Comp. Phys., 62, N°1, 40-65 (1986) [36] M. Israeli and A. Lin, Iterative numerical solutions and Boundary Conditions for the Parabolized Navier-Stokes. Comp. & Fluids 13, N°4, 397-410 (1985) [37] C.J. Chen and H.C.Chen, Finite Analytic Numerical Method for Unsteady Two-Dimensional Navier- Stokes Equations. J. Comp. Phys. 53, N°2, 209-226 (1984) t383 T. Hino, Numerical Simulation of a Viscous Flow around a Ship Model, Proc. Soc. Naval Arch. Japan, Spring Meeting (1987) [39] A. Masuko, Y. Shirose and S. Ishida Numerical Simulations of the Viscous Flow around Ships including Bilge Vortices, Proc. 1 7th. Symp. Naval Hydrodynamics, Den Haag (1988) [40] V.C. Patel, Some Aspect of Thick Three- Dimensional Boundary Layers, Proc. 14th. ONR Symp. Naval Hydro. Ann Arbor, 999-1040 (1982) [41] V.C. Patel, Ship Stern and Wake Flows: Status of Experiment and Theory, Proc. 17th. Symp. Naval Hydrodynamics, Denn Haag (1988) [42] L. Larsson, Boundary Layers on Ships, PhD Thesis, Univ. Chalmers (1975) [43] L. Lofdahl and L. Larsson, Turbulence Measurements Near the Stern of a Ship Model, J. Ship Research 28, N°3 (1984) f44] K. Wieghardt and J. Kux, Nomineller Nachstrom auf Grund von Windkanalversuchen, Jahrb. Der Schiffbautechnischen Gesellschaft (STGJ, 303-318, (1980) [45] K. Wieghardt, Zur Kinematik einer Nachlaufstromung, Z. Flugwiss. Weltraumforsch. 7, Heft 3 (1983) / - Fig.lc Fig. la - Physical Domain Fig.lb - Computational Box Fig.lc - Transverse plane in the Computational Domain. [46] H.C. Chen and V.C Patel, Calculation of Stern Flows by a Time Marching Solution of the Partially Patabolic Equations, Proc. 15 th. ONR Symp. Naval Hydrodynamics, Hamburg (1984) [47] H.C. Chen and V.C. Patel, Near-Wall Turbulence Models for Complex Flows including Separation, AIAA Paper 87-1300 [48] M. Wolfshtein, The Velocity and Temperature Distribution in One-Dimensional Flow with Turbulent Augmentation and Pressure Gradient, Int. J. Heat Mass Transfert 12, 301-318 (1969) [49] F.H. Harlow and J.R. Welsh, Numer-ical Calculation of Time Dependent Viscous Incompresible FLows of Fluid with a Free Surface, Phys. Fluids 8 2182-2189 (1965~. See also Los Alamos Rept. LA-3425 (1966) [50] J. Piquet, P. Queutey and M. Visonneau, COmputation of the Threedimensional Wake of a Shiplike Body, Proc. 11th. Int. Conf. Num. Meth. Fl. Dyn., Lect. Notes Phys. Springer-Verlag (1989) [51] J. A Meijerink and H.A. Van Der Vorst. An Iterative Solution Method for Linear Systems of which the coefficient Matrix is a Symmetric M-Matrix, Maths. Comp. 31, N°137, 148-162 (1977) [52] J.A. Meijerink and H.A. Van Der Vorst., Guidelines for the usage of Incomplete Decompositions in Solving Sets of Linear Equations as they occur in Practical Problems, J. Comp. Phys. 44 (1981) [53] T. Cebeci, K.C. Chang and K. Kaups, A Three Dimensional General Method for Threedimensional Laminar and Turbulent Boundary Layers on Ship Hulls, Proc. 12 th. ONR Symp. Naval Hydrodynamics, Washington (1978) t54] K.H. Mori & N. Ito, Near-Wake Computations by Solving the Vorticity Transport Equation in a Body- fitted Coordinate System, Proc 4th. Int. Conf. Num. Ship Hydro., Washington D.C., 512-528 (1985,) [55] L.E. Erikkson, Generation of Boundary Conforming grids around Wing-Body Configurations using Transfinite Interpolation, AIAA Journ., 20, N°10 1982. 305 NW tNC NE _ WC = r E C 5 _ _ , ~ . , SWI ISC h 15£ 1-' -1' - - -1
N t n ., (i) he'd 1._ __~ ~ ' ' : 1 _/ ' 1 1 1 ~ ,l,~ _~. _: ~ D In it T . ~ ~ _ ~ __l__,L__ _ it 1 {/ I.~/1 ~ Fig.2 - Colocated grid. a, location of velocity components and of pressure. cat \ ~1 n .5 .S .7 .3 .9 1.~ I.t 1.: 1.3 1.' 1.; 1.i 1.7 1.] t.9 2.0 2.1 2.2 2.3 2.' 2.5 !~' Fig.3 - Perspective view of the grid surface j = 1; 11 = 0. ~ \\ p F :~ Fig.4 - Longitudinal Pressure Distributions.c,, experiments , computations. ,,. . 1 I,. . 111 ,,. Fig.5 - Girthwise Pressure Distributions. , experiments; , computations. 306
Al . - lI, ·: i: - ~ al~ ~J 11 . 1 l l ~ l ~ l 0 .04 .o~ .12 .16 .20 .20 x or (o no. A. co N (O LD CO tO O o C' .n a; _ a.... ( t !~ ~i O . I Fig.S - Girthwise Pressure Distributions. Q. experiments; , computations. U W =2 o N ~0 ~1 O- o ~1 ~ : _ _ _ · 0 .04 .08 .12 .16 .20 .24 x tO I (~ 1 ~r p O- . . , · . 0 .04 .08 .12 .16 .20 .24 Fig.8a - Velocity components and Pressure at a given crosssection X/L = .944, for several depths Z = const., as a function of Y.a, experiments; , computations. 307
Fig.8b - W 1 1 . 1 1 1 1 1 0 .04 .08 .12 .16 .20 .24 0~ _- ~ . ._ ~o a, - Fig.8b - U ~ ~ I ~ ~ ~ ~ ~ ~ ~ ~ 0 . Oq . 08 . 1 2 . 1 6 . 20 . 24 0 . 04 . 08 . 1 2 . 1 6 20 . 24 Fig.8b - P ~ , , , , , , 1 0 .04 .08 .12 .16 .20 .24 0 .04 . 08 . 1 2 . 1 ~ .20 . 24 o ' a) o ~D (o rn r~ _~. Fig.8c - P . o- ~ o , o - ~: (.D- ~ ~. O o. N ~, (D ' _ - (~1: _ :: ~ t0 _ ~ O _ A A, e , mA (~0 ~ ~ A A A A o O - ~ _ _ .04 .08 .12 .16 .20 .24 Fig.8c - W w~ .. ll ~ ~ ~ ~ ~ ~ O ~0~ A~ 1 , ~ I · I 0 .04 .00 .12 .16 .20 .24 0 .04 .08 . i2 . i6 .20 .24 Fig.8b, c - Velocity components and Pressure as a function of Y. at given crosssections: XIL = .962 and XJL = 1.007, for several depths Z = const. , experiments; , computations. 308
ins - i~''2_22,'2 Fig.6 - Secondary velocities at several crosssections X = const.; left, experiments; right, computations. wit by G ~~ ~~ O X/L = 0.942 o,/~7 X/L = 0.955 ~n,.l .~y o- lo (D lo . At_ tool .~ Fig.7 - Isowake lines; left, experiments; right, computations. 309
~ x = -157mm - i /, /,, ~ /, /, ~ / 32^~" _ ~ , ~ , 1 _', ' ~ , i ; x = -53mm /' i., it,, /'','' /, A, / , - ., /, ~ , ,, , ~ ;. s . /,` ~ ,' /' ',2 ~ /'; ,, i, ~4 ,,"- ~; =_` ,, ~ _ ~,~ _ _ _ _ _ /~/ " /,' ~ /' 1 ! oc . I C : r . 0c ~ _ 1 G~ . ~ ~Jt Ot 1 ~ _~_ ot .( CO C~` r~ ~ .~ CO~ 'Dt I L t O . Q2 . 0~ 06 . 08 . 1 0 . 1 2 . 1 4 X X1L=.949 ,~~ - ~- = , , ~ -r- r O .02 .O] .06 .08 .10 . 1;? . 11 1 ',. X/L = .962 ~ '' , I r I I I G . 02 .01 . 06 . 013 . 1 Ct . 12 . 1'1 Fig. 10 - Contours of isolevels of the longitudinal component of the vorticitv fi~la in tr~ncv`~r~- crosssections. Left, experiments [443; right, computations. 3~0 , _~_ .. ~ ~v · _^ v~
o- x=l~ / . _~~ ~ ~ e ~5 = 1 1 1 . ... _ _ :'°~ / ~ '0 ~ ~ I WS ~- _~ \\ / / 0 .02 .~] .06 .~S .10 .12 .11 Fig. 10 - Contours of isolevels of the longitudinal component of the vorticity field in transverse crosssections. Left, experiments [441; right, computations. Fig.9 - Detailed view of the aft part of the HSVA tanker. ---, \(x) exact intersection between V and ~B; , present approximate "stair" representation of this line (note: the sting is omitted in the computations). ~Groupe de 1~1ode1 isat ion Nu~erique RSCETE - . ~ ~ . -- ~ ., , =, ; ' ', , -, ,_ ' - ~ .+ - - ^ - '; . _i __ ' - ' ' ' ' · ~ A;~ .~ ~ ' j ~ ~ , ; , A ~ y* ii ~ 't;; ' ~ ~ . ! '7n ~ - i i T. i · r 7 ~ ~ ~ ;; ~ , ., 7 7 ; ~ ~ ; ~ . ~ .._. .- ~ ' _ '~.; 7.~. 2' ~ i; ~ 7 ~ ~ ; . ;~i _ - ~ ~ +~ ~ ~ ~ ! j7~ 7~A j1 ~ _~_;7 ~ t-; - ~ ~ ~ · ~+ j ~ !~ .~ ~ 7 _ ^;, ~ .: 2 ~ , ,. - h,, _,, 7 ~, ~ ~ ~ ,~ ~ _ ! ~ ~ ~ ~ ~ 7 ; ~ T i 7 i i . _ ;~ - ~~ ~ i '' ~ ~; -t ~ ; , ~~~ , ~, ~ i - ~ -' ' ~ , ~ , $~7 ; j r ~' .;~J ~' 2 ~ 't ; ' ; 1 ~i ! ~'.'.~;'~~~ 72:.,,~,, .~' -=r-'_, j ,_ + ~ '';,, ~_ ~ f. + - ,: _. ~. " '~, ~ , +, . :', ., T.': - ~ ;; ~+~ _3,~ + - - - . - . . . - i i . . . : . _ . ? .: _ _ 7: :: ~ . L _ . 7 7 - ~ 7 : A i~ _ =- -~ +; ' -i '~ - - ; ~ ;'i ;- - '- - ' . ' i- ii7 .; - ~ - ,. Fig. 1 1 - Perspective view of the aft part of the HSVA tanker; skin-friction lines. 311
DISCUSSION by Y. Kodama You used a non-aligned grid system near the stern, which greatly reduces load for grid generation and improves grid orthogonality. However, I think a special treatment should be needed in order to satisfy conservation there. Did you do those treatments there? If not, what is your opinion about the effect of possible violation of conservation there on the computed results? Author's Reply Because of the use of an H-grid associate to the constraint x=x(~),the grid which has been built here is no more body-fitted in the intersection between the overhang and the plane of symmetry(y=0). This implies two consequences; (i) the exact geometry is not well represented, so that the no slip condition is not written at the right place in the plane of symmetry. (ii) the development of the boundary layer below the overhang is not correctly captured because of the lack of discretization points and a distribution law which is not in agreement with the strong variation of the velocity profiles in the plane y=0. Due to the used topology, this weakness is not easily surmountable because the grid lines which are supposed to describe the overhang symmetry line boundary layer come from the surface of the hull. It is believed that the correct representation of the body surface is more important than conserving momentum on a geometry which is incorrectly reproduced in the symmetry plane. The best remedy to this fundamental weakness is to relax the constraint x=x(~), working on more general (and less orthogonal) grids. 312
