Numerical evaluation of the centerplane potential of a centerplane source distribution

NUMERICAL EVALUATION

OF THE CENTERPLANE POTENTIAL OF

A CENTERPLANE SOURCE DISTRIBUTION

by

F. Noblesse

The preparation of this report was supported by a grant of the National Science Foundation (Grant No. ENG75-03974)

IIHR Report No. 179

Iowa Institute of Hydraulic Research The University of Iowa

Iowa City, Iowa

ABSTRACT

The paper is concerned with the development of

a

numerical procedure for evaluating the linearized velocity potential

4)(x,y,o),

on the plane z=0, for the steady free-surface gravity flow of a Uniform stream past a given distribution of sources on the same plane z=0 assumed to be vertical and parallel to the oncoming stream (this plane z=0 is

thus a plane of symmetry for the flow, which is referred to as the center- .

plane). The centerplane Havelock source potential i8 written in the form

G(x,y,ovp,V,o) =-1/r + Q(x',y')/r' + W(x',y') where r is the distance between the field point (x4y,o) and the source (u,v,o), x'

Ex-.t, y' E

y+\),

,2 2

r' =- (x + yl ) is the distance between the field pOint and the image

(p,-v,o) of the source with respect to the free surface, and W(X1,y') is the potential of the wave-like disturbance trailing behind the source. The strength of the image singularity Q(x',10) tends to +1 as r'40 and to

-1 as r140D. The numerical procedure is based on universal matrices which

may be computed Once and for all and stored on magnetic tape or disk.

ACKNOWLEDGEMENTS

The preparation of this report was supported by a grant of the National Science Foundation (Grant No.ENG75-03974). Additional funds for computer timewereprovided by the Graduate College of the University of Iowa. The author i8 grateful to both institutions for their support.

I also wish to thank Drs. Landweber, Dagan and Tatinclaux for their advice in the development of this numerical procedure.

Fig. 2 Integration contours 20

Fig. 3 Strength of the image singularity 21

Fig. 4 Wave integrals 22

LIST OF CONTENTS

Page

INTRODUCTION

THE CENTERPLANE HAVELOCK SOURCE POTENTIAL

THE CENTERPLANE POTENTIAL 6

OUTLINE OF FUTURE DEVELOPMENTS 11

Appendix 1. Expression of the near-field disturbance in 13

terms of the complex exponential integral

Appendix 2. Numerical evaluation of the strength of the image 14 singularity, Q(x',y')

Appendix 3. Numerical evaluation of the wave integrals 16

REFERENCES 19

FIGURES 20

LIST OF FIGURES

Fig. 1 Definition sketch for the centerplane Havelock 20

NUMERICAL EVALUATION OF THE CENTEAPLANE POTENTIAL OF A CENTERPLANE SOURCE DISTRIBUTION

I. INTRODUCTION

The success enjoyed by Gadd [1] in predicting the "residuary resistance" of a number of ship hulls by the method of Guilloton [2], as well as studies by Dagan [3] and Noblesse [4] rationalizing Guilloton'S method and suggesting Possible modifications, have called for the develop-ment of a numerical procedure for evaluating the flow around a ship hull

by a Guilloton-type method. The motivation for developing such a numerical procedure has recently been strengthened by further theoretical

develop-ments, of Noblesse and Dagan [5] pointing to the possibility of obtaining

a complete,second7order thin-ship approximation by means of centerplane

source distributions alone.

From the computational point of view, the essential feature

of Guilloton'S and related methods is that they involve source distributions

on the Ship centerplane alone,

_and

furthermore that most quantities of practical interest, e.g,, wave resistance, hydrodynamic lift and moment, sinkage and trim, velocity distribution on the hull and free-surface profile along the hull, can be obtained from flow variables evaluated at the centerplane. This means that the computational problem in these methods of evaluating the (three-dimensional) flow around a ship hull is two

dimensional in nature, in that functions of two variables only need be

evaluated,

as

will be shown below.

The velocity field in the centerplane, taken as the plane z=0,

can be obtained by numerical differentiation of the centerplane Veldcity potential 0(x,y,0), and other quantities of interest such as the Vertical

2

and longitudinal zapping functions (displacements) can likewise be

easily derived from x,y,0). The basic computational Problem therefore consists in evaluating the centerplane velocity potential 0(X,y,0),

simply denoted by

0(x,y)

