Experience in computing the wavemaking of source/sink models

1.

Udr

2 L,v-v

L

edø-*f

Lab.

y.

Scheepsbouwkunde

Technische

_HO9OQreport

Ship 190

Delft

September 1975

National

Physical

Laboratory

Ship Division

EXPERIENCE IN COMPUTING THE

WAVEMAKING OF SOURCE/SINK

MODELS

by R G Standing

Department of Industry

27 JÄ

76

ARCHIEF

SIVITR-7526

No extracts from this report may 1e reproduced. without the prior written consent of the Director, National Physical Laboratory.

The source must be acknowledged.

Approved on behalf of Director, NFL, by Mr J A H Paff ett, Superintendent of Ship Division

EERIENCE IN COMPUTING THE WAVEMAKING OF SOURCE/SINK MODELS by R.G. Standing Ship Division

SU4ARY

In the linear theory of ship wavemaking the hull is normally replaced

by an array of fluid sources and sinks. There is usually no particular

difficulty in computing the free wave pattern downstream, but the

disturbance close to the hull is more complicated, involving both single

and double integrals for each source. The purpose of this report is to

put on record the experience gained in writing computer programs to

evaluate these terms. The programs are used to compute wall and bottom

corrections and contributions to the bow wave of a simple ship model. NFL Report Ship 190

Introduction

General approach ₃

2.1

Principal value integrai 6

2.2

Integrals with respect to e io

2.3

Alternative forms of the integrals 12

Point source in a steady stream of infinite depth 13

3.1

Theory 13

3.2

Computing method when the integrand has 16

exponential decay

3.3

Computing method when there is no exponential decay

20

J-t. Point source in a steady stream of finite depth

23

4.1 Theory

23

4.2

Computing method

₂₁

!.3

The effect of restricted water on flow past a model

28

Uniform vertical panel of sources in a steady

30

stream of infinite depth

5.1

Theory

30

5.2

Computing method for

i(u,)

when < 0 32

5.3

Computing method for

i(w,)

when = O 314

5.14 _{Bow wave calculations for a simple ship model 31}

Conclusion 39

Notation 140

Page

Appendices 1. Computation of sine and cosine integrals

The sign of A when

- /2 ' li.Î

Analytic evaluation of integral (21)

)4 _{Analytic evaluation of integral (32)} )49

Asymptotic behaviour of _ZL when R is large 50

The theory of ship wavemaking often combines the use of linear wave theory

with a source and sink model of the hull In thin ship theory the sources

lie on the ship's centre-plane, while in other models (for example

reference i) they lie on the hull surface. The disturbance consists of

local and free wave components, of which the latter carry away wave energy

and are associated with wave drag. These components can be expressed in

ternis of single and double integrals, derived by Lunde (2) for point

sources moving in deep and shallow water, in channels and open water. The

double integral normally contains the local (or transient) disturbance and

sometimes part of the free wave term, while the single integral contains

the remainder of the free wave. The two integrals interfere upstream to

leave only the rapidly decaying local term, and. they reinforce downstream.

The single integral is often quite straightforward to compute,

especially if the water is bounded either side by channel walls, as in

reference (3), when the single integral becomes an infinite sum. More

serious difficulties arise when computing the double integral. Wigley _»4)

had already overcome some of these problems by 1949. He published a short

table of the double integral for a point source moving in deep open water,

and. referred to an earlier unpublished table of the single integral.

Guilloton also, in a series of notable papers, particularly reference (5),

tabulated surface wave elevations and sub-surface pressures due to moving

linearised 'wedges', equivalent to panels of distributed sources. He used

a graphical rather than numerical integration procedure.

Recently investigations have required further calculations of these

integrals. Similar programs have also been developed in the context of

hovercraft wavemaking. This report describes the techniques used. Many

of these are common to several programs and some are described in earlier

-2

some of these calculations many years ago, but was limited by the low speed

and capacity of early computers. Simple methods have been used wherever

possible because these are easily programmed and modified. Thus for

example the trapezium rule is often used because it is simple and

readily adapted to cope with singularities and local rapid oscillations.

The programs have been tested over a reasonably wide range of

conditions and checked against published data. Refined integration

procedures (e.g. halving the step length) and checks on continuity suggest

errors which are generally less than 1-2% and rarely more than 5%, but no

error analyses have been carried out and accuracy cannot be guaranteed in

all cases.

Some of the most frequently recurring problems are described in

section 2. Techniques are suggested for eliminating poles, truncating

infinite ranges of integration, and coping with rapid oscillations.

Specific applications and further details are contained in the following

sections.

Ship hulls are often modelled by arrays of point sources. Expressions

are derived in section 3 for the velocity components and surface elevation

induced by a point source in a steady stream of infinite depth. These

expressions are transformed by contour integral methods before evaluation.

Shallow water terms complicate the integrals considerably and. prevent

the use of contour integral methods. The computer program

described in section 4 is used to calculate the effect of a tank's walls

and bottom on the surface wave and axial velocity close to simple 2- and

1f-source ship models. The walls affect both the local and free wave terms

but the bottom affects the local term long before the free wave, in depths

considerably greater than the conventional half transverse wavelength limit.

speeds. _{This program, supplied by NPL Division of Numerical Analysis and}

Computing (DNAC), has many other applications in describing the local wave field

of ship models, but is quite slow to run and becomes inefficient at great distances.

The velocity field around a point source is singular at the source

itself. _{It may be necessary to model certain parts of the hull surface by}

panels of distributed sources instead. In other situations, particularly

at low craft speeds, a simple hull form may be modelled better by a small

number of distributed source panels than by large arrays of point sources.

Section 5 describes methods for use with vertical panels of sources in a

deep steady stream.

The resulting program has been used to calculate the surface wave pattern

around a simple hull consisting of wedges with vertical walls joined

back-to-back. _{Thin ship theory models this by a vertical panel of sources at}

the bow and one of sinks at the stern. _{The program is quite fast and}

accurate both near to and far from the model.

A program, using similar techniques, has been developed to compute

diffracted regular wave loads on fixed bodies of arbitrary shape. It uses

Hess and Smith's method (i) and. is documented fully in reference (6).

Further similar programs to compute the wavemaking of hovercraft are

described in a sequel to this report, reference (ï).

2. GENERAL APPROACH

The disturbance typically takes the form

