Irrotational flow about ship forms

(1)

IRROTATIONAL FLOW ABOUT SHIP FORMS

by

L. Landweber and M. Macagno

Sponsored by Office of Naval Research Contract Nonr - 1611( 07)

IIHR Report No. 123

Iowa Institute of Hydraulic Research

The University of Iowa

Iowa City, Iowa

December 1969

This document has been approved for public release and sale; its distribution is unlimited.

(2)

-IRROTATIONAL FLOW ABOUT SHIP FORMS

Introduction

A method of computing the potential flow about ship forms would, in spite of the neglect of viscous effects, be valuable in the preliminary design of a ship, or in investigating means of improving the performance of

an existing ship. If an efficient procedure for performing such calculations were available, one could determine, without recourse to model tests, whether the streamlines along the forebody are such that bilge vortices would be generated and what would be the effect on these streamlines of various

modifications of the bow. If separation at the stern is not severe and the bow not too blunt near the free surface, useful results could be obtained for the wavemaking of a ship form on the assumptions of irrotational flow and the linearized free-surface condition, but with the exact boundary

condition on the hull. The last condition is significant since it would enable the effects of local modifications of form, especially at the bow,

to be studied.

A well-known method of calculating the potential flow about a

three-dimensional form is that of Hess and Smith [l]*. In common with the Hess-Smith approach, the method to be described herein determines a distribution

of sources on the surface of a given body by solving the basic integral equation of potential theory for such a Neumann problem [2]. The methods

differ in the treatment of the singularity of the kernel of the integral equation, the selection of an iteration formula for solving the integral equation, the quadrature formula used to reduce the integral equation to a set of linear equations, and the procedure employed to calculate the velocity distribution along the hull once the source distribution has been determined.

Although the procedure to be described has been available as a

computer program for several years, and early results with a body of

revo-lution (a spheroid) and a three-dimensional form (an ellipsoid) showed very

good agreement with the exact solutions, publication has been delayed be-cause anomalous results were obtained when the method was applied to a

(3)

-2-mathematical form with parabolic lines and sections, having sharp edges

at the bow, stern and keel. Presently an attempt is being made to apply the method to a Series-60 model [3]; but since this form is serving as a vehicle for development of a procedure for fitting a mathematical equation to an arbitrary ship form, the hull coordinates and direction cosines of the normals to the hull, required as input in the potential flow program,

are not yet available. Consequently it has been decided to present the method and the computer program without further delay, since others may be more efficient in obtaining the required input for their ship forms.

Statement of Problem

Our problem is to develop means of computing the flow about a

ship form, including free-surface effects. We shall assume the fluid to be

incompressible and inviscid, and the flow irrotational. We shall suppose that a ship form has been prescribed, and that its draft and trim are known. We shall assume that the surface disturbance is sufficiently small that the

boundary condition on the free surface may be linearized. The boundary con-dition on the hull will be satisfied exactly.

An obvious criticism of these assumptions is that, in neglecting the effects of viscosity and sinkage and trim, and employing the approximate

linearized form of the free-surface boundary condition, in comparison with the elegant, classical, thin-ship theory, only one of several equally

im-portant corrections will have been made. Our view is that the solution of the present problem can be made the basis for incorporating additional

corrections by iterative techniques. For example, the resulting pressure distribution can be used to calculate the equilibrium trim and sinkage of the ship, which can then be applied to obtain a second approximation for the

flow about it.

Let f(x, y, z) = 0 be the equation of the hull surface S, with x in the direction of the stream U, and z positive upwards, with the origin

in the undisturbed level of the free surface. Denote the velocity potential

by

(I) = Ux + ct, (1)

(4)

Laplace equation

v2q)

,

0 ₍₂₎

The boundary condition on S is then

4

U

ax

aN BN

where N denotes distance in the direction of the outward normal to S. The

free-surface boundary condition will be taken in the linearized form

3

cb

2

k0'

= 0- _{k0 =} z = 0

ax-

U

