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.

Numerical Evaluation of the Complete Wave-Resistance Green's Function Using Bessho's Approach H. T. Wang Naval Research Laboratory Washington, USA J. C. W. Rogers Polytechnic Institute of New York New York, USA Abstract This paper presents a numerical method for calculating the total Green's function for the wave resistance case of a source in steady translation below the free surface. Starting with a representation of this function in the complex plane given by Bessho a series of transformations of variables are used to reduce it to three real single integrals. The integrands are regular and are entirely in testes of elemen- tary functions. Two of the integrals are even in the direc- tion of flow while the remaining integral is odd. The even and odd integrals may also be conveniently expressed in terms of the near and far field components. While the method is applicable nearly everywhere (except for the well known difficulties at the limiting cases when the source is located near the free surface or the field point is close to the source), computer time varies in the computation domain. A computer code has been written to implement the above method. Sample calculated results are given in several forms to show the accuracy and computer time requirements of the code. A series of line and contour plots are given to show typical shapes of the integrands at dif- ferent locations in the computation domain as well as to exhibit the relative behavior of the various component integrals. 1. Introduction For a number of decades, interest in the important Green's function for a submerged nonoscillating source mov- ing at constant forward speed in fluid of infinite depth has largely centered on its far field characteristics. This was both due to the greater simplicity of the far field evaluation as well as its direct applicability to finding the hull resis- tance component due to wavemaking. The far field wave pattern required only the evaluation of a single integral with a regular integrand over a one-dimensional wavenumber space while the complete Green's function involves the addi- tional calculation of a double integral with a singular integrand over a two-dimensional wavenumber space. The applicability of using only the far field analysis to obtain the wave resistance was aided by the pioneering work of Michell [1] who showed that reasonable estimates of the source strengths modeling the hull surface for the case of thin ships could be obtained directly from ship geometry, without need to ascertain the near field mutual influence of the sources. The landmark paper by Eggers, Sharma, and Ward [2] presents a comprehensive survey of different methods of using the single-integral far field Green's func- tion to obtain the wave resistance. In recent years, with the availability of ever larger and faster computers, interest has been enlarged to include the near field terms of the Green's function as well as a more accurate calculation of its far field behavior. An evaluation of the near field terms would give a more accurate determi- nation of the source strengths for hull forms which do not conform to thin ship theory as well as a more detailed defin- ition of the flow field and pressure forces on or near the hull. In the far field case, the use of modern day remote sensing technology makes it of interest to assess the ship wake for wavelengths which are significantly shorter than those applicable to the wave resistance problem. Efforts at rendering the initial double integral represen- tation of the Green's function amenable to numerical calcula- tion usually involve expressing it as a series of single integrals. Noblesse [3] gives several alternate single integral representations. A popular form is to express the Green's function as two single integrals consisting of a near field part N which is even in the flow direction x, and a wave distur- bance part W which is defined only downstream of the source. The near field part N has an integrand which is in terms of the higher order derived exponential integral func- tion. Noblesse [4] and Euvrard [5] have conducted detailed studies of the behavior of N and W. especially at limiting regions of the computation domain. In terms of actual numerical implementation, Newman [6] has developed a procedure for calculating N in terms of extensive sets of tabulated coefficients of Chebyshev or ordi- nary polynomials. The coefficients take on four different sets of values, depending on the radial distance R from the source. In the case of the wave disturbance part W. Nob- lesse [7] and Newman [8] have developed procedures for the specialized case of the vertical xz centerplane while Baar and Price [9] implement a more general calculation for the entire domain with the exception of the region near the free sur- face. In these studies, the computation region is again divided into several regions depending on the values of x/z or x/v y~ + z ~ . In many cases, the solution is in terms of 133