/2 ( _a(k,O) -Pv

I

klzI+ikw

o Jo k_k* e dk 11/2 +

I

b(0)e_k*jZj+ik*wd6,

je

o I (i)

-1-where PV denotes the principal value of the integral, k and O are

wave ni.miber and direction, O O _'ff12, w

= X COS

O

+ y Sifl

O, k j

the real positive root of

k* = K sec2 O tant k*d

o

where such a root exists, d = water depth, K = glu2, g = acceleration

due to gravity and U = undisturbed upstream flow speed relative to

Cartesian coordinates moving with the steady wave pattern, with

x measured

positive upstream and z positive upwards. The functions a and b

depend on the quantity to be evaluated but are generally analytic functions

of k and O.

The root k* is most conveniently found by iteration. The

Newton-Raphson method gives a rapidly convergent sequence of approximations

a(n

i, 2, 3...) to

k*d using

a = a for a

>1.1

o o

= /3(a -i) for i < a 1.1

o a = n+ i a

_on

T -a

n a ao(i_Tn2)_i'

where a

Kdsec2O, T

= tant a , and no root exists when a < 1.

o o n n o

The problem now is to deal with the singularity at k = k*, the

infinite range of the principal value integral and rapid oscillations as O

approaches

_/2

If z is reasonably large the last two items present no

problems because the exponential term makes the integrand decay rapidly as

k or O increases. The only remaining problem, the singularity at

k = k*, can be ignored if integration steps are chosen suitably to step

symmetrically over it. The integration process is then straightforward,

and any improvements are merely ways in which the process can be made more

longer an exponential decay factor and a straightforward integration is

often impracticable.

Most of the integrals which have been evaluated

successfully with

z = O

assamie deep water.

They are then amenable to

contour integral methods, as explained below, to separate the local and

free wave ternas and transform the double integral into single integrals of

simple functions.

The following notes give guidance on overcoming sorne of these problems.

(i)

Divergent and slowly convergent integrals

Experience has shown the value of examining the asymptotic form

_{of the}

integrand over an infinite range of integration and

near singularities.

The aim is to find out whether the integral is convergent or divergent,

where its range can be truncated safely and how to remove the singularities

and speed up the integration process.

As an

example of the kind of problem which may arise, the integral

is of course infinite.

A simple program might truncate the range

i

when the integrand is small, say at t = 100, giving an entirely erroneous

result log

100 = L.6.

This kind of behaviour may be hidden by the

complicated form of the integral, as for example in Lunde's (2) ecjuation 15.114

and its various alternative forms discussed on page 214.

The following example shows how the integration process can be speeded

up by subtracting a suitably chosen function.

In the integral

dt

(t2 +

sint)-Ji

the integrand behaves as

l/t2 as t --

,

and so must be integrated up to,

-6

first and integrated separately, thus

Sifl

t

1+

dt,

i

similar accuracy is achieved with the range truncated at t = 14 and a

substantial saving in computer time.

Similarly when evaluating

r'

logt cos t dt

Jo

i: is best to subtract out the singularity first, leaving

J'

(cos t-1)log t dt - 1.

In this form the integral is well-behaved at t = O and can be evaluated

straightforwardly.

(ii) Integration method

The trapezium rule was used in almost all the programs described

because of its simplicity. Simpson's rule might have been more efficient

and almost as simple. More sophisticated methods for integrating

oscillatory and exponential functions were not thought worth the extra

programming problems.

2.1 Principal Value Integral

(i) Contour integral methods for deep water

Contour integral methods can sometimes be used to eliminate the pole at

k = k*. This technique was used by Wigley (14), Jinnaka (8), Yeung (9),

Adee

(io)

and Standing (ii), among others. It is most useful when the

water is of infinite depth and a(k,O) is simple. A typical integration

contour is shown in figure 1. The principal value integral is the

around the complete closed contour, which excludes any pole, is zero.

Thus an alternative form of the k-integral consists of a residue term

coming from integrating round the small semi-circular arc around k =

and an integral along the sloping line L passing through k = O. The

line L and large circular arc must be placed in the correct quadrant so

that the integrand is exponentially small on the arc, and in the limit its

contribution to the contour integral is zero. The residue term can usually

be combined with the single integral to form the free wave disturbance,

while the remaining integral along the sloping line represents the local

term. A suitable choice of L (see note ii below) can greatly simplify

and speed up the integration process.

Sections 3.2 and 5.2 contain examples of the use of this method.

(ii) Best choice of contour to reduce oscillations

and eliminate the singularity

Wigley () chose the imaginary k axis as his line L. This choice

transforms the exponential decay term into

an

oscillatory function, and the

oscillatory functions into exponential-type terms. Thus the original form

of the integral is time-consuming to evaluate if z is too small compared

with w; Wigley's form is expensive to evaluate if z is too large. As

a compromise the program should choose one or other form depending on source

depth.

It may be preferable, however, to choose the line L on which the integrand

decays without oscillations. This is simply k = u/(lz -iw), where u

is real and positive. The integration is then straightforward unless

w/zl is small, when the line L approaches the pole at k = k*.

Yeung

₍₉₎

and Adee (io) used this method, expanding the integrand in series

form near the singularity. To avoid this Standing (ii) chose a line L

which never reaches the pole, but on which the integrand. oscillates at most

h(k) + h(k*)log 1 (2) k1 k_k*

-8-k = u(1 +

ic)/(IzI -

i)

where E: exp(-j/zj)

Here the functions in the integrand are slightly more complicated but the

numerical process considerably more straightforward.

(iii) Subtract out the pole when the water is shallow

The functions a(k,O) and b(0) are usually complicated when the water

is shallow. An infinite number of poles on the imaginary k axis must now be

avoided, as well as the single pole k = k* on the positive real axis.

These poles k. = ip.. correspond to the real roots p. of

p. =

K sec2ütanp.

a.

J o J

These complications make a contour integral transformation generally not

worthwhile, though there are exceptions (e.g. Garrison and Chow (12) and

Hogben and Standing

(6)

for wave diffraction problems,and reference (7)

for two-dimensional problems). In these special cases the integral along

the imaginary axis is zero and the k-integral part of equation (i) can be

replaced by an infinite series, each term of which corresponds to a residue

at a pole k.= ip..

In other cases the only practicable approach seems to be the evaluation of

equation (i) as it stands, usually a slow and inefficient process. First

the pole can be removed by subtracting a suitable function whose integral

is known. In Monacella's method (13) the principal value integral is

written

roe h(k) h(k)_h(k*)

pvI

dk =

making the first integrand well-behaved, over its entire range. Note that

both k* and k1 must be finite and that the original integrand h(k)/(k_k*)

is integrated from k = k1 to infinity. A convenient choice is often

k1 = 2k* to make the logarithmic term zero.

If an exactly odd number of steps covers the range k = O to 2k* it

is not strictly necessary to subtract h(k*)/(k_k*) as its contribution at

symmetrically placed points either side of k = k* cancel (see for example

reference i1). But care must be taken not to accumulate too large a

contribution on one side of the pole before subtracting the similarly large

contribution from the other side, as errors may become excessive.

(iv) Truncation of infinite k range

The range of integration must be truncated at some finite value of k.

If there is an exponential decay term the integration process can be stopped

at kizi = 5, say, where exp (-kIzi) < 0.01. But if z is small the

integrand will oscillate many times before this point is reached. Rounding

errors and run time may become excessive, and a different approach may be

appropriate (see notes y and vi below), or an earlier truncation point

needed (see notes i and ii of section 2.2)

(y) At large distances omit the local term

The free wave can usually be written as a single integral, and it may

be adequate to compute only this term several wavelengths downstream of the

source, although this too will prove a time-consuming process at very great

distances. Note that the free wave is not normally identical to the single

integral in equation (i), but has the same form with different limits of

integration.

(vi) Write k-integral as a standard function when z O

If z = O the k-integral of sines, cosines and algebraic functions

lo

-=

_{- Ci(u)cosu-}

_{si(u)sin u,}

where u > O. These can be calculated quickly using Abraniowitz ana.

Stegunts (15) rational approximations

5.2.36

to

5.2.39

for u - 1 and

series expansions

(5.2.1)4

and

5.2.16)

of Ci and si for u < 1 (see

appendix i).

Section

5.3

contains an example of the use of this method.

(vii) Reverse order of double integration at supercritical speeds

The order in which the two stages of the double integration are carried

out is quite arbitrary. It is usual to integrate with respect to k first

at subcritical speeds 1J2/gd < 1. At supercritical speeds however the

pole k = k* does not exist over the range O < e < e with O cos1 V'gd/U2,

and there are problems as e approaches e0. In this case it is easier to

reverse the order of integration and integrate with respect to O first.

This is a principal value integral with pole 0* = cos'

/KO/tanh

kd in

the range O . 0* .

71/2,

and is treated as described above in note liii).

This method is used in section

)4.2.2

of reference (7).

2.2 Integrals with respect to O

(i) Rapid oscillations: vary step length and stretch

range to infinity

The integrand oscillates more and more rapidly as O approaches ¶12.

and shorter integration steps should be used here if the exponential decay Two particularly useful functions

1(u) ₌ = g(u) = j are -ut J sirit e I

dt

dt

(3)

Jo

dt

= t+u Ci(u)sinu- si(u)cos i

cost

1+t2 o u, -ut te dt = u o ..Io l+t2

oscillation cycle into roughly the saine number of equal steps (say 10-20

steps per cycle), but with shorter steps in regions of slow phase variation.

The aim was to associate roughly the saine error with each step. Reference

(1), section 3.2.2, describes this procedure in more detail.

A transformation t = tan O or t = sec O before integration makes

the problem less acute. The tail end of the range is stretched out to

