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.

An Interactive Approach for Calculating Ship Boundary Layers and Wakes for Nonzero Froude Number Y. Tahara, F. Stern (Iowa Institute of Hydraulic Research, The University of Iowa, USA) B. Rosen (South Bay Simulations Inc., USA) ABSTRACT An interactive approach is set form for calculating ship boundary layers and wakes for nonzero Froude number. The Reynolds-averaged Navier-Stokes equa- tions are solved using a small domain with edge condi- tions matched with those from a source-doublet-Dawson method solved using the displacement body. An overview is given of both the viscous- and inviscid-flow methods, including their treatments of the free-surface boundary conditions and interaction procedures. Results are presented for the Wigley hull, including comparisons for zero and nonzero Froude number and with available experimental data and the inviscid-flow results, which validate the overall approach and enable an evaluation of the wave-boundary layer and wake interaction. NOMENCLATURE A¢,B¢,etc. Aij,Bij,aj,bjk i bj CD,CP,CU,Cnb finite-analytic coefficients (nb = NE, NW, SE, etc.) 2 friction coefficient (= 2~W/pUo) pressure coefficient residuary-resistance coefficient (= 2R/pSU2o) Froude number (= U<J~) conjugate metric tensor in.general curvilinear coordinates t turbulent kinetic energy characteristic (ship) length normal unit vector piezometric pressure residuary resistance Reynolds number (= UoLtv) wetted surface area boundaries of the solution domain Cf Cp CR Fr . . glJ k L n p R coefficients in transport equations influence coefficients geometric tensor Re S Sb,Se,etc. S¢,S t U,V,W UX~VY,WZ Uc UO Us uu,vv,etc. x,y,z x,r,0 x+,y+,z+ a* £ vt (,q,( p * Tij,Iij NEW Subscripts p INTRODUCTION characteristic (freestream) velocity wall-shear velocity (= Reynolds stresses Cartesian coordinates cylindrical polar coordinates dimensionless distances (= U:xtv, etc.) displacement thickness rate of turbulent energy dissipation free-surface elevation dipole strength kinematic viscosity eddy viscosity body-fitted coordinates density source strength time increment fluid- and external-stress tensors wall-shear stress transport quantities (U. V, W. k, £~; velocity potential edge value freestream or zero Fr value inviscid-flow value The interaction between the wavemaking of a ship and its boundary layer and wake is a classic and important problem in ship hydrodynamics. Initially, the interest was primarily with viscous effects on wave resistance and propulsive performance due to the lack of Reynolds number (Re) similarity in model tests. More recently, also of interest are the wave-bounda~y layer and wake interaction effects on the details of ship wakes and wave patterns due to the advent of satellite remote sens- ing. The present study is central to the aforementioned problems, i.e., it concerns the development of an interac- tive approach for calculating ship boundary layers and wakes for nonzero Froude number (Fr). Thus, both the effects of wavemaking on the boundary layer and wake and, vice versa, the effects of the boundary layer and wake on wavemaking are included in the theory, although the focus here is somewhat more on the former. source functions time; arclength in tangential direction tangent unit vector velocity components in cylindrical polar coordinates velocity components in Cartesian coordinates wake centerline velocity 699

Historically, inviscid-flow methods have been used to calculate wavemaking and viscous-flow methods the boundary layer and wake, in both cases, without accounting for the interaction. Recent work on wave- making has focused on the solution of the so-called Neumann-Kelvin problem using both Rankine- and Havelock-source approaches. Methods implementing these approaches were recently competitively evaluated and ranked by comparing their results with towing-tank experimental data [11. In general, the methods under- predicted the amplitude of the divergent bow waves, were lacking in high wave-number detail in the vicinity of the bow-wave cusp line, and overpredicted the ampli- tudes of the waves close to the stern. These difficulties were primarily attributed to nonlinear and viscous effects. The methods using the Havelock-source approach generally outperformed those using the Rankine-source approach, except with regard to the near- field results (i.e., within one beam length of the model) for which one of the the lager methods [2] was found to be far superior. Considerable effort has been put forth in the development of viscous-flow methods for ship boundary layers and wakes. Initially, three-dimensional integral and differential boundary-layer equation methods were developed; however, these were found to be inapplicable near the stern and in the wake. More recently, efforts have been directed towards the development of Navier- Stokes (NS) and Reynolds-averaged Navier-Stokes (RANS) equation methods; hereafter both of these will simply be referred to as RANS equation methods. At present, the status of these methods is such that practical ship geometries can be considered, including complexi- ties such as appendages and propellers. Comparisons with experimental data indicate that many features of the flow are adequately simulated; however, turbulence modeling and grid generation appear to be pacesetting Issues with regard to future developments (see, e.g., the review by Patel [3] and the Proceedings of the 5th International Conference on Numerical Ship Hydrodynamics [43~. Relatively little work has been done on the inter- action between wavemaking and boundary layer and wake. Most studies have focused separately on either the effects of viscosity on wavemaking or the effects of wavemaking (i.e., waves) on the boundary layer and wake. Professor Landweber and his students have both demonstrated experimentally the dependence of wave resistance on viscosity and shown computationally that by including the effects of viscosity in inviscid-flow cal- culations of wave resistance beKer agreement with exper- imental data is obtained (most recently, [5~. Such effects have been confirmed by others, including other more detailed aspects of the flow field such as surface- pressure distributions and wave profiles and patterns [61. Most studies concerning the effects of waves on boundary layer and wake have been of an approximate nature utilizing integral methods and assuming small- crossflow conditions (see Stern [7] for a more complete review, including references). In [7,8], experiment and theory are combined to study the fundamental aspects of the problem utilizing a unique, simple model and computational geometry, which enabled the isolation and identification of certain important features of the wave- induced effects. In particular, the variations of the wave- induced piezometric-pressure gradients are shown to cause acceleration and deceleration phases of the streamwise velocity component and alternating direction of the crossflow, which results in large oscillations of the displacement thickness and wall-shear stress as com- pared to the no-wave condition. For the relatively simple geometry studied, f~rst-order boundary-layer calculations with a symmetry-condition approximation for the free- surface boundary conditions were shown to be satisfac- tory; however, extensions of the computational approach for practical geometries were not successful [91. Miyata et al. [101 and Hino [11] have pursued a comprehensive approach to the present problem in which the NS equations (sub-grid scale and Reynolds aver- aged, respectively) are solved using a large domain with approximate free-surface boundary conditions. In both cases, the basic algorithms closely follow those of MAC [12] and SUMMAC [131. However, [10] uses a time- dependent free-surface conforming grid, whereas [111 uses a fixed grid which does not conform to the free sur- face. The results from both approaches are promising, but, thus far, have had difficulties in accurately resolving the boundary-layer and wake regions and, in the case of [10], have been limited to low Re. The present interactive approach is also compre- hensive. Two of the leading inviscid- [2] and viscous- flow [14] methods are modified and extended for inter- active calculations for ship boundary layers and wakes for nonzero Fr. The interaction procedures are based on extensions of those developed by one of the authors for zero Fr [151. The work of [7,8,15] is precursory to the present study. Also, it should be mentioned that the present study is part of a large project concerning free- surface effects on boundary layers and wakes. Some of the related studies under this project will be referenced later. In the following, an overview is given of both the viscous- and inviscid-flow methods, with particular emphasis on their treatments of the free-surface boundary conditions and the interaction procedures. Results are presented for the Wigley hull, including comparisons for zero and nonzero Fr and with available experimental data and inviscid-flow results, which validate the overall approach and enable an evaluation of the wave-boundary layer and wake interaction. In the presentation of the computational methods and results and discussions to follow, variables are either defined in the text or in the NOMENCLATURE and are nondimensionalized using the ship length L, freestream velocity UO, and fluid density p. COMPUTATIONAL METHODS Consider the flow past a ship-like body, moving steadily at velocity UO, and intersecting the free surface of an incompressible viscous fluid. As depicted in figure 1, the flow field can be divided into four regions in each of which different or no approximations can be made to the governing RANS equations: region 1 is the inviscid flow; region 2 is the bow flow; region 3 is the thin boundary layer; and region 4 is the thick boundary layer and wake. The resulting equations for regions 1 and 3 and their interaction (or lack of one) are well known. Relatively little is known about region 2. Recent exper- iments concerning scale effects on near-field wave patterns have indicated a Re dependency for the bow wave both in amplitude and divergence angle [161; however, this aspect of the problem is deferred for later study. Herein, we are primarily concerned with the flow 700

in region 4 and its interaction with that in region 1. As discussed earlier, the description of the flow in region 4 requires the solution of the complete RANS equations (or, in the absence of flow reversal, the so-called partially-parabolic RANS equations, however, this simplification will not be considered here). There are two possible approaches to the solution of the RANS equations: a global approach, in which one set of governing equations appropriate for both the inviscid- and viscous-flow regions are solved using a large solution domain so as to capture the viscous-invis- cid interaction; and an interactive approach, in which dif- ferent sets of governing equations are used for each region and the complete solution obtained through the use of an interaction law, i.e., patching or matching conditions. Both approaches are depicted in figure 1. The former approach is somewhat more rigorous because it does not rely on the patching conditions that usually involve approximations. Nonetheless, for a variety of reasons, both types of approaches are of inter- est. In [151, both approaches were evaluated for zero Fr by comparing interactive and large-domain solutions for axisymmetric and simple three-dimensional bodies using the same numerical techniques and algorithms and turbu- lence model. It is shown that both approaches yield sat- isfactory results, although the interaction solutions appear to be computationally more efficient. As men- tioned earlier, the present study utilizes the interactive approach. This takes advantage of the latest develop- ments in both the inviscid- and viscous-flow technolo- gies; however, a large-domain solution for the present problem is also of interest and a comparative evaluation as was done previously for zero Fr is planned for study under the present project for nonzero Fr. Viscous-Inviscid Interaction Referring to figure 1, there are two primary dif- ferences between the interactive and large-domain approaches with regard to the solution of the RANS equations: (1) the size of the solution domain, i.e., the placement of the outer boundary SO; and (2) the bound- ary (i.e., edge) conditions specified thereon. For the large-domain solution, uniform-flow and wave-radiation conditions are appropriate, whereas the interaction solu- tion requires the specification of the match boundary (i.e., SO) as well as an interaction law, and also a method for calculating the inviscid flow. In the present study, solutions were obtained with the match boundary at about 2b, where ~ is the boundary-layer and wake thickness. The interaction law is based on the concept of displacement thickness 8*. A three-dimensional 5* for a thick boundary layer and wake can be defined unambiguously by the two require- ments that it be a stream surface of the inviscid flow continued from outside the boundary layer and wake and that the inviscid-flow discharge between this surface and any stream surface exterior to the boundary layer and wake be equal to the actual discharge between the body and wake centerplane and the latter stream surface. A method for implementing this definition for practical geometries is presently under development [171; how- ever, in lieu of this, an approximate definition is used in which two-dimensional definitions for a*, i.e. ~ = `: (MU-) dr (1) for the keelplane and waterplane at each station are con- nected by a second order polynomial. In summary, the inviscid-flow solution is obtained for the displacement body 5*. This solution then provides the boundary conditions for the viscous- flow solution, i.e. U(So) = Up(So) = Ue W(SO) = Wp(SO) = We p(So) = pp(So) = pe (2) Because S* and Vp (SO) are not known a priori, an initial guess must be provided and the complete solution obtained by iteratively updating the viscous- and inviscid-flow solutions until the patching conditions (1) and (2) are satisfied. Viscous Flow The viscous flow is calculated using the large- domain method of Patel et al. [14], modified and extended for interactive calculations and to include free- surface boundary conditions. The details of the basic method are provided by [141. Herein, an overview is given as an aid in understanding the present modifications and extensions. Eguations and Coordinate System The RANS equations are written in the physical domain using cylindrical coordinates (x,r,8) as follows: au + ~ a ~rV' + ~ aW = 0 (3) Dt = ~ ~ (¢ + Hi) ~ ear Ad) ~ r am (u-w) UV+ReV2U (4) DDV W2 = - ad `~y -~9 +~ - r are (v-w) - r (A - ww) + Re (V2V-r~ at -rim) (5) DW VW ~ a - D +-= - ~ (uw) ~ ar (vw) r at (6 + ww) - r (vw) + Re (V2W + ~ DO ~ r=) (6) with 701

DD` =~+Uaa +vaa +_ a and v2=~+~+ ~ a + ~ a2 Closure of the RANS equations is attained through the use of the standard k-£ turbulence model without modifications for free-surface effects. The lim- ited experimental data available for surface-piercing bodies [18] indicate that, near a free surface, the normal component of turbulence is damped and the longitudinal and transverse components are increased. This effect has also been observed in open-channel flow [19] and in recent measurements for free-surface effects on the wake of a submerged flat plate [20] and a plane jet [211. Such a turbulence structure cannot, in fact, be simulated with an isotropic eddy viscosity turbulence model like the present one; however, this aspect of the problem is deferred for later study. In the standard k-£ turbulence model, each Reynolds stress is related to the corresponding mean rate of strain by the isotropic eddy viscosity vt as follows: au av 1 au aw -uv=vt(~r + ~X) -Uw=vt(r as + ax) 1av aw w -Vw=vt(r all + fir ~ r ) - UU =Vt (24) - 3-k - VV =Vt (2 ear) 3 k - ww = vt (r aa~ 2 Vr ) 3 k (7) Vt is defined in terms of the turbulent kinetic energy k and its rate of dissipation £ by k2 Vt = Cot £ where Cp is a model constant and k and £ are governed by the modeled transport equations Dk_ a (1 ak fit ax Rkax) + r Or(Rkraark) +~aa~ (Rika~ )+ G - £ (9) Do a 1 as 1 a 1 as Dt = aX (R aX' + r ar (R£ r ar) +~a~ (Ready + C£1 k G - C£2 k (10) G is the turbulence generation term G=Vtt2~(aaXU)+(aavr)+(-aa~w3+vr-) + (aU+av'2+ ~ aU+aaxw' +(rDa~v+aDw- r-~ }(11) The effective Reynolds number Rib is defined as Rip Re ~, (12) in which ~ = k for the k-equation (9) and ~ = £ for the £- equation (10~. The model constants are: Cal = .09, C£1 = 1.44, C£2 = 1.92, t7U = <'V = oW = ok = 1, C7£ = 1.3. The governing equations (3) through (12) are transformed into nonorthogonal curvilinear coordinates such that the computational domain forms a simple rect- angular parallelepiped with equal grid spacing. The transformation is a partial one since it involves the coor- dinates only and not the velocity components (U,V,W). The transformation is accomplished through the use of the expression for the divergence and "chain-rule" defi- nitions of the gradient and Laplacian operators, which relate the orthogonal curvilinear coordinates x1 = (x,r.,8) to the nonorthogonal curvilinear coordinates 41= (t,q,(~. In this manner, the governing equations (3) through (12) can be rewritten in the following form of the continuity and convec~ave-transport equations am (bin + b2V + b3W) + as (b2U + b2V + bow + Ha: (b3U + b3V + b3W) = 0 (13) all a2¢ + g22 a2¢ + g33 02~¢ = 2At a: (8, + 2B~ ad) + 2C~ am + If am + so (14) Discretization and Velocity-Pressure Coupling The convective-transport equations (14) are reduced to algebraic form through the use of a revised and simplified version of the finite-analytic method. In this method, equations (14) are linearized in each local rectangular numerical element, AK, = Aq = A: = 1, by evaluating the coefficients and source functions at the interior node P and transformed again into a normalized form by a simple coordinate stretching. An analytic solution is derived by decomposing the normalized equation into one- and two-dimensional partial-differen- tial equations. The solution to the former is readily obtained. The solution to the latter is obtained by the method of separation of variables with specified bound- ary functions. As a result, a twelve-point finite-analytic formula for unsteady, three-dimensional, elliptic equa- tions is obtained in the form 702

1 8 (P= R {) Cnb~nb l+CpECU+CD+ - ~ 1 + Cp(CU¢U + CD¢D +-~ - S)1 (15) It is seen that up depends on all eight neighboring nodal values in the crossplane as well as the values at the upstream and downstream nodes MU and ~D, and the values at the previous time step ~ I. For large values of the cell Re, equation (15) reduces to the partially- parabolic formulation which was used previously in other applications. Since equations (15) are implicit, both in space and time, at the current crossplane of calculation, their assembly for all elements results in a set of simultaneous algebraic equations. If the pressure field is known, these equations can be solved by the method of lines. However, since the pressure field is unknown, it must be determined such that the continuity equation is also satisfied. The coupling of the velocity and pressure fields is accomplished through the use of a two-step iterative procedure involving the continuity equation based on the SIMPLER algorithm. In the first step, the solution to the momentum equations for a guessed pressure field is cor- rected at each crossplane such that continuity is satisfied. However, in general, the corrected velocities are no longer a consistent solution to the momentum equations for the guessed p. Thus, the pressure field must also be corrected. In the second step, the pressure field is updated again through the use of the continuity equation. This is done after a complete solution to the velocity field has been obtained for all crossplanes. Repeated global iterations are thus required in order to obtain a converged solution. The procedure is facilitated through the use of a staggered grid. Both the pressure-correction and pres- sure equations are derived in a similar manner by substi- tuting equation (15) for (U,V,W) into the the discretized form of the continuity equation (13) and representing the pressure-gradient terms by finite differences. Solution Domain and Boundary Conditions The solution domain is shown in figure 1. In terms of the notation of figure 1, the boundary condi- tions on each of the boundaries are as follows. On the inlet plane Si, the initial conditions for ~ are specified from simplellat-plate and the inviscid-flow solutions. On the body surface Sb, a two-point wall-function approach is used. On the symmetry plane Sk, the conditions imposed are 8(U,V,p,k,£~/a~ = W = 0. On the exit plane Se, axial diffusion is negligible so that the exit conditions used are a2~/aX2 = 0, and a zero-gradient condition is used for 19. On the outer boundary SO, the edge conditions are specified according to (2), i.e., (U,W,6) = (Ue,We, be) and a (k,£~/ar = 0, where (Ue,We,pe) are obtained from the inviscid-flow solution evaluated at the match boundary SO. On the free-surface S,~(or simply 11), there are two boundary conditions, i.e. and V n = 0 * tijnj = tiinj (16) (17) where n is the unit normal vector to the free surface and tij and All are the fluid- and external-stress tensors respectively, the latter, for convenience, including surface tension. The kinematic boundary condition expresses the requirement that ~ is a stream surface and the dynamic boundary condition that the normal and tangential stresses are continuous across it. Note that ~ itself is unknown and must be determined as part of the solution. In addition, boundary conditions are required for the turbulence parameters, k and £; however, at present, these are not well established. In the present study, the following approxima- tions were made in employing (16) and (17~: (a) the external stress and surface tension were neglected; (b) the normal viscous stress and both the normal and tan- gential Reynolds stresses were neglected; (c) the curva- ture of the free surface was assumed small and the tan- gential gradients of the normal velocity components were neglected in the tangential stresses; and (d) the wave ele- vation was assumed small such that both (16) and (17) were represented by first-order Taylor series expansions about the mean wave-elevation surface (i.e., the water- plane Sw). Subject to these approximations, (16) and (17) reduce to the following: (UX11 X + Vylly - WZ) ~ = 0 (18) Sw 19(Sw)=ll/Fr2-~] afi I (19) az ~ low Ia~v,k,£' =0 US ~ NEW (20) where Cartesian coordinates (x,y9z) have been used in (18) and (19~. Conditions (18) through (20) were implemented numerically as follows. The kinematic condition (18) was used to solve for the unknown free- surface elevation ~ by expressing the derivatives in finite-difference form and 11 in terms of its difference from an assumed (or previous) value. A backward dif- ference was used for the x-derivative, a central difference for the y-derivative, and the inviscid-flow lip was used as an initial guess. The dynamic conditions, (19) and (20), were used in conjunction with the solution for ~ in solving the pressure and momentum and turbulence model equations, respectively. Backward differences were used for the z- and ~derivatives. Inviscid Flow The inviscid flow is calculated using the method of Rosen [2], i.e., the SPLASH computer code. The method is an extended version of the basic panel method of Maskew [22,231 originally developed for the predic- tion of subsonic aerodynamic flows about arbitrary con- figurations, modified to include the presence of a free surface and gravity waves both for submerged and surface-piercing bodies. As is the case with the basic 703

method, lifting surfaces and their associated wake treatments as well as wall boundaries are included; however, the present overview and calculations are for nonlifting unbounded flow (see [24] for SPLASH results for lifting flow). The details of the basic method are provided by [22,231. Herein, an overview is given as an aid in understanding the extensions for the inclusion of the free surface and gravity waves and the present interaction calculations. The flow is assumed irrotational such that the governing differential equation is the Laplace equation V20 = 0 (21) where 0 is the external perturbation velocity potential, i.e. Vp = UOx + Vo (22) A solution for ~ may be obtained by defining also an internal perturbation potential ¢~ and applying Green's theorem to both the inner and outer regions and combining the resulting expressions to obtain ¢= - J FOUR-~) + R~} dS (23) sb where RpQ is the distance from the surface point Q to the field point P and ~ = ¢~ - O and c, = D(¢ - o~/anQ are the dipole and source strengths, respectively. In [22], the nature of solutions to (23) is investigated for two dif- ferent specifications for hi, i.e., ¢~ = 0 and UOx. In both cases, (23) is solved for the surface potential (i.e., ¢(Sb)) by representing the body by flat quadrilateral panels over which ~ and ~ are assumed constant and uti- lizing the far-field a ~ O and body D¢/an = - UOnx boundary conditions. The zero internal perturbation potential formulation Cot = 0) is shown to produce "results of comparable accuracy to those from higher- order methods for the same density of control points." In this case, the velocity normal to the external surface vn is Vn = UOnx + 80/0n = UOnx + ~ (24) and, the velocity tangent to the external surface V' is V~=Uotx+a~l~t=uotx-~/8t (25) where tx is the x-component of a tangent vector and t is arclength in a tangential direction. For solid surfaces, Vn is usually zero, but it may be a specified nonzero value to simulate body motion, boundary-layer growth, inflow and outflow, control-surface deflection, etc. Hence, in the basic method, (24) is used to evaluate the source strengths directly. The corresponding doublet strengths are then given by solution of the discretized form of (23~. Values of V' are subsequently computed using (25) with a central difference for the t-derivative. It should be recognized that the so-called zero internal perturbation formulation is, in fact, equivalent to methods based on Green's third formula applied directly to the external perturbation potential (e.g., [251~. In the SPLASH code, the internal zero-perturba- tion boundary condition is satisfied not only inside the submerged portion of the configuration, but also on the "other side" of a finite portion of the free surface. Both are represented by source-doublet singularity panels and flow leakage from one side of the free-surface to the other, at the free-surface outer boundary is assumed to be negligible. This assumption is valid if the outer boundary of the free surface is sufficiently far from the configuration, and if wave disturbances are eliminated before reaching the free-surface outer boundary. In this case, the discretized form of (23) is (i = Hi, Aij ale + Hi, Bij Hi = 0 (26) Sb + Sw Sb + Sw The free-surface shape is determined by representing the undisturbed free surface by panels, whereupon free-surface boundary conditions linearized with respect to zero Fr are imposed [261. The zero Fr velocities, UO, VO, and WO, are obtained by first considering all free-surface panels as solid and fixed (in contrast to a traditional approach which employs the double panel or image model). The nonzero Fr velocities are then expressed as small increments to those for zero Fr. The velocities tangent and normal to a free-surface panel are, respectively Ux ~ UO + AU Vy ~ VO + TV (27) and since WO = 0 for a free-surface panel. Through Bernoulli's equation, the pressure on free-surface panels is a function of local velocity, and is approximated by retaining only first-order incremental velocity terms Vn = Wz ~ WO + AW ~ AW (28) pa = ~ ~ 1 (U2 + Vy + W2) ~ ~2 1-(U20+V2)t-{Uo~U+Vo~V } ~ 2 ~ ~ - (U20 + V20 ~ ~ - { Uo (Ux - Uo) +Vo (Vy~ Vo) } (29) Free-surface boundary conditions are linearized in a similar manner, retaining only h~rst-order incremental velocity and surface-elevation terms. The kinematic free- surface boundary condition (18) is approximated by Wz = Vn ~ UO 11X + VO lly ~ (UO + VO)1/2 also (30) where the subscript so denotes differentiation along a zero Fr streamline. The dynamic free-surface boundary condition (19), after differentiation along so, and substi- tuting for 1lsO from (30), becomes aSO Fr2 BUS + V2~1/2 (31) 704

A five-point backward difference is used in the ~ and ~ directions, and the free-surface grid metrics, are used to compute the pressure gradient COOP + vO8P an ax ay asO- (U2+v2~ll2 = uO (aP a; + ap all) + V (aP at + ap all) at ax ~ ax at by Do, by - (32) (U2 + V2)1/2 The pressure-gradient algorithm is structured to permit the use of any blocked free-surface grid arrangement. Also, using less than a f~ve-point backward difference tends to dampen wave amplitudes. This wave-damping mechanism is employed on panels near the outer boundary of the finite free-surface model, so that wave disturbances are eliminated before reaching the free- surface outer boundary. At this point, a sufficient number of linear dependencies have been established to permit the elimi- nation of the unknown free-surface source strengths in (26), i.e., (24) relates source strength to panel normal velocity, (31) relates free-surface panel normal velocity to streamwise pressure gradient, (32) with backward differences relates streamwise pressure gradient to free- surface pressures, (29) relates free-surface pressure to free-surface panel tangential velocities, (25) relates panel tangential velocities to the local surface gradient of dou- blet strength, and central differences relate the local surface gradient of doublet strength to doublet strengths. Hence, free-surface source strengths can be expressed as a linear combination of free-surface doublet strengths, A. <,j-aj +2 . bjk ink w Substituting for Hi from (33) into (26) yields (i - ~ Aij ~j + ~ Bij ~j . Sb + Sw sb +v sw Bij (aj +2 bjk ink) w With free-surface source strengths eliminated, and source strengths on the solid body evaluated directly, solution of (34) yields the corresponding dou- blet strengths. The free-surface source strengths are then given by (33), and (24) and (25) are used to compute the resulting velocities on both body and free-surface panels. Pressures on free-surface panels are given by (29~. A similar linearized formula is used for pressures acting on body panels, and configuration forces and moments are obtained by panel pressure integration. For interactive calculations, the SPLASH code calculates the inviscid free-surface flow about the equivalent displacement body resulting from the previous viscous calculation. For this purpose, the equivalent displacement body is treated as a solid fixed surface. The inviscid flow velocities required for the next viscous flow calculation, at off-body points on the viscous grid outer boundary SO, are obtained using the computed source-doublet solution and velocity influence coefficients. A sub-panel velocity influence-coefficient algorithm was developed which utilizes a bilinear varia- tion of source and doublet strength across each panel. The continuous variation of source and doublet strength on each panel, and across panel edges, enhances the accuracy of off-body velocity calculations at points close to any body andlor free-surface panels. WIGLEY HULL GEOMETRY AND EXPERIMENTAL INFORMATION The Wigley parabolic hull was selected for the initial calculations since the geometry is relatively simple and it has been used in many previous computational and experimental studies. In particular, it is one of the two hulls, the other being the Series 60 CB = .6 ship model, selected by the Cooperative Experimental Program (CEP) of the Resistance and Flow Committee of the International Towing Tank Conference [27] for which extensive global (total, wave pattern, and viscous resis- tance, mean sinkage and trim, and wave profiles on the hull) and local (hull pressure and wall-shear-stress distri- butions and velocity and turbulence fields) measurements were reported. It was for these same reasons that the Wigley hull was selected as the first test case of the basic viscous-flow method [141, including comparisons with some of the zero Fr data of the CEP. Herein, compar- isons are made for zero Fr with this same data and for nonzero Fr with the appropriate data of the CEP. As will be shown later, the nonzero Fr data is not as complete or of the same quality as that for zero Fr, which was the motivation for a related experimental study for the Series 60 CB - .6 ship model [28] for which calculations and comparisons are in progress. However, the comparisons are still useful in order to validate the present interactive approach and display the shortcomings of both the computations and experiments. The coordinates of the Wigley hull are given by y = 2~4x(1 - x)~1 - (z/d)2} (35) where B = .1 and d = .0625. Waterplane and typical crossplane views are shown in figure 2. RESULTS In the following, first, the computational grids (figures 2 and 3) and conditions are described. Then, some example results are presented and discussed for zero Fr, followed by those for nonzero Fr, including wherever possible comparisons with available experi- mental data, and, in the latter case, with inviscid-flow results. The convergence history of the pressure is shown in figure 4. Figure 5 provides a comparison of the large-domain and interactive solutions. The free-sur- face perspective view and contours, wave profile, and surface-pressure profiles and contours are shown in fig- ures 6 through 10, respectively. The axial-velocity con- tours, crossplane-velocity vectors, and pressure, axial- vorticity, and turbulent kinetic energy contours for sev- eral representative stations are shown in figures 11 through 13. Lastly, the velocity, pressure, and turbulent 705

kinetic energy profiles for similar stations are shown in figures 14 through 16. On the figures and in the discus- sions, the terminology "interactive" refers to results from both the interactive viscous and displacement-body inviscid solutions. When the distinction is not obvious it will be made. The terminology "inviscid" or "bare- body" refers to the noninteractive inviscid solution. Computational Grids and Conditions The viscous-flow computational grid was obtained using the technique of generating body-fitted coordinates through the solution of elliptic partial differ- ential equations. Because of the simplicity of the present geometry, it is possible to specify the axial fit and cir- cumferential f3 control functions as, respectively, only functions of I, and (; however, in order to accurately satisfy the body-surface boundary condition and resolve the viscous flow, f2 = f2~t,ll,(~. Partial views of the grids used in the calculations are shown in figures 2a,b for a longitudinal plane and typical body and wake cross- planes, respectively. Initially, a large-domain grid was generated. Subsequently, a small-domain grid was obtained by simply deleting that portion of the large- domain grid that lay beyond about r > .2. The outer boundary for the small-domain grid is shown by the dashed line in figure 2. For the large-domain grid, the inlet, exit, and outer boundaries are located at x = (.296,4.524) and r = 1, respectively. The first grid point off the body surface is located in the range 90 < y+ < 250. 50 axial, 30 radial, and 15 circumferential grid points were used. As already indicated, the small- domain grid was similar, except 21 radial grid points were used. In summary, the total number of grid points for the large- and small-domain calculations are 22,500 and 15,750, respectively. The inviscid-flow displacement-body and free- surface panelization is shown in figure 3. 423 panels are distributed over the displacement body and 546 over the free surface for a total number of 969 panels. The panel- ization covers an area corresponding to 1 ship length upstream of the bow, 1.5 ship lengths in the transverse direction, and 3 ship lengths downstream of the stern. This panel arrangement was judged optimum based on panelization dependency tests E161. The conditions for the calculations are as follows: L= l;Uo= l;Re=4.5x 106;Fr=0and.316; end on the inlet plane the average values for ~ and Us are .0033 and .0455, respectively. These conditions were selected to correspond as closely as possible to those of the experiments of the CEP with which comparisons will be made [5,29,301. Initially, large-domain calculations were per- formed for zero Fr. A zero-pressure initial condition was used and the values for the time a`, pressure ap, and transport quantity Ad (where ~ = k and £) underrelaxation factors and total number of global iterations were .05 and 200, respectively. Next, small- domain calculations were performed, first for zero Fr, and then for nonzero Fr. For zero Fr, the interaction calculations were started with a zero-pressure initial condition and freestream edge conditions (Ue = 1,We= Pe = 01. After 200 global iterations, the edge conditions were updated using the latest values of displacement thickness. Subsequently, the edge conditions were updated every 200 global iterations until convergence was achieved, which took three updates. For nonzero Fr, the calculations were started with the zero Fr solution as the initial condition and with nonzero Fr edge condi- tions obtained utilizing the zero Fr displacement body. This solution converged in 200 global iterations. Most of the results to be presented are for this case; however, some limited results will be shown in which the nonzero Fr edge conditions were obtained using an updated nonzero Fr displacement body. The values for at, up, and At (where ~ = k and £) used for the small-doma~n calculations were the same as those for the large-domain calculations; however, for nonzero Fr, in addition, a value of .01 was used for Al (where ~ = U) for grid nodes near the outer boundary. The a~9/3z term in (19) was found to have a small influence and was neglected in many of the calculations; however, this may be due, in part, to the present grid resolution. The calculations were performed on the Naval Research Laboratory CRAY XMP-24 supercomputer. The CPU time required for the calculations was about 17 minutes for 200 global iterations for the viscous-flow code and 1 minute for the inviscid-flow code. Extensive grid dependency and convergence checks were not earned out since these had been done previously both for the basic viscous-flow method [14] and for other applications. However, some calculations were performed using both coarser and finer grids. These converged, respectively, more rapidly and slower than the present solution. Qualitatively the solutions were very similar to the present one, but with reduced and somewhat increased resolution, respectively. The convergence criterion was that the change in solution be less than about .05% for all variables. Usually the solu- tions were carried out at least 50 global iterations beyond meeting this criterion. Figure 4 provides the conver- gence history for the pressure and is typical of the results for all the variables. In figure 4, the abscissa is the global iteration number it and the ordinate is the residual R(it), which is defined as follows: imax imax R(it)= Ad, ~ptit-l) - pkit) i/ At, Pith ~ (36) i=1 i=1 where ill and imax are the total number of iterations and grid points, respectively. Referring to figure 4, global iterations 1 - 200 correspond to the final iterations of the zero Fr solution and global iterations 200 - 400 to those for the nonzero Fr solution. Zero Fr Figure 5 provides a comparison of the zero Fr large-domain and interactive solutions and experimental data. The two solutions are nearly identical and show good agreement with the data, which validates the pre- sent interactive approach. The agreement with the data for the large-domain case is, of course, not surprising since this was already established in [14] for a similar grid and conditions, i.e., the present zero Fr solution is essentially the same as that of [141. Some additional aspects of the zero Fr solution are displayed in figures 11 through 15 for later comparison with the nonzero Fr solution. Reference [14] provides detailed discussion of the zero Fr solution, including comparisons with the available experimental data. In summary, there is a downward flow on the forebody and an upward flow on the afterbody in response to the external-flow pressure gradients. The boundary layer and wake remain thin and attached and the viscous-inviscid interaction is weak; 706

however, on the forebody, the boundary layer is rela- tively thicker near the keel than the waterplane, whereas the reverse holds true on the afterbody and in the near wake. The stern vortex is very weak. In the intermedi- ate and far wake, the flow becomes axisymmetric. As indicated in figures 5 and 14 through 16, the agreement between the calculations and data is quite good; however, there are some important differences, which are primarily attributed to the deficiencies of the standard k-e turbulence model with wall functions. In particular, the axial velocity and turbulent kinetic energy are overpredicted near the stern and there is a more rapid recovery in the wake. Nonzero Fr Figure 5 also includes nonzero Fr results for comparison. On the waterplane, the surface and wake centerplane pressure displays very dramatic differences, the wall-shear velocity shows similar trends, but with reduced magnitude, and the wake centerplane velocity indicates faster recovery in the intermediate and far wake. As will be shown later, the first closely follows the wave profile, the second is due to an increase in boundary-layer thickness near the waterplane for the nonzero Fr case, and the third can be explained by the wave-induced pressure gradients. On the keel, all three of these quantities are nearly the same as for zero Fr. The free-surface perspective views (figure 6) and contours (figure 7) vividly display the complex wave pattern consisting of both diverging and transverse wave systems. The bow and stern wave systems are seen to initiate with crests and the shoulder systems initiate with troughs, which conforms to the usual pattern described for this type of hull form. Very apparent is the reduced amplitude of the stern waves for the interactive as com- pared to the inviscid solution. Also, the diverging wave system is more pronounced and at a smaller angle with respect to the centerplane. Note that the axial and trans- verse wave-induced pressure gradients can be discerned from these figures, but with an appropriate phase shift, i.e., increasing and decreasing wave elevations imply, respectively, adverse and favorable gradients. The wave profile along the hull is shown in figure 8, which, in this case, includes experimental data for comparison. On the forebody, the two solutions are nearly identical and underpredict the amplitude of the bow-wave crest and the first trough. On the afterbody, the interactive solution indicates larger values than the inviscid solution, with the data in between the two. The wave profile for the nonzero Fr displacement body (figure 3b) is also shown in figure 8. The differences are minimal on the forebody, whereas, they are significant on the afterbody and depart from the data. It appears that the present simple definition ( 1 ) is insufficient for "wavy" displacement bodies. The surface-pressure profiles (figure 9) show similar tendencies as just discussed with regard to the wave profile. On the forebody, the two solutions are nearly identical, but, in this case, in very close agreement with the data. The pressure on the forebody shown by the dashed line is that obtained from the inviscid dis- placement-body solution. On the afterbody, here again, the interactive solution indicates larger values than the inviscid solution, with the data in between the two. The wave-induced effects are seen to diminish with increasing depth and the agreement between the two solutions and the data on the afterbody shows improvement. The surface-pressure contours (figure 10) graphically display the differences between the two solutions and the data. Note that the axial and vertical surface-pressure gradients can be discerned from these figures, i.e., increasing and decreasing pressure imply, respectively, adverse and favorable gradients. The larger wave elevation and pressure on the afterbody for the interactive solution results in the closed contours near the stern displayed in figure 10b. As already mentioned, the viscous-inviscid interaction is weak for the Wigley hull, which is the reason that the inviscid and viscous pressure distributions are quite similar. However, it appears that the interaction is greater for nonzero as compared to zero Fr. Figures 11 through 13 show the detailed results for several representative stations, i.e., x = .506, .904, and 1.112, although the discussion to follow is based on the complete results at all stations. Note that for zero Fr the upper boundary shown is the waterplane, whereas for nonzero Fr, it is the predicted free surface. Also, the axial-velocity, -vorticity, and turbulent kinetic energy contours are not shown for the inviscid solution since, in the former case, their values are all very close to 1 and, in the latter two cases, they are, of course, zero. Solid curves indicate clockwise vorticity. On the forebody (figure 11), the boundary layer is thin such that many aspects of the solutions are simi- lar; however, there are some important differences. The nonzero Fr pressure fields show local and global effects of the free surface, i.e., near the free surface, regions of high and low pressure coincide with wave crests and troughs, respectively, and at larger depths, the contours are parallel to the free surface. Also, for nonzero Fr, the crossplane-velocity vectors are considerably larger, especially for the interactive solution. The inviscid solu- tion clearly lacks detail near the hull surface. The extent of the axial vorticity is increased for nonzero Fr and is locally influenced by the free surface. In both cases, as expected, the direction of rotation is mostly anticlock- w~se. On the afterbody (figure 12), almost all aspects of the solutions show significant differences. The boundary layer is thicker near the waterplane for nonzero as compared to zero Fr. This behavior begins at x ~ .825, which coincides with a region of adverse axial wave-induced pressure gradient (see figure 7~. The differences for the pressure field and axial-vorticity contours are similar as described for the forebody; however, in the case of the crossplane-velocity vectors, there is an additional difference that near the free surface the interactive solution displays downward flow. This is consistent with the fact that the free-surface elevation is above the waterplane and the pressure is generally higher near the free surface than it is a larger depths, i.e., 11 ~ O and ap/az ~ 0. Note that, as expected, in both cases, the direction of rotation for the axial-vorticity is mostly clockwise. The turbulent kinetic energy contours are nearly the same for both Fr. In the wake (figure 13), the solutions continue to show significant differences. Initially, the low-velocity region diffuses somewhat and covers a larger depthwise region; then, for x > 1.2, recovers quite rapidly. A similar behavior was noted earlier for the wake centerline velocity for x ~ 1.2, both of which, as already mentioned, are consistent with the wave pattern. The zero Fr pressure field is nearly axisymmetric and fully 707

recovered by the exit plane. The nonzero Fr pressure field continues to show free-surface effects, i.e., the contours are parallel to the free surface, but also fully recovered by the exit plane. Note the considerably larger wave elevation near the wake centerplane for the inviscid as compared to the interactive solution, which was pointed out earlier with regard to figures 6 and 7. Here again, the crossplane-velocity vectors are larger for nonzero as compared to zero Fr, especially near the wake centerplane for the interactive solution. The interactive and inviscid solutions display differences near the free surface, which appear to be consistent with the differences in their predicted wave patterns. The zero Fr axial vorticity decays fairly rapidly, whereas, for nonzero Fr, the decay is slow with a layer of nonzero vorticity persisting near the free surface all the way to the exit plane. The turbulent kinetic energy contours are similar for both Fr, but recover faster for the nonzero case. Figures 14 through 16 show the velocity, pressure, and turbulent kinetic energy profiles for similar stations as for figures 11 through 13, i.e., x = .5, .9, and 1.1. Also, included are both zero and non zero Fr experimental data. At the largest two depths, z = .05 and .0625, data for both Fr are available, whereas, at the waterplane, z = 0, only zero Fr data are available. At the intermediate depths, data are available for both Fr, but for different z values. Since the interest here is primarily nonzero Fr and the zero Fr data and comparisons were already displayed in [14], only nonzero Fr data are shown for z = .0125, .025, and .0375. For zero Fr, a corrected pressure is also shown which includes a constant (= -.03) reference-pressure correction as described in [141. Turbulent kinetic energy data are only available for zero Fr. At x = .5, consistent with previous discussions the differences between the two solutions are quite small and the agreement with the zero Fr data is good. However, the nonzero Fr data show some unexpected differences. In particular, the axial-velocity profile has a laminar appearance and the boundary-layer thickness is relatively large, the vertical velocity is upward, and the pressure shows considerable scatter. It is pointed out in [5] that the pressure-measurement error was appreciable. At x = .9 and 1.1, here again, consistent with previous discussions the differences between the two solutions are significant and the agreement between the zero Fr solution and data is good, except for the aforementioned discrepancies. The nonzero Fr solution shows larger axial velocities than the measurements for the inner part of the profiles. Here again, the measured profiles have a laminar appearance and the boundary layer is thick. However, no doubt, a part of the difference is due to the calculations, i.e., as is the case for zero Fr, due to deficiencies of the k-e turbulence model an overprediction of the velocity near the wall and wake centerplane is expected. The transverse velocity is small and with similar trends for both calculations and measurements. The calculations indicate downward vertical velocities near the free surface and upward values for the midgirth region and near the keel. The agreement with the data near the keel is satisfactory, but in the midgirth region and near the free surface the data display greater upward flow than the calculations. In the wake, the nonzero Fr data show surprisingly small vertical velocities near the wake centerplane. Here again, the nonzero Fr pressure data shows considerable scatter and is difficult to compare with the calculations. Consistent with earlier discussions the turbulent kinetic energy profiles are nearly the same for both Fr. Lastly, Table 1 provides a comparison of the cal- culated pressure-resistance coefficient and experimental values of the residuary-resistance (i.e., total - frictional) coefficient. The experimental values cover a range of Re, including the present value, and clearly show a dependency on Re. Interestingly, the inviscid result compares well with the data at the highest Re, whereas the interactive result is close to that that the data implies at the present Re. WAVE-BOUNDARY LAYER AND WAKE INTERACTION The comparisons of the zero and nonzero Fr interactive and inviscid-flow results with experimental data enables an evaluation of the wave-boundary layer and wake interaction. Very significant differences are observed between the zero and nonzero Fr interactive results due to the presence of the free surface and gravity waves. In fact, the flow field is completely altered. Most of the differences were explicable in terms of the differences between the zero and nonzero Fr surface- pressure distributions and, in the latter case, the addi- tional pressure gradients at the free surface associated with the wave pattern. The viscous-inviscid interaction appears to be greater for nonzero as compared to zero Fr. It should be mentioned that other factors undoubtedly have important influences, e.g., wave-induced separation, which are not included in the present theory. The interactive and inviscid nonzero Fr solutions also indicate very significant differences. The inviscid solution clearly lacks "real-fluid effects." The viscous flow close to the hull and wake centerplane is clearly not accurately resolved. The interactive solution shows an increased response to pressure gradients as compared to the ~nv~sc~d solution, especially in regions of low velocity. Also, the inviscid solution overpredicts the pressure recovery at the stern and the stern-wave amplitudes. CONCLUDING REMARKS The present worlc demonstrates for the first time the feasibility of an interactive approach for calculating ship boundary layers and wakes for nonzero Fr. The results presented for the Wigley hull are very encouraging. In fact, in many respects, the present results appear to be superior to the only other solutions of this type available, i.e., [10,111. This is true both with regard to the resolution of the boundary-layer and wake regions and the wave field. Furthermore, it appears that the present interactive approach is considerably more computationally efficient than the large-domain approaches of [10,111. This is consistent with the previous finding for zero Fr [151. However, a complete evaluation of the present method was not possible. In the former case, due to the limited available experimental data. As mentioned earlier, a related experimental study for the Series 60 CB = .6 ship model [28] was recently completed for which extensive measurements were made at both low and high Fr for which calculations and comparisons are in progress. In the latter case, due to the considerable differences in numerical techniques and algorithms and turbulence models between the present methods and those of 708

[10,11]. As mentioned earlier, the pursuit of a large- 10. domain approach to the present problem is also of Interest and will enable such an evaluation. Finally, some of the issues that need to be addressed while further developing and validating the present approach are as follows: further assessment of the most appropriate free-surface boundary conditions; improved definition and construction of displacement bodies; the inclusion and resolution of the bow-flow region; extensions for lifting flow; and the ever present problem of grid generation and turbulence modeling. Also, of interest is the inclusion of nonlinear effects in the inviscid-flow code. ACKNOWLEDGEMENTS This research was sponsored by the Office of Naval Research under Contract N00014-88-K-0113 under the administration of Dr. E.P. Rood whose sup- port and helpful technical discussions are greatly appre- c~ated. REFERENCES 1. Lindenmuth, W., Ratcliffe, T.J., and Reed, A.M., "Comparative Accuracy of Numerical Kelvin Wake Code Predicitions - "Wake-Off"," DTRC/SHD- 1260-1, 1988. 2. Rosen, B.,"SPLASH Free-Surface Code: Theoretical/Numerical Formulation," South Bay Simulations Inc., Babylon, NY, 1989 (proprietary report). 3. Patel, V.C., "Ship Stern and Wake Flows: Status of Experiment and Theory," Proc. 17th Office of Naval Research S ymposium on Naval Hydrodynamics, The Hague, 1988, pp. 217-240. Proc. 5th International Conference on Numerical Ship Hydrodynamics, Hiroshima, 1989. 5. Shahshahan, A., "Effects of Viscosity on Wavemaking Resistance of a Ship Model," Ph.D. Thesis, The University of Iowa, Iowa City, IA, 1985. 6. Ikehata, M. and Tahara, Y., "Influence of Boundary Layer and Wake on Free Surface Flow around a Ship Model," J. Society of Naval Architects of Japan, Vol. 161, 1987, pp. 49-57 (in Japanese). Stern, F., "Effects of Waves on the Boundary Layer of a Surface-Piercing Body," J. of Ship Research, Vol. 30, No. 4, 1986, pp. 256-274. 8. Stern, F., Hwang, W.S., and Jaw, S.Y., "Effects of Waves on the Boundary Layer of a Surface- Piercing Flat Plate: Experiment and Theory," J. of Ship Research, Vol. 33, No. 1, 1989, pp. 63-80. 9. Stern, F., "Influence of Waves on the Boundary Layer of a Surface-Piercing Body," Proc.4th International Conference on Numerical Ship Hydrodynamics, Washington, D.C., 1985, pp. 383-406. 709 Miyata, H., Sato, T., and Baba, N., "Difference Solution of a Viscous Flow with Free-Surface Wave about an Advancing Ship," J. of Computational Physics, Vol. 72, No. 2, 1987, pp. 393-421. 11. Hino, T., "Computation of a Free Surface Flow around an Advancing Ship by the Navier-Stokes Equations," Proc. 5th International Conference on Numerical Ship Hydrodynamics, Hiroshima, 1989. 12. Harlow, F.H. and Welch, J.E., "Numerical Calculation of Time-Dependent Viscous Flow of a Fluid with Free Surface," The Physics of Fluids, Vol. 8, 1965, pp. 2182-2189. 13. Chan, R.K.C. and Street, R.L., "A Computer Study of Finite-Amplitude Water Waves," J. of Computational Physics, Vol. 6, 1970, pp. 68-94. 14. Patel, V.C., Chen, H.C. and Ju, S., "Ship Stern and Wake Flows: Solutions of the Fully-Elliptic Reynolds-Averaged Navier-Stokes Equations and Comparisons with Experiments," Iowa Institute of Hydraulic Research, The University of Iowa, IIHR Report No. 323, 1988; also J. of Computational Physics, Vol. 88, No. 2, June 1990, pp. 305-336. 15. Stern, F., Yoo, S.Y. and Patel, V.C., "Interactive and Large Domain Solutions of Higher-Order Viscous-Flow Equations," AIAA Journal, Vol. 26, No. 9, 1988, pp. 1052- 1060. 16. Longo, J., "Scale Effects on Near-Field Wave Patterns," M.S. Thesis, The University of Iowa, Iowa City, IA, 1990. 17. Black, R., "Definition of Three-Dimensional Displacement Thickness Appropriate for Ship Boundary Layers and Wakes," M.S. Thesis, The University of Iowa, Iowa City, IA, expected 1991. 18. Hotta, T. and Hatano, S., "Turbulence Measurements in the Wake of a Tanker Model on and under the Free Surface," Fall Meeting of the Society of Naval Architects of Japan, 1983. 19. Rodi, W., "Turbulence Model and Their Application in Hydraulics," Presented at the IAHR Section on Fundamentals of Division II: Experimental and Mathematical Fluid Dynamics, 1980. 20. Swean, T.F. and Peltzer, R.D., "Free Surface Effects on the Wake of a Flat Plate," NRL Memo Report 5426, Naval Research Laboratory, Washington D.C., 1984. 21. Ramberg, S.E., Swean, T.F., and Plesniak, M.W., "Turbulence Near a Free Surface in a Plane Jet," NRL Memo Report 6367, Naval Research Laboratory, Washington D.C., 1989. 22. Maskew, B., "Prediction of Subsonic Aerodynamic Characteristics: A Case for Low- Order Panel Methods," Journal of Aircraft, Vol. 19, No. 2, 1982, pp. 157-163.

23. Maskew, B., "A Computer Program for Calculating the Non-Linear Aerodynamic Characteristics of Arbitrary Configurations," NASA CR- 166476, 1982. 24. Boppe, C.W., Rosen, B.S., Laiosa, J.P., and Chance, B., Jr., ''Stars & Stripes '87: Computational Flow Simulations for Hydrodynamic Design," The Eighth Chesapeake Sailing Yacht Symposium, Annapolis, MD., 1987. 25. Stern, F., "Comparison of Computational and Experimental Unsteady Cavitation on a Pitching Foil, ASME J. Fluids Eng., Vol. 111, 1989, pp. 290-299. 26. Dawson, C.W., "A Practical Computer Method for Solving Ship-Wave Problems," Proc. 2nd International Conference on Numerical Ship Hvdrodvnamics, Berkeley, CA., 1977, pp. 30-38. 27. "Report of the Resistance and Flow Committee," Proc. 18th Int. Towing Tank Conf., Kobe, Japan, 1987, pp. 47-92. 28. Toda, Y., Stern, F., and Longo, J., "Mean-Flow Measurements in the Boundary Layer and Walce and Wave Field of a Series 60 CB = .6 Ship Model for Froude Numbers .16 and .316," Iowa Institute of Hydraulic Research, The University of Iowa, IIHR Report No. xxx, 1990 (in preparation). 29. Sarda, O.P., "Turbulent Flow Past Ship Hulls - An Experimental and Computational Study," Ph.D. Thesis, The University of Iowa, Iowa City, IA., 1986. 30. Kajatani, H., Miyata, H., Ikehata, M., Tanaka, H., Adachi, H., Namimatsu, M., and Ogiwara, S. "The Summary of the Cooperative Experiment on Wigley Parabolic Model in Japan," Proc. 2nd DTNSRDC Workshop on Ship Wave-Resistance Computations, 1983, pp. 5-35. Table 1. Residuary-Resistance Coefficients L(m) T(°C) UO(m/s) Fr Re CR ExperimentIHI 6 12.8 2.423 0.316 11.9 x 106 1.803 x 10-3 Expenment SRI 4 10.6 1.978 0.316 6.14 1.998 Experiment UT 2.5 17.3 1.564 0.316 3.6 1.866 Inviscid ~0.316 -- 1.79 Interactive ~0.316 4.5 x 106 1.92 710

DIVERGING WAVE 80W WAVE '~ST E R N WAVE REGION 2. BODY SURfACE, Sb / / TRANSVERSE / / Bow FLOW yt ~ \ /| / WAVES ~ ~ Uo ~ /4~MMETRY PLANE, Sk x S- >- ~\l~ ---------I (INTERACTION} REGION 3 ~ ~ ~ 1 ~ S. \ THIN BOUNDARY LAYER \ / o REGION 4. ~tINTERACTtoN, THICK BOUNDARY LAYER ~ WAKE REGION 1. INVISCID FLOW INLET BOUNDARY, S. (LARGE DOMAIN: EXIT BOUNDARY, Se OUTER BOUNDARY, S (LARGE DOMAIN) f (a) (x,y)plane (a) Fr = 0 O. (b) Fr=.316 Figure 3. Displacement bodies. a' REGION 1. INVISCiD FLOW ARC I' \ So ,/ _ {lNTERACTlON) / OUTER BOUNDARY, So - ,/ (LARGE DOMAIN) SYMMETRY PLANE, Si REWATER PLANE, -W (b) (y,z) plane Figure 1. Definition sketch of flow-field regions and solution domains. ° VATERPLANE o _ . ~ . - o to ~- Fr=o - F`r=0.316 %~% - %% _ ) 50 100 150 200 250 300 350 400 44 to Global Iteration Number Figure 4. Convergence history. (a) - Fr=0.318 -- i'r=O:luLeracLive O ['r=O:Large-Doulain O t:xp.-l'stmuff & J. _9~ _ _ _ _ ~_ _ _ _ 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 X (a) longitudinalplane 1=7 X= 0.~ff 1=~D X= 1.002 l 0.6 0.8 1.0 .2 - F.r=0,316 -- F`r=O:lnteractive O Fr=O:Laree-Domain O Esp.-Ylatmuff dc J. _~ , , , D.4 0.5 0.6 0.7 0.8 0.9 1.0 (C) - FY=O.SIG -- F`r=O:lDteractive . O Fr=O:Large-Domain O Exp.-Sarda - I , T , 0.00 0.05 0.10 0.15 0.20 0.00 0.05 0.10 0.15 0.20 Y Y (b) body and wake crossplanes Figure 2. Computational grid. 1.50 1.75 2.00 X Figure 5. Comparison of interactive and large-domain solutions for the watelplane: (a) surface and wake centerplane pressure, (b) wall-shear velocity; and (c) wake centeIplane velocity. 711

(a) interactive ~111111111111111111111 (b) inviscid Figure 6. Free-surface perspective view. ..~ 0.2 0.4 O.B 0.8 1.0 X (a) z/d=.04 -- Present cal.- Inviscid - Present cal.-VjBCOUB O Bare body - Inviscid O Exp.-IHI "o~ODOOODOaC -- Present cal.- Invincid - Present cal.- Viscous O Bare body - loviscid O Esp. - IHI [ ~& l l ).0 0.2 0.4 0.6 0.8 X (b) zld = .92 Figure 9. Surface-pressure profiles -0.2 0.0 O.Z 0.4 0.6 0.8 1.0 1.2 1.4 1.6 X (a) interactive - or~ ~, o = 0.001 " X (b) inviscid Figure 7. Free-surface contours. l - Disp.body - O_Fr - - Disp. body - Non_O_Fr O Bare body O Esp. ~ =~__= =~; I ~ 1 ~ 0.0 0.2 0.4 0.6 O.B 1.0 X Figure 8. Wave profile. ·c, ~0 1 ~ 0.0 0.2 0.4 0.6 X (a) expe~iment .. 0.0 o.z ~ . ~ . : : . . . . 0.8 1.0 0.4 0.6 0.8 X (c) inviscid Figure 10. Surface-pressure contours. .015 1.0

Interactive (a) O_ o U~ o o - (C) O o (d (e) o . o o C~ O o o - U] o o ... ,. o - o O ~ Inviscid Fr = 0 Fr = .316 Fr = .316 F.: Do.soo (0.900 ,. l l 0.00 0.05 0.10 (b o ,f- ' .. o_o:oi5 ..- o-o,olo , I ~ U.05 0.10 0.15 0.20 . , ," . ,' ~ , ~ ,.' .' ' , ~' 0.1 . r , , O.00 0.05 0.10 0.15 0.20 \~J/~ 0.00 0.05 0.10 . . 0.00 0.05 0.10 `~de /' ,- .,.' ' b,_'o,,05o., / . b, 0 02S. /-'''' '- · - -o~o;OZO 0.00 0.05 0.10 0 15 0;20 | -- .-n n2n ~-0,030 . . . 0_o.076 . . , 1 0.1 r- , ', , , 0.00 0.050.10 0.15 0.20 0.000.05 0.10 O.15 0.20 y 0.00 0.05 y Figure 11. Comparison of solutions at x = .506: (a) axial-velocity contours; (b) pressure contours; (c) crossplane-velocity vectors; (d) axial-vorticity contours; and (e) turbulent kinetic energy contours. 713

(o) 1 all (b) ~ o ~ . - ,/ ~ O a a O . - ~. . O . . . . 0.000.05 0.10 0.15 O.ZO o.eo 0.05 0.10 0.15 O.ZO (C) an_ _ (d - o 0 . . n an n no n In . . . O.UO 0.05 0.10 Y Disco Or = 0 ~ = .316 ~ = .316 ..... . . . . O.00 0.()5 0.10 0.15 O.3O / o O_ o 6- ~ O_ I 1 0.00 0.05 0.10 0.00 0.05 0.10 0.15 O.ZO , .. . . . . . . 0.00 0.05 0.10 0.15 0.20 0.00 0.05 0.7 O.l5 O.ZO ~0 Y Fugue 12. Carson of spurns ~ x = .9~: (a) ~i~-vel=1~ cone; (b) Is cares; ~) s~l=~-ve~1~ vow; (d) =1~-v~c1~ congas; ad (e) opulent anemic entry congas. 714

(a) o o o ~- L~ o o o - O_ _ (b) o o o o o - o in - o Interactive _. Inviscid Fr= 0 Fr = .316 Fr = .316 7J Al= J o 1 1 0_ ~ 0 t1() 0.05 0.10 0.00 0.05 0.10 _._ _-~005 it// _ / / i 15 / .0-10 I I ~I 0.00 0.05 0.10 0.15 0.20 0.00 0.05 0.10 0.15 0.20 (C) o o . U: _ . . _ (d) (e) 0.00 0.05 0.1 0 (). 15 C).''0 - 1 0.00 0.05 0.10 y . ~1 . .. 0.00 0.05 0.10 0.15 0.20 0.00 0.05 0.10 ().15 0. 'U y 1 0.00 0.05 0.10 y Figure 13. Comparison of solutions for x = .1.112: (a) axial-velocity contours; (b) pressure contours; (c) crossplane-velocity vectors; (d) axial-vorticity contours; and (e) turbulent kinetic energy contours.

l - - CAL. [Fr=0.316] -- CAL. [Fr=O] EXP.-Sarda [Fr=O] EXP.-Ali S. [Fr=0.313] o o _ ~Z= 0.0000 o _ ~~~~~ ~Z= 0.0125 0 ~ o o o o c ,~-Z= 0.0250 Q':;~ z=0.0375 i~£~ z= 0.0500 _- ~ Z=0.625 0.00 0.01 O.OZ 0.03 0.04 0.05 y o u~ 0 CAL. [Fr=0.31~;] -- CAL. [Fr=O] O EXP.-Sarda [Fr=O] EXP.-Ali S. [Fr=0.313] ~_ -~- Z= 0.0000 0.0250 Z= 0.03~5 Z= 0.0500 ~X: Z= 0.0625 CAL. [Fr= 0. 316] -- CAL. [Fr=O] O EXP.-Sarda [Fr=O] ~ EXP.-Ali S. [Fr=0.313] 0 o _ ~Z= 0.0000o ~c ~z= 0.01250 . ~o ~~-°~~~~~~ ~Z= 0.02500 e ~Z = 0 . 0 3 7 50 __ . ~:A Z= 0~0500 0 Z=0.0625 0 cz 1 1 1 1 O 0.00 0.01 0.02 *0.03 0.04 0.05 y CAL. [Fr= 0 . 316] - - CAL. [Fr=O] O EXP.-Sarda [Fr=O] X EXP.-Sarda:CORRECTED o ~ EXP.-Ali S. [Fr=0.313] o o - o 0 ~0 0 .~_ Z= 0.0000 o .~ ~ ,,< Z= 0.0125 o o _~__ ~Z= 0.0250 o 0 0 0 ~ ~_0_ ~Z~ 0.0375 0 0 ~Z= 0.0500 0 >~ 0~__, --~w.~- 0 ~Z= 0.0625 o ~ ~ ~ ~ ,~, . O ~ ~=';7941~-W-' ~ ~0 O 1 1 1 1 0_ 0.00 0.01 O.OZ .0.03 0.04 0.05 y CAL. [Fr=0.316] - - CAL. [Fr=O] O EXP.-Sarda [Fr=O] ~' Z= 0.0000 Z= 0.0125 ~~ Z= 0.0250 Z= 0.0375 ~ Z= 0.0500 ~ Z= 0.0625 1 1 1 1 0.00 0.01 O.OZ *0.03 0.04 C) 05 Figure 14. Velocity, pressure, and turbulent kinetic energy proD~les at x = .5. 716 0.00 0.01 0 02 0.03 0.04 0.05 y

- - In - - CAL.tFr-0.316] -- CAL. [Fr-O] O FXP.-Sarda [l`'r=O] EXP.-Ali S. [Fr=0.313] ~, V9~-Z=O.OOOO - 1, _ Z=0.0125 ~Z=0.0250 :;~ Z= 1).0375 ~, :~Z= 0.0500 _ ~ Z= 0.0625 O- 1 1 1 1 0.00 0.01 O OZ *0.1)3 0.04 0.05 oF 1 o o 0 of O F o o CAL. [Fr=0.316] -- CAL. [Fi=O] O EXP.-Sarda [Fr=O] X EXP.-Sarda:CORRECTED EXP.-Ali S. [Fr=0.313] ~Z= 0.0000 ~; o - o o o o c~ ~~~~~ Z= 0.0125 ~-Q-^ ^~ Z= 0.0250 ~Z=0.0375 ~Z= 0.0500 - z= 0.0625 (}~ 1 1 1 1 0.00 0.01 0.02 *0.03 0.04 0.05 y - CAL. [Fr=0.316] -- CAL. [Fr=O] O EXP.-Sarda [F'r=O] EXP.-Ali S. [Fr=0.3 13] 1\ _ ~a~ Z= 0.0000 Z= 0.0125 Z= 0.0250 ~9~ Z= 0.0375 Z= 0.0500 9~ ~ Z= 0.0625 ,~ ~ 0.00 0.01 0.02 *0.03 0.04 0.05 y CAL. [F'r=0.316] -- CAL. [Fr=O] O EXP.-Sarda [Fr=O] Z= 0.0000 - CAL. [Fr=0.316] -- CAL. [Fr=O] O EXP.-Sal da LFI =()] EXP.-Ali S. [Fr=n.313] . ~Z= 0.0000 __ _ = Z= 0.0125 __~ Z= 0.0250 Z= 0.03~5 - ,~4Z= 0.0500 Z= 0.0625 . ~5L=1 1 O.00 0.01 0 02 0.03 0.04 0 05 * y ~~ Z= 0.0125 -Z= 0.0250 ~Z=0.0375 = 0.0500 Z= 0.0625 1 1 1 0.00 0.01 0.02 .0.03 0.04 0.05 y Figure 15. Velocity, pressure, and turbulent kinetic energy proD~les at x = .9.

o . o o -- ar o o( - o o o - o o _ - N o o C'AL. [Fr= 0 . 3 1 6] - CAL. [Fr= 0. 31 6] -- CAL. tFr=O] -- CAL. [Fr=O] ': EXP.-Sarda [Fr=O] O EXP.-Sarda [Fr=O] | i\ EXP.-Ali S. [Fr=0.313] | ~A EXP.-Ali S. [Fr= .313] | ~)Z= o.oooo I r =c Z= .0000 ~O. Z=O.OIZ5 1 °~ Z= .0125 1 ° =oo250 1 0~:z= .0250 1 °~ :Q ~, z= 0.0375 1 °~=~\ z= 0375 1 ° ~ ~:~ ~--bo Z = 0 . 0 5 0 0 ~° ~ = / ~ ~ Z = . 0 5 0 0 ~ _ ~ ~ Z= 0.0625 o- ~=~ ~ Z= 0.0625 o 0.00 0.01 0.02 0.03 0.04 0.()5 0.00 0.01 0.02 0.03 0.04 0.05 * * Y Y o - 0 ~ o o o o o o P~ o o o 1 ~4 o o o o o o o o o o o o 1 1 1'1 1 ° 0.00 0.01 0.02 *0.03 0.04 0.05 y - CAL. [F'r=0.316] -- CAL [Fr=O] O EXP.-Sarda [Fr=O] X EXP.-Sarda:CORRECTED ~ EXP.-Ali S. [Fr=0.313] ~ Z= 0.0000 Z= 0.0125 ,:^,~ ~ Z= 0.0250 _ _ e _ _ _ _ _ _ _ ~ Z= 0~0375 ~Xz= 0 0500 A Z= 0.0625 - CAL. LFr=0.3 1~3] - - CA L. lFr= 0] O EXP.-Sarda tPr=()] EXP.-Ali S. [Fr=0.313] Z= 0.0000 ~= Z= 0.0125 ~= Z= 0.0250 Z= 0.03~5 ~3~Z= 0.0500 ~_ _Q ~ ~ - A Z= 0.0625 . . . . 0.00 0.01 0.02 *0.03 0.04 0.05 y CAL. tFr=0.316] -- CAL [Fi=O] O EXP.-Sarda [Fr=O] ~Z= 0.0000 Z= 0.0125 Z= 0.0250 Z- 0.0375 ~~ Z= 0.0500 - 11~ - O 0625 Z- . 0.00 0.01 O.OZ *0.03 O.04 0.05 y Figure 16. Velocity, pressure, and turbulent kinetic energy profiles at x = 1.1. 718

DISCUSSION Kuniharu Nakatake Kyushu University, Japan I appreciate your interesting paper. In your calculation, the flow field and pressure distribution around wave surface are not taken into account. If you distribute source panel on the calculated wave surface and determine its strength from the kinematic condition on the wave surface, you may obtain more plausible results. This is possible in the linearized framework. AUTHORS' REPLY We agree that a higher-order treatment of the free-surface boundary conditions would be preferable. We hope to make progress in this area in the future. DISCUSSION Hoyle Raven Maritime Research Institute Netherlands, The Netherlands This valuable paper addresses the difficult problem of prescribing free-surface boundary conditions inside the viscous domain. The authors solve the wave elevation from integration of the kinematic condition, and thus from the velocities at the undisturbed free surface. Prescribing boundary conditions for these velocities here is, therefore, critical. If I understand it correctly, in (20) not only the stresses but also the associated strain rates have been neglected, which are of leading order in the wake elevation. This is not a definite solution of course. What alternatives for the free-surface condition do you consider? AUTHORS' REPLY The approximations (a) through (d) utilized in deriving equations (18) through (20) have been clearly stated in the text, and therefore do not require further explanation. An alternative treatment of the free- surface boundary conditions within the present overall framework is to retain more terms in the approximations (a) through (c) and higher- order terms in approximation (d). DISCUSSION William B. Morgan David Taylor Research Center, USA In the authors' presentation of their problem they stated that they were using an "Unsteady RANS Code using a Kit turbulence model." I believe this statement is inconsistent in that a RANS Code with a kin model cannot be unsteady. I can understand attempting to use this code in a Quasi-steady" sense, but I believe it is not applicable to the unsteady problem. Would the authors please comment? AUTHORS' REPLY Although not discussed in the present paper, we recognize the limitations of the kit turbulence model with regard to simulating unsteady flow. Issues concerning this point were discussed in one of our earlier papers on another topic [31]. In our presentation, we only wished to emphasize the general capability of the IIHR basic viscous- flow method for unsteady flow notwithstanding the limitations of current turbulence models for such applications. [31] Stern, F., Kim, H.T., Patel, V.C., and Chen, H.C., "A Viscous-Flow Approach to the Computation of Propeller-Hull Interaction,. Journal of Ship Research, Vol. 32, No. 4, December 1988, pp. 246-262. DISCUSSION Kazu-hiro Mori Hiroshima University, Japan 1. It must be important to be of the same order in approximations when the viscous effects are taken into account in the free-surface computation. From this standpoint of view, the use of the displacement thickness method is not consistent where the thickness is calculated exactly. This may be crucial when the 3-D separation is dominant. My suggestion is that the viscous flow should be taken into account as the double model flow directly in the inviscid computation. 2. According to our computation and experiment, the separation of the flow at stern is much affected by the bow wave system. This means that we should not expect precise discussions on the interaction between the viscosity and the free surface by the iterative procedure as done in the present study. AUTHORS' REPLY We are not clear as to the precise meaning of your questions; however, we would like to point out that viscous effects have been included directly in the inviscid-flow computation through the use of the displacement body. As discussed in the text, the limitations of this approach are not yet known. The present results are encouraging, but a complete evaluation requires further validation through comparisons with experimental data and a large~omain approach. Work along these lines is in progress. We thank both the oral and written discussers of our paper for their pertinent remarks. 719