Image Sink _ (X,Y, - Z) Free Surfaced 0 Free Surface (x,y,z) ~ / Field Pod Y ~ ~ x - .(X,Y,Z) Current U Source Point Figure 1. Definition of Coordinate System and Flow Configuration Go=km| dd| dk or ~o ~ O exp [k(z +Z) + ik (xX) cos ~ + ik(yY) sin 8] (2) k cos2 ~1incise In the above equation, all length variables have been made dimensionless by multiplication by the factor g /U2 where g is the gravity constant. Go represents an inverse double Fourier transform over the two-dimensional k wavenumber space. The parameter ~( > 0) is added to ensure that the radiation condition is satisfied at infinity. Also, following Ursell, the terms (xX), (yY), and (z + Z) appearing in Eq. (2) are simply replaced by x,y, and z for the sake of convenience. 2.2 Concise Statement of Bessho's Approach The usual way of simplifying the double integral expression for Go is to perform an initial integration over k or 8, thus reducing the double integral to a single integral. The integrand, however, is not entirely in terms of elemen- tary functions but contains the higher order derived exponen- tial integral function E~(u) defined by ~ e-& E~(u)=i d'A (3) u ~ 2. a series of functions, and there is often a delicate balance between the regions of convergence of near field and asymp- totic expansions. Our paper presents a numerical implementation of a representation of the Green's function developed by Bessho [10]. By means of a series of ingenious transformations, he succeeds in reducing the entire Green's function to a single integral along a curved path in the complex plane. The integrand is regular and is entirely in terms of elementary functions. This remarkable representation has been noted by previous investigators [4,6] and has prompted Ursell [11] to give a more complete derivation, including the justification of an important intermediate step. Our work differs from previous numerical implementations in several aspects. Our work is entirely in terms of integrals as opposed to the use of expansions in series. Our implementation is for nearly the entire computation domain and for the complete Green's function as opposed to specialized domains or particular parts of the Green's function. The paper starts with a concise statement of the critical points of Bessho's contribution. Then, a more complete description of the derivation, as given by Ursell, is given in outline form. One reason for giving this derivation is to point out the starting point for our numerical work, which occurs before the final single integral representation is reached. A detailed description is then given of the transfor- mations required to convert the original complex integrals to three real integrals, two of which are even in x (Ge) and one is odd in x (Go). The simple relationship between Ge and Go and the commonly used N and W is pointed out. Pecu- liar features of each of the three integrands are discussed, and the procedures used for their integration are pointed out. A special limiting process is used to obtain Ge and Go on the axis directly downstream of a source at the free surface. The variation of computer time requirements and accuracy . . . Of the calculated results in the computation domain are dis- cussed. Numerical results are presented to check on the accuracy of our procedure and to illustrate in graphical form the behavior of the integrands and their integrated values. The paper concludes with a summary of the principal find- ings. 2. Derivation of Bessho's Single Integral Representation 1 Initial Representation in Double Integral Form In this work, we will consistently follow the notation used by Ursell [11]. Figure 1 shows the coordinate system used in our work, where x corresponds to the direction of the current flow U. y is the horizontal direction perpendicu- lar to x, and z is the vertical direction positive downwards. For a stationary nonoscillating source located at (X, Y,Z) the Green's function G(x,y,z; X, Y,Z) at field point (x,y,z) is given by G(x,y,z) = Rid R + Go (1) where Ri = [(xX)2 + (y_ y)2 + (z_Z)2]~/2 R = [(xX)2 + (y _y)2 + (z+z)2]~/2 In Bessho's approach, the initial integration to reduce the double integral to a single integral is not performed at the outset. Instead, after a series of transformations of the double integral, including changes of variables, deformation of the paths of integration in the complex plane, expressing Go as derivatives of more convergent integrals, and inter- changing order of integration, the double integral is finally converted to single integral form by means of the following crucial equality I ~ exp (imP) dm (4) _= m Q Pi exp (i P Q) (P > 0, Im Q > 0) Sari exp (i PQ) (P < 0, Im Q < 0) 0 (P Im Q < 0) 134

where P is real and Q is complex. Thus, the order of integration is reduced entirely in terms of elementary func- tions. The following section gives a summary of Ursell's detailed derivation and justification of Bessho's approach. 2.3 Summary of Ursell's Derivation Using the equality 1 cos ~ d ~ = 1 0 cos ~ d rewrite Go in Eq. (2) as the following sum of A and B Go (x,y,z) = A(x,y,z) + B(x,y,z) (6) where A and B are now integrated from ~r/2 to +'r/2 over 8. By introducing the following two sets of new variables m = k cos 8, tanh vsin ~ (7) y = p sin cY, z = p cos cY, x = p sinh,B, R = p cash ~ (8) and writing A and B as A = Ao + i at Al, B = Bo + i aX Be (9) the following expressions are obtained for Ao, Bo, A i, and Be, Ao =| dm | dv exp[m p cash (vice) +imx] (lea) Bo =| dm | dv exp[m p cash (vicY)imx] (lob) A 1 = 1 km I dm I ~ do exp[m p cash (vion) +imx] 7r `-o O -= mcoshvis (lOc) Be= slims dml" do exp[mpcosh(vin)imp] 7r To 0 - ~ mcash v +ie (led) The reason for introducing the new integrals Al and Be is that they are more strongly convergent than the original integrals A and B. The double integrals Ao and Bo are simplified by Ursell by making the change of variables vicev, and using equality (9.6.24) of Ref. [12] to convert the integrand in terms of the Bessel function Ko Ao + Bo =lo dm KO(mp) cos me 2 2 = , = Vp2 + x2 R Actually, the integral Ao + Bo can be more directly simplified (without the need to use Bessel functions) by not- ing that the order of integration can be interchanged since it is absolutely convergent, carrying out the simple integral with respect to m, and making the change of variables w = ev in the resulting single integral to arrive at the same result given above. By making similar shifts in the paths of integration of Al (v = for + 2 Hi + w, w real) and Be (v = ion 2 Pi + w, w real), the sum A ~ + B ~ takes the form Al + B.` = km | dm | dw 7r ~o ~ co exp (imp sinh w + imp sinh A) + (single integrals m~ s~nh(w + fix)ze due to residues at poles crossed by the shifts) (12) By interchanging the order of integration, Bessho notes that the integral with respect to m is of the form given in Eq. (4) and thus the double integral is converted to a single integral entirely in terms of elementary functions. Ursell points out that the double integral in Eq. (12) does not satisfy the sufficient (but not necessary) condition of absolute convergence in order to justify the interchange of order of integration. His proof of the validity of the inter- change basically consists of extending Bessho's idea of gen- erating derivative functions, shown in Eq. (9), to generate the following still more strongly convergent integrals A2 and B2 A2 = 1 lime dm5 ~ dv exp[mpcosh(yict) +imx] To 0 -= (1i)(mcoshvie) (13a) B2= 11iml dam do exp[mpcosh(vin)imp] 7r To O -= (l+im)(mcoshv +ie) (13b) which are related to A ~ and B ~ as follows Al = (1 at) A2, BY = (1 at) B2 (14) Ursell proceeds to operate on A2 and B2 in a manner similar to Bessho's transformations of Al and Be. That is, once again the paths of integration are shifted for A2 (v = for + 2 Hi + w, w real) and B2(v = ice 2 pi + w, w real), and the resultant residue evaluated at the pole of the integrand crossed by this shift (for some but not all m) where the pole V is defined by 135