in the sequel, due to a given source distribution, of strength T(x,y), On the centerplane z=0. The present report i8 then

concerned with the numerical evaluation of' the centerplane potential 0(x,y), given by

0 (X,Y) = (Jc,y;11,v)T(p,v)dpdv

Where G(x,y;p,v) denotes the cehterplane Havelock source potential G(x,y,0;p,v,0) and is a given area in the lower half plane V<0 where.

the source distribution

T(p,v)

is defined. The OpreSsion of G(X,y;p,V), as well as the precise meaning of every symbol appearing in the above

equation, is given in the., sequel.

II. THE CENTERPLANE HAVELOCK SOURCE POTENTIAL

Let us consider the steady free-surface gravity flow of a uniform stream of Velocity U past a unit source located at a point S below the

free surface, see Fig. 1. The vertical plane passing through the source S and parallel to the oncoming stream U is a plane of symmetry for the flow. This plane is taken as the plane Z=0 and, in reference to the thin-ship'

theory, is referred to as the "centerplane". The horizontal undisturbed free surface is taken as the plane Y=0 with the Y axis directed upwards. The X axis is taken parallel to the oncoming uniform stream U and pointing

downstream. The (dimensional) coordinatesof the source S are denoted by

{XS' YS<0,0} and the point P where the flow is observed has coordinates

{X,Y<0,0}. The coordinates of the image S' of the source S with respect

to the free surface are {XS'

3

Variables are Made dimensionless'

with

respect to the acceleration of gravity g and the velocity U taken as reference quantities. Thus,

dimensionless coordinates x,y,p,v are defined as

= (g/U2){X,Y,X81Y } (2.1)

The linearized velocity potential, on the centerplane z=0, for

the flow, described above is given by the "centerplanen Havelock source

potential G(x,y,0;0,v,0), Which is

simply

denoted by G(x,y;p,v) and is written in dimensionless form as

G(x,Y;11,v) = -1(x-P)2+(y-v)21-1/2 + R(xp,Y+v) (2.2)

where R(xp,y+v) is a function of the two variables xp and y+v defined -and regular - in the lower half space y+v<0. We have

. k[(y+v)+i(X-p)cose] , 2 2,-1/2 4

f

dkdeir/2(co p R(x-p,y+v) = [(x-p +(y+v) j + 0 TO k cos2e-1.

Where the convention is made that only the real part of the right-hand

side is to be taken. Equation (2.3) can easily be derived from eq. (13.36)

p. 485 in Wehausen [6] (it must be remembered that the equation given by Wehausen holds for the case of a uniform stream flowing past a sink in

the negative x direction).

Let us introduce the complex number C defined

by

= (y+v)sec20 + i(x-p)sece (2.4)

Eq. (2.3) can then be written as

R(x-p,y+v) [(x_u)4.(171.v)2]-1/2.2 + 4/11./2I(C)sec2ede

+4i

7412e sec2ede (2.5)

7 0 0

2 e[(y+v)+1(N-0-coselsec2e 2

where the change of variable k=y sec28 was performed in the double 1 i ex. eYC 'I(C) = t, ---- dy (2.6) 0 y-1

The Cauchy-principal-value integral I(C) can be written in terms of the complex exponential integral El(C), see Appendix I. It yields

I(c) = e

B1(c) +

ihsgn(x-p)e ,(2.7)

where sgn(x-p)=( -p)/lx-pl and integral, with I(C) defined as

ro, e-t

E1(C) =c -t

dt

as defined

in

eq. (5.1.1) P. 228 in Abramowitz and Stegun [7].

By using eq. (2..7), eq. (2-5) can be written as

R(x-p,y+v) = N(x7P/Y+v) + W(x11.,y+v) 2.9)

where N(..x-p,y+v) and W(x-p4.114.4 are defined as (larg0<w)

(7/2 C

w(x=p,y+v)

= 81H(x-p) e sec2Ode

'

In eq. (2;11), H(xp) is the Heaviside step function defined as li(x-p)=01

for x<u and H(x-p)=1 for x>1.1. Only the real parts of the right-hand

sides of eqs. (2.10) and (2.11) are to be taken, in accordance with the

convention adopted in eq. (2.3). It is easily seen that N(p-x,y+v) =

N(x-p,y+v) so that this term represents a disturbance symmetric upstrearn and downstream from the source. Furthermore, N(x-p,y+v) represents a

nonoscillatory "near-field disturbance", as shown below. The term W(x-p,y+v),

on the other hand, is zero upstream from the source and represents an oscillatory wave-like disturbance. In particular, far downstream i from the

(2.8)

(2.11)

N(x-p,y+v) = -P) +(Y+v)2 2,-1/2J + --J4 (7r/2 c

e-E1 (c)sec2.8de

source, W(x-p,y+v) behaves like

- 8(ir/2)1/2 ey+v(x-)-1/2sin(x-p+7/4) (x-p>>1)

as can readily be obtained by the method of stationary phase. Hereafter,

the terms N(x-p,y+v) and W(x-p,y+v)are referred to as the "near-field"

and "wave" terms, respectively.

A good picture of the behavior of the near-field term N(x-p,y+v)

is obtained by writing it as follows:

N(x-11,y+v) = Q(x-p,y+vTnx-1.1)2+(Y+v) J

where.Q(X-pty+v) is defined as

= 1 + !=[(x-p)2+(y+v)2]1/2f7r/lie{eC 2

El(C) }sec ede

The function Q(x-p,y+v) has been evaluated numerically, see Appendix 2, and is represented in Fig. 3 for 0<x-p<3, 0>y+v>-3. It is seen in Fig. 3 that Q(x-p,y+v)+1 as r'E[(x-p)2+(Y+v)21 ÷0. On the other hand, Q(x-p,y+v)÷-1 as r'->00, as one readily obtains by substituting

eE1

() in eq. (2.13) by

the first term, 1/C,

in

the asymptotic expansion of _{eE1-(C) for CI}>>l, see eq. (5.1.51) p. 231 in Abramowitz and Stegun [7]. Keeping two more terms

in this asymptotic expansion, i.e., taking

eE1

(0% 1/C-1A2 +2/C3, yields

2 X-

1 6 y7fy, _2-y'/r',

rl 1-y'/r' ) i77-2Er"r' _(1-y7r1)21

where x'Ex-p, y'Ey+vi and r'=(x'2+YI)1/2 as already defined above.

In summary, by using (2.2), (2.9) and (2.12), the centerplane

Havelock Source potential G(x,y;p,v) can be written in the form

G(x,y;p,v) =

-.r +

x

)r

+ W(x',y') (2..15)

where rEpx-u)2.4-(17.1.v) 2'1/2, _{and Q(x',y') and W(x',y') are giVen by (2.13)}

and (2..11), respectively, _{Thus, the velocity potential of a submerged}

(r ' >>1) (2.14) (2..12)

6

source in an oncoming uniform stream is written as the sum of three

terms: (i) the potential -1/r of the source in an infinite domain,

(ii) the near-field disturbance potential Q(x1iy1)/r' which may be interpreted as the potential of an image sink of strength 12(x!,y1) in

an infinite domain, and (iii) thepotential W(x',10) of the wave-like disturbance trailing behind the source. It is convenient to refer to

the function Q(x%y!) as the "strength of the image Singularity". Sindel

Q4-1 as r'-)40 (i.e., as the field point is very far from the source)

and Q4+1

as r'--)-0

(i.e.,

as the source and the field

point

approach each

other and the free surface),. the "Image singularity" of a unit source Is seen to become a unit source or sink as r'4-0D or r1÷0, respectively.

III. THE CENTEPPLANE POTENTIAL

We start by defining notations which are used in several instances

in the present Section. Let a continuous function f(x,y) be given. We

pq

define four functions f (x,y;h,k), (p,q=0,1) in terms of the values of f(x,y) at the points x±h, y±k by the following expressions

4f00 (X,y;h,k) = f(x+h,y+k) + f(x-h,y+k) + f(X-h,y-k) + f(x+h,y-k) 4hf10 '(x,y;h,k) = [f(x+h,y+k)-f(x-h,y+k)] + [f(x+h,y-k)-f(x-h,Y-k)] 4hkf11 (x,y;h,k) = [f(x+h,y+k)+f(x=h,Y=k)] - [f(X-h,y+k)+f(x+h;y-k)] (3.1) 4kf01(x,y;h.,k) = [f(x+h,y+k)-f(x+h,y7k)] + [f(x-h,y+k)-f(x-h,yk)]

7

--It is readily seen that f00 (X,y;h:i.k) represents an average value of

the function f(x,y) over the rectangle x4h<x,<x+h, y-k<y<y+k, while

f10

and

f01

are average slopes in the x and y directions, respectively,

and

f11

gives an average value of the second partial derivative f,. xy We now proceed with the proper subject of the present section, which is the numerical evaluation of the centerplane potential

0(x,y)

defined in eq. (1.1). Let the centerplane Green function G(x,y;p,v) be

written in the form G(x-p,y±v). The centerplane potential

0(x;y) is

then given by

0_(c,y) = if

G(x-p,y±v)T(p,v)dpdv

a

0(x,y)

is to be evaluated

at

the nodes of a rectangular numerical mesh x <x <...<x

<...<xM+1

and y1=0>y2>...>yn>...>17N4.1. The value of

0(x,y)

1 2 at the node x=xm,

y=yn

is given by M N

.0(x,yn)

= E E. i=1 j=l ijmn Where

0_.

is defined by ijmn X'+1 Yi 0

=1

-

dpif

dv' G(xt -p',Y -±v') T (111 ijmn _Yj+1 n 1

w"introducethequantitiesx.,y.,h.andk_defined as

1 1

(x.+x.1

)/2, y. = _)/2 1 1 1+ 3 7+ (xi.+1-2Y/2, = (Yj'''Yj+1) /2 (3.2) (3.3) (3.4)

_

With

pg = x.4-p

and v

=T+v,

eq.

(3.4)

then becomes

1 3

hi fki

0. ₎

'

du G(x -;,.-11,Y _{7±i7.±9) T6-c-.4.11,7.1-v)}

ijmn h. -IC_ m 1 _{3 1} 3

1 3

In each rectangle -h.<p<h , -k.<v<k., the source density function

3 T _Y,A-V) is approximated by

(3.5)

1 1

T(Tc+pi-17.+v) =

Z E Tlx167".;h.,k.)PPvg

1 3

p=0 q=0

1 3 1 3

We introduce the following notations

. = x -x., y . = y

-y., y

y

-mi M 1 n3 n 3 nj n

We also define the quantities

h. k. 2 ,2,-1/2 p q J -P V pq. Mi nj i h. k. mi YA] 3 )c.

Npq (xmij ;h.,IS-)n 3 h.1dPI 3dvisilx .-P,Y!.+v)p9vcik. . mi n3

3

h. k.

Wpqmi j13

:h.,k.)

= I

ldpl 3dvW(x .-11,10.+v)p

n h. k. mi n3

a. 3

where N(x-p,y' +v) and W(x -p,10 .+v) have been defined in eqs. (2.10)

mi nj mi n3.

and (2.11). Eqs. (2.2), (2.9) and (3.7)-(3-11) show that eq. (3.6) becOmes

1 1

4. E E

ijmn p=0 q=0 1 3

1.

3 ijmn

Pq

-where yb.. is defined by 13mn

Pq

O. = -S (x ,y ;11.,k.) + N (x ,10 ;h.,k.) + W (x ;h.,k.) (3.13)

ijmn

pq mi

nj 1 3

pq mi

nj 1

pq mi

nj 1

Closed-form analytical expressions for the "source terms"

(x,y

;h.,k.) can easily be Obtained- The "neat-field terms"

pq mi

nj 1 j

N, (x,y' Jh.,k.) can be expressed in terms of the "source terms"

pq MI

nj 1 3

(x;h ,k.)

by

writing the near-field disturbance N(x --P,Y';+v) as

pq mi

nj 3 mi nj

in eq. (2.12), and approximating the function Q(xmi-p,y;ii+v) by

1 1

r rSr s

Q(xi43

=

r0 s0

E E (-1) Q (x

mi

,y'nj1

01.,k.)44 v

m==

in each rectangle -h.<p<ita., -k <v<k.. _{The coefficients Qr}S_{(x ,y'} ; ,k )

1-7 mi nj j

(3.7)

F

where the four coefficients T (x,,Y.;h:',k.) are 'related to the Values of

3.

the source density function

T(x.±hyj ±Xj1 1 - ) at the corners of the rectangle

by means of eqs. (3-1)

(3.8)

(3.12)

are ielated to the Values of the function Q(X ±h.,y' ±k.) at the

mi 1 nj

corners of the rectangle by means of eqe. (3.1). The "near-field terms"

N (x .,Y1.;h.,k..) then become

pq mi

n3 1 3

N (X,10 ;h.,k.) =

pq mi

nj 1 3

obtains

I = cos3efe(y+v)sec28 + i(x-p)secelh,k

00 -h -k

I10 =

3

cosOfe(y+v)sec28 + i(xp)sec8 (p-i

cos8}

h,k

-h,-k -,

3 (y+y)sec28 + i(x-p)sec8 h,k

101₀₁ cos Afe (v-cos20)1_-h,-k

I11

cos38{e(y+v)sec28 +)sc8(p-i

cos8)(v-,cos-8}-n,-k

0 h k

where the arguments ofI pq(0;x,y,h,k) have been omitted for shortness..

We introduce the quantity Q (x,y) = Qc (x,y) + i (x,y) defined as

R(x,y) = Qc (x,y) + i.Rs (x,y)

rw/2

e

sec28 + ix

sec°cosP8d8 (3-20)

p

where y<0, x>0 and p>1. Equations (3.17), (3.18) and (3.19) then give

9

1 1

rrs

E (-1) Q (x,y' ;h.,k.)

s

(x ,y'

r=0 s=0

i

1 g p+r,q+s ml nj j

These "near-field terms" can thus Simply be expressed in terms of the

"strength of the image singularity" Q(x1,y1) and the nine quantities

h yk ,-4 .

S(X,17;,h,k) = r dpl dVL(p-x)2 +(v-Y)2J PaV ;

h

which can be written in closed-form analytical expressions.

For lxi

_{m -}

1>h.,

eqs. (2.11) and (3-11) yield

W (x,Y1 ;11,,k,) = 8i 171(x .)f7/2I (8;x ,y. ,h.,k.)

sec2ea

pq mi

nj 1 g. mi 0 pq mi

j

3

where the quantity Ipq(8;x,y,h,k) is defined by

I rh -1(P7x)seceppdv rk e(v+y)sec28vcidv

(;x,y,,

_{1h e} I

pq8

hk) =

Lk 1

I It is emphasized that eq. (3.17) is not valid for ix .1<h.. One easily

mi 1 (3.15) (a,13=0,1,2) (3.16) (3-17) (3.18) (3.19)

W (x= 32h.k

H(x. ) 0c11(x 00 mi n3 1 3 1 3 mi 1 ml n3 1 3 W = H(x .)(-c01+ 0s ) 10 1 mi 1 211 11 W = 32h k_ H(x .)(0c10: -

cc)

01 3 mi 1 - ' W = 32h.k. H(k .)(0c00+ 0s +0c _{-0s411 )} 11 mi 1 2 10 3 01 10

where the notations defined

in

eqs. (3.1) have beenused. Thus, the.

"wave terms" W (x,y' ;h.,k ) are expressed simply in terms of the four

pq mi

nj 1

integrals Rc1, QS2, 0c3 and Qs4. These integrals are referred to as the "wave integrals", and are examined in Appendix 3.

In summary, the centerplane potential

4)(x,yn)

is given by

eqs. (3.12), (3.13), (3.15) and (3.21). These equations involve

the following basic quantities: the "strength of the image singularity"

Q(x-11.,.y+v) defined in eq.. (2..13), the nine "source terms" Sa6(x,y;h,k)

defined in eq. (3.16) and the four "wave integrals"- p (x,y) defined in

eq. (3.20). The numerical evaluation of the functions Q(x-1.1,y+v) and

2 (X,y) is discussed in Appendices 2 and 3, respectively.

A similar numerical procedure has previously been developed by Yeung [8] for the calculation of sinkage and trim. The present procedure follows closely Yeung with respect to the treatment of the "Wave disturbance"

(the wave integrals

c(x)

defined in (3.20) are identical to

Yeung's

c,s

J integrals), but the "near field disturbance" is treated somewhat

P+1

II

IV. OUTLINE OF FUTURE DEVELOPMENTS

The present stage of development of the numerical procedure is essentially as described in this report, with the strength of the image

singularity, Q(x-p,y+v), and the wave integrals Qci(x-p,y+v), S2s2, Rc3

and Os4 evaluated at a sufficient number of points (x-p, y+v) for these functions to be defined with accuracy 10-6 by means of a bicubic spline

interpolator routine. The computed values of Q and of the four wave

integrals are stored on a magnetic tape.

One way, probably the shortest way, from this stage on,Would be to select

a

constant-size rectangular numerical mesh, i.e., hi E h

for i = 1,2,...,M and k. E k for j = 1,2,...,N, and .form the eight matrices

Npq (xmi ,y'j01,k) and W

(x ;h,k),(p,q = 0,1) defined in eqs. (3.15)

n pq mi nj

and -(3.21). For given values of h and k, these matrices are universal,

and thus need be constructed once and for all, then stored on a magnetic tape (if desired, the four matrices S (x ,y ;h,k) can also be evaluated

pq mi nj

and stored on a tape). With these matrices formed, the evaluation of the potential (1)(xm,yn) at any mode

(x ,y )

of the numerical mesh amounts to

m n

simple additions to be performed by the computer, as shown by eqs. (3.3),

(3.12) and (3.13)..

In order to cover the range of Troude numbers of practical interest,

say 0.1 < Fr < 1, several valtes of h and k must be considered.. It would

seem sufficient to consider four cases, and thus four sets of matrices:

(i) h

= 0.025, k = 0.0125 for 0 < x' < 5, 0 > y' > -1, (ii) h = 0.05, k = 0.025 for 0 < < 10, 0 > y' > -2, (iii) h = 0.1, k = 0.05 for

0 < x' < 20, 0 > y' > -4, and (iv) h = k = 0.1 for smaller Froude numbers. This represents fairly large matrices, containing some 16,000 elements each.

12

The size of these matrices appears to be too large for direct storage in, . the memory of the computer (for a given Froude nttber, hence for given and k, eight or twelve of these matrices should be Stored in the memory Of

[ the computer, as discussed above). For permanent record, these matrices';

would be stored on Magnetic tape, and temporary storage on disk could be

used for the purpose of: numerical calculations.

Although the procedure deSCribed

above

Should be feasible, the access time from the disk to the computer might be significant. In addition,

this procedure is somewhat rigid due to the fixed-size uniform numerical'

mesh. It may thus be worthwhile to develop asymptotic expansions and

rational approximations for the near-field disturbance

N(x-1,y+\)

and the four wave integrals. Figure 3 and eqs. (2.12) and (2.14) suggest that the

function N(x',y') should be written as

l/r' + f(17'/ ') + r'f1(171/r1, r') for r' small

N(X'ilru) =

71/r.' + E a(y'/r')/r for r' large

k=i

where a (y'/r') are the functions appearing in the asymptotic expansion of N(x',y') for large r', as in eq. (2.14), and fo(y7r1) and

f1(y7r1)

I

denote functions of y'/r'. A preliminary analysis shows that four or five additional terms in the asymptotic expansion (2.14) would yield an approximation to the near-field disturbance N(x',y') with errors less than

-6

10 for rl> 20 or even 15. Asymptotic expansions for the wave integrals

13

Appendix 1. Expression of_the_near=field disturbance in terms of the

complex exponential integral

Let the complex number C, defined in eq. (2.4), be written in

the polar form c=aeia, where a>0 and 7/2 <a<37/2 since y+v<0 and x-p>0. _

Following Hess and Smith [9] and others, see Yeung [8], the Cauchy

pincipal value integral 1(C), defined in eq.. (2.6), can be expressed

as the sum Of an oscillatory term and a nonoecillatory term by means of

an appropriate contour integration in the complex y plane. Let the

p

complex variable y be written in the polar form

y=peit

. We then have .

i(a+tp)

yC=ape . The exponential.

teraLelC

In-the integral (2.6) can be

rendered monotonic decreasing by choosing 4)=7-a, which yields yC=-ap.

Performing the contour integrations shown in Fig. 2 then yields

=ap e (.3.7-,-a). T

T=1

Y 17e +

f°

pe

e

dp = 0 -1

Where the - and + signs correspond to the upper and lower contours (case x->0 and x-p<O), respectively. This readily gives

cc

,ap

= ef dp

0

p+e

Equation (2.7) is finally obtained by successively performing the change of variables X=ap and t=A+C in the above integral.

14

Appendix 2. Numerical evaluation of the strength of the image singularity

For simplicity, we use the notations x'=x-p,

r=y+v

already

introduced in Section II. The complex_number C defined in eq. (2.4)

.

then becomes C=y'sec20+1 x'sece, where y'<0. We have eCE1(E) = eCF1(C).

see Abramowitz and Stegun [7] p. 229, and eq. (2.13) then readily shows

that Q(-x',y') = Q(x',y'). The numerical evaluation of Q(x',y') may thus

be restricted to x'>0, y'<0.

As 6 increases from 0 to V2, the point c moves along the

parabola C =(177)02)C2 from C

=10+

ix' to infinity, in the second quadrant

of the complex plane C = Cr+iCi = y'sec26+i x'sece. The rate at which

moves along this parabola, given by dC/d6, depends on y'+ix'. In particular, dc/d6=0-0 as y'

Performing the change of variable t=C-:tan 6, with c an arbitrary

constant, in eq. (2.13) yields

Q(x'iY') = 1 +

f

()}

dt

c 0 '

_

where C. =_

-(t2

+c2) ) ia(t+c2)1/2 with

x'/c

and 13=-,y'./c2. The

efficiency of the numerical evaluation of the above integral was greatly improved by using the asymptotic expansion of eEl(C), see eq. (5.1.51) p. 231 in Abramowitz and Stegun [7]. The first three terms in this

asymp-totic expansion were retained and Eq. (A2:.1) was Written as

+

ft°

Re{eE1

(c)} dt tr c 0

4rtr.

2

, +

} dt ) C j to

A closed-form analytical expression was developed for the last integral

(the asymptotic expansion (2.14) was obtained

from

this expression by

setting t0=0). The value of to was so chosen as to yield Q(x'iy') with

(A2.1)

(A[2.2)

- : 15

an error less than 10-6.

The first integral in eq. (A2.2) was evaluated numerically.

-7

Simpson's integration

rule

was used; the specified accuracy was 5.l0. The integraMA.El(C) was evaluated by means of the rational approximation developed by Hershey [10], supplemented by the asymptotic expansion (5.1.51) in Abramowitz and Stegun,, for large values of

Id.

The values of the

exponential integral obtained by using Hershey's rational approximation were checked against the values given by Abramowitz and Stegun [7].

The specified accuracy for the inteqrand

eCE(0

was 10-8. For the

numerical integration of the first integral

in

eq. (A2.2), it was

found

convenient to select the constant c as c=r'/(5+r'), which remedies the difficulty Mentioned above, i.e., that the rate at Which the point C

2

moves along the parabola Cr=(171/x12

)Ci vanishes as C40, and yields a

relatively uniform value for the upper limit of integration to (t varies

o

between 4 and 20, the higher values being attained for snaller values of

y' and r'). The computing time on the IBM 360/65 at the computer center of the University of Iowa varied between 0.1 sec to 1 sec per value of Q(x',y') depending on x' and y', the computing time increasing fairly

rapidly as r'-)-0. This represented a cost varying between less than

Appendix 3. Numerical evaluation of the wave integrals

For the purpose of numerical integration, the change of

variable sinh t = tan e was performed in the wave integral (3.20) whicli

was then written as

-(y)

eix

- 9

where t.

Y .

cosh . [(log

E)/y]23--1 1/(p+1) cosh [(1/C ) I, Finally, 16 2 icosh t ix(cosh t-1) _e p+1 co$h t

The integrand in the above integral is the product of two

Monotonic decreasing functions, namely --exp(y cosh2t) and SechP-Ilt,

and the oscillatory function exp[ix(cosh t-1)]. Each of these three

Unc-tions depends on one parameter only, i.e., y < 0, p > 1 and x > 0,

respectively. It is readily seen that exp(y cosh2t) < ey if t > ty

Similarly, the dt sech t < e if t > t P P th

k oscillation of the oscillatory

-function exp[ix(cosh t-1)] occurs at t = tx,k = cosh 1(1 + 2k7r/x), k =

V(p+1)

It is then easily seen that t < t if y e log c.

Y P -P

This shows that - except for y = 0 - the function exp(y cosh2t) may be considered the fastest decreasing function. Next, by comparing t

,k and

x

t, it is found that exp[ix(co$h t-1)] undergoes at least 5 oscillations

-7 -7

(k = 5) 'before exp(y cosh2t) decreases to 10 (e -= 10 y if x > 107r/

[4(-y)-3/4=1],

in

which case the limits of integration were selected as

t k ,= 0,1,2,..., For-x < ift/[4(-1)4-11, the fixed integration step

At = t

Y/5 (with eY

= 10-7) was used. For the case y = 0, the integrand

undergoes 5 oscillations at. least before p+1 O

-6

if x > .125. Thus, for y = 0, the limits of integration were selected

astk = 0,1,2,....

x,k'

suMmary, two cases were considered: (i) for x < 107/[4(-y

integration

(ii) for x > 10*/[4(-y)-2-1], the limits of integration were Selected as

tk1

x, k = In both cases, Simpson's integration rule was used

and the integration was pursued until the contribution

of

the last inte-gration step-was less than 10-7. The upper limit of integration was

found to vary between about 1 and 4 (except for x = y = 0 Where it is

about 9). The computing time required to evaluate the four wave integrals

varied between less than 0.04 sec

tp

to 0.25 sec, representing a cost varying between 0.4

up

to 2.

The numerical calculations were Checked against the following

results:

QC 10) == .1, (0,0).= 2/3i Rs (0,Y) = Rs (0,y) = 0,

1(0, Rc3 2 4.

QC

(0'

y) - -Y_ _eY/2

1 2 : [K1(1.) -

K0(1)1

Rc3(0,y) = y e (1 - K (-Y-) + K (--Y-) ]

3 1 2 3 2 2

X

Qc (x,0) = --x Y (t)dt - Y (x)]

1 2 o 1

Jo

Where. Y , Y

, K

and Ki are the usual Bessel functions, see Abramowitz

o 1 o

and Stegun [7].

The behavior of the four wave integrals _{c1, Rs2, Rc3 and Qs4} as functions of x-p for y + v = 0, and as functions of y + for X-p = constant, is depicted on Figs. 4.

For large values of x, we may obtain an approximation to these wave integrals by means of the method of stationary phase, see Erdelyi

[11]. It is Convenient to perform the change of variable t = tan

e in

(3.20).

This yields 00 eyt2

et

(xiy) = 1 dt 1 tp

(t-471)2 (t-1)

18

which is in. the form of eq. (4) p. (48) in Erdelyi.

We easily obtain

0 (x,y)

"TA dY .i(x+71/4)

{1

i

1/4

- -

p - 2y

P . 2 2x

3 3 + 24p_+ 16p2 + 16y(1-4P+4Y)

...)

8

2

19

REFERENCES

Gadd, G.E., 1973, "Wave resistance calculations

by

Guilloton's Method", Trans, RINA, Vol. 115, P-

317-Guilloton, R., 1964, "L'gtude thgorigue du bateau en fluide parfait",

ATMA, Vol. 64, P. 537.

Dagan, G., 1975, "A method of computing nonlinear wave resistance of thin ships by coordinate straining", to appear in the J. Ship. Res. Noblesse, F., 1975, "A perturbation analysis of the wavemaking of a ship, with an interpretation of Guilloton's method", to appear in

the J. Ship Res.

Noblesse, F. and Dagan.G., 1975, "Nonlinear ship waves theories by continuous mapping", submitted for publication.

Wehausen, J.V. and Laitone, E.V., 1960, "Surface waves", Handbuch der Physik, Springer-Verlag, Berlin,

Abramowitz, M. and Stegun, I.A., 1964, "Handbook of mathematical functions", Dover, New York.

Yeung, R.W., 1972, "Sinkage and trim in first-order thin-ship theory",

J.S.R., Vol. 16, p. 47.

Hess, J.S.,and Smith, A.M.O., 1967, "Calculation of potential flow about arbitrarybodies", Progress in .Aeronautical Sciences", Vol. 8, Pergamon Press, New York.

Hershey, A.V., 1959, "Computing programs for the complex exponential integral," NPG Report No. 1646, U.S. Naval Proving Ground, Dahlgren, Va. Erdelyi, A. 1956, "Asymptotic expansions", Dover, New York.

Fig. I Definition sketch for the

centerplane Havelock source potential

20 V

y

SI(Xs 1Ys

10

S101.9-11,0

S(Xs,Y 8,0)

S(1L,V,O)

P(X,Y,o)

P(x,y,o)

Fig. 2 Integration contours

case x >IL

1,47171.5"---Assts-2,

Artara.

/00000/AtArAfirvAiAgrAirArAtivavAr

L

iprzeiroArmmi

-Amwamw.w.Aw

Y,,,Artivialaril:1211rA

PWA IA 1r

22

I I

4

5

6

7

8

9

10 11 12 13

Behavior of the wave integrals as functions of xiu. for y +v=0

Behavior of the wave integrals as functions of y+v for

constant