infinity and the phase increases only linearly or quadratically with t,

and so at a more uniform rate (see references 8, 11,

14).

(ii) Stationary phase approach to truncate integration range

Following the above transformation to variable t the range must be

truncated at some finite point. This is straightforward if the exponential

decay is fast enough, for the integration continues until the oscillation

amplitude is small. Even if there is no exponential decay, as in

references 8 and 11, the algebraic decay of other terms may be fast enough,

though care is needed (see note (i) of section 2.1).

Eut if the decay rate is slow the integration must be stopped before

this point is reached. A kind of stationary phase argument can be invoked

to identify the parts of the range which contribute significantly to the

integral. A few oscillation cycles are included on either side of a point

where the phase gradient is zero. This technique was suggested by Tuck,

Collins and Wells (16) who also applied an artificial fade factor to the

iritegrand to make it decay rapidly on either side of each stationary phase

point.

If no fade factor is applied care must be taken to choose a suitable

truncation point in the wave cycle. Huang and Wong

(i4)

suggested

inter-grating as far as a maximum or minimum of the integrand. This represents

a 'mean' value of the integral about which it oscillates on further

+

12

-stationary phase of cos (t) and the integral is truncated at t = t1,

r rt

h(t)cos(t)dt

i

h(t)cos(t)dt

+ + - ç-

-h(t)1

d h(t) sinip(t) - sinì(t) - '(t)J dt

If coslp(ti) is a maximum and h(t1)/4'(t1) is finite, then

ij(t1) = mi

where n is an integer, and the second term involving sini(t1) is zero.

Contributions to the third term from successive half cycles are alternately

positive and negative and normally decrease in magnitude. They may be

regarded as terms in an alternating series. The difference between the

exact and partial sums of such a series is no greater than the first

neglected term. The error in truncating the integral at t = t1 can

therefore be estimated by integrating as far as the next zero of sini(t).

In practice it was found best to compare integrals up to three successive

zeros of sin-ip(t), stopping when the total difference was small.

Reference (7), section 3.2.2, describes one application of this

technique.

2.3 Alternative Forms of the Integrals

The integrals which follow can be written in real or complex form, the

former being more useful for computations, the latter more concise. To

simplify comparisons the following identities will be useful. writing