m + it for B2 In either case, Eq. (17a) is expressed as the following sum cash V(m) = m it forA2 (15) of two single integrals By making the further change of variables from m to V, as given in Eq. (15), in the single integral resulting from the residue evaluation, the following expressions are obtained for A2 and B2, where each is given in terms of a single integral and a double integral A2 = km| dm I°° dw 7r~-o 0 1 im -= exp (imp sinh w + imp sinh I) 16 mi sinh(w + ICY)it ( a) 2i 1 1= dV exp(pcoshVcosh(VicY)+ixcoshV) 2 ~` 1i cash V BY = ~ km I dm I dw ~ -so 0 1 + im -= exp (imp sinh wimp sinh A) (16c) m + i sinh (w + in) + it 2i 11(~ Or ) dV P( P cash Vcosh (V-ZCY)ix cash V) (16d) where cY is taken to be0 in Eqs. (16b) and (led). Ursell shows that the double integrals (16a) and (16c) are absolutely convergent, and thus it is permissible to inter- change the order of integration. By writing m for m in Eq. (16c) the two double integrals can be combined into a single form, with order of integration interchanged 1 1= dWI" dm exp(impsinhw+impsinh~g 7r - ~ - ~ (1im) [mi sinh (w + ice)] = ~ I dw l(w,O (17a) The presence of the extra factor (1im) in the denomina- tor makes the equalities given by Eq. (4) not directly appli- cable for evaluating I(w,,8). Instead, Ursell evaluates it by means of an elaborate contour integration, accounting for the residues due to the poles which are enclosed for O<w<,B and ,B<w< as. Actually, Eq. (17a) can be put in a form to which Eq. (4) is applicable, by rewriting the denominator, as follows (1im)[mi sinh(w + de)] sinh(w + i(Y) + 1 ~ mi sit | (w + icy) - m +l 136 1 1 I(w, ,8)dw 7r 00 2i 1 it exp [p (singsinhw) sinh (w + )] dw (18a) + 2i | P~[P+( Dh§( + ~ )] dw (lab) The component single integrals (16b), (16d), (18a), and (lab) which together define A2 + B2 form the starting point for our numerical analysis and implementation. Ursell proceeds to essentially sum these component integrals and perform the derivatives indicated in Eqs. (14) and (9) to 16b obtain Bessho's representation of Go = A + B as a single ( ) integral in the complex plane. To begin the compacting process, these four integrals are all expressed in terms of the variable u by letting V = u + 2 7ri in Eq. (16b), V = u2 Hi in Eq. (16d), and w = uin in Eqs. (18a) and (lab), resulting in A2+B2=2i t1-~+l~i+ I-~-l~i] 1 + h exp [(p sinh (uin)x) sinh u] + 2i | em+ 1 + sinh u exp [ - (p sinh (uin)x)] (19) In carrying out the differentiations a/8x indicated in Eqs. (14) and (9) it is important to recognize that the limits of integration involve the variables cr and ,6. From the transformation Eq. (8), it is seen that ax = 0 1 = p cash ~ ant _ aa~ = ~ (20) Performing first the operations indicated in Eq. (14) on A2 + B2 to obtain A l + B l and then the transformations indicated by Eq. (9), using the result given in Eq. (11), the following complex single integral representation for Go is obtained Go=A+B= R +i a (Al +51)=2 tI oo- Phi + I -00+ 1 pi ~ sinhu exp[(p sinh (uice)x) sinhu]du (21) The above integration paths are not the same as those used by Bessho, shown in Fig. 2. In order to obtain the integrals in the Bessho form, make the substitutions x = X, X = p sinh b, u = v, and note that Go(cY) = Go(ax), resulting in

:0' 2 i] ['25 (is, in) O \ L0~- 2 iJ -1 ~ Figure 2. Bessho's Integration Paths at- 1 Hi co+ 1 Hi Go (-X,y,z) = -2 50 + Ib+i<x (no' - ~ i: sinhvexp [(p sinh(vin)X) sinh v]dv (22) 3. Computational Implementation of Bessho's Method In this chapter we will present our computational implementation of Bessho's method. As mentioned previ- ously, our analysis starts from the four component complex single integrals given in Eqs. (16b), (led), (18a), and (lab). Through a series of transformations, we simplify and con- dense these complex integrals to three real integrals, two of which are even in x and the other is odd in x. Taken together, these integrals represent the entire Green's func- tion. Reference [13] describes the FORTRAN 77 computer program which evaluates these integrals. 3.1 Simplification of Complex Integrals We start by making the change of variables ~ = iV for part of the integration of Eqs. (16b) and (16d) and then rewrite each integral as two integrals with real limits of integration, resulting in the following expression for A2 + B2 A2 + B2 = Eq. (18a) + Eq. (lab) 2 1 2 ~ exp[p cos; cos({a) + ix cost] d ~ (16b.1) + 2i5 ° expEpcoshVcosh(Vit) + ixcoshV] dV `~6b.2' + 4 | coshVsin~x coshV) or +25 2 ~ exp[pcos~cos(~+a) ixcos(] d; (16d.1 2i 5 ° exp[pcoshVcosh(Vin) ixcoshV] dV (16d.2 In the above, use has been made of the identity cash (iV) = cos V relating trigonometric and hyperbolic cosines. Carrying out the differentiations indicated in Eq. (14), the following integrals are obtained for A 1 + B 1 Al +s1 =(1 a )(A2 +s2) = 2i 5 exp[p(sinh,Bsinhw)sinh(w + in)] dw (23a) 21 2 exp[pcos~cos(;or)+ixcos(]d; (23b) + 2i 5 exp[pcoshVcosh(V - in) + ixcoshV] dV (23c) or + 2 5 2 exp[pcos~cos({ + a) - ixcos(] do (23d) o + 2i 5 exp[pcoshVcosh(Vin)ixcoshV] dV (23e) The above equation is a result of the following considera- tions. Recalling from Eq. (20) that 8,6/8x = 1/R, the derivatives 8/ax of the variable limits of integration appear- ing in Eqs. (18a) and (lab) cancel each other. Recalling from Eq. (8) that x = p sinh ,B the operation (18/8x) applied to the integrands results in Eq. (lab) going to zero and a removal of the denominator in the other five integrals. Proceeding similarly, using Eqs. (9) and (11), we obtain Go in the form Go = R + i ax (A 1 + B1) = G1 + G2 + G3 + G4 = 2 5 sinh (w + in) exp[p(sinh,Bsinhw)sinh(w + in)] dw (24a) +25 ~ cosfexp[pcostcos((~x)+ixcos(]d; (24b) _ _ct + 2 5 2 cosf exp[pcos~cos({ + cY)ixcos{] d{(24c) exp[p coshV cosh(Vice)] dV (24d) In the above equation, the two integrals in terms of the vari- able V have been combined into a single integral by using the well known identity e:XCoshv _ e-`XC°shv = 2i sin (x cash P)- 137

