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.

Nonlinear Ship Waves Y.-H. Kim (David Taylor Research Center, USA) T. Lucas (University of North Carolina-Charlotte, USA) ABsTRAcr A boundary element method is presented for solving a nonlinear free surface flow problem for a ship moving a constant speed on a calm sea. The method of solution is based on distribution of simple Rankine type~llr) sin- gularities on the body and on the true location of the free surface, which must be obtained by iteration. The key is the iterative algorithm for determining the free surface and wave resistance using a new one parameter family of up- stream finite difference methods A verification of numeri- cal modeling is made using Wigley hull and validity of computer program is examined carefully by comparing the details of wave profiles and wave-making resistance with Series 60, CB =0.60 model. =RoDucrIoN Many of problems facing engineers involve such difficulties as a nonlinear governing equation and nonlinear boundary conditions at known or unknown bound- aries.The determination of surface ship waves is a non- linear problem wherein the governing equation is linear but there is a nonlinear free surface condition imposed at an unknown boundary. In practice the solution is approxi- mated using numerical techniques, analytical techniques and combinations of both. Foremost among the analytical techniques is the systematic method of perturbations in terms of a small parameter. In the past, problems of this type were generally treated after the boundary condition on the free surface was linearized. Limited studies on nonlinear boundary condi- tions have been conducted for simplified two-dimensional cases. With the recent advent of supercomputers and numerous new techniques in computational fluid dynamics there is growing interest in solving the three-dimensional nonlinear free surface problems by various schemes. Dommermuth and YueL11 reported their study on the non- linear three-dimensional gravity waves created by a moving surface pressure distribution using a higher-order spectral method. During the 5th International Conference on Numerical Ship Hydrodynamic a significant amount of work was been reported on this subject including Jensen[2], Musker[3l, and Bai[41. SSPA's series of reports by XiaL5l, Nit6], and KimE71 also deal with this topic. 439 Many different numerical schemes have been pro- posed. The method proposed here is a boundary element method with Rankine type source as the kernel function. The drawback to this approach is, as well known, the necessity of proper numerical implementation of the radi- ation condition required since the potential flow model contains no viscosity. The upstream finite difference oper- ator first introduced by DawsonE8] successfully satisfies the radiation condition by producing an inherent local damping of the free surface waves by unwinding. Since then numerous researchers have tried to improve Dawson's method and to extend it to the nonlinear free surface problem. SWIFT (Ship Wave Inviscid Flow Theory) is a computer code developed by the authors and S.H. KimE9] based on the linearized free surface con- dition. SWIFT uses a higher order panel method with either constant or linearly varying singularity strength distribution across a curved panel . The stability and accu- racy of the computational results have been examined and found to be quite robust. In the present paper the SWIFT approach is extended to the nonlinear free surface problem and, though devel- oped independently, is very similar to Nit6] and Kiml7] of SSPA. The convergence of the method has been major focal point among researchers who have been working on the nonlinear problems. XiaL5] experienced severe con- vergence problems during iteration. NiE6] used a a higher order method and test results were converged if an under relaxation coefficient was used to modify the wave eleva- tion and velocity distribution for the next iteration. KimE7] obtained convergent solutions successfully by satisfying the kinematic and dynamic boundary conditions simulta- neously without a relaxation factor. They all used the three point finite difference operator instead of Dawson's four point operator. Musker[3] extensively studied stability and accuracy of a nonlinear code, and experienced the usual difficulties in convergence with the four point scheme. In this paper we introduce a new one parameter family of four point schemes which can give an arbitrarily large damping. The four point method of Dawsont8] and three point method of Xia[5] are special cases of this approach. We have experienced greatly enhanced stability for all tested ship hullforms over a ~broad range of speeds. Wigley mathematical hull is used to verify the numerical modeling of the problem. The measured wave profiles and wave-making resistance of Series60, CB=0.60 are used to validate the developed computer program.