= x cose + y sine, (71 f71

It!2

1:

i

I

id@

e dO = 2Re ( e dO = 2Re e Jo dt.

f/2

= 1 j

/2

= 2 cos(x cos O)sin(y sin e)dO,

r72

rr/2

i

I

Im i e dO = - Re I je dO

/2

= 2 _{sin(x cose)cos(y sinO)dO.}

Further useful identities are

'7F

J(t)

= -

cos(t sine - nO)dO, . .. (8)

o

where J _{is the Bessel function of the first kind of order n (reference}

n

15, equation 9.1.21), and the source potential

1 2

(y)

1= -

dO

_{dk ecos(kxcosO)cos(kysinO).}

(x2+y2+z2) 7F

o

(9)

Other identities come from integrating or differentiating these formulae

with respect to

x, y or

z.

3. POINT SOURCE IN A STEADY STREAM OF INFINITE DEPTh

3.1 Theory

In linear theory the Green's function or potential _{at any point}

(x, y, z) due to a source of unit strength at

_{(, r,}

_{) in a uniform}

stream of infinite depth moving at constant speed U is given by

Wehausen and Laitone (ii), equation 13.36,

7F 7Í 1W (ji) Im

J

eW

_{d8 = - Re} J ie dO o 'o (6) (Î)

wI e re i 4K 71/2 r y, z) =

- - + -

+

_2

dO PV dk r r1 ir o o o

k(z+)

k(x-)cose]cosk(y-n)sinOJ} / (k cos2O-K0) {e cos

K (z+)sec2O

+4K

e sec2O o

o!

j.0

sin[K0(x_)secëjcos[k0(y_rDsinosec2do,

r 2 = (x-)2 + (y-r)2 + (z-)2, o

(io)

r12 = (x-)2 + (y-ri)2 + (z+)2, K = _glU2,

and PV

denotes the principal value of the integral. Cartesian coordinates

have been chosen with x measured positive upstream, z positive upwards

from the mean free surface.

The velocity of fluid particles is

u = (-u +

/z)

where the relative velocity

r R 4K 2 -- +

dOPvjkdk

By

Bz,

r3

ir J r JA(k,B) kcos2O-K o

and.

K (z+C)sec2ü

o

+ K

2

J

E(0) e

secO dO,

o

o

R = (x-c, y-n, z-c),

o

= (x-c, y-n,

A(k,o) =

{os8

_{sin[k(x_)cosecos[k(y_n)sinoJ,}

- sin® cos[(x_)cosoisin[k(y_n)sinoJ,

cos[k(x-)cosecosfk(y-n)sino},

IB(0) = {cosO cos[K(x_)sece]cos[k(y_n)sjnosec2e1

- sinO sin[Ko(x_)seco1sin[Ko(y_n)sinesec2o1,

sin[k(x_)sece1cos[K(y_n)sinesec2o1}

The ternis involving r and r1 represent a source and image sink in

an unbounded fluid. They may be combined with the double integral term by

means of the identity

1 1 2 - - - = çir/2 {e r r ir o -o o

[(X-

cosocosrk(y-n)sino1 dk

but care must be taken with the modulus signs if the expression is written

in terms of hyperbolic functions.

The surface wave elevation Z comes from the steady Bernoulli ecjuation

2gZ + (q u)2 + q) 2 _{+ q)2}

x

y

Z

at the free surface. When linearised this gives

u

z

= q)/3x on

Z ₌

g

-kz-

_elzl}

16

-The source and sink terms cancel, leaving

I2 gZ(x,y) lK - Í

cosOdOPV'kdk

U T

_J0

_Jo {e'sin[k(x_)coscos1k(y_n)sinej}/(}çcos2e

- K)

K sec2O

+K2j

sec36 e o o cosLK(x-)secO1 cosllc(y_n)sinesec2olde.

it stands, but it is usually more efficient to remove the pole at k = K0sec2O

and reduce the number of integrand oscillations. This is done by writing

the principal value integral as part of an integral around a contour in the

complex k-plane. A typical contour is shown in figure 1. The sloping

line L is arbitrary. Standing's (ii) line L,

(12)

The complex form of this equation is simply

gZ ( x,y) 2K ik ek

fI

= Re joe cosßdO PV U I kcos2O-K o o

K sec2O (i+)

o

+2K2

o

sec3ee

dO },

(13)

where = (x-)cosO + (y-n)sine.

3.2 Computing method when the integrand has exponential decay

and writing P = = P

+iP.,

r i 1+is P. = P - iP. r i = p 2 +

p2

= r i

i+2

where P =

,P.

-r 1+c2

i+2

, ₍₁₅₎

results in a new expression for Z _{which includes the residue term at the}

pole in the single integral, so that

gZ(x,y) 2K

=

_{Re { s}

j sgn()cosOde

U

/2

iu(ucos2O+K0)e1

k = - u(1 + ic)/(

+ i),

(ia)

will be chosen. It includes Adee's (io) and Yeung's (9) as the special

case = O. The reason for choosing non-zero c will be explained below.

Noting that 1. ik {i sgn() ikc;fl Re 1ie = Re e J.

p12 J

_o

u2cosO+2K uP cos2e+K2p2

_or

_o du + where

sgn()

= = 1 -1 O if if if w > O, < O,

w=O,

(see appendix 2)

18

-r/2

_{K sec23(+iw)}

+ 2K2

[i- sgn )]sec3O e ° d o 2K

J2 sgn()

cos6e

u(uPlcos2OfK0P2)e_Udu /2

rl2

Jo u2cose+2Kupcos2e+K2p2

K sec2G

+ 4K 2 sec3O e cos(K u sec2O)dO,

o

tan1

P.coscu + P sincu, i r (P 2-P.2)sincu + 2F P.coscu, r i

ri

x-- with - /2 e' . /2.

t yn

Similar expressions were used by Standing (ii) when evaluating velocity

components at the free surface,

I 2 2 2

uC(u,O) edu

Bx 3y z

z0

dO I

-I

-T/2 i0

n

u2cos0+2K uP cos2 21n12)

e'

_{K sec2B} + 4K 2

_{B(0)sec0 e}

_de, o

i_L

_{"I ¿} where C(u,O) = {Kp2(uP1cos2O-I-K0P2)sgn()cos0, Klpt2(uPicos2O+K0P2)si(w)sinO, - u(uPgcos2O+K P4)cos2ø} o p ₌ (p 2 - P.2)cosEu - 2F P. siricu r i

ri

= r -

ri

p.)coscu + (p.3 - 3F 2P.)sincui r i

(16)

0f

where Pl = P2 = O' =

B(e) =

{cosecos(Ksec2e), sinOcos(Ksec26),

sin(K sec2O)},

o

an alternative form of E(0) in equation (ii).

The choice of a _{is arbitrary, but if} a < i the sine and cosine

functions of au in P1, P2, P3 and P1.f oscillate at most once before

the integrand has decayed to a negligible size. Yeung (9) and Adee (io)

chose a = O. This eliminates all oscillatory behaviour, but as

becomes smaller the line L in figure 1 approaches the singularity at

k = K sec2ü. The integrand. develops a sharp peak and trough either side

of this point, and Adee and Yeung had to expand their integrands in series

form before integrating. With Standing's (ii) choice of non-zero a,

a =

_exp(jj/),

(i8)

the line L makes an angle of at least

45°

with the real k axis and so

never approaches the singularity. In return for a small additional

complication in the algebra,the singularity is removed and. the integrand

oscillates at most once before decaying. The integration is therefore

quick and easy.

A program RGSGTS1, written in BASIC, calculates the three terms in

equation (17'), and was used in the Guilloton transfoimation calculations of

reference (ii). The wave elevation gZ/U is simply the first of these

terms, 3q/Bx. The trapezium rule is used throughout. The infinite u

integral is truncated at u = 5.5 where the exponential term e

Steps u = 0.1 are small enough for practical purposes.

Asymptotic values of these infinite integrals, which omit all but the last

term in the denominator, are used when K0 sec2ep > 30. The ranges

- /2 - e e' and O e 1'/2 are treated separately because the

chosen for the double integral, tSO rr/100 for the oscillatory single

integral, the exact values chosen being divisors of the range in question.

These steps are sufficiently small in many cases of interest, but if is

very small or large the oscillations of the single integrand become too

rapid and a finer division is needed.

3.3 Computing method when there is no exponential decay

The depth decay term becomes unity when calculating the surface

elevation due to a source on the free surface itself. The technique now

is to replace the k-integral by standard functions (c.f. Lamb, reference

18, 213). Setting = O in equation (13) gives

gZ(x,y) U gZ(x,y) U 20 -2K 1 ik = ° cosûdO PV cos26-K Jo _o ¶ _{J_1T/2} k -+ 2K 2

J

sec36 eiK w sec2edOO o

Replacing this by a contour integral as before, where the line L in

figure 1 is now the imaginary axis k = iu for u real and positive, gives

= _sgn()cosOd0 Re { 2k

.ITJ/2

cos2e+K e

f

iu(iu

K 2+u2cosO

o o 0' -iK wsec2ê

+LK2

_secOe

3 0 _dO /2 du

where 2K -Uk)! - o

_s()cos36d@

u2 e ir J_o K 2+u2cos_o re! + K 2 ¡ sec3ûcos(K sec2e)de o

_12

o

Writing y = u cos2O/K makes

O' =

ynl

sgn()

= + 1 - 1 O - -K vlwlsec2ü r ₂ o f u e

v2e

J cos3O du

= Ksec3ê

K +u cos O o o do du - ir!2 . O' . Tri2 1 y2 dv

provided -ir!2 < O < ir!2 and w O. Using Abramowitz and. Stegun's (15)

equation 5.2.12 makes this expression

secO

K sec3O f(K lisec2e),

l-o o

u)!

where f(u) = Ci(u) sin u - si(u)cos u. Thus

gZ(x,y) 2K o U -

_J

sgn()

secO -ir !2 - K sec3O f(K _llsec2e)]de o _o

+ bic

2 _{sec3Ocos(K wsec2O)dO.}

o o

-_ll

(19) (20) if if if w > < = O, O, O. (see appendix 2)

Although the first term can be integrated separately it is better to retain

it in the integrand. The second term on its om behaves badly at O = ±

and needs the first term to regularise it. The singularity in the first

terni at O = O', where = O, is best removed by integrating the term

analytically over a finite range containing this point. In this respect

note that b rb secO secO sgn(w) = dO for ''T/2 <(a, b)< 'T/2 wi a w 1tan b dt = where t = tanO tan a (x-)+(y-n)t (x-g) + (y-n)tan b = log (y-n) 22 -(x-c) + (y-n)tan a

which is well behaved even if y = ri

provided x

and tan a, tan b are

finite. The function f is evaluated quickly and easily using

approxinia-tions in reference (15), quoted in appendix 1.

The last integral presents the most serious difficulties. The

oscillations of the integrand become more rapid and unbounded as O

approaches -r/2. The transformation t = tan O makes the integrand

better behaved, but it still oscillates unboundedly. No computations have

been attempted. It may be possible to ignore oscillations beyond a certain

point (see note i of section 2.2). If not it may be enough to find the

asymptotic formula only (as for exemple Ursell (19)) 'by the method of

1i. POINT SOURCE IN A STEADY STREAM OF FINITE DEPTH

)4 i Theory

In linear theory the Green's function or potential at the point

(x, y, z) due to a source of unit strength at

(, n,

) in a uniform

stream of depth d moving at constant speed U is given by Kostyukov (20)

equations 2.177 to 2.119,

1

y, z) =

- - - - + -

de v dk

r r ir

o 2 o

{e'cosh

k(z+d)cosh k(+d)cos[k(x-)cosej

cos(y-n)sin(kcos2e + K) -

K}/

{ (kcos2ü - K tanh kd) cosh kd}

+ K

J/2

dO {cosh k*(z+d)cosh k*(+d)

sin Lk*(x_)cosó cos [k*(y_n)sin

(cosZ800sh2k*d

_{- Kd)}

₍₂₁₎

where K g/U2, k* is the real positive root of

- K sec2O tanh k*d = O o r02 =

(x-)2

e o

+ (y-)2 + (z-)2, r22

=

(x-)2 + (y-)2 +

=

cos1

Vk for K0 d < 1 = O

for Kod-1

The principal value PV is required over the whole range O e ir!2

214

-at supercritical speeds Kd

< 1. Again Cartesian coordinates have been

chosen with x measured positive upstream, z positive upwards from the

mean free surface.

Kostyukov's expression, equation (21), differs from Lunde's (2)

equation 15.11, Wehausen

and

Laitone's (17) equation 13.37 and. Wehausen's

(21) equation 3.51b. The potential is not unique

and

the term

14 (ff12

o

-

aOPVJ

-K dk

Ttj

(kcos2S-Ktanh kd)cosh kd

represents an arbitrary additive constant which makes the integrand

well-behaved as k O and k - oe so that the integral (21) is finite.

Lunde

omitted this term, so that his integrand behaves as 1/k as k -3 O.

Wehausen included this term, but without cosh kd in the denominator.

(This is a misprint, acknowledged in a private consrnmication.) His

integrand is well-behaved at k = O but behaves as 1/k as k -

oe

In

both cases the integral is infinite. Reference (ii) also contains a

misprint: moving the opening square brackets to just inside the integral

'Tr/2

+ 14K seC2O

ek*+

fl*(x)co5eJcos*(yn)sinode

oj

o

(22)

where k*

Ksec2O.

This is an alternative form of equation (io), obtained

using the identity

signs recovers equation (21) except for a finite constant.

Note that as d - oe _{tanh kd - 1, hr2 - O, cosh kd} ekd,

cosh k(z+d)

ek),

cosh k(+d)

1 2 Tl"2 foe so that i kcos2O+K0 y, z)

- -

+

-r o

f

o dO PV o

i

kcos2O-K J

I2

k(z+)

2

f

e

äi

e

cos[k(x-)coscosrk(y-n)sin.

ri

The velocity of fluid particles in shallow water is

u

=

(- U +

/Bz),

where from equation (21)

R

(/x,

/3z)

= +

r3

r23

o

k

(kcos2O+K )cosh k(+d)

+ -

dO PV

dk k A1(k,e) e

d

o

Tr

-

(kcos2ü-K tanh kd)cosh

o o

o

+ 1K

{

k*B1(0)cosh k*(+d)

o I

2000sh2k*d

- K d

} de,

cos o o

and

equation

R =

(x-c,

y-ri,

o

=

(x-e, y-n,

(k,O)

=

cosh k(z+d)sin[k(x-)cosOjcos[k(y-n)sinO]cose,

- cosh

k(z+d)cos[k(x-)coséJsin(y-ri)sinOJsine,

sinh k(z+d)cos[k(x-)cosOjcos[k(y-n)sine] }

B1(o)

=

{cosh k*(z+d)cos[k*(x_)cosO cos[k*(y_n)sineJ

cosO,

- cosh

k*(z+d)sin[k*(x_)cosOJsin[k*(y_n)sinOjsinO,

sink k*(z+d)sin*(x_)cosOJcos[k*(y_n)sineJ }

The free surface elevation

Z

comes from the linearised Bernoulli

so that 26 -U Z

= -

/3x on

z =

O, g

-

I cosOdO PV dk

k

e_kd

I

kcos2ü+K

7rJ

J I kcos2O-K

tanh

kd

cosh k(+d)sin(x-)cosOJcos[k(y-n)sinO1

k*cosücosh k*d + 4K o e o cosj*(x_)cose]cos[k*(y_n)sineJde, where k* - K sec2etanh

k*d.

= O. o

The first two terms combine, using equation (9), to give

1

-

cosüdO

dk

k ecosh k(+d)sìnk(x-)cosèJ

Ti

'o

cos[k(y-n)sinej,

and

using the identity

i + tanh kd = e sech kd

gives finally the wave elevation,

Ti

/2

foe ksech kd gZ(x,y) o K J cosOdO PV J dk

U 7F (kcos2ü-K tanh ka)

o

cosh k(+d)sin[k(x-)cosO]cosc(y-n)sine

k*cosOcosh k*d 1K

cosh k*(+d)cosLk*(x_)cose

(cos2Ocosh2k*d_K a) e o o [k*(y_n)sine]de. (cos26cosh2k*d_K d) o cosh k*(+d) (2).)

gZ x,y) (x-e) (x-e)

k sech kd (kcos2B-Ktanh kd.) k*cose cosh k*d -o + 2K +

f

cosh

k*(+d)e*W

dG } o I (cos2Gcosh2k*d_K d)

_j_/2

-.- o o

where

= (x-)cosG+(y-)sine

and the single integral is now over the

range _/2 to _/2 excluding -G to O

o o

1.2 Computing method

In an unpublished note dated September

1970

G.F. Miller of NFL DNAC

carried out a similar analysis to obtain equation (2)), and separated it

further into local and free wave components by means of Wigley's ()4)

contour integral transformation. Under Miller's direction B.T. Hinde

wrote a KDF9 USERCODE procedure DELINT to evaluate these components when

Kd >

1. This procedure was later supplied to Ship Division.

DBLINT uses an integration subroutine INTSPT, based on Newton-Cotes

five point quadrature rules. The integration interval is selected

automatically according to the local behaviour of the integrand. The

pole at k = k* is removed by a variant of Monacella's method (13). The

k

single and double integrands decay roughly as e and e , but

oscillate as cos kX and cos k*X where X2 = (x-)2 + (y-n)2. The

integration process is only efficient when X/ is of moderate size, and

run times of an hour or more were experienced when

Ix/j

exceeded 10.

The accuracy of such calculations comes into doubt because of the accumulation

(25)

The complex form of equation (2lt) is

gZ(x,y) =

Re12K

I2

cos@dOPV

c U

t1T

cosh k(+d) ie

- 28

of round-off errors. Wigley's (14) contour integral transformation might

be expected to give a more suitable expression for evaluation when

is large. This contour includes the imaginary axis k = Ip however, and

the integrand has an infinite number of poles at k = ii. where p are

the real roots of

p.cos2O - K tan p.d O

J o

The resulting expression for Z is complicated and no attempt was made to

evaluate it.

14.3 The effect of restricted water on flow past a model

Blockage corrections are normally applied to measurements made in

towing tanks to eliminate tank wall and bottom effects. Subroutine DBLINT

was used to investigate these effects theoretically and justify the use of

wind tunnel blockage factors. Corrections were calculated to the axial

velocity around a simple 'ship' model consisting of point sources and sinks

on the tank centre-plane y = O, the tank having width w and depth d.

The correction due to finite tank depth, but with laterally unbounded (or

open) water, is computed first. On this is superimposed a wall correction

in which the walls are replaced by source and sink images at laterál points

y = ±w, ±2w,

±3w,

etc., the wall correction being the flow field associated

with these images in water of depth d. Comparisons were made between

results at zero and non-zero Froude number F = TJ/v' where L = model

length. At F = O the upper surface is a rigid lid, and tank walls, top

and bottom can be replaced by a doubly infinite array of images in the

transverse and vertical directions.

The first comparisons were made using a model consisting of a point

opposite sink at (0, 0, - 0.1L). Figure 2(a) shows profiles of the axial

c J x

velocity (equal to wave elevation gZ/U2) on the model centre-plane

U

y = O in deep open water. The profiles at F = 0, 0.223 and 0.316 have

completely different forms. Thus at F = 0.316 the wavelength is about

0.6L, reduced to 0.3L at F = 0.223. At F = O the free surface

n n

wave effects disappear.

Although the open water profiles are dissimilar the bottom and wall

corrections are markedly similar in form. Figures 2(b) and (c) show

corrections due to a tank bottom at d = L with additional corrections due

to one or two image pairs when the tank width w = L and 2L. These

images represent a substantial part of the wall correction, but it was

impossible to compute corrections due to images at greater distances

because computer run-time became excessive. Indeed the second image

corrections may be inaccurate because of the number of computer operations

required. The effect of further images may perhaps be estimated by

comparison with corresponding zero Froude number results which show the full

wall correction.

Figures 2(d) and (e) show corrections for shallower water d = 0.5L

when w = L and 2L. Figures 2(f) and (g) show the effect of changing

the source depth, with source and sink now at (L, 0, - 0.05L), (o, o, - 0.05L).

Deep open water profiles are shown in figure 2(f), corrections in figure 2(g)

when d = L, w = L. Finally figures 2(h) and (i) show the effect of

splitting the source and sink in two, thus simulating a slightly more

realistic ship model. Half strength sources are placed at (L, 0, - 0.1L),

(0.75L, 0, - 0.1L) and half strength sinks at (0.25L, 0, - 0.1L),

(o, o, - 0.1L). Again the corrections in figure 2(i) are similar but the

30

-In most cases the corrections increase only slightly with Fraude

number, the differences being least in the deepest and widest tank, w = 2L,

d = L. The difference is more noticeable on reducing either w or d,

especially at the higher Fraude number Fn = 0.316 when the corrections

also lose their fore and aft symmetry. Changing the source and sink depth

or splitting them in two had little effect on the corrections.

The conclusion is that although deep open water profiles change

radically, the corrections due to tank walls and bottom are relatively

insensitive to Fraude number. Wind tunnel corrections (F_n = o) seem to

be adequate up to about F = 0.3 for typical models in typical towing tanks.

These results are discussed more fully and compared with other blockage

estimates in an unpublished note entitled 'A short note on tank blockage'

by G.E. Gadd of NPL Ship Division.

Note that the finite width correction comes almost equally from the

local and free wave terms. The finite depth correction is entirely a local

effect, the water being too deep (gd/U2

= 5,

10 and 20) for the free wave

to reach the bottom. It is a common feature of such problems (e.g.

two-dimensional hovercraft wavemaking in reference 7 and section

5.)

that the

local term feels the bottom long before the free wave, in water depths

considerably greater than the usual half transverse wavelength.

5. UNIFORM VERTICAL PANEL 0F SOURCES IN A STEADY STREAM

0F INFINITE DEPTH

5.1

Theory

Although it is usually convenient to model ships by arrays of point

sources it may occasionally be more efficient to use distributed sources

and integrate analytically over each source panel before attempting

numerical work.

Again Cartesian coordinates are chosen moving with the source sheet

upstream. The undisturbed flow speed is U upstream. A wedge with side

walls defined by y = ± (x-1) tan o. extends downstream from x = to

x = and from the free surface down to z = ?. This wedge is modelled

according to thin-ship theory by a panel of sources with uniform density

U tan a/2rr on the rectangle x - z of the centre-plane

y = O. The water is assumed to be infinitely deep. The wave elevation

Z at the point (x,y) comes from integrating equation (13) over the source

panel, where gZ(x,y) tana U2 K , çç

kW1i±W2)

1-e )(e Re

J

e pV k(kcos2@-K)

z

L z _Jo

K

sec2e _{f iK w1sec2Q iK} 2sec2O

o ti o o I 1-e

)e

-e

)dO i-z!2 where K = g/U2, = (x-1)cosO + y sine, W2 = (x-2)cosU + y sine,

PV denotes the principal value of the integral.

This equation can be split into four components,

gZ(x,y) tana -

t'"

o)

_-

) Tr

-

o) + I(2, ç) U2 I ¿ K /2

roer

i

ek

= Re

_2H

ri _J

LkK

k(kcos28-K O O O

K (i+)sec2O

1

+1

e dO J

-e'd1ç (26)

32

-the first term in -the double integral being an arbitrary addition to remove

the singularity at k = O. A different term is required when ç O.

This is given in equation (31). The factor k can be removed entirely

from the denominator using

i cos28 k(kcos2t3-K ) K kcos26-K k o o o Then since

ikII

cik

f

2+ç2 Re ( (1-e )e

-r = log I

H0

kJ

ç2

(see appendix

3),

the principal value integral can be rewritten

dk

5.2

Computing method for i(,ç) when ç < o

The numerical integration process is speeded up by means of a contour

integral transformation of the principal value. This transformation

removes the pole at k = K sec2O and ensures that the integrand decays to

a negligible size within one oscillation. A suitable contour is shomm in

figure 1, where the line L is

k = -

u(1+iE)/(ç+i),

and once again, for reasons discussed in section

3.2,

a suitable choice of

E is

}.

(21)

(28)

E = _exp(lI!ç). e .1 pij j

'-i

-I

LiÇIWI

e''

-+

0 kKO k(kcos2O-Ko =

i

I

2+ç2 'Ï

'(+Re

( PV

ekD

log1

2K ç2 _J o kcos2O-K

Thus where again I

Re4PV.

dic l. kcos2û-K o o ir K sec2O -o sin(K Iwsec2û)

= e

o' K o

cos2û _{f (ucos2û+K P )coscu-K P.sincu}

+

or

oi

e_Udu,

_{... (29)} K

u2cosO+2uK P cos2û+K 2fP2

O

O -

or

o

-2 =

pl

i+2

- P + iF., r P = , P. = r

_i+2

i 2

Equation (29) is a simplified version of equations

₍₅₎

and

₍₉₎

of reference (ii), but unfortunately the second term was omitted from equation (5),

reappearing correctly in equation (6). Putting equation (29) into i(,)

gives finally i

12

= _{- - cos2ûdû} I

eu

du Tr JO

(ucos2û+K P )cossu-K P.sin su

_or

oi

u2cosO+2uK P cosû+K 2Jp2

-

or

o -2n

lo[

2 -

2+2

dû

-r°' J K sec26 o - 2 e sin(K sec2e)de, o /2 where e

_{= tan'[-(x-)/yI,-1T/2 .}

e' . /2, = (x-)cosG+y sine, s

= exp(/),

P P. and

jJ

are defined above.

r i

The infinite integrals now contain no singularities and decay within one

wavelength. The integration process is simple and fast.

A program RGSZWD1, in BASIC, evaluates equation (30) using the

trapezium rule. This represents the finite draught correction to the

program described in the next section 5.3. Step lengths u = 0.5 between

u = O

and u

= 4.5 were found to be adequate, also se = ¶/10 between

O = /2 and O = /2 in the first two terms of equation (30). Smaller

steps are needed if the panel is of shallow draught. Ari asymptotic formula

replaces the infinite integral when Kipi > 6 cose. For convenience the

last term of equation (30), the free wave correction, is combined with the

free wave part of equation (33).

5.3 Computing method for

i(w,)

when = O

Because the additional arbitrary term in equation (26) becomes infinite,

a different formula is needed when = 0. A suitable choice is

K 'rr/2 .00 _, ikiwl) f i-e

I(,0)

= Re J dO PV j k(kcos2e-K ) ir _J o o iK sec2O

+iJ

e dej.. - 3 -/2

Following Jinnaka

(8)

this integral is transformed by means of a contour

integration to remove the pole at k

Ksec2O.

A suitable contour, shown

in figure 1, has the line L along the imaginary k axis, k = iu where u

is real and positive. This time k is retained in the denominator. The

principal value integral in equation (31)

( CO

ikw(

-uIwl (1-e ) e -1 du

.RePV

dk = K

-k(kcos2O-K ) ° u2cos1O+K 2 u o o o o

+ - sin(Kjsec2O)

(32) o

= -

y - log

T -

g(T) + tr Sin

T},

(see appendix L)

K

o

where T =

_Klwjsec2e

and

g(T)

is defined by Abramowitz and Stegun's

(15)

equation

5.2.1,

g(i) = - Ci(T)cos

T -

si(T)sin T.

Thus

Because there is now no exponential decay terni as e -- ± in2 it is

convenient to stretch the range to infinity to make the oscillations behave

better. Jinnaka's

(8)

transformation t = tanO makes

t,

1 dt dt

i(,o)

- - {y + log

T + g(T)}

+ 2 sin 't

i+t2 i+t2

(33)

= ZL + Z, say,

i(, o)

= -

- y - log

T -

g(T)+T(1-sgn())sin T

where

sgn()

= i if w >

-O,

(see appendix 2)

36

-where _ZL is the local term, Z the free wave term,

t' =

tane'

=

-T =

Rt-t'

/(1+t2)/(1tt2)

R2 = K

[x_)2 + y2J.

Jinnaka's (8) method for computing wave profiles and contour maps of

and _ZL formed the basis of a more recent program RGSZWD1, written in

BASIC. This program also computes finite draught corrections (see section

5.2). Because t' becomes infinite when y = O an alternative expression

is used for T

T = Rsin8' - tcosO'

where tanG'

= -

(x-)/jy

and /2 . O' Tr/2.

The trapezium rule is used throughout.

When evaluating _ZL an asymptotic formula is used for R - 5,

cos2O'

ZL - y - log 2R

R

(see appendix 5). Although {y + log T + g(r)}/(1+t2) has no singularities

and tends to zero as t t' and t -'- its convergence is slow because

of the logarithmic term, and a long range of integration is needed.

Despite its singularity at t = t' it was found more convenient to evaluate

i i g(T)

ZL =

-y-log2R--j

12dl for R<5

-(see appendix 5 (i)). The last term is quite rapidly convergent and the

range is truncated at

t = ± (3 + 1/R) for R > 0.125,

A polynomial approximation to g(T) is used when T

IS

large and a series

approximation when T

IS

small (see appendix i). Step lengths

ät = 0.25 are chosen for t-t' > 0.5, cSt = 0.1 for 0.5 > It-t' > 0.1

and an approximation to .t'+o.i g(_T) dt

(1)

if .t'-ò.i 1+t2

over the innermost region (see appendix

6).

With hindsight, because of

these complications, it might have proved more efficient and accurate to

remove the singularity by subtraction over a finite range surrounding

tt'

before integration.

When evaluating Z the range is truncated at the lower end at

t = - 2 for R > i,

t = - 2/R for i > R >

0.33,

t =

- 6

for R <

0.33.

The upper limit is t', but the range is truncated at t = 5 when t' > 5.

Step lengths are

= 0.05 for R < 10,

0.5/R for R > 10,

so that the integration becomes slow when R 20. As noted previously the

finite depth free-wave part of equation

(30)

is calculated at this stage.

K (i+t2)

This involves multiplying sin T in equation

(33)

by (1-e °

5.) _{Bow wave calculations for a simple ship model}

Experiments were carried out in NPL No. 2 tank to examine the

non-linear effects of craft speed and bow entrance angle. The model, shown

in plan view in figure

3,

consisted of two wedges joined back-to-back

38

-o o

entrance angles of 5 and 10 . The model was towed in both directions at

two speeds, nominally 1.22 and 1.83 m/s, and the wave patterns were

measured by automated and manually operated probes. Reference (li)

compares the test results with the predictions of thin-ship theory and a

transformation method due to Guilloton. Thin-ship theory describes this

double-wedge model by a uniform panel of sources representing the forward

wedge and a uniform panel of sinks representing the rear wedge, both panels

being on the model centre-plane. To simplify the calculations the curved

mid-section was ignored, effectively replaced by a mean straight section

without wavemaking.

Sample computations are shown in figures L and 5. Figure u(a) is a

contour map of wave heights when U = 1 .22 m/s and the model travels with

o

its fine (5 ) end forward, and figure 1#(b) when U = 1.83 m/s with the blunt

(10°) end forward. In the first case the forward panel covers a section of

the x-axis O > Kx > - 26. the rear panel - 30.1 > Kx > - 43.2, and the

draught is K

= -

4.36. Only the bow wave is shown,roughly bounded by

the so-called Kelvin angle at 19° to the x-axis. In the second case the

forward panel covers O > Kx > - 5.86, the rear panel - 1.5 > Kx > - 19.22,

and the draught is K = - 1.94. Here the shoulder wave disturbs the how

wave pattern within one or two wavelengths. In both cases wave height

contours are marked at intervals of 0.02 and 0.04

of K0Z

respectively.

Figures 5(a) and (b) show how the corresponding wave profiles on y = O

are made up. They show the separate contributions of Z, ZL, the finite

draught correction terms and the total sum of all these. In both cases the

model is ciuite deep-draughted and the correction to the free wave term is

small. The correction is almost entirely in the local term. The effect

of this term is a general lowering of the water level around the bow and

the local infinite draught term

ZL.

Upstream where

x > O, Z

= O

and

the correction makes the local term die away more rapidly.

Further calculations, comments and experimental comparisons are

presented in reference (ii).

6.

CONCLUSION

This report describes programs which compute the linearised wavemaking

of point and distributed sources in deep and shallow water.

Single and

double integrals are evaluated to give both local and free wave contributions

to the disturbance.

The Ship Division programs use methods which are simple

to code, of quite general application and easily modified when difficulties

arise.

These difficulties include singularities, rapid oscillations and

infinite ranges of integration.

The programs have been tested for accuracy

over a reasonable range of conditions.

Checks against published data and

refinements of the integration process suggest errors which are generally

less than 1-2% and rarely more than 5%.

Specific problems discussed include the wavemaking of point sources in

deep and shallow water, and of sources distributed over panels.

Sample

results include the computed wave pattern around the bow of a 'mathematical'

double wedge model.

Estimates of blockage corrections to the wave profile

around simple two and four source models show that the zero Froude number

(double body) correction is quite satisfactory over a range of operating

speeds used in routine tank tests.

Programs to predict the wavemaking of hovercraft are described in a

further report, reference (7).

NOTATION

(U _t-1

Y+lou+j

at,

Water depth,

Froude number U/v'gL,

n

í'(u) = Ci(u) sin u - si(u)cos u,

g(u) =-Ci(u)cos u - si(u)sin u,

g = acceleration due to gravity,

i(w,) =

part of wave elevation due to a finite draught source panel:

Im = imaginary part of expression see equation 26,

K

= glU2,

o

k = wave number,

= real positive root of k* = K0sec2ütanh k*d,

L = model length,

P = P

+iP.

=

r i

= r +

P]. = P.

COS EU +

P sin eu,

i r

P2 = (P 2P.2)

S±fl EU +

2F P. cas eu,

r i

ri

P3 =

(p 2p2)

COS Cu +

2F P. sin eu,

r i

ri

PLf _{= (Pr3_3PrPi2)COS}

Cu +

i3_r2Pi)S

Cu,

R,r = general radii,

R = (x-e, y-n,

-o

= (x-e, y-n,

= (x-e, y-n,

Re = real part of expression,

t

4O

-Ci(u) =

d =

r 2 =

(x-)2 + (y_n)2 + (z-2),

O r12 (x-)2 + (y--n)2 + (z+)2, r22 = (x-)2 + (y-n)2 + (z+2d+)2,

sgn()

= + i if w > 0, = - i if w < 0, = O

if w=0,

ru IT Sin -t si(u)

= --+

dt, 2 t o

t, u, y = variables of' integration,

U = model or source speed,

w

= tank width,

x, y, z = Cartesian coordinates, x positive upstream, z positive

upwards from the undisturbed free surface,

Z = wave elevation

ZL, Z = local and free wave parts of the wave elevation due to an

infinite draught source panel: see equation 33,

a = wedge angle,

y = Euler's constant 0.511216,

Sk, t, 6u, = integration step lengths,

s =

arbitrary quantity chosen to be exp(w/),

= Cartesian coordinates of source point,

e = wave angle,

e = cos i/K d. for K d < 1, 0 - e . 1T/2,

o o o = O

for Kd>1,

o

e'

= tan I

x-I-îi]2

e'

ynJ

Ktanhkd ,oe*r/2,

0* =

cosl/

k

-= jth real positive root of u

= Ksec2O tan

t K

IIsec20,

o

= velocity potential in disturbed flow, so that velocity

components are

(-U+/x,

4/y,

B/z),

= phase function,

HESS, J.L. and SMITH, A.M.O. 'Calculation of potential flow about

arbitrary bodies'. Progress in Aeronautical Sciences, vol. 8, 1967,

1-1 38

LUNDE, J.K. 'On the linearised theory of wave resistance for

displacement ships in steady and accelerated motion'. Trans.

SNANE, vol. 59, 1951, 2)4-76

EVEREST, J.T. and HOGBEN, N. 'A theoretical and experimental study

of the wavemaking of hovercraft of arbitrary planform and angle of yaw.'

Trans. RINA, vol. 111, 1969, 3)43-365

)4. WIGLEY, W.C.S. 'L'etat actuel des calculs de resistance de vagues'.

BUll. Assoc. Tech. Mar. Aero, vol. b8, 19)49, 533-587

GUILLOTON, R. 'The waves generated by a moving body'. Trans. RINA,

vol. 102, 1960, 157-173

HOGBEN, N. and STANDING, R.G. 'Wave loads on large bodies'.

Proc. mt. Symposium on the dynamics of marine vehicles and structures

in waves, University College London, 197)4

STANDING, R.G. 'Experience in computing the wavemaking of hovercraft'.

NFL Ship Division report No. 191, 1975

JINNAKA, T. 'Wave patterns and ship waves". Soc. Nay. Archit.

Japan, 60th Anniv. Series, vol. 2, 1957, 83-9)4

YEUNG, R.W. 'Sinkage and trim in first-order thin ship theory'.

J. Ship Res., vol. 16, 1972, pp 1t7-59

ADEE, B.H. 'Calculation of the streii1ines about a ship assuming a

linearised free surface boundary condition'. J. Ship Res, vol. 17,

1973, pp i)4o-i)46

STANDING, R.G. 'Phase and amplitude discrepancies in the surface wave

due to a wedge-ended hull form'. J. Fluid Mech., vol. 62, 197)4,

GARRISON, C.J. and CHOW, P.Y. 'Wave forces on submerged bodies'.

J. ASCE Waterways and Harbors Div., vol. 98, 1912, pp 375-392

MONACELLA, V.J. 'The disturbance due to a slender ship oscillating

in waves in a fluid of finite depth'. J. Ship Res., vol. 10,

1966, pp 2L2-252

lit. HUANG, T.T. and WONG, K.K. 'Disturbance induced by a pressure

distribution moving over a free surface'. J. Ship Res., vol. 1)1-,

1970, pp 195-203

ABRAMOWITZ, M. and STEGUN, l.A. 'Handbook of Mathematical Functions'.

Dover, 1968

TUCK, E.O., COLLINS, J.I. and WELLS, W.H. 'On ship waves and their

spectra'. J. Ship Res., vol. 15, 1971, 11-21

WEHAUSEN, J.V. and LAITONE, E.V. 'Surface Waves'. Encyclopedia of

Physics, vol. 9, 1960, Springer

LAMB, H. 'Hydrodynamics'. Cambridge University Press, 6th edn.,

1932

URSELL, F. 'On Kelvin's ship wave pattern'. J. Fluid Mech,

vol. 8, 1960, 1t18-1t31

KOSTYUKOV, A.A. 'Theory of ship waves and wave resistance'.

E.C.I., Iowa, 1968

11. WEHAUSEN, J.V. 'The wave resistance of ships'. Advances in

Computation of sine and cosine integrals

Sine and. cosine integrals and related f and g functions are

defined by Abrainowitz and Stegun (15) equations (5.2.1) to (5.2.7),

ii

11sint

si(u)

= --+

dt,

They are calculated rapidly using series expansions of Ci and si

(reference (15), equations 5.2.1)4, 5.2.16) when O < u < i,

n 2n

(-i)u

Ci(u)

=y+logu+

n=i 2n (2n)' n 2n+1 (-i) u Si(u) = - ir/2 +

n0 (2n+i)(2n+i)

where Euler's constant y = 0.57721566. These series are rapidly convergent

and are truncated when the denominator is small. When i - u < rational

approximations (reference (15), equations (5.2.38), _(5.2.39)),

i ' u8+a1u6+a2u+agu2+a4 1(u)

= -

, u u8+b1u6+b2u+b3u2+bz 2 o t cos t-i Ci(u)

= y+logu+

dt, t o f(u) g(u) = = Ci(u)sin u - si(u)cos u, - Ci(u)cos u - si(u)sin u.

- 46

i

u8+clu6+c2u+c3u2+cL

g(u)

=

-u2

u8+d1u6+d2u+d3u2+dj'

where

a1 =

_36.027264

b1 =

_40.021433

a2

= 265.187033

b2

= 322.624911

a3

= 335.677320

b3

= 570.236280

a

38.102495

b1

= 157.105423

c =

_42.242855

_al

=

_48.196926

C2

= 302.757865

= 482.485984

= 352.O18498

d3

=1114.978885

cz. =

21.821899

= 449.690326

that

Appendix 2

The sign of w when 'rn2 'E e' 'E 'T/2

Assume that (y-n) - 0. Let (x-e) = - R sine', (y-n) = R cose', so

= (x-)cose + (y-ri)sine = R sin(o-e').

X-e' =

tan1

Note that this depends on having (y-n) - 0, w being of opposite sign when

(y-n) < o.

the range - _{/2 E e /2.} Then

y-n

- ,

this point occurring once in

w < O and sgn(w) = - i for e < e',

-

-Appendix 3

Analytic evaluation of integral (27)

I

I -t e dt

edt

-nç-r]1w1)

= um

. r+O I t

.-

t -r

_-r(+iwl)

=

{i[-rJ

-EER

_Eit_r(+iI)I+Ei[-R(+iwI)1}

where E1 is defined by Abramowitz and Stegun's (15) equation 5.1.1.

Further equations 5.1.11 and 5.1.52 show that

E1(R)--O

as R--,

E1(r)'-y-logr+O(r)

as

r--O,

so that e (1-e )e R ikl

ii0

k

¿+iIwI

}

log =

Re{1o[

= 2 dic ikiwl 1g

(1-e

)e

-k o = 11m OJ

R-R r (1 eikiwl dic 1g k

Appendix

4

lwIsec2S.

o o o

This integral may

be

written

+ _log(R2+1) -

iog(r+1

g(T),

K2

o

where Abramowitz

and

Stegun's (15) equations

5.1.1,

5.2.13 have been used

and

g(T) =

-

Ci(r)cos T - si(r)sin r.

Further equations 5.1.11

and 5.1.52 show that

E1(R)-O as R+,

for finite non-zero T ₌

KIwsec26.

In the limit T -- O the integral

becomes zero because

g(T) " - y - log T O (T) as T -- O

(reference (15), equations 5.2.7'

and

5.2.16).

o"

Analytic

evaluation of integral (32)

-u w e -1 du u 1 I

f

-TV e

-1

dv

u2cosO+K 2

o =

K2

o o

[

v+i jv

for finite non-zero K sec2O. Here u

= y

K sec2O.

and

r

K

R l (

VT

i » -Tv ve

-

lina

¶ -

I-Tv

e -1' I dv +

dv

-K 2 o

r

Lv

v2+j

K 2_{o o} v2+1 = l _{1(r) - E1(RT)+log r - log R}

K

2 o R-+ Taking limits o r - O, R --u w e

-1

E1(r)-y-logr+O(r)

as

r-O.

gives finally

du

1 =

--+lOgT+g(r)

u

K2

o u2cos4e+K 2 o

so that

-

50 -Appendix 5

Asymptotic behaviour of _ZL when R is large

i at i ,2 - - (y + log r) = - y - -

log[RJsin(o_Ot)Jsec2do

lT i+t2 if /2 i

f/2

= _y

--

j

[log R - 2 log cos O + ioglsin(e-O')I]dO.

/2

Now

_n2

log cosO dO = 2 log cosO dO =

-

Tf log 2

/2

from reference (15) equation 4.3.1lt5. Similarly

fr/2

iogsin(O-O')dO

= log cosO dO = - r log 2,

-71/2 so that co dt

--

(y+logT)

- -y-log2R

Ir i+t2 exactly co dt

When R is large the main contribution to - - g(T)

TI i+t2

-co

comes from a small band of width 0(1/H) around t = t'. Writing

T = (t-t')R makes T = Tj

4±

[

+ t1j21/(1+t12) 1'

_ITI

+ 0(1/H), co co 1 dt 1. dT

--

g(r)

--

g(T)

nR -co -co

2

g(T)dT

rrR(1+t'2)

_J0

i cos28'

cos2O'

- y - log 2R - for large R

R

R2 = (x-)2 + y2, tanG'

- (x-)/y.

,

R(i+t'2) - R

where reference (15) equation

1

5.2.30 has been used. Thus finally

dt

ZL

=

--¶

(y+logr+g(T))

-

52

-Appendix 6 Approximation to integral

(31)

t'

+'

g(t)dt g(t)dT 1t'-r -r

i+(t'+r)2

where

T =

t-t', t

= RITI

/1+(tt+T)2/(1+t'2)

t 'T

RITID

+

o(T2fl

2

for small T. If also

RITI «

1,

g(t) - - log T

+ rrt/2 +

0(t2),

so that the above integral is approximately

1

f"

i 2t'T t'T TRITI

-

I i _- y log RITI + dT r j 1+t'2 1+t'2 i+t'2 2

-r

2 lT

(i +t '

2)

i:

1TRT

y+logRT

dT 2

2r

if =

y- i +logRr--Rr

iï(1+t'2)

This approximation is valid for

r «

i

and Rr « 1.

k=k*

CONTOUR OF INTEGRATION

IN THE COMPLEX