In the following, we express the above complex integrals as the sum of real and imaginary parts, and con- sider only the real part. By means of a series of transforma- tions, we reduce the resulting integrals to compact form. 3.2 Conversion to Real Integrals Consider first G1 given by Eq. (24a). By first expand- ing sinh (w + ice) and then making the change of variables cosh w dw = dw', Go becomes the following G ~ = 2 5 [cos ~ ~ cos(x sin ~ ~ yw A) + sin ~ sin (x sin ~ ~ yw ~w2 + 1)] e -xcosa w + zw2 dw (25) Let us now consider G2 and G3 given by Eqs. (24b) and (24c). By using the well known relationship e+'XC°st = cos (x cos 0+ i sin (x cos 0, the real part of these two integrals can be conveniently written out. By making the change of variables ~ = if' in G3, G2 + G3 can be written as the following single integral _x G2 + G3 = 2 1 2~ +a cos ~ cos (x cos Me -PCostcos(~-a)d; (26) By making the change of variables ~ = if' + cx/2, using various trigonometric identities for the cosine terms, keeping only the even functions of ~ (since the odd functions give a zero integral), and (very similar to the final transformation indicated for Eq. (25)) making the change of variables cos (d; = dw, the integrand of G2 + G3 takes on a form resembling that of Go given in Eq. (25) a G2 + G3 = 41 2 [COS (X COS (X COS a! it) . CX sm cos(xsin~2 w)w :2 sin(xcos 2A) 0` P (cosa+1 - 2w2) sin(x sin 2 w)]e 2 dw (27) Finally, let us consider G4 given by Eq. (24d). By making use of the identity cosh (Vice) = cosh V cos ax - i sinh V sin cz the real part of G4 is written as 00 G4 = 4 1 cosh V sin (x cosh V) cos (p cosh V sinh V sin cx) e -P cosh2 vcos a dV (28) By making the change of variables cosh V dV = dw and recalling Eq. (8), G4 becomes G4= 4lo sin(x~)cos(yw~)e-z(w +~)dw (2 By making the further change of variables dw = sec28 d8, we can transform G4 to the single integral form given in Eq. (13.3b) of Wehausen and Laitone [14] G4 = 4 1 sin(xsec28cos8)cos(ysec28 sin8)e -z see dsec249db (30) In our numerical computations we find it somewhat more convenient to make the change of variables sec28 sin ~ = u and use the following alternate form proposed in [2] G4 =410 sin(xs)cos(yu)e-s Z Z~ du (31) where s(u) = [ 2 ~ = see 8. Adding the component integrals given by Eqs. (25), (27), and (31), the resulting form for Go in real form is given by Go = GX + Go + Go = Ge + Go = 2 | [cos cx ~~ cos(x sin ax yw A) + sin cx sin (x sin cx ~ yw ~)]e -xcos`Yw+pcos<YW2dw (32a) cos + 41 2 [cos 2cos(xcos 2 ~)cos(x sin '2Y w) . ax wsm ,~ sin (x cos 2A) sin (x sin 2w)] ~ (cosa+1 - 2w2) e ~ dw (32b) + 4 1 sin (xs) cos (yu)e -s Z 25 - du (32c) where s(u) is defined in Eq. (31), and p and ax in Eq. (8). Since Go given by Eq. (32c) corresponds to the single integral part of the total Green's function given in [14], Ge, the sum of GX and Go given respectively by Eqs. (32a) and (32b), corresponds to the double integral part. Equation (32c) shows that Go is odd about x = 0 and has an infinite upper limit of integration while Eqs. (32a) and (32b) show that the Ge is even about x = 0 and has upper limits of integration which depend on the field point. 138

The integrals Go and Ge are also simply related to the commonly used near field and wave disturbance components N and W [4,6,9]. Comparison of the expressions for Go and W shows that W(x, y, z) = Go(X, y, z) + Go(~X ~ Eye Z) (33) Since W + N = Go + Ge, N is related to Ge and Go as fol- lows N(x, y, z) = Ge(X' y, z)Go( AX I, y, Z) (34) 3.3 Behavior of the Integrands It is of interest to investigate the behavior of the three Integrands in different regions of the x, y, z computation domain in order to derive effective integration schemes. While the behavior of the integrand for the odd integral Go is relatively well known, we do not believe that the Integrands of the even integrals OX and GC have been previ- ously considered. Figures 3 to 5 respectively show the Integrands for Gx, Gc, and Go at various (x,y,z) locations. In these figures, the horizontal and vertical coordinates have been normalized so that the independent variable lies between 0 and 1, and the integrand lies between -1 and 1. To ~ _ Z . 6 ¢ llJ Z . 2 L o 111 - . 2 J <t > - . 6 \ 1.0 T\ _: \ - . 6 ant'/ -1.0 ~` ~ 1 , 1 1 1 1 1 0 . 0 . 2 . 4 . 6 . 8 I NDEPENDENT VAR I ABLE U Figure 5. Integrands for Go i\ Am\\ /\ x ~ 1, y ~ 1 . \ ~ x= 10, yin l. \\ ~ ^=20. ye 20. - l- - -\\ ~ ~\- \ jl \ \ / / ~ 1- 1 \1 Amp 1 ~ . Z-1 1 11 1 1 , 1 O ''J' ~ . O O 1~ \ Z .6 - \ ~ - \ z .2 '\ ~ _ ~ /, L D - . 2 _ J _ '3: > _ _ _ - . 6 _ _ -1 . O 1 1 1 1 1 1 1 1 Figure 3. Integrands for Ox / \ .8 1.0 ~ ~ or: 'A ~ \ 1 x - 1. y 1. Z 1 1 x ~20. ye ~ . 1. Z ~ . 1 x 1, A= 1, Z=20 \ / I \ J I .2 .4 .6 1 These figures show that the behavior of the Integrands varies greatly, depending on (x,y,z). Thus, it is not likely that any single integration approach will be optimum throughout the entire computation domain for a given integral. It is this great divergence in behavior which has led previous investigators to consider only specialized com- ponents of Go and to divide the computation domain into several regions [4-9]. Our principal effort has been to deter- mine minimum integration domains for the two integrals which have theoretically infinite upper limits of integration: Go and (as x /p m) Gx Figure 3 shows that Gx has a smooth integrand for small or moderate values of x /p, but at larger values, the integrand seems to concentrate at thin strips at the left and right ends of the domain of integration. This is due to the fact that the argument of the exponential function in Eq. (32a), f(w) = x cos cue w + p cos car w2 is zero at the ends of the integration interval, w = 0 and w = x/p, and takes on negative values in between. Noting that f(w) is symmetric about the midpoint of the interval, w = x/2p, the following integration scheme is adopted for Gx. First, the limits of integration lying on either side of x/2p are determined by sewing f(w) = A, where A may be chosen to be suitably large, say, 20 (i.e., e-A = 2 x 10-9), resulting in x ~ tow = +~ -/(x/2p)2 - 4A/(pcosa) (35) 2p 2p For cases where the term under the square root sign is nega- tive, Aw0. Integration for Gx is thus confined to the two strips Aw as follows Gx= 50 g(w)e f (w) dw+| P g(w)e f (w) dW = 5 P [g(w) + g( - w + x/ p)]ef (W) dw (36) .e i.0 where g(w) are the trigonometric functions multiplying ef (w) INDEPENDENT VARIABLE U in Eq. (32a) and the last equality results from the symmetry Figure 4. Integrands for Gc of f (w) about x /2p. 139

The error incurred by using a finite upper limit UM in the integral for Go, Eq. (32c), may be estimated as follows. By neglecting the trigonometric factors, which have magni- tudes less or equal to 1, and changing the variable of integration from u to s, defined in Eq. (31), the error ~ is bounded by ~ S 4 1s e-s Z s ds (37) where SM = 2 + I. By making the change of variables t = s2, the above may be integrated by parts to give z `~t + z Issue Z ds () ds (38) where d ( s ) = d (52 _ 1)-~/2 c 0. Therefore, it is ds u ds necessary to consider only the first term in deriving esti- mates of the maximum error incurred by using finite values of UM. By considering ~ and z to be given parameters, Eq. (38) may be solved iteratively for the required upper limit UM. Table 1 gives the calculated values of UM for ~ = (10-3, 10-4, 10-6) and z = (10., 1., 0.1, 0.01). Table 1 shows the manner in which UM increases with decreasing values of ~ and z. Table 1Variation of the Upper Limit of Integration UM in Go with Error Criterion and Source Submergence z ~ 1 Z=10 1 z=10 1 z-0~1 1 z=0.01 10-3 0.01 6.2 77. 880. 10-4 0.085 8.4 99. 1100. 10-6 0.60 13.0 140. 1500. 3.4 Limiting Behavior as pO Use of the above error estimate appears to indicate that UM ~ °° as z0. Yet it is well known [4,8] that Go(x, y = 0, z = 0) is given by Go(x' 0, 0) = 2 W(X) = -2~7rYI(x), x > 0 (39) where Ye is the Bessel function of the second kind of order one [12]. The existence of the limit shows that a higher order error analysis, accounting for the oscillatory tri- gonometric functions in Eq. (32c), would lead to a more efficient procedure for the integration of Go. In order to obtain the complete function Go for (x, 0, 0) it is necessary to calculate Gx and Gc defined in Eqs. (32a) and (32b) as y0 and z0. Since Go exhi- bits singular behavior in the neighborhood of this axis, it is necessary to use care in taking this limit. By starting with Eq. (29) for Go, and neglecting terms which are well behaved near the axis p = 0 and which are of the order of 140 1/ fix ~ as fix ~ no, we have derived in [13] the following expression giving the first order far field behavior in the neighborhood of the axis (x, 0, 0) GO(x,y,z)~sgn(x)2:e-xZ/4(Y +Z)sin( ~ Y 2 + ~ + ~)+2~(sinx+sgn(x)cosx)e-Z (40) This equation generalizes, for the case of large x, the expression given in [5] for the wave disturbance W for the specialized case z = 0. The above equation confirms the findings in [5] that the limit z = 0, y0 is singular, and that (provided z ~ O) the limit y = an, zO is indepen- dent of a. Since the near field component N is well behaved, Eq. (34) shows that Ge must show a similar behavior. We find it convenient to compute the values of Gx and Gc for y = 0, z0. In this case, ~ = tan-~(y/z) = 0. The integral Gc may be calculated in a straightforward manner by setting z = 0 and then evaluating the well behaved integrand between the finite limits O and 1. In the case of Gx, the limiting process must be applied with-some care. In particular, taking care to preserve the previously noted symmetry of the exponential factor about the midpoint of the integration interval x/2z, the integral Gx has the fol- lowing limit km Gx (x, 0, z) = km 2 5 w e -xw+zw2 dw zo zO 0 ~/~ = 1im25x/2z ~ w + x/zw ~ e-xwdw Z-° ° L~ ;(x/z -w)2 + 1 ~ 50 ~~e_Xwdw+2 (41) Thus, Go may be conveniently calculated on the axis (x, 0, 0) by adding Go, Go, and Gx. 3.5 Computer Time Requirements Computer time depends on various factors, of which the values of x, y, and z, the error criterion I, the type of integral, and the integration rule are the most important. The simple trapezoidal rule and the higher order Simpson's rule were used to integrate each of the three integrals over a range of x, y, z, and c. The integration starts with an initial partition of the given range into two intervals, and then suc- cessively doubling the number of intervals until the integrated values from two successive iterations agree to within the specified error e. We have found that on average, the use of Simpson's rule for Ge(=GX + Gc) results in a computer time which is 0.3 of the time using the trapezoidal rule, while the corresponding ratio for Go is 1.5 [13]. These trends may be due to the fact that the (usually) smoother integrands of Ge benefit from the parabolic fit through the points used in Simpson's rule, while the (often) oscillatory integrands- of

Go seem to be better approximated by the straight line fit of the trapezoidal rule. We realize that, based on the diverse behavior of the integrands shown in Figures 3 to 5, the use of different integration rules (for the same integral) based on x, y, and z would result in even greater savings of computer time. We have not, however, pursued such detailed refine- ments. Figure 6 shows the average CPU time per field point for Gx, Gc, and Go, respectively, using the Hewlett Packard (HP) 9000, Model 550 minicomputer for ~ = 10-4 for vari- ous intervals in the region 0.1 c x, y, z c 40. It is well known that for values of x, y, z near the origin (0, 0, 0) and the axis (x, 0, 0) the behavior of Go is singular. For a given x-interval, the figure shows two sets of numbers. The upper set gives the average CPU time for a fixed value of y (indi- cated on the vertical axis) and a uniform xz grid of field points over the interval indicated on the horizontal axis. The lower set of numbers gives corresponding results for a fixed value of z and a uniform xy grid. Several runs were repeated using c= 10-6 and these suggest that average CPU time is approximately tripled compared to the ~ = 10-4 results. Also, some runs have been made on the Cray XMP/24 mainframe computer and the corresponding CPU times are typically 60 times smaller. O.. 1.0 _ N >I 10 _ 40 z: (.02, .02; 1.0) (.01, .0s; - ) ( 01, .10; .10) (.02, .10; -) (01, .25; .02) (.01, .25; - ) x o 1 y (.02, .05; _) 1 0 (.05,.30; -) 10(.20,1.3; -)40 1 . (.65, .20; 2.3 _ (.20,1.0; -) (.35,.20;.20' (.20, .40; - ) (.20, .45; .02) _ (.10, .20; - ) (.10,.10; 1.5) (.05, .25; - ) (.10, .10; .15) (.02, .20; - ) (.02, .25; .02) _ (.02, .20; - ) (01, .25; .01) (.01, .30; .01) (.10, .45; .01) Figure 6. Approximate Computer Times per Field Point for ~ = 10-4 Most of the computer time trends can be directly inferred from the expressions contained in Eq. (32) or shown in Figs. 3-5. One trend is that the overall sum of the CPU times for all three integrals tends to be a minimum around (1,1,1), which may be inferred from the smooth behavior of all the integrands at this point. Another trend is the increase in CPU time with increasing x since this leads to increasing oscillatory behavior of the integrand. With increasing z, computer times decrease for Go (larger effect of the exponential decay factor) and GX (smaller values of x/p) but tend to increase somewhat for Go (more confine- ment of the exponential factor to a narrow region near w = cos(cx/2)). It may be noted that in the analysis given in [9] for W. which is closely related to Go' the calculation is limited to values of y /z less than tan (86.4°) = 15.9. 4. Numerical Results We first present various checks which have been made on the accuracy of our formulation. We then present con- tour plots of the various integrals to show their overall behavior. Finally, we give line plots of these integrals to show the finer details of their individual and relative behavior. 4.1 Numerical Checks Newman [6] gives accurately calculated benchmark values of N(x,y,z) for the field points (R. 0,0), (O,R, O), and (O,O,R) where R = 1, 4, 10. The last two sets are most convenient for us to calculate and we shall discuss these first. For both of these cases x = 0, where Go = 0, and Eq. (34) then shows that N _ Ge. In addition, since x = 0, then x /p = 0, the integral GX O. reducing the cal- culation to N= Gc. For an ~ = 0.0001, our calculated values [13] agree with those given in [6] to at least four decimal places. The set (R. 0,0) provides a more thorough check on the accuracy of our formulation since all three of the integrals Gx' Go, and Go would be involved. Unfortunately, our upper limit of integration for Go ~ as z0. We perform this check in two ways. One way is for us to calcu- late the three integrals and form N according to Eq. (34) for small values of z. For an ~ = 0.001 and z = 0.01, we obtain agreement to at least two decimal places with the benchmark results [13]. While this method is computation- ally expensive and inefficient, it does show the stability of our formulation as z0. The second and more convenient method is to calculate N(x, O. O) with the aid of Eqs. (39) and (41) for Go (x, O. O) and GX (x, O. O), respectively. Using this approach, our results agree with those in [6] to six decimal places for ~ = 10-6. 4.2 Contour Plots Figures 7a to d respectively show contour plots for the component integrals Gx, Gc, Go' and N for the case of z = 0.1 over a 101 x 41 rectangular grid with 0 s x s 50, and 0 s y s 20. To illustrate the effect of z, Figures 8a and b respectively show Ge form = 0.1, and 1.0 for the same rectangular grid. Figure 7a shows that the largest waves Of Ox are con- fined in a triangular region which does not extend to the x axis. Figure 7b shows that the wave pattern for Gc is rela- tively simple, consisting of waves whose crests are nearly parallel to the y axis. Figure 7c shows the familiar far field wave patterns confined to the Kelvin sector. Figure 7d shows that the near field integral N has contour lines which are nearly circular and asymptotically decay as 2/R. Figures 8a and b show the well known disappearance of the clutter corresponding to the short waves as z Increases. 141

/ ~# ~^ a = Figure ah. Conlour Plol of at, z = a. ^ Z = ~ Figure ad. Contour Plot of ad, z = 0.1 Figure 8a. Conlour Plot of at, z = 0.1 ~ ~^ z = 142

4.3 Detailed Line Plots Figures 9a and b respectively show line plots of the component integrals Gx, Gc' and Go at y = 0 and 4 for z = 0.1 and 0 s x c 20. To obtain a view of the behavior of the integrals over a wider range of x, including upstream values, Figures 10a and b respectively show plots of Go at y = 0 and 20 for z = 0.1 and100 c x c 100. 3 . ~ . O ; . O i.0 N\/'; 3.0 G x G c - Go - \ / \\ / \ / 0.0 4.0 8 . 0 A X I A L D I S T A N C E X Figure 9a. Line Plots of Go, Go, Go; y = 0 5.0 _ Z 3.0 _ o _ Z 1 . 0 _ 11 U) _ Z -: .0 _ 3: _ G x G c G o ~_~_~\ ~ ~ ~ yN 16.0 20.0 Figure 9b. Line Plots of Go, Gc' Go; y = 4 5. 0 _ Z 3.0 _ o L) z 11 u) z I1J TIC 1.0 _ I -1 .0 _ -3.0 _ AXIAL DISTANCE X Figure lea. Upstream and Downstream Line Plots of Go, y = 0 z z u) -1 . C LL c: 3.0 1.0 -3.0 _ -5 .0 -60 .0 -20 .0 20 .0 AXIAL DISTANCE X 60.0 100 .0 Figure lob. Upstream and Downstream Line Plots of Go, y = 20 Figure 9 shows that the relative behavior between the component integrals varies with the location y. Figure 9a shows that on the x axis, y = 0, the even integral com- ponent GC resembles the odd integral Go, while the com- ponent GX serves as a correction factor indicating the differ- ences between Ge and Go. At the off axis location y = 4, Figure 9b shows that the reverse is true, i.e., now GX resemble Go' while GC tends to serve as the correction fac- tor. Figure 10 shows the expected trend that Go 0 far upstream of the source. This is a verification of our approach in that the separate calculations of Ge and Go do combine to give the required annulment far upstream. This figure shows that the local disturbance is nonoscillatory, decays rapidly near x = 0, and then shows a slow but per- ceptible decay far upstream. 5. Conclusions A numerical procedure has been developed to calculate the complete wave resistance Green's function Go for the case of a nonoscillating source translating below the free surface. The numerical implementation is based on the unique work of Bessho, who succeeds in representing Go as a single integral in the complex plane. We have recast this integral as three real single integrals Ox' Gc' and Go, where GX + GC corresponds to the double or even integral Ge and the third integral to the single or odd integral in the usual representation of Go. By simple rearrangement, we can also express our results in terms of the more physical near field and wavelike components N and W. Computer time depends on a number of factors, principally the submergence of the source z, the required accuracy a, and the horizontal dis- tances x and y from the source. A number of checks have been performed on the accu- racy of our numerical analysis and computer code develop- ment. For field points located on each of the three coordi- nate axes, our calculated values for N agree well with accu- rately calculated benchmark values. Our comparison for the points on the x axis is facilitated by using specialized limit expressions for y = 0, zO for two of our three integrals. 143

Also, the calculated results show the required mutual annul- ment of Ge and Go far upstream. A number of line plots are given to shown the often dramatic differences in the behavior of the integrands in different regions of the compu- tation domain. Contour and line plots show the relative and detailed behavior of the various integrals. 6. Acknowledgments This work was conducted as part of a research program in free surface and marine hydrodynamics supported by the Naval Research Laboratory and by Code 12 of the Office of Naval Research. The second author performed this work under the sponsorship of the U.S. Navy-ASEE Summer Faculty Research Program. 7. References 1. Michell, J.H., "The Wave Resistance of a Ship," Phi- losophical Magazine, Vol. 45, No. 272, pp. 106-123, January 1898. Eggers, K.W.H., Sharma, S.D., and Ward, L.W., "An Assessment of Some Experimental Methods for Determining the Wavemaking Characteristics of a Ship Form," Transactions of the Society of Naval Archi- tects and Marine Engineers, Vol. 75, pp. 112-157, November 1967. Noblesse, F., "Alternative integral representations for the Green function of the theory of ship wave resis- tance," Journal of Engineering Mathematics, Vol. 15, No. 4, pp. 241-265, October 1981. 4. Noblesse, F., "The Fundamental Solution in the 13. Theory of Steady Motion of a Ship," Journal of Ship Research, Vol. 21, No. 2, pp. 82-88, June 1977. 5. Euvrard, D., "Les mille et une faceties de la fonction de Green du probleme de la resistance de vagues," Ecole Nationale Supe'rieme de Techniques Avancees Report 144, June 1983. 6. Newman, J.N., "Evaluation of the Wave-Resistance Green Function: Part 1The Double Integral," Jour- nal of Ship Research, Vol. 31, No. 2, pp. 79-90, June 1987. Noblesse, F., "The Steady Wave Potential of a Unit Source, at the Centerplane," Journal of Ship Research, Vol. 22, No. 2, pp. 80-88, June 1978. 8. Newman, J.N., "Evaluation of the Wave-Resistance Green Function: Part 2The Single Integral on the Centerplane," Journal of Ship Research, Vol. 31, No. 3, pp. 145-150, September 1987. 9. Baar, J.J.M. and Price, W.G., "Evaluation of the Wavelike Disturbance in the Kelvin Wave Source Potential," Journal of Ship Research, Vol. 32, No. 1, pp. 44-53, March 1988. 10. Bessho, M., "On the Fundamental Function in the Theory of the Wave-Making Resistance of Ships," Memoirs of the Defense Academy, Japan, Vol. IV, No. 2, pp. 99-119, 1964. 11. Ursell, F., "Mathematical Note on the Fundamental Solution (Kelvin Source) in Ship Hydrodynamics," IMA Journal of Applied Mathematics, Vol. 32, pp. 335-351, 1984. 12. Abramowitz, M. and Stegun, I.A., Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, U.S. Government Printing Office, Washington, D.C., June 1964. Wang, H.T. and Rogers, J.C.W., "Calculation of the Odd and Even Integral Components of the Wave Resis- tance Green's Function," Naval Research Laboratory Report 6411 (in press). 14. Wehausen, J.V. and Laitone, E.V., "Surface Waves," Encyclopedia of Physics, Vol. 9, Springer-Verlag, Ber- lin, pp. 446-778, 1960. 144