MATHEMATICAL FORMULATION We consider a ship which moves in the positive x- direction with constant forward speed UO in Figure 1. Let Oxyz be a Cartesian coordinate system with Oz opposing the direction of gravity and z=0 coincides the undisturbed free surface, the positive x-axis in the direction of the ship's forward velocity. The fluid is assumed to be inviscid and incompressible and it's motion is irrotational such that the velocity field of the fluid Vcan be defined as v~x,y,z) = V¢(x,y,z) `1' where ¢(x,y,z) is the velocity potential and satisfies the Laplace equation v2¢ = o (2) in the fluid domain D and the boundary condition An =Uon (3) on the body surface SO where n=(nX,ny,nz) denotes the outward unit normal vector on the boundary. On the free surface z = ((x,y), where ~ is the free surface elevation, the kinematic boundary condition and dynamic boundary conditions can be given respectively as ~x~x+¢y~y-¢z = 0 ~4' go + ~ (V) ~ V0- Uo2) = 0 (5) where g denotes the gravitational acceleration. Combining both dynamic and kinematic conditions Eqs. (4) and (5) becomes V¢eV:2(V¢~23 + gig = 0 (6) Finally, energy consideration requires that velocity poten- tial approaches the uniform onset flow potential and that there be no waves far upstream of the ship, and that waves always travel downstream. The problem described in equations (1) thru (5) is nonlinear, since the free surface boundary conditions(4) and (5) themselves are nonlinear and should be satisfied on the true free surface which is unknown and should be obtained as a part of solution. In their exact form the problem is difficult to solve. In order to be able to treat the problem, the equations are approximated by ones which are more tractable. Following linearizations are similar to NiL6] and Kim(~71. First we define D~o,() = Add + tY(Y ~ (z D2~,,() = go + 2(Vt V) - Uo2) We denote ¢° and (° the velocity potential and the wave profiles, respectively, obtained from the previous step and introduce two small parameters 6o and b: with respect to the previous solution. Take Taylor series expansion of D~(c,,() and discard higher order terms greater than Act and S(, then the kinematic free surface condition becomes D~o,() ~ D~°,(°) + ha D,(o,(°~6c; + a3` D~(c,°,(~ =40~0 +~0~0_40 + 6¢xCx + bPy~y ~ Adz + ¢°~( + ¢°~` +~°xz(°x + (yz(Y ~ (zz)~( (7) We redefine a new velocity potential ~ = 0° + ~ and the linearized kinematic boundary condition can be expressed as follows: D~o,() ~ (x~x + (Y(Y - (z + {x6(x + (y6(y + (¢xz (x + tyz By ~ (zz )~( = 0 ~8' Similarly the dynamic boundary condition can be linearized as follows; D2(c',() ~ 6` (1 + g(¢xtxz + (Y¢YZ + (zfzz)) + ~ - ~ {Uo2 + ¢° + ¢° + ¢° - 2~°¢X + ¢°¢y + LIZ LIZ)} = 0 (9) Now the new linearized boundary value problem is for- mulated. The singularity distribution and wave profile are continuously updated until the solution converges and satisfies original Eqs. (4) and (5~. Once the singularity strength distributions are determined, the wave resistance can be computed by Rw = J~J. pnx4S SO (10) where SO is the wetted ship hull surface and the fluid pressure p is given by the Bernoulli equation p = PI ~v¢~2 u2~ pgz NUMERICAL SCHEME (1 1) The method used here is a boundary element method in which simple Rankine sources are distributed across each panel. The boundary-value problem formulated above then reduces to a determination of an unknown singularity distribution over the boundary surface of the fluid domain. Once the singularity distribution is determined, the hydro- dynamic quantities of interest, the velocity and the pres- sure, are also determined. 440

The present method uses a panel scheme to achieve a numerical solution to the problem of Eqs. (1) thru (9~. Panel schemes proceed by first dividing the boundary surface into panels. Source distributions are assigned to each panel. These distributions are expressed in terms of unknown singularity parameters Hi associated with the panel and neighboring panels. A finite set of control points (equal in number to the number of singularity parameters) is selected at which the boundary conditions are imposed. The construction of each network requires numerical development in these areas: A. Surface geometry defi- nition; B. Singularity strength definition; and C. Control point selection and boundary condition specification. Essential features of the computational scheme in each of these areas are as follows A. Geometry input for a network is assumed to be a grid of corner point coordinates partitioning the network surface into panels. Panel surface is obtained by fitting a paraboloid to corner points in an immediate neighborhood by the method of least squares. B. Discrete values of singularity strength at certain points on each network are assigned as singularity pa- rameters. Singularity splines are constructed for each network by fitting a linear distribution on each panel of the network by the method of least squares. C. Certain standard points on each network are assigned as control points. These points include panel center points as well as network edge points. Boundary conditions involving the specification of potential or veloc- ity are applied at panel center points for the purpose of controlling local properties of the flow. Curved Panel Approximation We use curved panels to patch the ship and free surface. Let S be a true surface element bounded by four corner points four corner points P~,P2,P3, and P4 as shown in Fig. 2. We assume that S may be represented in the approximate form All, rig = (0 + (~ t + (~ rat + 2 (~ ~2 + (~ ~q + 2 (~2 where ~ (,q,() are orthogonal coordinates local to S. The six coefficients ~ (0, Hi, All, Aid, (an, (~ ~ are obtained by requiring that the approximate surface given by Eq. (12) pass through its four corner points exactly and through its neighboring points approximately in a least square sense. That is, for each panel the expression R = -~W,`(~,31~) - (,`) 2 ,` (13) is minimized with respect to the six unknown coefficients appearing in Eq. (12~. The summation in Eq. (13) ranges over the N neighboring grid points. The choice of the N points and the weight Wk has been made in an attempt to .. . . .. minimize any 1rregu arltles t net may appear in t ne paraboloid approximation on the true surface. In the present formulation N=16 and Wk=1 or 108 are used depending on the location of the panel in the surface panel network. Here we used the term network to represent the collection of quadrilaterals. The choice of these features is required by the argument given by Hess[10] that certain computing methods which use flat panels do not obtain increased accuracy over the basic zero order method~constant source on flat panels) even though they use higher order source distributions. For purposes of computational efficiency Eq. (12) can be reduced by an appropriate coordinate transformation to the canonical form ~=at2+bll2 (t,11)~I; (14) where ~ is the quadrilateral formed by the projection of the corner points on the (mu) plane. In the case of flat panel, a and b are both zero. Singularity Strength Definition The true singularity distribution will be approximated by a truncated Taylor's series on each panel. Such a rep- resentation is valid on any interior part of the network providing the paneling there is sufficiently fine. We will consider a linearly varying singularity distribution on source panels. There may be an advantages in using even higher order distribution, but as pointed out by Hess[10], it would be necessary to consider a higher order panel geometry definition for the sake of consistency. Specifically, we assume that the singularity strength ~ at a point ~ A, 11, ~ ~ on a panel S is given by o(t,11) = co + Cl`t + still (15) The coefficients in Eq. (15) are determined by using a method of weighted least squares over the panel and up to eight of its neighboring panels. This method requires that the form of Eq. (15) gives exactly the singularity value at its centroid and approximate values in the least square sense at centroids of the neighboring panels. Induced Velocity Potential The velocity potential at P=(x,y,z) induced by a singularity distribution on S is given by ~ = - 4 J| -ds (16) where r = ~ (~-x)2+~11-y)2+~-z)2 ~ I/2 and ~ is given in Eq. (15~. We are now in a position to evaluate the integral Eq. (16) using the relations of Eqs. (13), (14), and (15~. The ~ can be expanded in closed form by using a com- bination of closed form calculations and recursive rela- tionships. Details are fully included in Johnsonfll] and Kim et al.~91. Derivatives of velocity component along z- direction in Eqs. (8) and (9) are new and their final closed forms are summarized in the Appendix. METHOD OF SOLUTION As usual we discretize the fluid boundary into a finite number of panels. Because of symmetry, half of the ship as well as the free surface domain is used as an input. A typical free surface domain is truncated at a half ship length(L=2.0) ahead of the ship and is extended more than one ship length aft of the ship and a half ship length or more to the side. Since the convective term in the free surface condition involves x and y directional derivatives,

we use an algebraically generated waterline-fitted coordi- nate system, independent of the double model streamline coordinate, to create the panel network. As shown in Fig. 3 the transverse lines are straight and vertical to the center line of the ship. Our numerical experiments indicate that the solution becomes more stable and accurate when the aspect ratio of the free surface panel away from the ship closes to unity. The same finding was discussed by Musker(1989~. During iteration, the panels are rearranged such that the projection on the horizontal plane z=0 is unchanged, i.e. x and y coordinates are fixed and the z-coordinate is allowed to move up and down. The method of solution requires that the linearized free surface conditions (8) and (9) are satisfied for a given 0° and (° which were obtained from the previous step. The iteration starts from the solution which satisfies the well-known linearized free surface condition. In each iteration the free surface panels as well as body surface panels are adjusted according to the newly computed wavy surface and boundary con- ditions are applied to the updated boundaries. The itera- tion continues until the solution satisfies the exact free surface conditions. In practice, we provide a certain convergence criteria in the computation that requires the iteration to continue until the residual error satisfies the provided tolerance . This part will be discussed more extensively with numerical results later. One Parameter Family of Advection Methods to Enhance Convergence In this section we derive a one parameter family of advection methods to enhance convergence for the non- linear free surface problem. A number of researchers have been concerned with the question of convergence of various forms of the nonlinear method. Xia had severe convergence problems with his approach and recom- mended the three point finite difference operator instead of Dawson's more accurate four point operator. Ni and Kim have developed more powerful methods. Recently Musker has completed a detailed stability and accuracy test of a nonlinear code, and experienced the usual difficulties in convergence with the four point method. Musker used two approaches to enhance convergence: he used a spline scheme involving four upwind points and he elevated the free surface panels. His nonlinear iterations, unlike those of Ni and Kim, are all on the flat surface. In this paper we introduce a new one parameter family of four point advection schemes which can give an arbitrarily large damping. The four point method of Dawson and the three point method are special cases of this approach, and without the use of elevation of the panels (which we have not examined) we have experienced gready enhanced stability. We recommend the inclusion of this technique into existing codes. For most programs this would require only a minor modification of a few lines of code, as will be seen below. The four point method is: Slopei = -( A4fi + B4fi ~ + C4fi-2 + D4fi-3 ~ (17) where for an uniform mesh spacing of size 1, A4 = 1.6671 1, B4 = -2.51-1, C4 = 1-1, D4 = -0.1671-1. The related three point scheme is of the form Slopei = -( A3fi + B3fi-l + C3fi-2 ~ (18) where for an uniform mesh spacing of size 1, A3 = 1.5l-l, B3 = -2.01-1, C3 = 0.51-1. In considering all possible four point schemes for approximating f ' one might like to require that they be exact for quadratics. These three restrictions on a formula of the type Slopei = -( ANt + BNfi-l + CNfi-2 + DN0-3 ~ (19) leave only one degree of freedom and since both the three point and four point methods satisfy a quadratic exactly, the general derivation of this general fam~ly is very simple: AN = ~ 1 - Qmul ~ A4 + Qmul A3 BN = ~ 1 - Qmul ~ B4 + QmulB3 CN = ~ 1 - Qmul ~ C4 + Qmul C3 DN= (1-Qmul ~D4 (20) where Qmu~ is a parameter to be chosen greater than or equal to zero. The case Qmu~ = 0 gives the familiar four point formula while Qmu~ = 1 returns the more stable three point formula. We have found that Qmul = 1 is often fine for lower speeds but for higher speeds Qmul = 3 (or even sometimes Qmul= 5 ~ eliminates most of the problems with convergence for reasonably regular com- putational grids and with little loss of accuracy. To analyze the effect of the one parameter method and to get some feel for the meaning of the parameter Qmul, it is instructive to examine the coefficient of fi"' which is the leading term in the Taylor expansion of the error for both the three point and four point formula. (Recall that the basis of the four point formula was not that the fi"' term was eliminated but rather that the fi"" term was required to have a zero coefficient.) Here we will take an uniform spacing of size 1. The coefficient of fi"' in the three point method is of the magnitude 12/3 while in the four point method it is 12/6, leading again to the well known conclusion that the superior damping of the three point method will frequently enhance conver- gence over that of the four point method, an important finding of Xia for example. Thus the one parameter Qmu~ method has a coefficient for fi"' of the magnitude 12( (1 - QN,UI)/6 + Q~ /3 ~ = 12( 1 + Ql~lUI)/6 giving again 12/6 and 12/3 for Qmu~ = 0 or 1 as before. When denser meshes are used, or for large speeds where f"' may be smaller, the us~e of larger values of Qmul such as Q~ul = 3 or 5 gives the (larger) coefficients of fi"' of magnitude 2/312 or 12 giving in these cases twice or triple the damping of the 3 point method. We now examine briefly the spline scheme consid- ered by Musker, which for an uniform mesh has the formf refer to (201) with A4 = 1.5551-l, B4 = -2.1771-1, C4= 0.6891-l, D4 = -0.06671-~. This method does not quite make the quadratic term zero, so is not included in the above family, but it is very close with a value of Qmul of about 0.6. Thus it has much of the improved stability associated with the three point scheme, giving a coeffi- cient of f" of the magnitude 0.2712, closer to the 0.3312 of the 3 point method than the 0.16612 of the 4 point method. Thus we would be led to conjecture that the 3 point scheme would be more stable than the spline scheme, at negligible loss of accuracy. RESULT AND DISCUSSION To facilitate comparisons, the following non-dimen- sionalization has been made: 442

Cw~wave-making resistance) - RWI(1I2 pU2S) ((wave profile) = ~ /(L/2) = ~(U2/~2g)) (used in Figures) Although the higher order panel method gives better results, we used linearly varying source strength distri- bution across curved panel for the ship hull and constant source strength across flat panel for free surface during computation because of the following reasons: (1) The higher order panel method requires a much longer computing time, normally 3 or 4 times longer with the same number of panels. (2) The finite difference scheme used for the convective terms in the free surface condition requires such small panel sizes that the advantages of higher order panel scheme cannot be fully utilized over the free surface domain. Numerical experiments have indicated that the differences in results are negligible between using higher order panel and flat panel on the free surface domain. On the other hand, a Neumann condition imposed on the hull surface does not require any numerical differentiation and the higher order panel method gives better results than flat panel result with even fewer panels. Furthermore, the wave-making resistance is small quantity and to get this quantity accurately the integration of the pressure distri- bution over the entire wetted hull surface has to be carried out precisely. It is therefore required that the pressure distribution as well as panel surface be precisely approx- ~mated. Qmul Test with Wigley Hull Wigley parabolic hull was used to verify the pro- posed one parameter family of advection methods. The following computational conditions were used: Wigley model L=2.0,L/B=10.0,D/B=0.625 Bow = +1.0, Stern =-1.0 Computation domain x= -2.0, +2.0 (S SPA recommended) y= -0.75 Number of panels 132 on ship (23 Stax7VVls) 400 on free surface (41xl 1) Froude numbers 0.316, 0.408 Qmul's 0, 1,3,5 Computer Cray X-MP/24 We used the computational domain that XiaL5] and other SSPA reports recommended and also matched the panels on the ship hull and free surface. We selected higher Froude numbers because several researchers have experi- enced convergence problems particularly at higher speeds. Table 1 shows the computed wave making resistance using four different Qmul's: 0, 1, 3, and 5. The solutions diverged when Qmul = 0 and 1, and converged when Qmul=3 and 5. The experimental values of wave-making resistance coefficients were reported at the 17th I1lCE12] and Cw ranges between 1.7x10-3 and 2.3x10-3 at Fn=0.408. Qmul=0 is the four point scheme and the solu- tion diverged after 2 iterations. At Qmul=l, the three point scheme, the solution diverged after 5 iterations. Qmul=3 required only 4 iterations and Qmul=5 required just 2 iterations to get a converged solution. The magnitude of the wave making resistance is slightly smaller than that of Qmul=3. It is likely that the larger damping by using Qmul=5 caused the difference. Table 1 Qmul Test (Wigley Hull at EN = 0.408) Cw(Experiment) = 2.0 linear 2.248 2.304 2.317 2.231 . 1 1st | Iter 1.562 1.890 2.009 _ 1.965 2nd Iter 1.863 2.067 1.976 _ 1.939 3rd Iter Divrg 1.683 1.982 1.939 . 4th Iter 1.726 1.978 1.940 5th Iter 1.899 1.978 1.940 Cwx103 Non- linear _ Divrg Convg Convg The two small parameters 8o and 6( used when linearizing the free surface conditions were continuously monitored during the computation. In Table 2 EITor, the summation of the square root of bt squared terms, is tabu- lated at each iteration for different Qmul's and in Table 3, that of 8~. At Qmul-l, Error oscillates during the itera- tions until the solution diverges. With Qmul=3 or 5, the magnitude of Error decreases steadily as the iteration continues. In practice, the convergence rate after several iterations is rather slow so we assumed the solution to be converged when the ratio of the current Error to the first iteration Error is smaller than 1% for both Act and by. Table 2 ~ Convergence Test (Wigley Hull at EN = 0.408) N errorl= ~)2, i Qmul o 3 5 1st Iter 8.248 4.199 . 2.955 l 2.859 2nd Ite r 4.550 . 4.493 . 1.857 0.695 3rd Iter 152.5 5.488 0.507 0.298 l N = 400 4th Iter Divrg 2.683 0.333 0.033 error1 x1 n2 5th Iter 4.975 0.036 ._ 1 Non- I linear I 1 1 Divrg | Convg | 3 Table 3 ~ Convergence Test (Wigley Hull at EN - 0.408) N error2= ~)2, J ~ 1 l ..4 443 1st 1 Iter | 1.10 1 1 0.62 0.43 3 2nd Iter 1.30 1 1.1 0 o~o9 1 ~3 3rd Iter Diverg 1 .00 0.062 0.012 N = 400 1 1 4th I Iter 1 1 0.58 I 0.017 1 l 1 5th I Iter 1 1 0.84 0.009 0.001 1 | Non | linear | Diverg | Convrg | Convrg

Fig. 4 depicts the wave profiles at the hull surface during iterations including the wave profile obtained from the linearized free surface conditions. As the iteration progresses the bow wave grows and the first trough becomes deeper with the phase shifting slightly towards the bow. Only insignificant differences are observed aft during the iterations.Figure 5 shows the wave profile comparisons between the experimental, linear and non- linear computations. The nonlinear solution is now much closer to the experimental. Series60, CB=0.60 The Series60 with block coefficient 0.60 has ex- tensive model experimental results, including wave profiles and wave-pattern resistance. Experimental values used here were determined by Kim and Jenkins[131. Experimental data only for model fixed were considered in comparisons. Wave profiles were marked at every station along hull with a grease pencil during the run and were read after each run. Wave resistances were obtained from longitudinal wave profile measurements. The following computational conditions are used: Series60 model L=2.0, L/B=7.5, D/B=0.4 Model is fixed Bow = +1.O, Stern=- 1.0 Computation domain x = (-3.5, 2.0) y = (-1.5, 0.0) Number of panels192 on ship(25 Sta x 9WLs) 696 on free surface(59 x 13) Froude numbers0.22,0.25,0.28,0.30,0.32,0.35 ComputerCray X-MP/24 Figure 6 shows the comparison of wave profiles along the hull between computed and measured. The wave profiles obtained from linearized free surface condition are also plotted together. The nonlinear solution significantly improves the wave profiles, in particular, at the bow and the first trough and after that the difference between linear and nonlinear results seems insignificant. Qmul = 1 is used for Fn = 0.22, 0.25, and 0.28, and Qmul = 3 for Fn = 0.30, 0.32, 0.35. Table 4. Wave Resistance Coefficient Comparison Series60, CB = 0.60 (Model Fixed) Fn~ 0.22 0.25 0.28 0.30 0.32 0.35 linear 0.369 0.664 1.581 2.144 1.815 2.148 1 st Iter 0.278 . . 0.398 1.123 1.548 1.401 . ~ 1.886 2nd Iter 0.271 0.357 1.049 1.474 1.333 1.800 3rd Iter 0.268 0.354 1.043 1.467 1.350 1.816 4th Iter 0.266 0.355 1.043 1.460 1.354 1.828 (CWX1 03) . Nonlin ear 0.255 0.354 1.045 1.460 1.357 1.836 Table 4. summarizes the wave-making resistance computation results at each iteration and in Figure 7 the linear and nonlinear results are compared with the experimental data. It is interesting to see that linear solution consistently predicts higher wave-making resistance throughout the Froude number range. Nonlinear results show very close agreement with the measured data. Even though all the runs converged successfully, at the lower Froude numbers 0.22 and 0.25 agreement on wave profiles as well as wavemaking resistance were rather poor. It is likely that the panel size may be too large at lower Froude number where wave length is shorter and is inadequate to resolve the flow phenomena accurately. At FN=0.22 the estimated characteristic wavelength is 0.608 which is less than 1/3 of the ship length. In usual practice between 10 and 12 points per wave length are considered to get acceptable numerical accuracy for wave making resistance. This means at least 33 evenly distributed panels along the ship hull have to be provided. No further attempt was made to increase the panel density at the lower Froude numbers, because it exceeds Cray X-MP memory restriction 3 million words. This particular run requires 2.6 million words computer memory~maximum usable capacity is 3 million+) and approximately 70 CPU seconds were used at each iteration. As seen in Table 4. the most sizable improvement was made at the first iteration and in subsequent iterations Improvement was small. In practice we may need one more iteration over the linear solution to have reasonable results. Convergence of the solution is always monitored carefully during computation. Table 5 shows one exam- ple. Table 5. Convergent Test series60, CB=0.60 at Fn = 0.32 error1 0.9147 0.1783 0.0513 0.0219 . 0.0153 error2 0.670 0.170 0.01 7 0.037 0.013 SOUrCe 2.9383 3.0646 3.1 1 57 3.1411 3.1 52 1 3.1521 _ ~ 1 .8 1 5 1 .401 1 .333 1 .350 1 .354 1 .357 Here errorl and error2 are defined in Tables 2 and 3 and NF NF is the total number of free surface control points. For example, if we average errorl by NF=696, average wave height increment ~hl at each free surface panel after the first iteration is 0.001314 and after the fifth iteration 6~5 becomes 0.000022 which are nondimensionalized by half of the ship length. If Series60 hull length is 200m, 6~1=131.4mm becomes 6~5=2.2mm after five iterations. Fig. 8 shows two pairs of the contour plots using the locally developed visualization program WAVEt143. The Froude number is 0.32. The nonlinear solution was obtained after 6 iterations. The first pair shows contours of the wave height and source strength distribution on the free 444

surface for the linear step, and the second pa* for the nonlinear solution. The wave height is plotted for the positive y-plane and the source strength for the negative y- plane for convenience. 0.03 uniform spacing is used between contour lines for both plots. Solid lines indicate positive wave height and positive source strength, dotted lines negative wave height and negative source strength, and dashed lines zero. We observe the following: (1) the source strengths on the body toward the bow of the ship~not plotted here) are large and positive, but at the free surface near the bow a locally concentrated sink dism- bution is noticed, with a strong and wide source distribu- tion following immediately downstream,, (2) at the stern there is only sink distribution on the ship surface and the free surface, (3) the humps and hollows of wave height are out of phase to those of source strength. Zeroes of wave height occurs near the mean of either source or sink distri- bution along x-axis. CONCLUSIONS 1. The introduction of a one parameter family of advection methods has significantly enhanced the con- vergence. This new scheme can control the magnitude of damping in the numerical solution, and proper choice of ~ ] may eliminate the numerical instabilities that several researchers have experienced when using either the three point or the four point method. 2. Contour plots similar to Fig.8 for example, are strongly recommended when studying such a complex numerical model. Using the contour plots we can more easily visualize the quality of computational results and choices of truncation regions and panel densities which often may be misjudged just looking the local flow charac- teristics such as the wave profile or streamline tracing along the ship hull or especially a scaler quantity such as wave resistance alone. We have found that numerical instability sometimes starts at the edge of a boundary with negligibly small magnitude, and it progressively grows and spreads as iteration continues. 3. The most significant improvement was made between the linear solution and the first iteration. It is suggested that at least for the first iteration an under relaxation coefficient either obtained from Eq. (22) or provided with small value such as 0.5 used by Nit6] be used. For practical purposes, 2 to 3 iterations may be sufficient. 4. It now apples- that the long search for programs to accurately estimate wave resistance is now drawing to a successful conclusion, and the next emphasis should be on the more difficult problems of simulating the wave height on the hull and more generally the wave heights on the free surface. Here the comparison with experimental results combined with flow visualization methods is essential. ACKNOWLEDGMENT The first author was funded by the Applied Hydromechanics Research(AHR) Program supported by the Office of Naval Research and administered by David Taylor Research Center(DTRC).The second author was supported in part by an IPA grant from DTRC and super- computer grants from the National Science Foundation (Grant Number ECS 8515174) and the North Carolina Supercomputer Center. The authors express their gratitude to Dr. Wen-Chin Lin and Dr. Michael Wilson of DTRC for their valuable suggestions and encouragement during the course of this work. The authors also thank Mr. Peter Chang and Mr. Steve Fisher for their kind support. REFERENCES 1. Dommermuth, D.G. and D.K.P. Yue, "The Nonlinear Three-Dimensional Waves Generated by a Moving Surface Disturbance," Proceedings of the 1 7th Symposium on Naval Hydrodynamics, the Hague, Netherlands, 1988 2. Jensen, G., V. Bertram, and H. Soding, "Ship Wave- Resistance Computation," Proceedings of Fifth International Conference on Numerical Ship Hydrodynamics. Hiroshima, Japan, Sept., 1989. 3. Musker, A.J., "Stability and Accuracy of a Non- Linear Model for the Wave Resistance Problem," Proceedings of Fifth International Conference on Numerical Ship Hydrodynamics. Hiroshima, Japan, Sept., 1989. 4. Bai, K.J., J.W. Kim, and Y.H. Kim, "Numerical Computations for a Nonlinear Free Surface Flow Problem," Proceedings of Fifth International Conference on Numerical Ship Hydrodynamics. Hiroshima, Japan, Sept., 1989. 5. Xia, F., "Numerical Calculations of Ship Flows, with Special Emphasis on the Free surface Potential Flow," Chalmers University of Technology, Goteborg, Sweden, 1986. 6. Ni, S.-Y., "Higher Order Panel Methods for Potential Flows with Linear or Non-linear Free Surface Boundary Condition," Chalmers University of Technology, Goteborg, 1987. 7. Kim, K.J., "Ship Flow Calculations and Resistance Minimization," Chalmers University of Technology, Goteborg, Seden, 1989. 8. Dawson, C.W., "A Practical Computer Method for Solving Ship Wave Problems," Proceedings of the 2nd International Conference on Numerical Ship Hydrodvnamics, Berkeley, Calif., 1977. 9. Kim, Y.H., S.H. Kim, and T.R. Lucas, "Advanced Panel Method for Ship Wave Inviscid Flow Theory," DTRC-89/029, 1989. 10. Hess, J.L., "Consistent Velocity and Potential Expansions for Higher Order Surface Singularity Method," MDCJ 6911, pp.l-30, June, 1975. 11. Johnson, P.T., "A General Panel Method for the Analysis and Design of Arbitrary Configuration in Incompressible Flow," NASA Contract Rpt 3079, Boeing Com. Airplane Co., Seattle, 1980. 12. Norrbin, N.H.(editor), "The Proceedings of the 17th International Towing Tank Conference," SSPA, Goteborg, 1984 445

13. Kim, Y.H. and D. Jenkins, "Trim and Sinkage Here Effects on Wave Resistance with Series60, CB=0.60," DTNSRDC/SPD-1013-01, 1981. 14. Robinson, C.R., "Ship Wave Color Graphics Using DI-3000." Computer Science Senior Project, U. of North Carolina at Charlotte, 1986. APPENDIX Denvatives of Velocity Components in Z-direction ¢,~z = || (t )~( ~ sec(( · n~d~d (yz = IJ (ll Y)~( ~ sec(( · n~d~d (21) (22) (yZ=|| (~;s ) sec((·n~d~dq-||~sec((·n~d~dll where cs((,ll) is given Eq.~15), r= (~E,-x)2+(ll-y)2+~- zy2~ 1/2, and ~ shown in Fig. 2. We assume that surface panel is sufficiently fine that the term sec(( · n) becomes zero. After algebraic manipu- lation, the following expressions are obtained: 0~z = ooK,~,~(l,l) + o~,~xK~,,~(l,l) + K<r,~2,11) + a,,fyK,~,~(l,l) + K,~,~1,2~) (24) ¢yz = c~oK<',,~l,l) + o`(xK~,,~l,l) + K<,y(2~1~) + c','tyK`~,(1, 1) + K<5),(1, 2~) (25) zz where = csoK=(l,l) + c,`(xK,~(l,l) + K,~(2,1~) + (s,,(yK=(l,l)+KO~(1,2~) (26) P4 K`,x(m,n)=-hH(m+l,n,5) +a[ H(m+3,n,5)- 5 h2H(m+3,n,7) +2xH(m+2,n,5~-lOxh3H(m+2,n,7) ] +b[ H(m+ 1 ,N+2,5~-Sh2H(m+ 1 ,n+2,7) +2yH(m+ 1 ,n+ 1,5~- lOyh2H(m+ 1 ,n+ 1,7) +c[ H(m+l,n,S)-Sh2H(m+l,n,7) ~ (27) K<7y(m,n)=-hH(m,n+ 1 ,5) +a[ H(m+2,n+1,5~-5h2H(m+2,n+1,7) +2xH(m+ 1 ,n+ 1,5~- lOxh2H(m+ 1 ,n+ 1,7) +b[ H(m,n+3,5)-Sh2H(m,n+3,7) +2yH(m,n+2,5)- 1 Oyh2H(m,n+2,7) +cE H(m,n+1,51-5h2H(m,n+1,7) ~ (28) K`7z(m,n)= h2H(m,n,S) -a[ 2hH(m+2,n,5~-Sh3H(m+2,n,7) +4hxH(m+l,n,5~-lOxh3H(m+l,n,7) -b[ 2hH(m,n+2,5~-Sh3H(m,n+2,7) +4hyH(m,n+1,5~-lOyh3H(m,n+1,7) -c[ 2hH(m,n,5~-Sh3H(m,n,7) ~ (29) 446 c=ax2+by2-zO h = z - Zo z0= axO + byO 11 (t X)m-l(~1 y~n~1 p = ~ (~_X)2 + (ll y'2 +h ~1/2 The H(m,n,k) integral has recursion relations and actual computation can be done very efficiently and accurately. --~ = =- ~' n ~ ~9 F4 1. Coordinate system. gS,y,') rf P3 ;~ P1 ~ ~s Flg. 2. Field poinVpanel geometry. Figure 3. Typical Panel Arrangement

o-s o~ c~ o.o o~ 0.1 n.e t -0.2- ~0.8 -0.6 -0.4 -0.2 0.0 0~ 0.4 0.6 0.8 z/(L/2) Fig. 4 . Wave Prof iles at Hull Surface Wigley Hull at Fn=O. 316 -1.0 - .8 -0 6 -0.4 -o~ o.o o~ z/(L/2) Fig. 5. Wave Profile Comparisons Wigley Hull at Fn=O. 32 co o X3 _ ~ Linear Theory O Nonlinear Theory _ 1=~ T .~ \ Lo \ o \ ~ _ _ ~6 - at Hull Surf ace a5 c~ 2.5' 0.5 Fig. 6. Wave Making Resistance Comparisons Series60, CB=0 . 60 . Linear · Theory / Nonlinear / · Theo ~/ ~ Experiment 77 ~ / ~ 1 ~ / . ~ _ ~ 1 L\ '/ C - Rw _ w 2pU2S n~ O.l -nn -n l -0.2 _n ~ O1 -0.1 _ . . , . . . . . . · 0.2 0.25 0.3 0.35 ,,~ Fn= U/(gL)1/2 -0.3 447 . ~ -0.4 -0.2 0.0 0.2 0.4 x/(L/2) Fn = 0.22 1 1 1 1 ~r I I T ~: . = 1 1 1 | 7` Linear Theory L I O Nonlinear lheory I I Experiment 1rith Nodel Fixed I L I -0.2 0.0 0.2 x/(L/2) Fn = 0.28 I ~linear lheory L I O Nonlinear 1heory 1 I Experiment 1rith Yodel Fixed - 1.0 - .8 -0.6 -0.4 -0.2 0.0 0.2 x/(L/2) Fn = 0.30 1.0 0.( O. 1.0 F7FI ~1 ~'o'll 1 1 1 11 Ll 11 0.6 0.8 1.0

0.3 0.2 0.1 Cal -0.0 -0.1 -0.2 -0.3 e _ - 1.0 - .8 -0.6 -0.4 -0.2 0.0 0. x/(L/2) Fn = 0.32 | ~Linear Theory _ O Nonlinear Theory E3rper~ment Tenth Yodel Fixed ~ V ~ 0.6 0.8 1.0 - 1.0 -0.8 Figure 7. Wave Profile Comparisons at Hull Furface Series 60, Cg = 0.60 Aim,\ ( \~ \ \ five U\ight _ _ ~linear Theory O Nonlinear Theory _ _ I Experiment with Nodel Flax _ /. -0.6 - .4 -0~ 0.0 0.2 0.4 0.6 x/(L/2) Fn = 0~35 my', \~ . _ . '> 1 1 1 1 1 1 T 7 We___ 7 1 Aft//; //~ /Z~/L~Jsource Strength / LINEAR SOLUTION `~W \N \\ \ C\ ~\ Wa~eight In, FIR ~ 111 _~ 29 . __ _ _-_ 0.8 1.0 ~ 1 3 . 0 art /; Yew' ale / source Strength '/~ NONLINEAR SOLUTION Figure 8. Contour Plots of Wave Height and Source Strength Series60, CB=0.60 at Fn=0.32 448 . ~

DISCUSSION (SESSION V, PAPER 5) John V. Wehausen University of California at Berkeley, USA The agreement between the authors' ~nonlinear. computations and measurement is impressive, but I wonder is one should not be suspicious of such good agreement. Since the computations are based on the irrotational flow of an inviscid fluid, shouldn't a measured wave resistance based upon residuary resistance, or even wave-pattern resistance, show at least some discrepancy with the computed resistance? The wave resistance is a rather sensitive test, for it is the difference between two large numbers and one expects to see some effect of viscosity in the stern region. I shall be curious to see the results for similar computations for Series 60, CB = 0.80. I also have a question concerning the numerical procedure. ~Convergence,. or lack of it, appears to depend upon the choice of Qua, but Q,~,, = 3 or 5 both lead to convergence, but to different values. Which should I believe? and why? Even if they were identical, can one give an estimate of the difference between the obtained value and the solution to the exact problem (when 1 _ 0 in some fashion)? DISCUSSION Ronald W. Young University of California at Berkeley, USA Your results show remarkably good agreement in terms of wave resistance and wave/profile when compared with experimental data. Jensen and Soding (1988, Jahrbuch der Schiffbautechnische Gesellschaft) have also obtained similarly good agreement for the same formulation but using discrete sources. They also have very good conveyance in their iterative scheme. I was wondering if it will be possible for you to compare your results with theirs. To my recollection, their stern wave profiles do not have as good an agreement with experiments as yours. Of course, for an identical model, and identical mathematical problem, there should be only one set of correct results. For a blunt bow ship, Jensen & Soding had difficulty obtaining convergent and accurate results in the bow area. You may like to try a similar test with your code. I want to congratulate you and your colleagues, including Dr. S. H. Kim, for the successful development of this code. DISCUSSION Hauime Marno University of California, Santa Barbara, USA The examples, which are employed in this paper, are Wigley hull and Series 60 hull. In these cases the effect of nonlinearity is comparatively small, and the advantage of the nonlinear computation may not be well demonstrated. The nonlinear effect is remarkable in the case of hull forms as was employed in the paper by Maruo & Oginara (1985), for which the nonlinear computation is indispensable. DISCUSSION Hoyle Raven Maritime Research Institute Netherlands, The Netherlands Your method shows a very good convergence for the cases studied. However, this convergence depends on the numerical (artificial) damping caused by the difference scheme. This damping does not vanish upon convergence but remains present in the result. For the rather mild cases shown, its influence does not seem to be large, but for full hull forms the required value of AN could be much higher and the damping could be more detrimental. Similarly, a reduction of the panel size is likely to ask for higher Q i. Therefore, I believe that other inherently more robust formulations of the iterative procedure remain something to be searched for. Could you comment on this? DISCUSSION Kuniharu Nakatake Kyushu University, Japan I congratulate you on your fine results. You tested a new family of upstream finite difference schemes in order to satisfy the radiation condition. As well known, the calculated results vary with choice of different scheme. I think it is not reasonable. We are using a phase- shift method which we call KU (Kyushu Univ) method.' The feature of this method is nonuse of difference scheme. According to our experience, it is applicable to very narrow computation domain. And a theoretical proof for it is to be presented by Dr. H. Seto (Mitsubishi Heavy Ind.) in Nov. 1990. I recommend you to test KU method too. [1] J. Ando and K. Nakatake, Tran. of West-Japan Society of Naval Architects, Vol. 75, 1987. DISCUSSION Dimitris Nakos Massachusetts Institute of Technology, USA The design of iteration schemes, based on linearization about the previous solution step, is the most critical issue in numerical solutions of nonlinear steady ship waves. Such schemes are desired to be convergent within a few iterations, due to the large computational effort involved at each solution step. It is fair to state that the part of the wave flow which is expected to resist convergence is associated mostly with the diverging wave system. The authors have demonstrated effectively that excessive numerical damping is able to Filter out. all short diverging waves and consequently accelerate the convergence (see e.g. Figure 8). The question that arises is the following: Is it proper to disregard essential features of the flow in the name of numerical convergence, in particular, in light of the fact that the relatively long transverse waves are mostly unaffected by nonlinearities? AUTHORS' REPLY We would like to express our thanks to all of the participants who have shown their interest in our paper, made comments, and expressed their critiques. Several questions have a common interest and we have grouped them without advance consent. Reply to Professor Wehausen: The experimental values of the wave-making resistance we used to make comparisons were obtained by the longitudinal wave profile measurement suggested by Eggers[l]. In his derivation the effect of the viscosity is neglected and the corresponding computer code has been reported by Reed[l]. Table A summarizes the linear, nonlinear and measured wave-making resistance for Series 60, CB = 0.60. At the lower Froude number, the discrepancy between the measurement and the nonlinear solution is significant and even at the higher Froude number the difference is still noticeable. 449

Reply to Professors Wehausen, Marno, Yeung, and Dr. Raven: They have shown common interest in the nonlinear effects on a ship with blunt bow and large block coefficient. We have selected Series 60, CB = 0.80 which has a half entrance angle of 43 degrees as an additional test case. Figures A and B show the linear and nonlinear wave profiles together with a series of wave profiles obtained during each iteration for Fn = 0.20 and 0.25 with Q,l,'., = 3 and 5 respectively. No comparisons were made with measurements because of lack of availability. As Prof. Marno pointed out, the nonlinear effects on Series 60, CB = 0.80 is more significant than those on Wigley hull or Series 60, CB = 0.60. One noticeable difference observed between ships during the iteration is that for Series 60, CB = 0. 80 the wave profile after the first iteration has the largest hump and hollow near bow region and from the second iteration the wave profiles oscillate within the linear and the first wave profile band until the solution converged. On the other hand, for Wigley hull or Series 60, CB =0.60, the wave profiles always progressively grow as the iteration number increases until the solution converged. For Fn = 0.30, we increased Qua from 5 to 7 to 9, but failed to get a converged solution. As the iterations continued, the wave amplitude grew continuously and after the third iteration the wave amplitude near the bow region exceeded the ship draft and the solution diverged. Without having measurements, it is difficult to evaluate the quality of our converged solutions for a fullform ship. Whether the Fn = 0.30 case is beyond the capability of our code would require further study, but as the A,,,, method is very new, and has now been shown to converge nicely for lower Froude numbers even for this fullform ship, with further research we would anticipate convergence even at Fn = 0.30. Reply to Prof. Young: Dr. Jensen compares his wavemaking resistance computations with Ogiwara's experimental results which were conducted at free sinkage and trim conditions. Our code can be extended to handle the free sinkage and trim conditions. There are no plots on the wave profile in Jensen's paper. Reply to Prof. Nakatake: The phase shift and the shift of collocation point may be different terminology for the same approach. A shifting also introduces numerical damping and is algorithmically similar to the use of upstream finite difference formulae. We had earlier implemented the shift of the collocation method in our computer code and tested it. In our numerical experiments, we shifted the collocation points forward between 10% and 50% of the characteristic length of each panel. Our numerical results showed that the shifting apparently removed the oscillation of wave profiles, in particular near the bow region, but decreased the wave amplitude proportionally as we increased the shifting. According to our numerical experience, we haven't found significant advantages to applying the shifting logic to our code. Reply to Prof. Wehausen, Dr. Nakos, and Dr. Raven: It is well known that the panel method using a Rankine type singularity with an upstream finite difference scheme to satisfy the free surface and radiation conditions provides good results near the ship [4]. But because of an inherent local damping introduced by unwinding, the wave amplitudes attenuate significantly as they propagate and the results get poorer far downstream. Our aim was to develop a quality computer program which is practical for ship designers in the shipbuilding business. The linear codes of the past made many simplifications. We have provided an approach to solve the original nonlinear problem that eliminates these, and gives results that compare well with those of experiment, with a modest use of computer time. We haven't examined the relations between panel size, Proude numbers, and A s carefully at this point, however we will be considering these questions in the future. Your comments and critiques are very much appreciated. References: [1] Eggers, K.W.H., cuber die Ermittlung des Wellenwilderstandes eines Schiffsmodells durch Analyse seines Wellensystem,. I.II Schiffstechnik 9, pp. 79 -84, 1962; Disk 85; 10, pp. 93-106, 1963. [2] Reed, A.M., "Documentation for a Series of Computer Programs for Analyzing Longitudinal Wave Cuts and Designing Bow Bulbs,. DTNSRDC/SPD-0820-01, June, 1979. [3] Sclavanous, P.D. and D.E. Nakos, Stability Analysis of Panel Methods for Free-Surface Flows with Forward Speed, Proc. 17th Symp. on Naval Hydrodynamics, The Hague, The Netherlands, 1988. [4] Lindenmuth, W.T., T.J. Ratcliffe, and A.M. Reed, Comparative Accuracy of Numerical Kelvin Wake Code Prediction-'Wake-Off'," DTRC/SHD-12601-01, May, 1988. 450

Table A Wave-Making Resistance Comparisons Series60, Cg=0.60(Model Fixed) En i Linear I Nonlinear 0.22 0.369 0.255 . 0.25 0.664 0.354 0.28 1.581 1.045 0.30 2.144 1.460 0.32 1.815 ~1.357 0.35 2.148 1.836 (CWX103) Experiment 0.176 0.230 1.011 1.375 1.316 1.780 Fig.A. Series60, CB=O.BO at Fn=0.200 -0.5 -1.0 __ _ Symbol SWIFT ~ Linear Have _ I E} - Iteration No.1 O Iteration No.2 O Iteration No.3 I ~ Iteration No.4 _ - - ~ - - Nonlinear l ~ ~ = 1--1 . _N = 1 1 ~ _= === =. . Q ~ ===~N ~ =~- _~ fit _ ==~= _ =-__ _ _ _ , ~ _ - 1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 x/(L/2) 451

Fig.B. Series60, CB=O.8O at Fn=0.250 .o 0.5 - cot ~ o.o _' ._ o at -0.5 -1.0 _ I S~nbolSVIFr l ~Linear Have _ l ~Iteration No. 1 O Iteration No.2 O Iteration No.3 l ~Iteration No.4 _ l | - - ~ - - Nonlinear . . 14~: ~ W~+ _0 _ _ -=== == I I I AL ===~- -r ~- === - - 1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 x/L/2) 452 0.6 0.8 1.0