For a source of unit strength at the point P(C,

n,

in the same uniform stream, the velocity potential (Ds which satisfies (2) and boundary condition

(4)

may be written in the form

1 1

= Ux - - R' + H(P, Q)

where Q is a point Q(x, y, z) below the free surface, R is the distance from

P to Q, R' is the distance from P'(E,

n, c),

the mirror image of P in the

free surface, to Q,

R = [(x _ + (y n)2 (z _ = [(x

(y

n)2

(z

02]i'

(6)

and H(P, Q) is regular harmonic in the lower half space z < 0 and given by

(3)

(5)

H(P, Q) = iT

-f27 k(z

+ c) +

ik[(x - C) cos e + (y TO sin 0]

kdkde

00

k - k0 sec 0

7/2

(z c)sec

- 4k

f

ek 0 + 20_{sin[k lx - E)sec e]}

cos[ko (v

n) tan e sec e]

0

sec2e de (7)

fwhere

denotes the "Cauchy principal part". The velocity potential of

(5)

and

(6)

_{also satisfies the "radiation condition", that waves are} pro-pagated downstream from the source [4].

Now let M(P) denote the strength of a source distribution on S. The velocity potential of this source distribution which satisfies the

free-surface condition

(4)

and the radiation condition, by (5), is given by

-C4)

)2]32-, E)2 +

(5)

-(P)

(1)(Q) = Ux -

M

dS +

I

M(P) H(P, Q) dS S+S'

where S' denotes the mirror image of S in the plane z = 0. Taking into account the discontinuity in the normal derivative of the potential at a

surface distribution, the boundary condition on the hull surface S then yields

U 9x 3H Q)

M(Q) = 27 BN 27

f

M(P)

[311-] dSp -

T

I M(P) (PTT dS

(8)

DN

S+S'

a Fredholm integral equation of the second kind. The development of a procedure for solving this integral equation numerically for a given ship

form is our principal objective.

An essential difficulty in the numerical solution of

(8)

is that both of the integrands are singular. A means of removing these

singular-ities will be described and justified in the following two sections.

Treatment of Double-Model Integral

At points of S where the normal is continuous, we have by Gauss's

flux theorem

9 (1

-

I

DN dSp = 27

S+S,

This enables us to write

1 1

I

im(p)ii

dSp - M(Q) S+S' S+S M(P) a R Q NQ '

Q=

We shall now show that the singularity of the left member of (10) when P coincides with Q is not present in the right member.

The direction cosines at the point Q of S are given by

fxfz

9n

= M = n =

Q DQ

'QD'QD

where f f f denote partial derivatives of f(x, y, z) with respect to

x y z x, y, z, and 4-D = [f 2 f 2 4. f 2]2 ( 7 )

(9)

[1I]

R dsP - M(Q) (10)

(6)

and

'a

M(P) -3T-T (IT) m(Q)

-er

{R1

--Q

When, P is near Q, we can write the _{Taylor expansion}

1..0'71.(x1y,z)(3cOf--(Y

- n)f

Of +

4[(x

-

c)2f

xx

+ (y - r1)2f + (z, _)2f _{.4-- 2(y}

_n)(z

_4)f -yy zz _yz " Since we obtain and

-5-dfz =

-1

=+(-x)

1 [ 1 x) T;--c

pH+

P {1:):

[(x Idfx

(y _

ofy

(z _ M( Q) [(x

- Of

C

+ ,2

( DP (13) 3rd order terms (1)4) Similarly we have

- x)f + (n y)f

(c -

z)f =

fu

_{3rd order tents}

or (.x - Of (y n)f_ + _{- Of =- -i[Cx - :02f} ₊ _{+ 3rd order} termsE n, xx (15) -Also, we have M(P) M(Q) + - x)m + (n Y)M +

(c -

z)M + Then we have

a.

_{_} 1 '{(xx

--3

-3-- [(x_D '

RQ

-+

Of + (y

rom + _ TIOf -y (t n

of

(12)1 9NQ - _(z -1 - - ₊ -+ -+ - n)fn -- (z

--

+ + f(x, y, z) = f(E, n, = 0 (x - + - _{+ (z -} _- ₊ + - _- ₊ _-

_4[(x

-+ + - (z -= - -

(7)

--6

Substituting the results in

(1)4),

(15), (16)

and (17) into (13), we observe that terms of the first and second order in (x - E), (y -

rh),

(z

-cancel, leaving antisymmetric terms of the third order in the numerator. Since the denominator B3 is also of the third order, the ratio is

indeter-minate as R approaches zero. The integral of this ratio over a small

area symmetric about Q, however, is zero. For this reason we propose to set the integrand of the right member of (10) equal to zero when P coincides

with Q.

Treatment of Wave Integral

As in the prior case, let us consider the integral of the trans-posed kernel, 3H(Q, P)/N. For points Q below the free surface (z < 0), H(Q, P) is a regular harmonic function of P(E,

n, 0

for 0. Hence,

by Gauss's flux theorem, the flux of H(Q. P) through the closed surface, consisting of S and the surface So of the plane z = 0 capping S, is zero; i.e.

I

DH(Q P)

dS

= -

DH, P)

N

3N P j

Then we can write the wave integral in

(8)

in the form

m(p) DH(P Q) DN dS =

HP)

91:11P-2--C-11 - M(Q) 3Np -Ts KIN P) c=0 M(Q) dSP 0 From

(6)

we have

,... _{27ek(z +} + ik[(x - E) cos 6 + (y - n) sin 0] KI(P Q) 1

f

k

-k0

sec20

3NQ

₀₀

[n + i(2, cos 6+ m sin 8)] dedk

+ 41(6 f7/2sec36 . eko(z + sec26

2, cos [ko(x - ) sec e]

cos [ko(y - n) tan 0 sec

e]

k2 (19) dS (18) c=0 -J _P S

(8)

-

1

-mQ tan 0 sin [ko(x - ) sec 8] sin

[ko(y -

1-1) tan 6 sec 0]

+ nQ sec 6 sin [ko(x -) sec 8] cos

[ko(y

-

n) tan 6 sec e]} de (2o)

By comparison of the double integral in (20) with that in the relation

-2 ai -T [(z +

02

(x - (Y - n)2]2 =

1 F12:k(z + C) + ik[(x - cos 6 + (y - n) sin 6]

0 0

k[n + i(2, cos 6 + n sin 8)] dkde

the integrands of which are asymptotically equal for very large values of k, we see that the former integral is singular only at the free surface

z = = 0 when P coincides with Q. As in the previous section, this sin-gularity is removed from the first integral on the right of (19) by

subtracting the transpose from the original integrand, although the problem of treating the last integral of (19), which is also singular, remains.

Next let us consider the second integral of (20). With the

substitution X = tan 6, the integral becomes

vl 2,2 eko(Z + C)(1 + A2){

2, COS

[k0

(X - C)V1

+ A2] COS [k (y - n)A1/1 + x2] Jo

0

- mQ A sin [k0 (x - + Az=] sin [k0 (v - 11)AV1 + A2

+ flQ

Ii

_{+ X2 sin [ko(x -} _{+ A2] cos}

[ko(y -

)All + A2 dA

Convergence problems arise only when

z = c =

0. For this case let us

expand the integrand of (21) in powers of 1/A. We have, with

a = k0 (x - E), =

ko(y

-cos(all + A2) = cos [a(A. + - ".)]

1 1

= cos aA cos -

T7-

"*)]

-

sin cLA sin [a _{- "')]}

(21) +

(9)

-Hence, for very large values of A, we have

0,2

All + A2 sin [ail + AL] sin [13A/1 + A2] z

A2

- -4--+ 4) sin aA sin

[(A2 + 4)]

aX

+

cos ca sin [(3(A2 + i")] - -8-

sin ca, cos [(3(A2 + 4)]

(23h)

2

(1 + A2) sin [ail + A2] cos [ISA/1 + A2]

z

(A2 -

1) cos aX sin [(A2 + 4)]

+

22 -cos aA cos [13(A2 + 4)] + -sin aX sin [(A2 + 4)]

(23c)

2 8

Considering sin [3(A2 + 4)] and cos [13(A2 + 4)] in the forms

sin [130,2 +

= sin (f3A2) cos

g-

co. (fix2) sin

cos [13(A2 4. 4)] = cos ((32) cos

- sin (13A2) sin

S-3--2 2

it is seen that, with z =

= 0, the asymptotic form of the integrand of

(21) is a linear combination of the terms

sin aX

_cos

sin (.2),

8 sin

A cos aX

_(A2),

cos

A2 sin aX scio:

(12A2)

(24)

But from the table of definite integrals by Grbbner and Hofreiter [5], we

have

J sin

sin

W )

2,

dA

- 0,

0

cos

0

Since the derivatives of the first member of (24) with respect to a and

R yield the second and third members, one is tempted to conclude from these

(25)

,2,

_a

(22a)

(1

cos (a1/1 + A2)

- 10,2) cos °LA -

sin ccA

Similarly

,2

_a

(22b)

(22c)

(1

) +

-2--sin (a/1 + A2)

z

- 4A2

sin aX

cos aA

cos (IsAV1 + x2)

. cos [s(A2 +

+

42-

sin [s(x2 + 35)]

sin (ISA/1 + A2)

sin D(A2 + 4)] +

6_

cos [(A2 + 4)]

(22d)

The the terms of the integrand of (21) become asymptotically

V1 + A2 cos [a/1 + A2] cos [13A/1 + A2]

A cos aX cos [( A2

a

+ 4)]

.

- Tsui aA cos [3(A2 + 4)]

(23a)

4)]

(10)

-9-derivations of (25) that the infinite integrals of the second and third

members are also zero. However, since the integrals resulting from these differentiations are not uniformly convergent, this conclusion may not be valid, as is also shown by the integration by parts

(L L L

j sin aX sin (13X2) dX = - -1-cos aL sin ISL2 +=-1= I X cos aX cos (13X2) dX

a a

0 0 0

(26)

Although the limit of the left member as L co is zero, the oscillation

of the first term on the right between ± 1/a indicates that the last integral

is indeterminate. The mean value of the last integral would, however, be

zero in the limit.

If and c were not zero, the asymptotic forms occurring in

(24)

would have been multiplied by ek0(z + c)X2. With this factor the integral

of the first member of (2)-) and its derivatives would be uniformly

conver-gent, and consequently the derivative of the integral would be equal to the

integral of the derivative. For example we would have

B I k_.(z +

dA2

.

kn(z + dx2

(I3X2) dX u

Ba e u sin aX sin = e X

with a determinate value for the right member, no matter how close

z + c

is to zero. Again integrating by parts, we have

fLeko(z + C)X2 .

sin aX sin

(W)

dx .

L

ek0(z 4-

dr-'2

_{cos aL sin}

_a,2

a

0

2

+ _{[k0(z +} c) sin (13X2) + (3 cos (x2)] Xek0(z + C)X2cos aXdX (27)

0

Now the oscillation of the first member on the right with increasing values of L is damped by the exponential factor and yields in the limit

k,(z

+ c)X2

e u sin aX sin (ISX2) dX = 2 [k0 (z + C) sin (13X2) + (3 cos

($2,2)]

a

0

Xek0(z + _)X2cos _aXdX

(28)

which, when z + C is very small, becomes

0 cos aX sin dX 0

f

-W2)

(11)

k0(z +

)A2_{sinaX sin} _(f.ix2) dx 213

ek0(z +

C)X2X cos aX cos (132,2) dA 0

a

0 ₍₂₉₎

Hence, since by (25) the limit of the integral on the left of (29), as

z + c

0, is zero, we see that

lim

I ek()z

C)X2X cos aX

cos

sin

(13X2)

z + c

0 0

= M.V.O.

f

X cos aX

sin

(0,2) dX = 0

cos

0

where M.V.O. denotes the "mean value of the oscillation" for large values of

L. Clearly the foregoing result applies to either sin (3X2) or cos (3A2),

as is indicated in (30). Similarly, by integrating the last integrals in (26) and (27) again by parts, we can show that

lim

f'

ek(z

ox2

z = 0 0

-10-X2 sin aX

sin

(13X2) dA cos L

J

X2

= M.V.O.

sin aX

sin

(13X2) dX = 0

cos

In the above analysis it was assumed in (25) that 13 = k0(Y - O.

If = 0,

a

0, the terms of (24) become

sin aX, X cos aX, A2 sin aX and we can show by integration by parts that

lim

_{rek0(z +}

lim

c)X2sin aXdX = cos aXdX

z + c

-> 0 z +

c ,

0

j' Xek0(

z + c)X2

0 0

z

0

f

x2ek0(z + Ox2 sin aXdX = 0 urn (32)

If a

= 0 also, then we have

ek0(z +

C)X2XdA 1

2k0

1z + CI

which indicates that the second of the limits in (32) does not exist.

( 30) (31) + + -+

(12)

Our conclusion is that the integral in (21) is determinate _except when P and Q coincide and are at the free surface. _{When P and Q are at the} free surface, but not coincident, the integral must be determined as the limiting value as z + C -> 0 through negative values. Finally, the sin-gularity when P and Q coincide is removed from the first integral on the

right of (19) by subtracting the transpose from the original integrand, although, as for the double integral of (20), the last integral of (19)

remains to be treated.

Let us now consider the last integral in (19),

I aH(Q, P)

ac

So

[27ekz - ik[(x - E) cos e + (y - n) sin e]

1 dS = P 7 J k - 1(0 sec28 c=0 S00 0 k2dOdkdS

7/2

-4k2

j Isec4O

ekozsec28sin

[ko(x E) sec 8] cos

[ko(y - 71)

tan 8 sec 8]

0

S0

0

dedSP ( 33)

Let

n =

n(x)

_{be the equation of the hull waterplane at} _c _{= 0.} _Take

the origin at the midship section and let St denote the half length of the

ship. Then, interchanging the order of integration in

₍₃₃₎

_{with the}

inte-gration over

S0 taken first, we are led to consider the integral,

F(k, 0) =

11(x)eik(

cos 8

+ n

sin

_0)dndE

-2-41(x)

This becomes

F(k, 8) = 2 csc e fzeik E cos _{sin [k i} _{(E) sin 8] dE}

-Z

or, integrating by parts and noting that

n(z) =

n(-0 = 0, we obtain

2i sec 8 _{cos 0}

F(k, 6) =

_q'CE)

_{cos [lc n}

_{CE) sin 6]}

₍₃₆₎

where n'(E) denotes the derivative of Ti with respect to E.

(3/4)

(35)

(13)

-12-Along the parallel middle body of a ship form, n'(E) = 0,

and near the bow and stern, n(E) is very small. This suggests the

approx-imation

e) 2i sec 9

2eik E _{cos en,(}

(37)

F(k, E) [1 - 4k2n(02 sin2e] dE

Additional terms may be taken in the expansion of cos [k n (E) sin 8] if

required for greater accuracy. If the parallel middle body extends over

the range a S E

s

b, (37) may be written

F(k, 6) = 2i sec 6 {fa .4. r}eik E cos e

b

[ne(E) - 41c2

n'(E)

n(E)2 sin2e] dE

or, introducing p = E a in the first integral, and p = b in the

+ a - b

second,

F(k, 0), 2i sec 6eika cos Of eik(2. +

a)p cos

--1

[fl'(E) -

4k2n'(E)

n(E)2 sin26] dp

1 ,

+ eikb cos 6

f

e

ikj -

b )p cos 0[n'() - 4k2 nt(E) n(E)2 sin20] dp (38) 0

One can now fit either Fourier series or polynomials in p to the functions n'(E) and n'(E) n(E)2, the choice depending upon the particular form. If

n'(E) becomes infinite at the bow and stern, as it will if the radii of curvature are not zero at the extremities, a suitable fit which can satisfy this condition can be obtained from a polynomial for n(E)2 of the form

n(E)2 = (1 - 1_12) p(p) (39)

where DM is a polynomial such that p(±1) O. For then

n'(E) =

[4 (1

112) p'(u)

Pp(u)]

cl4j/dE

n(E)

is seen to become infinite at p = ± 1. In this way F(k, 6) can be ex-pressed as a series of functions of k and 0, each of which is regular even

as z approaches zero.

e

0

(14)

-

-13-Although not evident from the form of the last integral in (33), by returning to the complex exponential form from which its integrand was

derived one finds that the function F(k, 0) _{of (34), with k = ko sec20,} applies to this integral as well. _{This will not be developed here in}

detail, nor will the analysis of the wave kernel be carried any further,

since the application to a particular _{case, on which, it has been seen,}

the nature of the subsequent analysis would depend, has not yet been per-formed.

Convergence of Iteration Formulas

Hereafter we shall consider only the _{case where the boundary}

condition on the plane z = 0 is

a(1) n

az z = 0

i.e., the case of "zero" Froude number. _{The integral equation}

₍₈₎

then reduces to M(Q) = F(Q) + M(P) K(P, Q) dSp (41) S+S' where 1 a K(P, Q) = ,Tr , F(Q) = - 2rriQ DNQ R)

Equation (9) shows that the homogeneous integral equation

f(Q) = A j _{K(Q, P) f(P) dSp}

S+S'

has the eigenfunction f(P) = 1 when A = - 1. _{Thus A = - 1 is an eigenvalue}

of the kernel K(Q. P), and hence _{also of its transpose K(P, Q).} Consider the inhomogeneous integral _equation

M(Q) = F(Q) + A

f

M(P) K(P, Q) dSp

()43)

S+S'

which reduces to (41) when A = 1. _{The theory of this integral}

equation states that M(Q), considered as _{a function of the complex variable}

A, is regular in the unit circle about A _{= 0, and has a simple pole at A}

= - 1. _Writing

(42)

(4o)

(15)

-M(Q) = F(Q) + XF1(Q) + A2F2(Q) + IXI < 1

and substituting

(44)

into (43), yields the relation

F1(Q) =

Fn(P) K(P, Q) dSP S+S' Put M(Q) = F(Q) + AF (Q) + _{+ XriFn(Q)} 1

Then, by

()45),

we obtain the iteration formula

M0+1(Q) = F(Q) + A Mn(P) K(P, Q) _dSP

S+S'

According to (44), however, the sequence of functions M1.1(Q) defined by

(47)

may not converge when A = 1.

We can eliminate the pole at X = - 1 by considering

(1 + A) M(Q) = F(Q) + X(F + F1) + X2(F, + F2) + "' , 1X1 <

X21 (48)

where A2 denotes the next eigenvalue of K(P, Q), arranged in the order of

Comparison with

(47)

shows that the sequences Mn and Mr'l are obtained from

the identical iteration formula, but these sequences differ because of the

change in the initial approximations,

M0(Q) = F(Q),

M(Q) =

V2)2,

(51)

Thus, when X = 1, M(Q) = F(Q). Alternatively, if we observe from

(44)

and

()9)

that 1 m _ 1 + A (mn + XMh -1) we obtain when X = 1 Me = _{(Mn + Mh-1)}

(52)

increasing absolute magnitude. Defining

1 [F(Q) + X(F + F1 ) + +

x(1+ Fn)]

(49)

M'

=1

n + then, then, by (45), Me(Q) = F(Q) + X I Mr'l(P) K(P, Q) dSp

(50)

S+S'

-(45)

(4 6

(47)

(16)

For if (54) is valid, then by (53)

then

-15--i.e.,

the arithmetic means of successive pairs of members of the sequence

M(Q) form a convergent sequence.

Let us consider the modification of the iteration formula (50),

2M"1 (Q) = Mu(Q) + F(Q)_n+ + A j M"(Q) K(P, Q) dS

S+S'

with M" = M given by (51). We have, by (50),

0 0

M" = .3i [M'

+F+X.IM'

KdS] =i(M' + M')

1 0 0 0 1

M" = 41₂ Pi (M1 + 141)

+F+Xfi(M'

+ M') KdS] = i- (M' + 2M' + MT)

0 1 ₀

1 0 1 2

We can now show by mathematical induction that

jmt]

n n + M'] + 2nF + X [M' + M' + + IV] KdS 0 1 1

ra

n n 2nF E (1 + 1)nF = 1 2 n Then, by (50),

2nF + X I [M' + (1

+ _{+ M']KdS = M' +} 1JM' + (2iM' + + m' 0 1 2 3 n+1 Thus we have 1 n' [Vim, rinl

if

n 111)] ml 1 2M" = (

{M'

+ [

p

1 Loj n+1 2n 0 0 2j 2 n n+lf or since

Pa/1-1-1

) J1+1, n+1 1 ,, [n+1 n+1 n+1 Lmo +

imi

r+21],\I 2 /1+1

as we wished to show to complete the proof by induction.

(53) (54) 2M" = n+1 2n [M' +0 tln M1 + But = n _1

[Iv

2n 0 n' (_{1M y} 1 4. 2

imt

2 4_ =

(17)

In order to investigate the convergence of the sequence M"

let us take N sufficiently large so that, for r > N, IMI:

-

MI < e/2,

where M is the limit of the sequence {M'}.. From

(54)

we have

M" - M =

[(M' -

M)

+n](M'

- M) + "'

+Nj(N

-

M) + +

j(M'

-

M)] 2n 0 1 1 n n Then But 1M" -1 2n [[N+1] or (n) + +

_-Ni]

[1 + + [fl + = 1 2 n

-16-Then, if p is an upper bound of -

MI,

i = 0, 1, 2, N, and n << 2N, we have

nrn

IM" -

ml < -P-- [1 + ,l] + {n2j + (NI ] + i- < -P-217 (N + 1) [Nil + n 2n p(N + 1) nN e < + E N! 2n 2

by taking n sufficiently large. Hence the sequence

{M}

converges to

M.

The alternative iteration formulas for M' or M" arise when relation

(10) is applied to remove the singularity of the kernel in the integral equation (41). We obtain

M(Q) = P(Q) - M(Q) + f [M(P) K(P,

Q) - M(Q)

K(Q, P)] aSp (55)

S+S'

and the iteration formulas

Mn+1 = F - Mn + I [Mil(p) K(F, Q) - M(Q) K(Q, p)] dSP

(56)

< [ 2n n I N+2j Mo' -r + + N+1 n n, rn

lj

] <

-1 ' n N+2 2n 2n+1

+ ."

-also

(18)

-2Mn+1 = F + I _[Mn(P) K(P, q) - M(Q) K(Q, P)] dS (57)

By (9), the first of these is seen to be of the form (50), the second of

the form (53). Both begin with the same initial approximation which, by

(51), is given by

M0(Q) = 4 F(Q)

= -

ii+TT ZQ (58)

Although there is no a-priori basis for preferring one iteration formula over the other, comparison of the numerical results with the known exact solution for the case of an ellipsoid has shown that the sequence given by

(57) converged much more rapidly than that obtained from (56).

At points where the normal to S is not continuous, the integral equation (8), and of course the above iteration formulae, are not valid. At such points we can either set M(Q) = 0, as can be justified, or round sharp edges with small, nonzero curvature and continue to use the

iter-ation formula (57).

Distribution of Velocity Potential on S

Once the source distribution M has been found, the velocity potential cl) can be computed from (7). For points Q on S we again

en-counter a singularity when P coincides with Q. This singularity may be

removed as follows:

Let N(P) be a source distribution on S + S' which makes the surface an equipotential of potential This distribution satisfies the homogeneous integral equation

IN(P) K(P, Q) dSp = - N(Q) (59)

S+S'

with the same kernel as in (42). In fact, N(P) is the eigenfunction of K(P, Q) associated with the eigenvalue A = - 1. This equation can be solved by means of the iteration formula

N1(Q) =

- j Nn(P) K(P, Q)

dSP

(19)

-18-which, by applying (9), may be written in the singularity-free form

Nn+1(Q) = Nn(Q) - [Nn(P) K(P, Q) - Nn(Q) K(Q, P)] dSP

(60)

S+S

Since the matrices occurring in (60) have already been obtained for the numerical evaluation of M(Q) from (57), the corresponding values of N(Q) can be obtained from (60) with little additional computer time. Since

the potential is constant in the interior of an equipotential surface, its value may conveniently be computed at the origin as

cpo =

f

[x2 y2 4.

Z2]li

dSp

N(p)

(61)

S+S'

We can now apply the solution of this Dirichlet problem to eliminate the singularity from the expression for the velocity potential

(7), by writing

cD(Q)

= ux -

[m(p) -

N(P) ds + 6

P

'0 Nil+

S+S'

+ IM(P) H(P, Q) dSp (62)

Here also we can justify setting the first integrand of (62) equal to zero when P coincides with Q, by the same argument as was used in equation (10).

Application to a Double Ship Model -- Zero Froude Number

Since the x-y and x-z planes are planes of symmetry, it is

necessary to determine the source distribution over only one-fourth of the

hull surface of the double model. Let us consider only points Q for y, z, < 0.

Denote by Sl, S2,

S3, S4

the parts of S + S' for which y, z > 0; y < 0,

z > 0; y < 0, z < 0; y > 0, z < 0, respectively. Put R =

[(x

- )2 4- (Y '1)2 (z c)2];5 1 R

= [(x

E)2 (y n)2 (z

02]4

2

53 = [(x -

E)2 (Y

n)2 + (z + d2]2

R4 = [(x -

E)2

(y - r1)2 + (z + 02]2

- + 1 - --+

(20)

-19-the distances from

Q(x, y, z)

to the congruent points P E P1(

n,

P2(E5 -n,

), p3(

-n,

-d,

P4(

- ).

At congruent points we have

M(P) =

M(P2) = M(P3) = M(P4)

and, denoting the direction cosines at P.

by Z., m., n., i = 1, 2,

3, 4,

we obtain the following relations:

kp E2.1 =2.2

=2.3

=ft,

4

mp = m1= - m2 = -

_{m3 =}

np

E

n1 =n2

=

-n3

= -

n4

If the values of the integrand of

(57)

at congruent points P are collected, the resulting integral would extend only over

S1 and is found to be of the form 1

mn+1(Q) =

P(Q)Tr-f

_LMn(P)

J(P, Q)

-

Mn(Q) J(Q, P)] ds

(63)

Si where

j(p,

Q) =

_[(E

_ _x)st,Q + _(n

_{- y)mQ +}

-

_z)nQ]

mn

mQn+nQ

c

nc

QQ

-

2

(-TT

+ R _P R4 2 3

and, since the integrand of

(63)

vanishes when P coincides with _Q, we may set

J(Q, Q)

= 0

₍₆₅₎

Similarly the integrations over S + S' in connection with the Dirichlet problem in (60),

(61)

and (62) can be expressed in terms of

integrals over Sl. Thus

(60)

becomes

N1(Q) = Nn(Q)

- [Nn(p) J(p, _ Nn(Q) _{P)] dS}

₍₆₆₎

Si

in which J(P, Q) and J(Q, Q) are given in

₍₆₄₎

and (65),

(61)

becomes

o

= - 4

f

_[x

171)) ,

as

2 + y2 +

P

(67)

1 + + + 1 1 1 l 1 2 3

(64)

2

(21)

-and (62), without the wave integral, assumes the form 4)(Q) = Ux - L(P, Q) [M(P) N(P) m(Q)] dS + cp M(Q)

(68)

N(Q) P 0 N(Q) Si where 1 1 1 1 L(P, Q) = FT-- + 17- + + 1 2 3 4 and L(Q, Q) = o

(69)

Application to an Ellipsoid

With the equation of the ellipsoid in the alternative forms

x2 v2 z2 + = a or -20-b

y = - x2 cos e, z =

/a2

-

x2 sin e

a

the direction cosines are

=

LILL

m = AialL z Q a-, 5 Q b2

'

n

Q=

c2 where A(Q) = z2-a b' Also we have

(C - x)2,Q + (7 - y)mai +

-

z)nQ = - A(Q)[(: -2

x)2

-1y)2

-z)11

_c2j

(71)

The right member of (71) is preferable to the first for numerical computations, especially when P is near Q, since all terms on the right are then of second order of smallness, and a loss of numerical accuracy would be expected if the same result were obtained from the sum of the first-order terms on the left.

In terms of x and e as the independent variables, the element of area dSp becomes

[b2c2x2

x2 ,) 2 dS = + (1 - + c2cos20)] dedx a a

4

+ - - + + Q

(22)

with 0 varying from 0 to 27 and x. from -ato a.

Exact solutions of the foregoing, integral equations for the

ellipsoid are expressible in terms of the Lame ellipsoidal coordinates (10, it,v), which are related to. the rectangular coordinates (x, y, z) by [6]

puv 2 (p2 _k2)(k2 112)(k2 v2)

X =

hk Y = k2(k2 - h2)

(72)

When p a, k = - b2,

/2_

_{C2, (72) is equivalent to the equation of} the ellipsoid

(70,

Here

<, < h2 112 1k2 _{< h2}

Solutions of the present problem can be -expressed. in terms of the

Lame functions of the first and SecOnd kind, En and Finn, _{For the source dis=}

tribution we find.

3 v R[(Q)

47ahk Pr-_-(aY (a2 - p2)(a2 - v2 47 L- = (73) and fbr the velocity potential,

Fl(a) [1 - -1 I Ux

F (a)

1

Here Fl denotes the derivative of Fl with respect to its argument, and

1 1 do 3a = 3a

J2

17(02 h2)(02 k - LFlqb, AY -

E(0, A)]

(75) a h -= arccos A = -a k

where F(0, A), E(0, A), are, the Legendre incomplete _{elliptic integrals Of} modulus A and amplitude 0 of first and second kinds..

EXact sOlutions of the equipotential problem are also of interest. If the total strength of the distribution over the ellipsoid _{it G, the source,} strength N(Q) is given

by

(n)

GA(Q)

N""

F1( 1 2 (1)2 h2)(112

h2)(h2

_v2) h (k2 - h2) (74) (76) -21-z = i/a2 -

h

= k2 p2 Fl(a) = ₌ =

--

(23)

-Then

-22-On and within the ellipsoid the potential has the constant value

ds G

,(A

,\

(I) =

0 2 I [(a2 s)(b, s)(c2 = T-[ r'P' A'

(77)

0

Calculations have been performed for the case a = 1, b = 0.25, c = 0.50, for which Fl(a) 1 -Fl(a) = 1.12659 1 ,

u

M(Q) = - 1.12059

17

xA(Q) (I) = 1.12659 Ux N(Q) = 1-GA(Q) o = 1.76984 G FORTRAN Program

In the FORTRAN program given in the Appendix, the Gauss 16-point quadrature formula was selected for purposes of illustration. Because of symmetry about the x-z and x-y planes, with this quadrature formula there are 16

i's but only 8 ps, so that the matricesand

_{jij, kl} L. ki contain

128 x 128 =

16384

elements. An advantage of the Gauss formula is that it gives finer intervals at the bow and stern than amidships, as is required

by the rapid changes in form of the transverse sections at a ship's

extremi-ties. Nevertheless, for a particular form, it may be desired to use a

quadrature formula of the Simpson-rule type, with a fine mesh near the bow

and stern and a coarser one over the remainder of the hull.

The input data of the FORTRAN program are the coordinates

(x., y.., z..), the direction cosines 1.., m.. n.., and the weighting

fact-13 ij 13 13, 13

ors of the selected quadrature formula, Ak and Al. The various do-loops in

the program perform the following operations:

Loop 2 computes the first approximations for the iteration formulas

(63)

and (66), (m0)ij = SDI (I, J) and (110)ij = SDO (I, J). Here is also

computedtheproductoftheareaelementE.j

=E (I, J) by the weighting

i

factors, A.A.E.. = F (I, J).

j lj

+

(24)

-23-Loop 9 computes the kernel matrix of the integral equations (63) and (66), Jii, la which, multiplied by the weighting factors, is denoted by

C (IJ, KL).

The program yields simultaneously the transpose of the matrix

with its corresponding weighting factors, and then computes directly the

diagonal terms of the matrix. Also computed in this loop is the matrix

Lij

,

kl

= P1 (Li, KL),

including the weighting factors and the element of

area,usedincomputingqp.from (68).

Loops

22-23 compute successive approximations to the source

distribution Mn = SDI from the iteration formula (63), for 10 iterations.

Loops

25-26 compute successive approximations to the source distribution for the equipotential problem, Nn = SDO, from the iteration

formula (66), for 20 iterations.

Loop 28 computes (1)0 = PHO from (67).

Loop 35 computes cl)ij = PHI (I, J) from the expression for the velocity potential in (68).

Results for an Ellipsoid

The program was tested by applying it to an ellipsoid with axes

a, b, c in the proportions

a: b

: c = 4 : 1 : 2

Input data are given in Table 1, in which 8 is the parametric angle of the

equation of the ellipsoid, (70). The values of x and 0 in the table corres-pond to the abscissas of the Gauss 16-point quadrature formula.

Values of the source distribution, M(Q), computed from the exact expression as well as from the discretized iteration formula with 10

iter-ations, are given in Table 2. The results are seen to agree to within four significant figures for 0 near zero, but to only three significant figures

for 6 near

7/2.

Repetition of the calculations employing double-precision arithmetic showed essentially the same results, indicating that the errors are due principally to the discretization of the equations, rather than

round-off errors by the computer.

Table 3 gives the results for the source distribution N(Q) when

(25)

-21-strength on the hull was not normalized in the iterations, an unknown constant factor is present in the values computed from the iteration formula. The ratios of the exact values to those computed with twenty iterations show

con-sistency to five significant figures. The exact value of the constant poten-tial on the hull,

cA0E = -1.76977,

and the value of the potential at the

ori-gin computed from the source distribution of the twentieth iteration, cp_OC

4.07197,

are in essentially the same ratio as that given in the Table,

qb

OE OC

= 0.434621.

Finally, the distribution of the velocity potential over the hull

is given in Table

4.

The exact values are seen to agree with those computed from the previously calculated source distributions, M(Q) and N(Q), to within at least five significant figures. This is remarkable in view of the

afore-mentioned result that the values of the source distribution M(Q) from the

tenth iteration have considerably larger errors. A reexamination of the values of M(Q) in Table 2 shows that the computed values tend to be larger

than the exact ones at small values of 6, and smaller than the exact ones at the large values of 6, a possible explanation for the unexpectedly good

agreement in the values of the velocity potential.

(26)

References

[i] J. L. Hess and A. M. 0. Smith, "Calculation of Potential Flow about

Arbitrary Bodies", Progress in Aeronautical Sciences, Vol. 8,

Per-gamon Press, New York,

1966.

[2]

0. D. Kellogg, Foundations of Potential Theory, Frederick Ungar Publishing:Company, New York,

1929.

[3 ] F. H. Todd, "Some Further Experiments on Single-Screw Merchant Ship

Forms Series

60",

Transactions

of

the Society

of

Naval Architecture

and Marine Engineering, Vol.

61, 1953.

[4]

J. V. Wehausen, "Surface Waves", Encyclopedia

of

_{Physics, edited by} S. Flugge, Vol. IX, Fluid Dynamics III, Springer Verlag, Berlin,

_1960.

[ 5 ] W. GrtSbner and N. Hofreiter, Integraltafel, Zweiter TeiZ, Bestimmte

Integrale, Springer-Verlag, Wien, New York,

1966.

E. W. Hobson, The Theory

of

Spherical and Ellipsoidal Harmonics,

Chelsea Publishing Company, New York,

_1955.

C. von Kerczek and E. 0. Tuck, "The Representation of Ship Hulls by Conformal Mapping Functions", Journal

of

_{Ship Research, vol.}

_13,

no.

4,

December

1969.

(27)

--

-26-TABLE I

INPUT DATA - COORDINATES AND DIRECTION COSINES

OF A 14 : 1 x = : 2 ELLIPSOID

0755404

1

0.149245

0.161994

0.04E716

0.279081

0.957566

0.071991

0.442341

0.148048

0.140244

0.296377

0.929355

0.220094

0.710450

0.123216

0.215898

0.331168

0.864286

0.378599

0.970557

0.092529

0.270360

0.38095E

0.746614

0.545380

1.1E6584

0.061402

0.303743

0.435264

0.566082

0.700069

1.359729

0.034320

0.320359

0.476405

0.346307

0.808152

1.483732

0.014244

0.325388

0.495196

0.149402

0.855E39

1.5541146

0.002727

0.327584

0.499237

0.028841

0.265985

x = 0.265631

1

0.149245

0.123779

0.0372214

0.399553

0.914131

0.058726

0.442341

0./13/23

0./07/60

0.421931

0.882226

0.208931

0.719450

0.094149

0.164962

0.405/74

0.810548

0.355059

0.070557

0.070701

0.206581

0.525650

0.5E6941

0.501791

1.186584

0.040917

0.232090

0.586995

0.509045

0.629532

1.359729

0.026224

0.244785

0.630642

0.305678

0.713338

1.483732

0.010884

0.249393

0.649709

0.130712

0.748808

1.554146

0.002084

0.250306

0.653E18

0.0251E6

0.756233

x = 0.944575

1

0.149245

0.081162

0.024407

0.587142

0.807201

0.060627

0.442341

0.074174

0.070255

0.012317

0.769333

0.182195

0.719450

0.061733

0.108'109

0.652909

0.6E9016

0.301822

0.070557

0.946359

0.135455

0.716912

0.562963

0.411229

1.1865P4

0.030764

0.152120

0.769949

0.401221

0.496186

1.359723

0.017195

0.160505

0.804046

0.2341E6

0.546504

1.4E3732

0.007137

0.163526

0.81E094

0.098390

0.566518

1.554145

0.001367

0.164125

0.221000

0.019004

0.570613

x = 0.0E9401

1

0.149245

0.035899

0.010796

0.864215

0.501707

0.037719

0.442341

0.032898

0.031079

0.877978

0.465217

0.110316

0.719450

0.027306

0.047E45

0.009792

0.397762

0.174239

0.9/0557

0.020505

0.059914

0.925044

0.306740

0.224055

1.186584

0.013607

0.067311

0.943850

0.207696

0.256855

1.359/29

0.007605

0.070993

0.954540

0.117401

0.273970

1.483732

0.003157

0.072330

0.05E656

0.048936

0.280329

1.554146

0.000604

0.012595

0.959480

0.009378

0.281599

(28)

-27-TABLE I, Continued

x = 0.095012

1

0.149245

C.246102

0.074009

0.024054

0.996897

0.074948

0.442341

0.22491u

0.213060

0.025683

0.972764

0.230372

9.719450

0./27/91

0.327995

0.'129045

0.915587

0.401071

0.970557

0.140571

0,410734

0.034092

0.807036

0.589518

1.166584

0.093283

0.461450

0.039995

0.628266

0.776970

1.359729

0.052139

0.4E6692

0.044815

0.3934E2

0.918240

1.483732

0.021640

0.495853

0.047138

0.171776

0.9E4006

1.554146

0.004144

0.497669

0.047648

0.033248

0.993311

x = 0.281604

1

0.149245

0.237216

0.071337

0.073764

0.994468

0.0/4766

0.442341

0.216704

0.205367

0.078753

0.970062

0.229732

0.7191459

0.120432

0.316152

0.088994

0.912339

0.399648

0.070557

0.135495

0.395903

0.104319

0.803100

0.586642

1.1E6584

0.089915

0.444782

0./22/55

0.624000

0.771763

1.359729

0.050256

0.469118

0.136646

0.390133

0.910542

1.4E3732

0.920859

0.477948

0.143600

0.170125

0.974893

1.554146

0.003994

0.479699

0.145126

0.032935

0.988865

x = 0.458017

1

0.140245

0.219765

0.060069

0.128E03

0.988279

0.0/4345

0.442341

0.200246

0.190259

0.137370

0.903859

0.228263

0.719450

0.167159

0.292604

0.154966

0.904908

0.396393

0.970357

0.125528

0.366/78

0.181103

0.794152

0.580107

1.166584

0.023300

0.412007

0.211201

0.614565

0.766052

1.559729

0.040559

0.434002

0.235364

0.382012

0.893341

1.483732

0.019324

0.442788

0.246E02

0.166045

0.054615

1.554146

0.003700

0.1414141110

0.249376

0.032234

0.967870

x = 0.617E76

1

0.149245

0.1943E4

0.058456

0.104329

0.978170

0.073541

0.442341

0.177649

0.100285

0.20u949

0.952619

0.225459

0.719450

0.147853

0.259060

0.232674

0.609234

0.390228

0.970557

0.111030

0.324417

0.270390

0.777425

0.567228

1.186524

0.073680

0.364476

0.312904

0.507176

0.732522

1.359729

0.041162

0.324413

0.346471

0.369421

0.862231

1.4E3732

0.017092

0.391649

0.362157

0.160293

0.918231

1.554146

0.003273

0.393023

0.335660

0.930982

0.930270

(29)

TABLE 2

SOURCE DISTRIBUTION ON ELLIPSOID, M(P)

Exact

-0.002157 -0.006615 -0.011548

-0.017422 -0.025020 -0.035820 -0.052639

-0.077478

-0.002303 -0.007060 -0.012315 -0.018553

-0.026571 -0.037827 -0.054805 -0.072712

-0.002504 -0.007978 -0.013893 -0.020860

-0.029600 -0.041757 -0.059072 -0.060757

-0.003056 -0.000352 -0.016236

-0.024241 -0.034153 -0.047126 -0.064272

-0.0E2931

-0.00358G -0.010951 -0.018935

-0.02600 -0.039022 -0.052625 -0.069027

-9.0'64619

-0.004018 -0.012250 -0.021101 -0.031067

-0.042710 -0.056538 -0.072084 -0.025576

-0.004226 -0.012874 -0.022131 -0.03246F

-0.044395 -0.052253 -0.073343

-0.085945

-C.004272 -0.013011 -0.022357 -0.032774

-0.044757 -0.058615 -0.073004

-0.066019

Computed - 10th Interaction

-0.002157 -0.006615 -0.011548 -0.017422

-0.025021 -0.035822 -0.052640 -0.077497

-0.002303 -0.007061 -0.012317 -0.018555

-0.026573 -0.037830 -0.054902 -0.0/8720

-0.002505 -0.007921 -0.013297 -0.020866

-0.029699 -0.041771 -0.050094 -0.080748

-0.003052 -0.009355 -0.016243 -0.024251

-0.034166 -0.047141 -0.004278 -0.082904

-0.003585 -0.010950 -0.018932 -0.028056

-n.039014 -0.052608 -0.068994 -0.084577

-0.004014 -0.012238 -0.021079 -0.031030

-0.042666 -0.056420 -0.072016 -0.035525

-0.004219 -0.012853 -0.022096 -0.032417

-0.044328 -0.058172 -0.073252 -0.085889

-0.004264 -0.012988 -0.022312 -0.032718

-0.044685 -0.058529 -0.073515 -0.065963

(30)

TABLE 3

EQUIPOTENTIAL SOURCE DISTRIBUTION, N(P)

Exact

0.161174

0.166604

0.179037

0.290225

0.235190

0.293647

0.395723

0.556070

0.172067

0.178038

0.130946

0.213227

0.209772

0.310305

0.412026

0.504926

G.194614

0.201188

0.215395

0.239733

0.279090

0.342549

0.444086

0.579006

0.22E432

0.235833

9.251724

0.278506

0.321054

0.325591

0.483181

0.595210

0.267979

C.270157

0.293559

0.322490

0.366E21

0.431709

0.51E927

0.697323

0.300277

0.308914

0.327145

0.35692

0.401492

0.463800

0.541007

0.610193

0.315639

0.324C3u

0.343125

0.373103

0.417321.

0.477866

.551375

0.610837

0.319261

0.326087

0.346620

0.576656

0.420730

9.4E0244

0.553333

0.617372

Computed - 20th Iteration

n.370837

0.3E3791

0.011939

0.460628

0.541154

0.576101

0.910506

1.279444

0.395946

0.409638

0.439338

0.090604

0.574691

0.713969

0.949536

1.299E19

0.447777

0.062904

0.495592

0.551591

0.542155

0.786157

1.021727

1.333590

0.5255E7

0.542617

0.579120

0.641015

0.732699

0.2E9493

1.111734

1.309499

0.016579

0.635305

0.575437

0.742002

0.844004

0.993282

1.1939'80

1.397360

0. 60'692

0.710765

0.752712

9.P21363

0.923777

1.067139

1.246255

1.413170

0.72669E

6.746936

0.76942C

0.1252549

0.960213

1.009505

1.268039

1.019261

"./34570

0.754E78

0.797521

0.806E30

0.96E050

1.19E355

1.273145

1.420490

Ratio of Exact to Computed Values

6.434622

0.434022

0.034621

0.434621

0.434629

0.434620

6.43461°,

0.434619

0.434622

0.434621

0.h34b20

0.434620

0.434619

C.434619

0.454622

6.434022

0.434621

0.434620

0.434619

0.434622

0.034621

0.434620

0.434020

0.434619

0.434022

0.434622

0.434621

0.430020

0.430620

0.434619

0.034619

0.434622

0.034E21

0.034620

0.034020

0.434619

9.434E22

0.436622

0.434622

0.434621

0.434620

0.434619

0.030619

0.434622

0.430621

0.434020

0.434620

0.450019

0.430619

(31)

TABLE

4

VELOCITY POTENTIAL ON SURFACE OF ELLIPSOID

Exact

0.107040

0.317252

0.515997

0.696093

0.851031

0.975211

1.064149

1.114649

Computed

0.107038

0.317247

0.515989

0.696082

0.851017

0.975195

1.064130

1.114623

0.107038

0.317247

0.515989

0.696082

0.851017

0.975195

1.064129

1.114628

0.107038

0.317246

0.515988

0.696081

0.851016

0.975193

1.064127

1.114636

0.107038

0.317246

0.515988

0.696081

0.851016

0.975195

1.064134

1.114643

0.107039

0.317248

0.515991

0.696086

0.851023

0.975204

1.064146

1.114650

0.107040

0.317251

0.515997

0.696093

0.851033

0.975215

1.064156

1.114654

0.107041

0.317254

0.516001

0.696098

0.851038

0.975221

1.064161

1.114655

0.107041

0.317254

0.516002

0.696100

0.851040

0.975223

1.064162

1.114656

(32)

-31-APPENDIX

FORTRAN IV(d) Program for aBM 360/ Program for Computing Source Distribution and Velocity Potential on &Ship Surface

DVIENSTON X(16),Y(16,8),Z(16,8),Ft(16,8),EM(16,8),EN.(16,8),A(.16),

1CAS(128),SD(128),SIA(128),SD3(128),SDE(16,8),SDO(16,8),SO1(16,8)

2P111(16,8),C(128,128,),P1(128,128),F(16,8),E(16,8)

P = 12.56637061

,REAn(5,.1) (X(I),111=1,16) 1 FORi.IAT(8F10.8)

READ(5,1) (A(1),I=1,16)

REA1(5,1)

(( Y(I,J),1=1,8 ),J=1,E)

READ(5,1)

('( Y(I,J),I.=9,16),J=1,8)

READ(5,1)

(( Z(1,J),1=1,8 ),1=1,E)

READ(5,1) (( Z(I,d),I=9,16),,1=1,8)

READ(5,1) ((EL(1,3),I=1,8 );J=1,8>

REA1(5,1) ((EL(1,J),1=9,16),.1=1,8)

READ(5,1) ((01(1,J0,1=1,8 ),..1=1,8) READ( 5,1) ( ( EM( I, 1=9,1( ),J=1,8 )

REA3(5,1) ((EN(1,J),1=1,8 ),J=1,8)

READ(5,1) C(EN:(1,J),1=9,16),J=1,F)

DO 2 I' = 1,16

X2 = X(1) * X(tY

DO 2

J = 1,8 31 = J + 8

Y2 = Y(1,J)

Y(II,J)

Z2 = Z(I,J)

,J) E(1,,J)

F(1,.1') = E(11,J)

* A(1)

A(J1)

SD1(1,1.1) - EL(II,J)/ P

2 SrY0 ( I ,Li) SCIRT(X2 + Y2. + Z2),

KL = 0 DO 9' K '1,16 50 9 L KL K + 1 S93(KL)

SDOCK,L)

DAS(KL) = SOl(K,L)

SD(KL) F,011(K.,1.)

lu

= 0 St1t.1 = 0..0 DO 8 11-4 1,16 & J

it

= IJ + 1

= X(K) - X(I)

.Y5 = Y(K,1) - Y(1,J)

Yi Y(K,L)

Y(I,J)

ZO = Z(K,L)

Z(I,J)

Z1 = Z(K,L)

XO2 F X0 *

Y02 = YO * YO

Y12 = Y1 * Y1

ZO2 = ZO

ZO Z12 Z.1 * Li * * = * = -= = = 1,8 = = = 1,8 -= -+ Z(I,J) * =

(33)

APPENDIX, Cotitinued

SQRT(X02 + YO2 # Z02)

R22 = SC,RT(X02 + Y12 4 Z02)

R33 = SQRT(X02 + Y12 + Z12)

R44 = SCLRT(X02 + 11.02 + Z12)

R1 = R11 * R11 * R11

R2 = R22 * R22 * R22

43 = R33 * R33 * R33

R44 * R44 * R44

V2 = (1.01R2) + (1.0/R3)

a

V3 = (1.0/R4) + (1.0/R3)

W2 = EM(I,J) * Y(K,L) *

v2

W3 = EN(I,J) * Z(K,L) * V3

PIJ = X0 * EL(I,J) + YO * EM (I'

+ ZO *

1F(K.EQ.I.AND.L.EQ.J) GO TO 10

V11 = (1.0/R11) +(1.0/R22) +(10/433) +(1.0/R44)

V1 = (1.0/R1) + (1.0/R2) +

(

1.0/R3) + (1.0 / R4)

P1(IJ,KL) = V11

BA = PIJ

* V1 - 2.0 *

(W2 + W3)

BE = BA * F(1,J)

ro

TO 11

10 BA = 0.00000000

BE = 0.00000000

P1(IJ,KL)' = BE

Ii

0(IJ,KL) = RA * F(K,L)

8 SUM = SUM + BE

9 C(KIL,KL)

C(KL,KL1 - SUM

Ll = 1

13 DO 22

IJ 1,128

SUM = 0.0

')0 21 KL = 1,128

21 SUM = SUM + sn(KL) * C(IJ,KL)

22 Sn%(IJ) = 0.125 * SUfl + CAS(J)

IJ = 0

no

23 1

- 1,16

00 23

J, = 1,E IJ = IJ + 1 sn(IJ)

= soAm)

23 SlE(I,J) = SOA(IJ)

LA = LA +

1

1F(L\.LE.10) GO TO 13

WRITE(6,102)((SDE(I,J),1=1,8),J=1,8)

102 FOR1\T(408 SOJRCE DISTRHUTION AFTER 10 IOTERATI0NS/C8F10'.,6))

WRITE(6,110)((SDE(1,J),1=9,16),J=1,8)

110 FORMT(1H///((2F10.6))

LA = 1 27 DO 25 IJ = 1,128,

SUM = 0.0

DO 24- KL F 1,128

24 SUM = SUM + S03(KL)

C(TJ,KL)

25 S0(1J) = S03(4J) - 0.25 * SUM

1J = 0

DO 26 ii= 1,16

DO 25 J = 1,8 IJ = TJ + 1

S93(IJ) = SO(1J)

-32-R11 =

EN(I,J) = -= *

(34)

-

-33-APPENDIX, Continued

26 SDO(I,J) = Sr)(IJ)

=

LA +

1

IF(LA.L:.20) CO TO 27

URITE(6,104)((S00(1,J),I= 1,8 ),J=1,6)

104 FOR1AT(4nH DIRICHLET PRORLEM AFTER 20 ITERATIONS

_/(EF10.6))

RHO

= 0.0

DO 28

I

= 1,16

DO 22

J

= 1,8

VV= SORT(X(I)X(I) + Y(I,J)Y(I,J) + Z(I,J)*Z(I,J))

22 Pi0 = PHO - 4.0Sn0(1,J)

F( ,u) *

P

/

VV

IJ = DO 35 I = 1,16 DO 35 J = 1,8

Id

=

IJ + 1

COLE = SlE(I,J) / SD0(1,0)

KL =

0

SUM = 0.0

DO 33 K = 1,16

DO 33 L = 1,E

KL =

KL + 1

33 SUM = SU1

_{+ P1(1J,KL)(SDA(KL)-SD5(KL)COEF) * F(K,L) *}

P

35 PHI(I,J) = X(I) - SUM + PHO

* COLE

WRITE(6,107) (X(I),I=1,2)

107 FORAAT(1H ///(8E10.6))

WRITF(6,110)((PHI(!,J),1=1,2),J=1,2)

WRITE(6,107) (X(I),I=9,16)

WRITE(6,110)((PHI(I,J),I=9,16),J=1,8)

CALL EXIT

PM

(35)

-Technical Library Building 131

Aberdeen Proving Ground, Maryland 21005

Defense Documentation Center (20) Cameron Station

Alexandria, Virginia 22314

Technical Library

Naval Ship Research and Development Center Annapolis Division

Annapolis, Maryland 21402

Professor Bruce Johnson Engineering Department Naval Academy Annapolis, Maryland 21402 Library Naval Academy Annapolis, Maryland 21402 Professor W. P. Graebel Department of Engineering Mechanics

The University of Michigan College of Engineering Ann Arbor, Michigan ₄₈₁₀₄ Professor W. R. Debler

Department of Engineering Mechanics University of Mechanics

Ann Arbor, Michigan 48108 Dr. Francis Ogilvie

Department of Naval Architecture and Marine Engineering

University of Michigan Ann Arbor, Michigan 48108 Professor S. D, Sharma

University of Michigan Ann Arbor, Michigan ₄₈₁₀₈

DISTRIBUTION LIST FOR TECHNICAL REPORTS ISSUED UNDER CONTRACT Nonr-1611(07) TASK 062-183

OFFICE OF NAVAL RESEARCH

Professor W. W. Willmarth

Department of Aerospar,e Engineering University of Michigan.

Ann Arbor, Michigan 48108

Professor Finn C. Michelsen Naval Architecture and Marine

Engineering

445 West Engineering Bldg. University of Michigan Ann Arbor, Michigan ₄₈₁₀₄ AFOSR (REM)

1400 Wilson Boulevard

Arlington, Iiirginia ₂₂₂₀₄

Dr. J. Menkes

Institute for Defense Analyses

400 Amy-Navy Drive

Arlington, Virginia ₂₂₂₀₄ Professor S. Corrsin

Mechanics Department

The Johns Hopkins University Baltimore, Maryland 20910

Professor O. M. Phillips The Johns Hopkins University Baltimore, Maryland 20910

Professor L. S. G. KovasznaY The Johns Hopkins University Baltimore, Maryland 20910

Librarian

Department of Naval Architecture University of California

Berkeley, California 94720

Professor P. Lieber

Department of Mechanical Engineering University of, California

Institute of Engineering Research Berkeley, California 94720 Professor M. Holt

Division of Aeronautical Sciences University of California

(36)

Professor J. V. Wehausen

Professor J. R. Paulling

Professor E. V. Laitone

Department of Mechanical Engineering University of California

School of Applied Mathematics Indiana University

Bloomington, Indiana 47401 Commander

Boston Naval Shipyard

Boston, Massachusetts 02129 Director

Office of Naval Research Branch Office

495 Summer Street

Boston, Massachusetts 02210

Professor M. S. Uberoi

Department of Aeronautical Engineering University of Colorado

Boulder, Colorado 80303

Naval Applied Science Laboratory Technical Library

Bldg. 1 Code 222

Flushing and Washington Avenues Brooklyn, New York 11251

Professor J. J. Foody

Chairman, Engineering Department State University of New York Maritime College

Bronx, New York 10465 Dr. Irving C. Statler, Head Applied Mechanics Department

Cornell Aeronautical Laboratory, Inc. P. O. Box 235

Buffalo, New York 14221

Dr. Alfred Ritter

Assistant Head, Applied Mechanics Dept. Cornell Aeronautical Laboratory, Inc. Buffalo, New York 14221

Professor G. Birkhoff Department of Mathematics Harvard University

Cambridge, Massachusetts 02138 Commanding Officer

NROTC Naval Administrative Unit Massachusetts Institute of Technology Cambridge, Massachusetts 02139 Professor N. Newman

Massachusetts Institute of Technology Cambridge, Massachusetts 02139 Professor A. H. Shapiro

Department of Mechanical Engineering Massachusetts Institute of Technology Cambridge, Massachusetts 02139

Professor C. C. Lin

Department of Mathematics

Massachusetts Institute of Technology Cambridge, Massachusetts 02139 Professor E. W. Merrill

Department of Mathematics

Massachusetts Institute of Technology Cambridge, Massachusetts 02139 Prafessor M. A. Abkowitz

Massachusetts Institute of Technology Cambridge, Massachusetts 02139 Professor G. H. Carrier

Department of Engineering and Applied Physics

Harvard University

Cambridge, Massachusetts 02139 Professor E. Mollo-Christensen Room 54-1722

Massachusetts Institute of Technology Cambridge, Massachusetts 02139 Professor A. T. Ippen

Department of Civil Engineering

Massachusetts Institute of Technology Cambridge, Massachusetts 02139

(37)

-2-Commander

Charleston Naval Shipyard U. S. Naval Base

Charleston, South Carolina

₂₉₄₀₈

A. R. Kuhlthau, Director

Research Laboratories for the Engineering Sciences

Thorton Hall, University of Virginia Charlottesville, Virginia

22903

Director

219

Dearborn Street

Chicago, Illinois

6o6o4

Library

Naval Weapons Center

China Lake, California

₉₃₅₅₇

Library MS

60-3.

NASA Lewis Research Center

21000

Brookpark Road Cleveland, Ohio

₄₄₁₃₅

Professor J. M. Burgers

Institute of Fluid Dynamics and Applied Mathematics

University of Maryland

College Park, Maryland

₂₀₇₄₂

Acquisition Director

NASA Scientific & Technical Information

P. 0..Box 33

₂₀₇₄₀

Professor Pad

Institute for Fluid Dynamics and Applied Mathematics University of Maryland

₂₀₇₄₀

Technical Library

Naval Weapons Laboratory

Dahlgren, Virginia. ₂₂₎₀₁₈

Computation & Analyses Laboratory Naval Weapons Laboratory

Dahlgren, Virginia 22)018

-3-Professor C. S. Wells LTV Research Center Dallas, Texas

₇₅₂₂₂

Dr. R. H. Kraichnan

Dublin, New Hampshire ₀₃₄₎₁₄ Commanding Officer

Army Research Office Box CM, Duke Station

Durham, North Carolina

₂₇₇₀₆

Professor A. Charnes

The Technological Institute Northwestern University Evanston, Illinois

60201

Dr. Martim H. Bloom

Polytechnic Institute of Brooklyn Graduate Center, Dept. of Aerospace Engineering and Applied Mechanics Farmingdale, New York

₁₁₇₃₅

Technical Documents Center Building

₃₁₅

U. S. Army Mobility Equipment Research and Development Center Fort Belvoir, Virginia

₂₂₀₆₀

Professor J. E. Cermak

College of Engineering Colorado State University Ft. Collins, Colorado

₈₀₅₂₁

Technical Library

Webb Institute of Naval Architecture Glen Cove, Long Island, New York

₁₁₅₄₂

Professor E. V. Lewis

Webb Institute of Naval Architecture Glen Cove, Long Island, New York

11542

Library MS

285

NASA, Langley Research Center Langley Station Hampton, Virginia

₂₃₃₆₅

Dr. B. N. Pridmore Brown Northrop Corporation NORAIR-Div, Hawthorne, California

₉₀₂₅₀

(38)

Dr. J. P. Breslin

Stevens Institute of Technology Davidson Laboratory

Hoboken, New Jersey Mr. D. Savitsky

Hoboken, New Jersey Mr. C. H. Henry

Hoboken, New Jersey 07030

Professor J. F. Kennedy, Director Iowa Institute of Hydraulic Research

State University. of Iowa

Iowa City, Iowa 52240

Professor L. Landweber

Iowa Institute of Hydraulic Research State University of Iowa

Iowa City, Iowa 52240

Professor E. L. Resler Graduate School of

Aeronautical Engineering Cornell University

Ithaca, New York 14851 Professor John Miles

% I.G.P.P.

University of California, San Diego La Jolla, California 92038

Director

Scripps Institution of Oceanography University of California

La Jolla, California 92037

Dr. B. Sternlicht

Mechanical Technology Incorporated 968 Albany-Shaker Road

Latham, New York 12110 Technical Libraiy

USN Underwater 141,eapons &

Mr. P. Eisenberg, President Research & Engineering Station

Hydronautics Newport, Rhode Island 02840

Pindell School. Road

Howard County

Laurel, Maryland :20810 Professor A. Ellis

University of California, San Diego

Department of Aerospace & Mech. Engrg. Sci. La Jolla, California 92037

-4-(2)

Mr. Alfonso Alcedan L. Director Laboratorio Nacional De Hydraulics Antigui Cameno A. Ancon

Casilla Jostal 682

_ Lima,

Pera

Crommander

Long Beaeh Naval Shipyard Long Beach, California 90802 Professor John Laufer

Department of Aerospace Engineering University Park

Los Angeles, California 90007

Professor J. Ripkin

St. Anthony Falls Hydraulic Lab. University of Minnesota

Minneapolis, Minnesota

₅₅₄₁₄

Professor J. M. Killen

St. Anthony Falls Hydraulic Lab. University of Minnesota

Minneapolis, Minnesota

55414

Lorenz G. Straub Library

St. Anthony Falls Hydraulic Lab. Mississippi River at 3rd Avenue SE. Minneapolis, Minnesota _55/114 Dr. E. Silberman

St. Anthony Fails Hydraulic Lab. University of Minnesota

Minneapolis, Minnesota ₅₅₁₄₁₄ Superintendent

Naval Postgraduate School

Library Code 0212

Monterey, California

99940

Professor A. B. Metzner University of Delaware Newark, New Jersey 19711

Technical Library

Underwater Sound Laboratory

Fort Tru:Mbull

New London, Conneeticut 06121

nMIM

07030

(39)

Professor J. J. Stoker

Institute of Mathematical Sciences New York University

251 Mercer Street

New York, New York 10003 Engineering Societies Library 315 East 47th Street

New York, New York 10017

Office of Naval Research New York Area Office 207 W. 24th Street

New York, New York 10011 Commanding Officer

Box 39

FPO New York, New York Professor H. Elrod

Department of Mechanical Engineering Columbia University

New York, New York 10027

Society of Naval Architects and Marine Engineering

74 Trinity Place

New York, New York 10006 Professor S. A. Piascek

Department of Engineering Mechanics University of Notre Dame

Notre Dame, Indiana 46556

United States Atomic Energy Commission Division of Technical Information

Extension P. O. Box 62

Oak Ridge, Tennessee 37830

Miss O. M. Leach, Librarian National Research Council Aeronautical Library Montreal Road

Ottawa 7, Canada Technical Library

Naval Ship Research and Development Center

Panaman City, Florida 32401

Library

Jet Propulsion Laboratory

California Institute of Technology 4800 Oak Grove Avenue

Pasadena, California 91109

Professor M. S. Plesset Engineering Division

California institute of Technology Pasadena, California 91109

Professor H. Liepmann Department of Aeronautics

California Institute of Technology Pasadena, California 91109

Technical Library

Naval Undersea Warfare Center 3202 E. Foothill Boulevard 09510 (25) Pasadena, California 91107

Dr. J. W. Hoyt

Naval Undersea Warfare Center 3202 E. Foothill Boulevard Pasadena, California 91107

Professor T. Y. WU

Department of Engineering

California Institute of Technology Pasadena, California 91109

Director

1030 E. Green Street

Professor A. Acosta

Department of Mechanical Engineering California Institute of Technology

Naval Ship Engineering Center Philadelphia Division

Technical Library

Philadelphia, Pennsylvania 19112 Technical Library (Code 2)49B) Philadelphia Naval Shipyard