Three dimensional green function of a vessel with forward speed in waves

TECHNISCHE UPIIIIVER&Mrg Laboratorlurro WC? Scheepshydromethaniss Maketweg 2 - 2628 CD DELFT

JAN ROGER HOFF

THREE-DIMENSIONAL GREEN

FUNCTION OF A VESSEL WITH

FORWARD SPEED IN WAVES

on

UNIVERSITETET I TRONDHEIM

NORGES TEKNISKE HOGSKOLE MTA-rapport 1990:71

DOKTOR INGENIORAVHANDLING 1990:25

INSTITUTT FOR MARIN HYDRODYNAMIKK TRONDHEIM

NT1141,0001990 TECHNISCHE UNIVERSITET Laboratorlum voor

Scheepshydromechanica

Archlel

Makelweg 2, 2628 CD Delft

Mt: els - 786873 - Fa 015 - 781836

Three-dimensional Green function of a vessel

with forward speed in waves

by

Jan Roger Hoff

Division of Marine Hydrodynamics

The Norwegian Institute of Technology

?.H321144433T

100v

migatecesi

illt#46mInbyttagetsto8 isicknA

.ftleta COIS6S spewtoieM

Abstract

The three-dimensional Green function of a vessel with forward speed in waves has

been calculated. This Green function corresponds to the velocity potential of a

pulsating point source with unit strength moving below a free surface.

The Green function may generally be formulated in terms of a set of double Fourier integrals involving contour integration. In order to facilitate numerical evaluation,

the contour integrals must be reformulated as Cauchy principal value integrals. To further improve the numerical performance of the calculation, the complex exponential integral has been introduced and the double integrals are reduced to 'single' integrals. Effective algorithms for evaluating the complex exponential

integral have been given.

Formulations for the Green function and its derivatives with respect to the three

field point coordinates in terms of single integrals, involving the modified complex exponential integral, have been given both as function of the ordinary integration

variable as well as for the transformations required to handle the singularities in the integration range. A fairly detailed description of the calculation algorithm

has been given.

The given formulation and calculation algorithm have been found to be effective as long as the source is not too close to the free surface. A set of benchmark results accurate to five decimals are given for the Green function and its derivatives with respect to the field point coordinates.

Although the main attention of this work has been on the case r > 0.25, the

given formulations have been found to give accurate results for the case r < 0.25, including the steady case.

Since the formulation involves numerical integration, the resulting computing

re-quirements are fairly large. In order to utilize the vector capabilities of the Cray

supercomputer, a fully vectorizable algorithm has been formulated.

Further, an asymptotic formulation by the method of stationary phase has been

given for the Green function and its derivatives. This asymptotic formulation may

be used inside a wedge behind the source, when the distance between the source

and field point is sufficiently large.

Also the feasibility of approximative evaluation by polynomials has been

inves-tigated. Formulations in terms of Chebychev and ordinary polynomials in four variables have been given. It was found that the time required in the evaluation

of the polynomials would be of the same order as that required by numerical

inte-gration. Further, the calculation of the coefficients in the polynomials_{was found}

to be extremely timeconsuming.

Acknowledgement

The author wishes to express his gratitude to his supervisor, Professor Odd Faltin-sen, whose support is deeply appreciated.

This work was carried out while the author was on leave from the Norwegian Marine Technology Research Institute (MARINTEK). The work was mainly founded by

a scholarship from the Norwegian Institute of Technology (NTH). The financial

support from MARINTEK is appreciated.

The Norwegian Shipowner's Association's Foundation at NTH together with

MAR-INTEK provided financial support for the participation at the 17th International

Symposium on Naval Hydrodynamics.

The Norwegian Research Council for Science and Humanities (NAVF) provided computing resources on the CRAY-XMP at The Computing Centre at the

Univer-sity of Trondheim (RUNT).

The author also want to express his gratitute to Professor Ohkusu and his stu-dent Hidetsugu Iwashita of Kyushu University, Japan, who provided the data for

Contents

Abstract

Acknowledgement

iv

Contents

Nomenclature

ix 1 2

Introduction

The hydrodynamic boundary value problem

1

6

2.1 Coordinate systems 6

2.2 Exact formulation 6

2.2.1 Condition in the fluid 6

2.2.2 Body boundary condition 7

2.2.3 Free-surface conditions 7

2.2.4 Additional conditions 9

2.3 Decomposition of the unsteady potential 9

3

The Green function

11

3.1 The dispersion relation 13

3.2 Formulation as a Cauchy principal value integral 17 3.3 Alternative formulations of the Green function 19

3.4 Nondimensionalization of the Green function 25

3.5 Expressing the Green function in terms of the exponential integral 26

3.5.1 General case 26

3.5.2 Field point on the track 33

3.5.3 The singularity at 8 = 7, general case 34

3.5.4 The singularity at 8 = 7, field point on the track 36

3.5.5 Behaviour close to 0 = r/2, general case 37

3.5.6 Behaviour close to 9 = r/2, field point on the track

Derivatives of the Green function

41

4.1 The x-derivative 41

4.1.1 General case 41

4.1.2 Field point on the track 43

4.1.3 The singularity at 0 = 7, general case. 43

4.1.4 The singularity at 0 = 7, field point on the track 45

...

. . . . . .

. .. ...

. . .

. ..

. . . . .

...

. .

...

. . . .

. ...

. . . . . . . . 401

vi

4.1.5 Behaviour close to 0 = r/2, general case 46 4.1.6 Behaviour close to 8 = r/2, field point on the track 47

4.2 The y-derivative 48

4.2.1 General case 48

4.2.2 Field point on the track 50

4.2.3 The singularity at 9 = 7 50

4.2.4 Behaviour close to 9 = v/2 51

4.3 The z-derivative 52

4.3.1 General case 52

4.3.2 Field point on the track 54

4.3.3 The singularity at 0 = 7, general case 55

4.3.4 The singularity at 0 = 7, field point on the track 57 4.3.5 Behaviour close to 9 = r/2, general case 58 4.3.6 Behaviour close to 8 = r/2, field point on the track 59

5

Numerical evaluation of the Green function

60

5.1 Numerical integration using Gauss quadrature 61

5.2 Numerical integration using an iterative integration scheme. 61

5.2.1 Non-adaptive methods 62

5.2.2 Adaptive integration methods 63

5.3 Approximation by exponential series 64

6 Numerical evaluation of the exponential integral

66

6.1 Local formulas 66

6.1.1 Ascending series 66

6.1.2 Rational function approximation 66

6.1.3 Continued-Fraction expansion 68

6.1.4 Expansion close to the imaginary axis 68

6.1.5 Expansion for large arguments 70

6.1.6 Application of the local formulas 70

6.2 Global formulas 73

6.2.1 The Hess & Smith formula 73

6.2.2 Approximation by exponential functions 73

6.3 Comparison of formulas 75

7 Source point very close to the free surface.

77

7.1 Integration in the complex plane. 81

7.1.1 Transformation of the integrands 81

7.1.2 The complex square root function 86

7.1.3 Singular points in the complex plane 87

7.1.4 Moving the branch-point to the origin 88

7.2 Applying a stretching transformation to the integrand 89

. . . . . . . . .

...

. . . . . . . . . . . . . . a . . . . .. . . . . . . . . . . . .

...

. . . . . . . . .

vu .8

Asymptotic formulation far away from the source

95

8.1 Asymptotic behaviour of the integrals ., . . 95

8.2 Stationary phase approximation . . . .

.

96

8.3 Asymptotic behaviour of the Green function.... , 98

8.3.1 General case ₉₈

8.3.2 Field point on the track 120

&A Asymptotic behaviour of the x-derivative . . . . 124

8.4.1 'General case , ,

,

. 124

8.4.2 Field point on the track

_...

. . 125,

4,5 Asymptotic behaviour of the y-derivitive . . 126.

.8.5.1 General case _126'

8.5.2 Field point on the track . . .

.. .

. 128

8.6 Asymptotic behaviour of the z-derivatfve 128'

8.6.1 General case 128,

8.6.2 Field point on the track

.

129

8.7 Application of the asymptotic formulation - 130

,W !Generated wave systems. ₁₄₃

9.1 Curves of constant phase . _. 143

9.2 Wave systems from each term .. da= m. s. n. foe Iv 152

'9.3! Calculated wave patterns.. . 4.1 4 44.I. ..1.4 V rr.i. ..t 1: .1 160

la Computations

₁₆₈

10.1 Forming a composite algorithm . . . . - . . 168

10.1.1 Partitioning into subintervals and application of the rules! 169

10.2 Coding ₁₇₄

10.3 Computations on the CRAY . , 176

10.3.1 Initial calculations

10.3.2 Forming a fully vectorizable algorithm 178

11 Verification of the calculations

183

'12 Applications

₂₀₈

12.1 A vertical strut in steady motion .

12.1.1 An elliptical strut in steady motion ,

.. 208

12.1.2 A parabolic strut in steady motion

... .

12.2 A vertical strut in unsteady motion . . . .

. ... .

. . '210

13 Approximative evaluation by polynomials!

₂₂₆

13.1 Evaluation using orthogonal polynomials . . . 227 13.1.1 Evaluation of the Chebychev expansions .,

, ,

13.1.2 Evaluation of the Chebychev coefficients . '230 13.1.3 Numerical implementation . . . . .., . 231 177

208 209

viii

13.2 Approximation by ordinary polynomials . . 232

13.3 Application of the polynomial approximations

... 234

14 Conclusion 236

References 238

A Integration paths in the complex plane

242

A.1 The first integrals of To, and lin 242

A.2 The second integrals of /0/ and 102 243

A.3 The integrals /11,112113311-347 -143 and 144 245

A.4 The integrals 121 and Z2 246

Rational function approximation to the exponential integral

248

C Approximation of a function by an exponential series

252

C.1 General algorithm

...

. 252

C.2 Numerical implementation 254

Details in the stationary phase formulation

255

Relation between derivatives in spherical and Cartesian

coordi-nates

258

. .

. .

. . .

ix,

Nomenclature

a, b: 'Lower and upper limit of integration respectively air, hi, At Coefficients in exponential series

Breadth of vessel

'CF Intersection of body surface and the mean position of the free :surface:

Crest number Draught of vessel

E(x)

Principal-value exponential integral 4, Eq, Erj, Modified complex exponential integrals

Ei(z) Complex exponential integral

EST Estimate of integral by numerical integration Froude number based on length of vessel Green function

Acceleration of gravity

h : in Time specification with h-hours:, m-minutes

:1 Imaginary unit, i =

a, hile Unit vectors along the z, y, z axes

Nondimensional wave number, K = U2k1g: Ko Fundamental wave number, ICE, = g 1 U2

Wave number vector

Cartesian wave number components 11, ks, k3,,k4 Singularities in integration range

Contours of integration in Green function

1 Lenght of vessel

Characteristic length

Unit normal vector out into the fluid Pressure

Athmospheric pressure

R,st,,01 Spherical coordinatesucorresponding to he Cartesian

coordi-nates

Distance between source and field point

fir _{Distance between image source with respect to the free surface}

and field point. Cr,

X

r, Polar coordinates in stationary phase formulation Mean wetted surface of vessel

Phase function

T,

Chebychev polynomial of degree

rime

TOL Required tolerance in numerical integration Forward speed

Transformation variable

V Fluid domain

Fluid velocity vector Group velocity vector

Local velocity of wetted surface of vessel

Velocity vector of the steady flow relative to the moving ref-erence frame

Weights in quadrature rules X0, Yo, Z0 Space fixed coordinate system

X, Y, Z Coordinate system fixed with respect to the vessel Abscissae in quadrature rules

x,y,z

Field point coordinates x nondimensionalized by Ko

Nondimensional coordinates of field point in local coordinate

system with origin at image point of source with respect to

the free surface

z, Z Generally used to denote complex quantity Complex conjugate of z

Wave direction

7 Euler's constant. Also used as integration limit in Green func-tion

Small positive quantity

Requested accuracy in numerical integration

ER, Er Estimated error for real and imaginary part of integral Surface elevation

Rigid body motion vector Time derivative of 71, i.e. velocity

Second time derivative of II, i.e. acceleration Transformation variable

xi

0, k Polar coordinates corresponding to the Cartesian components kr.,k,, of the wave number vector k

90 Stationary point

En,

C Source point coordinates

Density of water Source density

Wave form parameter, 7 = 1I g

re Critical value of wave form parameter, 7, = 1/4. <1:( Velocity potential in space-fixed coordinate system (I) Incident wave potential in space-fixed coordinate system

Velocity potential in space-fixed coordinate system due to

steady forward motion

Velocity potential in steady moving coordinate system due to steady forward motion

Unsteady potential in steady moving coordinate system

(Po Incident wave potential in steady moving coordinate system

= 1,6 Radiation potentials

507 Diffraction potential

Wt Subscript t denotes partial derivative with respect to time ,Pet Subscript tt denotes second partial derivative with respect to

time

Subscript r, y, z denotes partial derivative with respect to the field point coordinates x, y, z

Subscript xx denotes second partial derivative with respect to

The ray given by b = arctan(1/2r)

Nondimensional frequency of encounter, 1.1 Ig

Rotational rigid body motion vector

ca Frequency of encounter

Wave frequency

4

1

Introduction

In order to predict the motions of a ship in waves, the ship is regarded as a rigid

floating body with six degrees of freedom. The hydrodynamic forces and resulting

motions are then calculated from linear potential flow theory, assuming that the

fluid is inviscid and incompressible. Further, both the incident waves and the

oscillatory motions of the ship in the six modes of motion are assumed small. Thus, the velocity potential satisfies the Laplace equation and the boundary conditions are imposed on the mean position of the fluid boundaries. Still, after these severe

assumptions, the resulting hydrodynamic boundary value problem is difficult to

solve.

Traditionally, the geometrical shape of the vessel, namely its dominant length in relation to its breadth and draught, has motivated simplification of the

three-dimensional effects of the flow. This has led to slender-body theories and the

so-called strip theories which replaces the complicated three-dimensional problem by a summation of simplified two-dimensional problems for different cross-sections along the hull.

Since the pioneering work of Korvin-Kroukovsky and Jacobs [1959], strip theories

have been widely used in ship motion calculations. The basic assumption in a strip theory approach is that the motion of the fluid in a plane perpendicular to the longitudinal axis of the ship is identical to that of the two-dimensional fluid motion around an infinitely long cylinder with identical cross section to that of

the ship, thus ignoring the three-dimensionality of the flow. This simplification is questionable for low frequencies, Ogilvie and Tuck [1969]. Further, the effect of forward speed is difficult to include in a rational way. The effect of forward speed is only included approximately in the boundary condition on the mean wetted surface of the ship and the frequency of encounter. All wavesystems generated by the ship are not properly accounted for. The most successful and widely used strip theory is that of Salvesen, Tuck and Faltinsen [1970]. Excellent reviews of the different seakeeping theories have been given by Faltinsen [1979] and Maruo [1989] on the evolution of the different slender-body theories.

With respect to seakeeping calculations for traditional ships, strip theories have been found to give good agreement with model tests in the prediction of heave,

pitch and roll motions for high frequencies and moderate Froude numbers.

How-ever, strip theories will give inadequate predictions of local effects as pressure distribution. Strip theory approach has also found extensive use in other areas of

marine hydrodynamics. Especially in the calculation of motion characteristics of

2 Chapter 1. Introduction semi-submersible drilling platforms, where the strip theory has been proven to be very successful.

The limitation on frequencies has in principle been overcome by the 'Unified

the-ories'. The unified theories of Newman [1978] and Sclavounos [1981] takes into account wave interaction between cross sections along the hull. From a rational

point of view, unified theory represents an improved slender body theory.

For non-slender ships the application of the different slender body theories and strip

theories are questionable in the predictions of their seakeeping abilities. These vessels require the solution of the full three-dimensional problem. Still, viscous effects must be neglected.

In addition to its obvious value in the determination of seakeeping qualities of

non-slender ships, the development of three-dimensional theories are also essential to the validation of the different slender body approximations.

The analytic simplification of the full three-dimensional problem into approximate theories, was necessitated by the available computers at the time. The limitations in memory, the low processing speed and the secondary memory access time set se-vere restrictions on the size of the calculations which could be carried out. With the computer technology of today, especially with the introduction of the

supercom-puters, memory and processing time are not a limited resource, and are expected to be even less so in the future. The need for simplifying the three dimensional

problem is thus not so predominant from a computational point of view.

The development of the three-dimensional methods was highly motivated by the offshore activities in the North Sea, which required the accurate determination of the diffraction forces of the large gravity base structures being designed. This

re-sulted in the so-called three-dimensional singularity distribution methods. These

methods involves the distribution of pulsating singularities in terms of sources or dipoles with a known frequency, satisfying the linear free surface boundary

con-ditions. The body boundary condition is then used to determine the strength of

the singularities. Imposing the body boundary condition will result in a Fredholm integral equation of the second kind, and thus these methods are also called bound-ary integral equation methods. Examples of these methods are those of Faltinsen and Michelsen [1974] and Garrison [19781. These methods are however limited to zero forward speed.

The unsteady three-dimensional problem with forward speed, bears strong resem-blance with the wave resistance problem of a ship advancing with constant speed

in calm water. In the limiting case when the frequency of oscillation approaches

zero while the speed remains finite, the unsteady problem will approach the steady

3

frequency of oscillation remains finite, the problem will approach the zero forward

speed problem. Experience from the numerical schemes used in the solution of these two limiting problems may be applied in the unsteady case with forward

speed.

Several alternative numerical methods have been developed to solve the steady problem. The most widely used are the boundary integral equation methods or panel methods. These may further be divided into two groups depending on the form of the Green function being applied. These are the Rankine source method and the Kelvin-Havelock source method. The Rankine source method is based

on the fundamental solution 1IR, while the Kelvin-Havelock method is based on the more complex Green function satisfying the linearized free surface boundary condition. The Rankine source method was first used by Gadd [1976] and Dawson [1977]. The complexity of the Kelvin-Havelock source relative to the simplicity of the Rankine source, seemingly recommend the use of the Rankine source method. However, the Kelvin-Havelock method offers several major advantages with respect to numerical implementation. The radiation condition, namely that steady waves do not exist ahead of the ship, must be numerically imposed in the Rankine source method, while it is automatically satisfied in the Kelvin-Havelock source method. Also, the infinite free surface must be truncated in some way in the Rankine source

method. With respect to numerical calculations the Rankine source method will require the solution of a larger equation system, since singularities will have to be distributed on the free surface in addition to the body surface. The

Kelvin-Havelock source method will require the solution of a smaller equation system, at the expense of a significantly more difficult Green function, both mathematically and numerically.

The need for imposing the radiation condition and truncating the infinite free sur-face, represents two basic difficulties with the Rankine source method. In addition the discretization of the free surface into panels will introduce discretization errors in imposing the free surface condition.

The main advantage of the Rankine source method in addition to being

compu-tationally simple, is the possibility of including nonlinear terms in the free-surface boundary condition. However, the nonlinear terms are only of importance close to the body. This has motivated the partitioning of the free surface into two domains, an 'inner' domain where nonlinear terms in the free surface boundary condition are retained and where the Rankine source is applied, and an 'outer' domain where the Kelvin-Havelock source satisfying the linearized free surface boundary condition is

applied. The two solutions are then patched on a patching boundary partitioning

the two domains. This method has been used by e.g. Aanesland [1986].

The above also applies to the unsteady case, but the Kelvin-Havelock source is now replaced by the more complex oscillatory source. Interaction between the steady

4 Chapter 1. Introduction and unsteady flow may now also cause nonlinearities in the free surface boundary

condition. Away from the ship, the free surface boundary condition reduces to that satisfied by the oscillatory Green function. This motivates for a formulation

in terms of an inner and an outer domain also in the unsteady case, with the

Rankine source used in the inner domain.

Also the Green function approach will provide a reasonable approximation for many

practical purposes. In this case the Green function approach may be used only,

and there is no need for a distribution of Rankine sources on the free surface. This approach was used in the first fully three-dimensional formulation for the unsteady problem given by Chang [1977]. Also later formulations have used this simplified approach, Inglis and Price [198213], Guevel and Bough [1982].

Also the body boundary condition will now include interaction between the steady and the unsteady flow. Two formulations are possible. In the first formulation, the

steady motion problem is solved and used to form the body boundary condition for the unsteady motion problem. In the second formulation, the perturbation of the flow caused by the steady motion is neglected, and the unsteady problem

is solved directly. Inglis and Price [1982a,198213] have presented calculations and comparison of both formulations.

When the forward speed is low, the free-surface boundary condition may be sim-plified. Formulations for low forward speed has been given by e.g. Hermans and Brandsma [1988] for the wave resistance problem, and Zhao et al. [1988] for the unsteady case.

In the unsteady case, the application of the Rankine source method will introduce

discretization errors in imposing the free surface boundary condition, since the waves will now have to satisfy a 'discrete' dispersion relation. Such a discrete

dispersion relation has been studied by Nakos and Sclavounos [1988]. They have

also studied the effect of this discrete dispersion relation on the stability of the

panel method.

An additional difficulty in the forward speed case to that of the zero speed case with

a surface piercing body, is that the discretization of the hull is not only dictated by the gradients of the body geometry, but also by the shortest characteristic wavelength of the wave system associated with the wave source potential. This implies that very small panels may be required in the modelling of the hull, thus

increasing the computing time significantly.

An alternative approach is performing the calculations in the time-domain. The

time-domain approach offers possible advantages because the Green function ap-plied in the time-domain will retain its relative simple form regardless of the forward

5

the Green function requires the calculation of a convolution integral. Examples

of seakeeping calculations in the time-domain have been given by eg. King et aL [1988].

With respect to seakeeping calculations in the frequency domain, the advantages of the Green function approach over those of the Rankine source approach cannot be exploited in practical calculations, unless the Green function and its derivatives can

be evaluated numerically efficient and to sufficient accuracy. That is, the task of

evaluating the Green function and its derivatives must be reduced to a 'black box'

routine, which will return the calculated values as efficient as possible given the

forward speed and oscillation frequency, in addition to the field- and source-point coordinates.

The forming of such a 'black box' routine for the unsteady case, is the task of this work.

Such a 'black box' routine may then either be used directly in a computer program for three-dimensional ship motion calculations, or used in the formation of an

ap-proximate routine utilizing interpolation or series approximation. Such schemes have been given for the Kelvin-Havelock source by Newman [1987] in terms of

Chebychev polynomials, and Telste and Noblesse [1988] in terms of table interpo-lation.

6 Chapter 2. The hydrodynamic boundary value problem

v21 = 0,

(2.2)

2

The hydrodynamic boundary value problem

We will here state the hydrodynamic boundary value problem for a vessel with

forward speed in incident waves.

2.1

Coordinate systems

We will define three Cartesian coordinate systems, with la = (X0, Yo, Z0) fixed in

space, ./1-e- = (X, Y, Z) fixed with respect to the ship, and i = (x,y, z) moving in

steady translation with the mean forward velocity of the ship.

The space fixed system, Y0, provides the coordinate system for which the simplest

formulations of the free surface boundary conditions may be given. Whereas the

vessel-fixed coordinate system, X, provides the system in which the body boundary

conditions on the vessel's wetted surface may be derived in the most convenient

way. The steady-moving coordinate system, , provides an inertial reference frame in which the motions of the vessel are periodic in time.

We will take Zo = 0 as the plane of the undisturbed free surface with the Z0-axis

pointing vertically upwards. The X0-axis is taken positive in the direction of the

forward velocity of the vessel.

The three coordinate systems are related by the following conventions.

The vessel-fixed coordinate system is defined such that X =

in steady state

equilibrium. The steady-moving coordinate system is related to the space-fixed

system by the transformation

(X0-ut,Y0,z0),

(2.1) with U the mean forward velocity of the vessel.

2.2

Exact formulation

2.2.1

Condition in the fluid

The fluid is assumed ideal and incompressible, with constant density p. 'notational fluid is assumed in the whole fluid domain. Surface tension effects are neglected. The fluid velocity vector V(X0, t) may thus be defined as the positive gradient of a time-dependent velocity potential (I)(X.0, t) governed by the Laplace equation

2.2. Exact formulation 7 throughout the fluid domain V.

The fluid pressure p(X0, t) is given by the Bernoulli equation

p= p(

-+ V2 + 9Z0)+ pa. (2.3)

Here pc, is the atmospheric pressure which is assumed constant.

2.2.2

Body boundary condition

On the submerged part of the vessel's surface, 8, the normal velocity of the hull is

equal to that of the adjacent fluid. The appropriate boundary condition is thus

( 17, ) 0 on S. (2.4) where V the local velocity of the vessel's wetted surface. The unit normal vector if. is defined to point out into the fluid domain.

2.2.3

Free-surface conditions

The free surface is defined by its elevation Zo = C(Xo, Yo, 0. The kinematic boundary condition is expressed in terms of the substantial derivative DIDt

aiat+v V, in the form

Dt((

= 0 on Zo = C, (2.5)

or, equivalently

a(

OF a(

84, ac

al)

ax0 -o

+-to azo

Since the position of the free surface is unknown, an additional dynamic condition

has to be applied, namely that the pressure on the free surface is atmospheric.

From the Bernoulli equation it follows that

Ft IN74N2 gC =0 on Zo C (2.7)

The two conditions (2.6) and (2.7), may be combined into one in terms of F. This gives an alternative boundary condition for the velocity potential

(ku

Oz,

2V<I Vc1:t 4V.:11. V(V4?- V(1.) -= 0 _{on Zo = C (2.8)}

The velocity potential due to the steady forward motion of the vessel will be

denoted by

4)(10, t)= (2.9)

on Zo = C.. (2.6)

4(1).

Chapter 2. The hydrodynamic boundary value problem

The velocity vector of the steady flow relative to the moving reference frame is

= V(4) Uz).

(2.10)

Following Newman [19781, we will write the total velocity potential in the form

t)= 0(1,0 = 0(1)d-

( t ) (2.11)

Introducing (2.11) into (2.8) and neglecting second-order terms in y), we get the

free surface boundary condition in the moving reference system V(W2) gik + coft + 21k. Vyo,

+VV. V(ITV. Vw) 4V9:7- V(W2)+

=0

on z = (.

(2.12) The nonlinear free-surface condition for the steady potential in the moving refer-ence system is then

W. V(W2) .gbr = 0 on z = (2.13)

with the steady free-surface elevation given by the dynamic boundary condition

(2.7), as

[Uz

9 z=c

1 r

--29 [W2 U2]z=c_

In the boundary condition for the unsteady potential it is implicit that the steady

velocity potential q5 and the steady velocity field W are known.

The use of a nonlinear boundary condition requires the utilization of an iterative

solution procedure. Furthermore, the boundary condition (2.13) requires extreme computing capacity. Thus, the free-surface boundary condition for the steady po-tential will have to be linearized in some sense to facilitate a practically obtainable solution.

The forward speed problem may be linearized by restricting the body geometry so

that perturbations of the steady flow due to the presence of the vessel are small and may be neglected. This implies that the vessel is either thin or flat.

The free-surface boundary condition can be further simplified by assuming that the frequency of oscillation is high and that the forward speed is low. This removes the explicit speed dependence from the free-surface boundary condition and enables

a solution to be found in terms of a speed-independent velocity potential. This approach has formed the basis for many strip-theory formulations.

(2.14) 8

+

2.3. Decomposition of the unsteady potential 9

2.2.4

Additional conditions

Equations (2.4) and (2.8) are the principal boundary conditions of the problem.

Additional conditions must be imposed, such as that the energy flux of the waves

associated by the disturbance of the vessel is transmitted away from the vessel at infinity. This condition is termed the radiation condition of the problem. In addition a condition of zero flow through the sea floor must be imposed. In the

case of infinite water depth this condition may be formulated as

V-0

8,8 Zo

CC

(2.15)

2.3

Decomposition of the unsteady potential

The unsteady motions will be assumed small so that the unsteady potential can

be linearly decomposed into separate components due to the incident wave, the scattered wave due to the disturbance of the incident flow, and each of the six

rigid body motions. The rigid body motions of the ship are denoted by

= [711, 2, 713iTe" (2.16)

= [ni, n2, Q317. eiwt

.E [n4,5,716]T eat (2.17)

Here 771,772,77, are the surge, sway and heave motions, and 774,77506 are the roll, pitch and yaw motions, respectively. Thus, the unsteady velocity potential can be expressed as

{

,,,p = CI + W7 +

E

,,,,,iii

eit.

3.1

Here (pa is the incident wave potential in the steady moving reference system. In the fixed coordinate system the incident wave potential is given by

4)0(X0, Yo, ZO, t) = .5fCG evZ° ei(vX° cos r3+4'1° sin '3+6j°1)

Wo

In the steady moving coordinate system the incident wave potential takes the form t) _9(a euzes(ux coe0+Lovsini3+wi)

(2.20) with = 0 corresponding to head sea, and

LI

woz

= .

The frequency of encounter w is defined by

ww0(1 WOU

6 (2.18) (2.19) (2.21) (2.22) +

Chapter 2. The hydrodynamic boundary value problem

The components so,,j = 1...6, are radiation potentials due to the motions of

the ship in each of the six modes of motion and çc is the scattered or diffraction

3

The Green function

Far away from the vessel the perturbation of the steady flow due to the ship may be neglected. In this case no information of the steady potential is required, and the mathematical model simplifies considerably.

Thus, IT/ = III, and the free surface condition reduces to

wit

2U1 +

+ g(pz = 0. (3.1)

Together with the Laplace equation, the body boundary condition and the

radia-tion condiradia-tion, the equaradia-tion (3.1) forms a set of boundary value problems for the

unknown potentials (pi = 1 ... 7.

The source distribution technique will be applied in the solution of these boundary value problems. The application of the source distribution technique to the solution of boundary value problems in marine hydrodynamics have been described by many authors, e.g. Faltinsen and Michelsen [1975] and Garrison [1978].

It has been shown by Brard [1972], that when applying a singularity distribution method to a surface piercing body with forward speed, a line integral contribution must be included in the fundamental integral equation. Thus, the velocity potential

at any point (x, y, z) in the fluid is given by

Cb(x,Y, z) = f fs er(, 71, C)C(x, 11,z1 ,77,C)cis

9 CF C)G(., y, zI,77, On. 4, (3.2)

where the contour CF is the intersection of the body surface and the mean position of the free surface. o is the source density distribution on the mean position of the body surface, and G is the appropriate Green function.

We will take the time dependence of the periodic motion at frequency ca as et". For infinite water depth, the Green function may be found as, Wehausen and Laitone [1960]. 1 1

2gIrr°3ti

G(X,Y? Zle, 77,

=

+ [ [

f

, lc) dlc dB, (3.3) R1 irLr a Jo J-1 JCL JL2 where f (9,k) = w (0 k)'

he k)

11 (3.4)

Chapter 3. The Green function

with

and

R2.= 02 + (y y)2.+ Of% (3c8),

=.(z .01+ (y

+ (z (3.9)1

Thus, 1/R is the Rankine source located at ,(E0,C)) and 1/R1 is its image about

the mean free surface..

Further there are two singular points in the second integral of equation (3.3) given by k = ki and k = k2, and two singular points in the third integral given by k = le, and Is = ka.

The integration path in the k-range is from 0 to co. At the singular points, the

integration path is defined as ki + ie, k2 ie, Is, + ie and k.v + de as c The singularities on the contours are given by

h(0,k)= k cos [k(y n)sin 91 eacp{k((z + () i(z e) cos 011, .(3.5)1

W(tick)= gk (w +Uk cos.9)Ft, (3.6)1

if r = caLTIg < 1/4

--1 arccos(*)

if

r =

1 Ig 1/4.

fill 1

11 4r cos 9

Via; j

w '27- cos 9 I Nrik; [11 + N/1 4r cos 91 2r cos 9 X3.10) ;(311)

The fundamental singularity G(x,y, zIE, 77, 0. defined in (3.3) is a function of the

parameters r and (21g. The behaviour of this, function is similar to the Kelvin

singularity, but produces four free, waves.

The nature of the function changes when r = 1/4. At this value some of the terms In the function becomes infinite, and a solution is not obtainable.

In the limiting case when r 0 0 while w2/g remains finite, G approaches the

fundamental singularity for body oscillations at zero forward speed, with waves, radiated from the singularity in all directions.. When r 0 0 while U21 g remains finite, G approaches the Kelvin singularity fOr steady motion.

12

+0.

3.1. The dispersion relation 13

3.1

The dispersion relation

The wavenumber vector ic = (kr, Icy) identifies a wave component with the

fre-quency w provided it points to a curve in the (kr, ki,) plane defined by the dispersion relation

W(9', k) = 0. (3.12)

In the figures 3.1 and 3.2 are plotted the solutions of the dispersion .relation as

obtained by (3.10) and (3.11) in the polar coordinates (fik I co2, 0).

A more informative picture may be obtained by expressing the disperSion relAtion in terms of the components of the wave number vector, kr, ku.

Introducing

k = Vk + kr!, (3.13)

And

hr = k cc:!9, (3.14)1

inW(6, k)1 now gives,

gikx24-k12,--w+Uktr

= 01. ₀₁₅₎₁

For a given ca the dispersion relation defines three or two curves depending on the parameter r being less or greater than the critical value 7-, = 1/4, respectively.

The curves representing the solution of the dispersion relation may be written in terms of r as

=

(r2+ 2rcr +icr2)3

_(3.0)

with

(3.17)

T3.18)

Figure 3.3 shows the curves obtained for r = 0.0, 0.2, 0.5, 1.0, 1.5. Figure 3.4

shows the behaviour of the curves for r = 0.1, 0.2, 0.25,10.5, DO in the vicinity of' the origin. Each curve corresponds to a wave system. Due to the symmetry of the

flow, the curves are symmetric with respect to the kr-axis. At r =0.25 the closed

loop around the origin meets the curve which lies in the positive k2-plane, and the two wave systems merge for values of r above re.

In the limit as w approaches zero while U remains finite, i.e., the steady case, the dispersion relation takes the form

= kr kr (1.19)

ler U2

U2 9

14 Chapter 3. The Green function

10

I ' I

20 40 GO 80 100 120 140 160 180

Figure 3.1:

The solution curves of the dispersion relation for T = 0.1,0.2,

0.5,1.0,1.5 respectively in the polar coordinates gk/w2,0. The curves corresponds to the wavenumbers hi, k3 in (3.10).

r = 1.5

11,11111

0 20 40 GO ao icio 120 140 160 180

0

Figure 3.2:

The solution curves of the dispersion relation for r = 0.1,0.2,

0.5,1.0,1.5 respectively in the polar coordinates gk/w2,9. The curves corresponds to the wavenumbers k2, k4 in (3.11). 100 80 - = 0.1 r -= 0.2 Sc -r = 0.2 40 -r = 0.5 = 0.5 20 -2 0 0

3.1. The dispersion relation 15

i 1I

I

-5

-4

-3

-

-1 0 1 2 _ 3

k ,

Figure 3.3:

_{The solution curves of the dispersion relation for r = 0.0, 0.2,}

0.5,1.0,1.5 respectively. 10

r

1.5

8

-= 0.0

6

4

2

-= 0.0 = 1.5

16 Chapter 3. The Green function

0.10

-0.5 -0.4 -0.3 -0.2 -0.1 -0.0

0 . 1

0.2

0.3

0.4

0.5

,

Figure 3.4:

The solution curves of the dispersion relation for T = 0.0, 0.1,

0.2, 0.5,1.0, 1.5 respectively.

0 .08

0 .06

0 .04

0 .02

-r = 1.0 _{= 0.5} = r =-0.25 0.2 = 0.25 0 .00 I

' II

I I ' I I ' I

3.2. Formulation as a Cauchy principal value integral 17

Thus the closed curve around the origin disappears, while the two open curves

becomes symmetrical with respect to the k,, axis.

In the limiting case when U approaches zero while co remains finite, the dispersion relation will take the form

wz

(3.20) g

Thus, the open curves will disappear and only the closed curve around the origin will remain in the form of a circle.

The group-velocity vector of a wave component at some point on the curves is

defined by, Lighthill [1967]

awiT

VW(ki, ky; w).

g 101sx' Oki

awtk5;0

This vector is by definition normal to the curve at the specific point. Thus the wave-systems corresponding to the closed curves for r < r covers the entire free

surface, while the wave-systems corresponding to the open curves are confined to

a wedge behind the source. The half angle of this wedge is given by the angle between the normal to the corresponding curve at the point of inflection and the

ki axis.

With forward speed and r < 1/4, there are four wavetrains associated with the singularities ki,k2,k3 and k4. Three of the waves are traveling behind the source and one in advance. For T > 1/4 the upstream wavetrain disappears in front of

the source, leaving only wave trains to the rear.

3.2

_{Formulation as a Cauchy principal value integral}

The Green function defined by (3.3) is difficult to evaluate numerically because of

the contour integrals. In order to obtain a numerical solution, we must therefore

recast the contour integrals as Cauchy principal value integrals.

Using that in the range 7 < 6 < r/2 we have

gk (cv + Uk cos 0)2 = U2 cos2 0(k k1)(k k2), (3.22)

and that in the range 7r/2 < 9 < r we have

gk (co +Uk cos0)2 U2 cos29(k 1c3)(k k.,), (3.23) we may now rewrite the contour integrals in the form

Irt

= f (0,k) dk de g. [141-c f k2-E = J'Y 0

+f

f +f

f,(0,k)dk dB, (MO

(3.21)

-=

18 Chapter 3., The Green function and,

r

f

k)dk

ft,

=

-21 r fk'-` +f

₊

f3/4-c+

+j

f2(9,k).dk,d0. (2.25) tr J4 JO kt kt+t ict t+t

k) denotes integration around he singularity k; and

h(9,k)

fr(9 k) U2 cos' 9(k _hock k2yr%

h(9, k)

=

In the limit as e and lc, reduces to Cauchy principal value integrals plus

the contributions from the four singularities k1, k2, 1e3 and Ira.

That is,

r

k)dk de

0

re,

U2 cos28(lc

ks)(k 14)

2ig h(9 , lei) li(9,k2)d9, , U2 cos2G(k2

=

fr fcc he 1, k) dk do'

0

'

2.

iv

h(9,k3)i

h(0,k4)de. U2.cos2 9(k4 k,) With the use of the relationships

k2 It1

g.Jl

4r cos 0'

for 7 < 9 < r/2

kt

ka j

U2 cos2 9

for r/2 <O

70,

(3:28)

(3.29)

0.307

'The derivatiVes of the Green function may be obtained by differentiating the terms

In (3.31) with respect to either x,y or z.

rrn

the Green function takes the

GkYrziC111.0 = 1 1 22. TZ R1 7r k form

fr

+

[ J ) + h(9, k2)1 MB, k)dk de (3.31) gk (to + Uk cos 9)2 de r [h(9,10 _ MO, k41 49. J-T

VI 4r cos 9

4r cos (3.26) (3:27) =

_-

-f

, - * 0, = -y +

--

-[h(9,

-3.3. Alternative formulations of the Green function

However, the formulation of the Green function given by (3.31) is not satisfactory

with respect to numerical evaluation. In each of the double integrals, there are two singularities in the k-range. With 8 = ir/2 one of these will be at infinity. In the first of the double and single integrals, the integrand will become infinite for

0 -y.

The singularities in the integration range will preclude numerical integration by

the use of ordinary quadrature methods.

The integration of a function which has a simple singularity in the range of

inte-gration, may be handled by a method developed by Monacella [1966]. However, Monacella's method is not applicable in the case of an end-point singularity. More sophisticated numerical integration schemes capable of handling internal as

well as endpoint singularities of several types, are available through subroutine

libraries.

3.3

Alternative formulations of the Green function

The form of the Green function given by (3.31) was applied by Chang [1977].

However, even though the singularities in the integration range may be handled

numerically, the direct computation of the double integral is very timeconsuming,

since the oscillatory behaviour of the integrand will require small step length in the numerical integration and thereby a very large number of evaluations of the

integrand.

Alternative and more efficient formulations are therefore required. Since the

ex-ponential function is present in the integrand, a natural approach is to seek to reformulate the double integral in terms of a 'single' integral involving the expo-nential integral. Such formulations have been given by Inglis and Price [1980],

Guevel and Bougis [1982], and Wu and Eatock Taylor [1987].

Inglis and Price [19801 replaced the double integral terms in the principal-value formulation for 8 > 7 by single integrals involving the exponential integral. The

principal-value integrals were evaluated along a ray in the complex plane so that

the oscillatory behaviour of the integrands vanished. Thus arriving at the following form of the Green function

G(z,y, = 1 1 2g

r

1, gk + {Qi(w+ ) + Qi(w- ) +f w{(23(w+)+ Q3(w-) h(0,k)dk dB (w + Uk cos 9)2 Q 2(w+) Q2(w--)} d0

Q4(w)

Q 4(w-)1dB

20 Chapter 3. The Green function where Go(M, Ai' ;t)

+2i[h(9 ,

ki) + h(9,k2)1c19 fr - 4r cos 9 [h(9, k3) h(9 ,k4)1d9

f

+2ii

N/1 4r cos (3.32)

h(9, k) k exp [k{(z + () + i(x 0 cos OH cos{k(y 71) sin 9}, (3.33) w+ -= (x

cos + (y

7i) sin , (3.34)

w_ (x cos (y n) sin 0, (3.35) .70.0±) {i (2 r exp(q,±) sgn (iv* ) irV1 4r cos 9 +exp(qi±)Et(q,±)}, (3.36) 43 + CI ITV). (3.37)

Here the time dependence was taken as et

the motion of the source is in

the negative X0-direction.

Guevel and Bougie [19821 also used the modified exponential integral to reduce the double integrals. This formulation involves the time-dependent terms e"" and

They arrived at the following formulation

M' ;t) = Go(M, Mt; t) + C1 (At, Mt; t) + G2(M, Aft; t), (3.38)

with M and Mt being the field point and source point, respectively, and

[1

1mmii ImArti1 cos wt, (3.39)

Gi(M, Mt; t)

---

1 Re fe,,, fir IKI[C1(K,0 + G1(1 f le)]

ri

1 Jo l Vi + 4r cos 0

K2{G1(K2e) + G1(K2e)]1 do} _(3A0)

8/1 -1- 4r cos 8 J G2(M , M'; t) =

Re { '

ti r i ( K3[G2(1C3) + G2(K3e)] (3.41) 7r.C. - 110 1/4-r cos 0 1 K4[G2(K.) + G2(1Cieli)de V4r cos 9 1 + r°' + G1(K3e)] - K41G3( K4) + G3(K4e)] de in V47 cos 8 1 K4G1(K4E) + Gl(Kcie)Id81} fo! K3[G3(K3e) + G3(K3e)]

4r cos

Here, the coordinates are nondimensionalized by the characteristic length 1, and

thus

X =

Y = Or Z = z11,

(3.42)

-=

-

-8-

--kiaz

-and eK3[G1(Ka0

-G(M,

-3.3. Alternative formulations of the Green function 21

with

K = Icl, F=

Ii

_{r = iDF,}

yg

= (X

X') cos 0 + (Y Y') sin 9, (3.45)

= (X

X') cos 0 (Y Y') sin 9, (3.46)

(3.47) with X, Y, Z being the coordinates of the field point and X', Y', Z' the coordinates of the source point.

Further

Eice) = Eve),

inn(e) >

E1() = E1( + ie),

Im (e) = 0,

El(e) = E1 (e) 27ri, Im (e) < 0.

G1, G2 and G3 represent the following modified exponential integral functions

G1() = J.Ei(e),

0 < arg(e) < 27r, G2(e)

= JE1(e),

7r < arg(e) < 7r,

G3() = ee[Ei(e)

27rij, 0 < arg(e) < 27r,

1 cos 0, cos 74;., 1 (3.48) (3.49) K Z + Z' + if2, Z + Z' + , (3.50) (3.51) (3.52) 1 + 27 cos 0 .N./1 + 47 cos 2F7? cos2 1 + 2r cos 6 + -11 + 47 cos 0 K2 (3.53) 2 F,? cos2 1 27 cos 47 cos 9

K3 =

if T < 1/4, (3.54) 2F,,? cos2 0

K3=

1 27 cos

i-V4r cos 6 1

if r >1/4,

(3.55)

2F,? cos2 1 2r cos 0 + 4r cos 61 K4 K4

lit < 1/4,

if (3.56) 2F,? cos2 9

1 27 cos

+ i1/47 cos 0

= 2F,? cos2 T > 1/4. (3.57) (3.43) (3.44) (3.58) (3.59) (3.60) (3.61) (3.62) (3.63)

=

=

2r 9 9 9 9 9 9

22 Chapter 3. The Green function Wu and Eatock Taylor [1987], also introduced the complex exponential integral to reduce the double integral. They wrote the Green function as the following sum of integrals 1 1 R +(101 + 102) + (II + hi) (121 1.22) + (133 + 134) (143 + 144)7 with 1 fi klexP(k1X1) .101 = o i.V47 cos 0 1 Ep(kixi) dB 1 r ic2exP(k2X1) r Jo W47 cos 0 1Eq(k2x1) d9,

1 1 ki exp(klx2)

/az =

E (kix2)d9 0 i-V4r cos 8 P 1 r k2eXP(k2X2)

r Jo iV4r cos

9 1Eq(k2x2)c10 ,

1 p k, exp(kixi)

=

1 En (k,x,)dB, Ir 7' 4cer cos 0

where ii = ,f, a = 1 for i = 1,2, and 7' = 0, a = 1 for i = 3,4.

Further

(Z)

Ea)

= El(Z) when Re (Z) > 0,

(Z) E1(Z) when Re (Z) <0 and Im (Z) > 0, (Z) (Z)

2ri

when Re (Z) < 0 and Im (Z) <

Eq(Z)

Ei(Z)

2ri

when Re (Z) < 0 and Im (Z) > (Z) E1(Z)

when Re (Z) <0 and Im (Z) <0,

Er(Z)

E1(Z) when Im (Z) > 0,i = 1,3,4,

(Z) = E1(Z)

2ri

when Im (Z) < 0, i = 1, 3, 4,

Er,(Z) E1(Z) + 22ri when Im (Z) > 0,i = 2,

Er(Z)

El(Z)

when Im (Z) < 0, i = 2.

(3.64) (3.65) (3.68) (3.69) 0, (3.70) 0, (3.71) (3.72) (3.73) (3:74) (3.75) (3.76) 1 ± i-V42- cos9 1 2r cos ca. 9k2, gki = (3.77)

ki,i = 1,2,3,4 in I,, are given by

1 2r cos 0 ± -V1 4r cos9 v, k2, kt = _{2r2 cos2 0} (3.78) , 1 + 27 cos9± 1/1 + 4r cos 0 k4, x3 v, (3.79) 272 cos2 0 1c1 and k2 in 101 and 102 are given by

(3.66) (3.67)

--

-=

=

=

=

-+

=

+

-A formulation in terms of a single integral has also been given by Bessho [1977]. This is a genuine single integral expression in the sense that only elementary func-tions are involved in the expressions for the integrand. However, this formulation is more mathematically involved than the formulations above. Also, the integration will have to be performed in the complex plane.

He wrote the Green function as

G(P,Q) =

1

(1

-

)

1 4

r

+

gT(z -

y

-with

T(z,y,z)

1 at (3.85) 2rig + cosech t

x {exp[iii{x+ psinh(t + ib)}lp.(t)- exp[itilx + psinh(t +

r- cosh-1 (Ca g) _dt

+2rig

_-

_{cosech t}

x

[e[-iv{x

- psinh(t + ib)}iv(t)+ expr-iv'{x - psinh(t + 2:5)}1vV)],

where P and Q are the field point and source points, respectively, and R = PQ, = PQ where Q is an image point of Q with regard to the free surface.

Further (3.84) sinh-1 zlp, (3.86) ii(t) 1 g (., 4w ± 2 I 1 + cosech t cosh t, (3.87)

iil(t) f

_g for 0 < p < ca or Li < p,

111(t) 1 _-= -co + g_v(t) (1

±

1/1

-

cosech

t)

cosh t, (3.88)

2 _g

for p < -,1

or

- w < p < 0,

-p sin 8, (3.89)

(3.90) p cos 6.

3.3. Alternative formulations of the Green function where v w2/ g; and

Xs = (z

C)

i[-(x -

cos 9 + (y - sin 0], (3.80)

X2 =

(z C)

1{(n

e) cos 9- (y - 77) sin 8], (3.81)

X3 =

(3.82) X4 kl (3.83)

no =

p =

P =

z =

y =

=

+ +

+

+

+

1 ,

24 Chapter 3. The Green function

Ohkusu and Iwashita [1989], have developed a numerical evaluation scheme based

on Bessho's formulation, utilizing a numerical implementation of the method of

steepest descent.

Demanche [1981] applied an integral transform on the equations in the boundary value problem. This transformed a function f (r, Si into a double sequence fr1,11fn,2 9 with

LAO =

f.,2(u) =

p00 (2c

(r, 9) cos ne de Jn(ur)r dr n = 0 ... oo, (3.91)

71 0 0

f (r,O) sin nO dB J(ur)r dr n = 1 ... oo. (3.92)

1 aD r2r

7 )0 j1'0

He then obtained the following form of the Green function

1 1

G=

Vr2 + (2 N/r2 + (2 +

oz

+4i,5

z,

Z2 Jo(ur)

U Jo

z? + 1

4

+ 1)

2 du

ct (27+1

Z' t+1 1

+

J(ur) cos n91eu(z+C)

.Z? +

where Z1 and Z, are explicit functions of u through a and #, defined by

cv

a .

l (3.94) Uu 1 0 = r§ V' VIt:, (3.95) Z1

I Nil

(a + 0)2 - i(a + 0)

ifa-I-0< 1

1

1 i ki(a

+ 0)2 1 (a +0)] if a + p > 1 (3.96) V1 fl) _ifla 01 < 1 Z2 = { ikl(a 13)2 1 (a 0)] if a

S> 1

(3.97) [V(a

+ (a

if a/9<-1

r2 (x e)2 + (y 77)2. _(3.98)

A direct comparison of the efficiency and performance of the different formulations is difficult to carry out from the analytic formulations. However, some comments may be given.

The formulation of Inglis and Price [1980] will still require the evaluation of a double integral when T > 0.25. This integration is awkward to carry out numerically since the integration in the k-range extends to infinity.

(3.93)

f

3.4. Nondimensionalization of the Green function 25 The formulation of Wu and Eatock Taylor [1987] is identical to that of Guevel and Bougis [1982] when r < 0.25 if the complex conjugate of the term involving e-twi is taken in the formulation of Guevel and Bougis. The formulation of Guevel and Bougie requires subdividing the interval [0,7] into two subintervals [0, arccos(1/2r)] and [arccos(1/2r),7], with different formulations in the two domains, when r> 0.5. The formulation of Ohkusu and Iwashita [1989] requires integration in the complex plane. This approach appears more cumbersome, since one has to locate the curves

of steepest descent and then perform the numerical integration along this curve. The integration has to be performed beyond a point at which the contribution from the remaining part of the contour is less than a specified error limit. In this computation the contour must not cross the branch cuts of the integrand. This may require integrating along a curve of ascent rather than descent, in the first stages of the integration. The benefit of this formulation is that the integration

when both the field and source point are close to the free surface can still be carried out effectively.

The performance of the formulation of Demanche [1981] is difficult to access, since no results are given.

3.4

Nondimensionalization of the Green function

Before we proceed, we will introduce nondimensional quantities in the expression for the Green function.

We will introduce the dimensionless quantities

K =U2k1g

= Ucalg a =U2G1g

= gz1112 = gyIU2

=gz/U2

6

=g/U2

ij = gnIU2

=g/U2

or with

Ko = U2'

K = klKo

=

VgKo a = GlKo

= Koz = Koz

= = Kan = Ko(

Now we will have

R2 (U2)2

[(i

t.)2 (g ,)2 = 2 .02 )2] (3.99) (3.102) (3.103) + + + (3.100) (3.101j

26 Chapter 3. The Green function

Where

K cos[K(g fl sin Olexp K[(i + t) cos Oil K)

(Cl + K cos O)] The singularities in the integration range are now given by

and

VT}

47 cos Si

\TR; 27 cos

N/R; [1 + 47 cos 81

27 cos 8

j

Li-

Wir cos 61 11

2r cos 8

/Tc

i

nir.

iv(47 cos 1}

27 cos

43.106)

43.107)

(3.108)

3.5

Expressing the Green function in terms of the

expo-nential integral

We will now express the Green function in term, of the complex exponential in-tegral. The following formulation is essentially as that given by Wu and Eatock

Taylor [19971.

3.5.1

General case

First, we decompose the integrand F(0, K) in terms id its singularities.

With 9 < 14, the integrand will have no singularities. However, the denominator

will have two complex roots given by

(3.109)

(3a10) Noting that the Jacobian of the transformation is

Ok 89

on ,Ko

3.104), 89 89 we will have 1 1.

2 [fl

+:

+

f 4- .1 I

I F(9K)KolcIK de, 3.105) Ti

It""

/4,

3 C2: = + 9

J=

= G(x, y , ,77,C)

i(i

F(9, [K 8

3.5. Expressing the Green function in terms of the exponential integral 27 Thus we may decompose the denominator in terms of its complex roots, and con-sequently

K (f2 + K cos 9)2 = cos2 9(K Kic)(K K2c) (3.111)

when 0 < 8 < 7.

Along the integration paths Li and £2, we will have

K (ft + K cos 9)2 = cos? 9(K K 1)(K _K2),

4r cos

K

=

cos2 when -y < 0 < r /2, and

K (II K cos 9)2 =- cos2 9(K K3)(K _Ks),

i-V4r cos 1

K1c K2e (3.112) cos2 9 47- cos K3 K4 = COS2 when 7r/2 < 9 < 7r. We may then write

1 r K1c Ku

F(9, K)

i./4r cos

1 LK K1

K K2,

X exp{K[(i + (7) 0 cos 9])- cos[K( 77) sin (3.117)

for 0 < < 7; and

1 r IC,

F(8, K) K2 1

N/1 47 COS 9 LK K1 _{K K2i}

x exp{K[(i + 0

i(i" 0 cos 9]) cosIK( ij) sin 6] (3.118)

for 7 < 9 < 2r/2; and

F(0, K) 1

1K3

K4

47 COS 9 LK K3 K

x exp{K[(i + i(i" -0 cos 911 cos[K( fj)sint9] (3.119) for 7r/2 < <7r.

The last cosine term may be expressed as

1

cos[K( fj) sin 9] = -2- fexp[iK(g ij) sin 8] + exp[iK(i ij) sin 9]}. (3.120) (3.113) (3.114) (3.115) (3.116) 9 = 9 9 V1

28 Chapter 3. The Green function

Thus we may write

exp {K((i" + 0 cos eq cos[K(5 71) sin 0]

[expfK[ci+

+ i[(1 0 cos

+ sin Bill

+ exp{KRi +

+ i[(i

-e) cos 71) sin 901. (3.121)

Consequently, we may write F(0, K) as

KLc

F(9, K) =

i(K

Kic)1/4r cos 6' 1

x1.2. [exp +

+ i[(1

cos 9 + (g 61}

+ exp{K. ₊

_{+ =Hi}

6 (if sin01]}]

K2,

i(K K2c)V4r cos 9 1

x

{explif{(i

-e) cos + ( sin 9 )]}

+

expf K [(i+

+ i[(1

) cos (g f7) sin 8]] (3.122)

for 0 < 9 < 7; and

F(9, K)=

K

(K K1)1/1 4r cos 8

x [exp{K{Ci +

+ i((i

) cos 9 + (5 sin ell}

+

expIK

+

+ if(1

9

(g )sin

91111

K2

(K K2)V1 4r cos 8

x-2 {exp{K{(i+ + 0 cos 9 + (5 73) sin 81]}

+

expIK

+

+ i[(i

0 cos 71) sin01]}}

for 7 < 6 < /2; and

F(0,K) =

K3 (K K3),/1 4r cos 8 (3.123) = K

+

+ 9

1)

9

Ri

3.5. Expressing the Green function in terms of the exponential integral 29

x

[expini

+ + cos 61+ sin Oil

+ exp{K [(i +

+ i[(1

cos ell I] K4

for r/2 < 8< r.

We may then express the Green function as the following sum

1 1

G = +

+(z, + 112) (121 +122) + (133 + 134) (1.43 144)/ where

Ko p

Kic exp(Kx,)dK d8

r JO Jo

i(K Kic)V47 cos 8 1

Ko K2c exp(Kxl)dK de

r

i(K

K2)/4r cos

1'

Kor

K1c exP(KX2)dK de

r Jo Jo

i(K KicW4r cos e 1

Ko r

Kze exP(Kx2)dK de 7i JO JO i(K K2c)V4T COS 8 1' Ko K, exp(Kx,) c do

r

(K K1 )V1 47 cos 6'

Ko j1 ç

KiexP(Kx2)dK de r .1-r Jo _(K 47- cos K0 K2exP(Kxi)dK dB

r

JO (K

K2)/1

_{4r cos}

Ko fir r.

K2exp(Kx2)dK do 7r .1-r Jo (K K2)1/1 _{4.7- cos 6}

Ko r

K3exp(Kx3)dKdO

7r Ji Jo (K K3)1/1

4r cos Or

(K

4r cos e

x12- {explK {(1"+ + f) cos 6 + sin6d}

+exp{K[(i

+ + -0 cos 9

(§ )sin61]

11 (3.124)

Zn 102

Ill =

112 = ,r K0 I re K2 eXP(KX2)dIC dB 1-21 122 133 = (3.125) (3.126) (3.127) (3.128) (3.129) (3.130) (3.131) (3.132) 102

r

= = =

30 Chapter 3. The Green function

with

134 =

r

K3 eXP(KX4)dK de

1' fir' ..lo (IC Ka)v'1 4r cos O'

Ko r r K.e

xp(Kx3)dK

a

43 = KO

r

1

Jo K4 eXP(KX3)dK d9 (K K4)V1 47 cos 9' fir

r°

K4 eXP(KX4)dK de

r J

Jo (K

_KIWI

47 cos t9'

Xi =

+

+ it(i

-) cos + (g

)sin9],

X2 =

+ () +

[ e) cos 6' 771) sin 8],

XS = XI

X4 =

X2-1C1x,

f

It

Ko,

K1 exP(Kix7)do

_-du.

e'

vil 47 cos 0 J-00

e-U , R e-u _e-"

au =

f du + f du a2ri,

_Kix, IL JCR IL KiX1 e-U

au = Ei(Kixi)

a2ri.

f*C017 U (3.133) (3.134) (3.135) (3.136) (3.137) (3.138) (3.139)

In order to express the integrals in terms of the complex exponential integral we introduce the new integration variable u = x,(K

Ki),i,j = 1,...,4.

For the integral /II this yields

Ko f

r

Ki exP(KX0dK dB Tr (K K)NA 4r cos 6 (3.140) (3.141) (3.142) (3.143)

Further details on the integration in the complex u-plane are given in Appendix

A.

The last integral may be evaluated by applying Cauchy's integral theorem in the

form

In (Kai) < 0

a = 1 Im(KiXi) °

a

CR is given by an arch in the first quadrant.

In the limit as R co, the integral over CR will vanish. The first integral on the

right-hand side is recognized as the complex exponential integral, and we will have

U

=

3.5. Expressing the Green function in terms of the exponential integral 31 The other integrals may be treated in a similar way, and we will then have

= Ko_It f7 KicexP(KicX1)E,(Kicxi)de Jo iN/47 cos 8 1 Ko [7 K7cexP(K2ai )_{E, (K200) de}

It Jo W4r cos 0

1

102 =

Ko 17 KlceXP(KIcX2)Ep(Kicx2) dB it JO iN/47"COSB 1

Ko r K2

eXP(K2cX2)E,(K2cx2) de It Jo iN/47 cos 8 1 where

E(Z) = E,(Z)

El(Z) when Re Z >0; and

E(Z)=E1(Z)Im Z > 0

Ep(Z) = Ei(Z)

2ri

Im Z < 0

Ea) =

(Z) 27ri Im Z > 0

E,(Z) = El(Z)

Im Z < 0 when Re Z < O. Similarly, Ko K1 exP(KIXI)Eri(KIX1)d9, = it .17 \,/1 4r cos 6) Ko

j KIeXP(KIX2)E(K)d9,

= it - 4r cos 9 K0 K2 eXP(K2X1) Er2(K2X1) d8,

21 =

it v/1. 4r cos Kofi K2 exP(K2X2) ₍ 2k 2X2)

22

it J-7 4T COS 17 KO

_{K3exP(K3X3)E(K)dO}

7.f

₌

_r

; 4r cos 0 f'" K3 exp(K3x4) ( 1 ,Ta it 1/1 47 cos 9 A-,r3k41,3X4)%v, K0 fr K4 eXP(K4X3) Er4(K4X3) Elva, .4.43 = ir ; N/1 47 cos 9 K0

_{KlexP(K4X4)E(K)dO}

-'14 it _{- 47 cos} (3.146) (3.147) (3.148) (3.149) (3.150) (3.151) (3.152) (3.153) (3.154) (3.155) (3.144) (3.145)

x

= K0

-=

9

32 Chapter 3. The Green function where for j = 1,3,4; and for" = 2, and for -7 <8 < 7r/2, and and Er j(2) = El(Z) Im Z > 0 Er j(Z) E1(2) 2ri

Im Z <0

Er(Z) = E1(2) + 2ri

Im Z > 0

Er,(2) = E1(Z)

Im Z <0

K [1 2r cos 9 N/1 4, cos K2 2 cos2 K3 {1 2r cos 9 N/1. 4r cos 01 K4 J 2 cos219

for 7r/2 <9 < r.

It is also understood that the exponential function E1(Z) takes the value of E1(2+ i0) when Im Z = O.

By changing the integration range from [7r/2, 7r] to [0,7r/2] in the integrals /33 . 144

by the transformation 0 -= r 9, these integrals may be rewritten as (with instead of i)) K01; K3 exP(K3X3) L-, f Kr N Jn 1.33 = ir o VI + 4r cos 0.nr3v13X3) au, (3.160) Ko f i K3 exP(K3X4) Er 3(K3x4) de, 1.34 =

_r

0 1/1 + 4r cos 0 (3.161) Ko f i K4 exp(K4X3)Er4(K1)(3) a, 1,13 = ₇ 0 .V1 + 4r cos 8 (3.162) K0 ti K4 eXP(K4X4) Er4(1C4X4)d97 (3.163)

r Jo

V1 + 4r cos 0 with

K31

[1 + 27 cos 9 N/1 + 47- cos 81 K4 cos29 2 (3,156) (3.157) (3.158) (3.159) (3.164)

X3 =

COS 0 + 0], (3.165)

X4 =

i[(1 0 sin 91, (3.166) Or

X3 = X2,

(3.167)

X4 =

(3.168)

where the overbar denotes the complex conjugate.

The only point of introducing this transformation is to facilitate the simultaneous

evaluation of the integrals II, ... /44 in the range [7,

=

2

-'

+

-3.5. Expressing the Green function in terms of the exponential integral 33

3.5.2

Field point on the track

When the field point is on the track of the source, F(9, K) will simplify to

1 F(9, K) = KK1cK1c K2K, 2c] expfni }, \/1 47 cos 9 [ (3.169) for 0 <9 < -y;

F

1 r K1 _K2

(8, K) = ]

expi

ln

+ 0

- z ) cosMI,

V1 47 cos 9 1K K1

K K2

(3.170)

for y<9 < i; and

1 r K3 K4

F(19, K)

47 cos 8 LK K3

K K4

exp{ KIP

+

cos01}

(3.171)

for I <9 <

We may then write

1 1 G = 102 + 12 + 13 14, R R1 where 2K0

r Kic exp(Kx)dK de

101 =

_{r Jo Jo}

(3.173)

i(K

Kie)V47 cos 9

l'

2K0 fl r K2 eXP(KX)dK dg

102 =(3.174)

7r Jo Jo

i(K

_{K2c)V4r cos 8 1'}

with IC, and K2, given by (3.109) and (3.110), respectively.

2-i 14

2K0 0, ro

K, exP(Kx)dK Jo (IC K1 )V1 47 cos 9' 2K0

_{[XK2 eXP(KX)dK}

zr J.y Jo _{(K K2)V1 47 cos 8'} with Ky and K2 given by (3.158).

Again we may transform the integrals 13 and 14 to the interval [0, win

2K0 f

f.

K3exp(K)dK de

7r Jo Jo _{(K K3)/1 + 47 cos 9'} 2K0 lc°

KaexP(K)dK dB

r

Jo Jo _{(K K4)%/1-f- 4/- cos} (3.172) (3.175) (3.176) (3.177) (3.178) r

-

-

+

-

--

-

-VI-

--

-=

-

-Thus

With E(Z), Eq(Z) and Er,(Z),i = 1, ... 4 given by equations (3.146), (3.147),

(3.156) and (3.157).

3.5.3

The singularity at 0 = -y, general case

For r> 0.25 the singular behaviour of the integrand at 9 = 7 precludes the use of

ordinary integration methods for the integrals rot In as given by (3.144) (3.151).

By noting that the singularity at 9 = 7 is of the square root type, the singularity

may be removed by a substitution of the variable of integration. We introduce the change of variable

de =

u du

2711

)2'

= 0 = V4r 1 t9 = -y u = 0 (3.187) (3.188) Introducing the new integration

To, =

x= (i +) - i(i -

COS O. variable is

x(K

= 1, 2K0

r

_Klcexp(Klcx) Ep(IC icX) (3.179) , 4, now gives (3.180)

r Jo

i.V4r cos 6 1

/In =

2K0

r Jo

r

K2cexP(K2cX)E(K)

g2a

de, (3.181)

i.s/4r cos 0 1 =

21(0 0 Kt

exP(KiX)E(K)

ri

IX,

da (3.182) 7r

J

-1

In

4r cos 0

/2 =

2K0 0 K2 exP(K2X)Er2(K2X) a, (3.183)

XI =

ji71 /K13e7CP4(71(C°318) Er3(K3) d9' 2rK0 (3.184)

r 10

Vi + 4-r cos 0

14 =

2K0 ti K4 eXP(KI)

Er4(K 4)dB.. (3.185)

r io

V1 + 4r cos 0

34 Chapter 3. The Green function

with K3 and If 4, given by (3.164). X is here given by

U = V4r cos 0 1. (3.186)

=

i

-3.5. Expressing the Green function in terms of the exponential integral 35 Then the integrals /01 and /02 may be written as

with

Ko f

K exp(K lexi)E,(K ',xi) du

Jo

2TiN/1 (4,.1)2

Ku - r/7777 K2cexP(K2cX0E,(K2A1)du,

_

r Jo

_{2ri1/1 41-t:1)2}

Ko 1`11--1 ALIC eXP(K1cX2)Ep(K IcX2) du

7 JO

_(ti)2

Ko --2cIC eXP(K2cX2)Eq(K2cX2) du,

r Jo

_2r7V1 1(1 u2) iu 2 (3.191) (31,21-,±:1 )2 2 -1(1

u2)+iu

K2c = (3.192) 2 (1.1±1)217. i /1.2 + 1 2 + 1)2 (14

47-XI = (i +

+

[(i

+ ),-- (if n)\11 1 7 (3.193)

X2 =

± + i { .112+1 1 /142 + 1)21 . (3.194)

(i

)47

(i

in

\

(i 4T

Similarly, we may use the substitution

u =

4r cos 6' (3.195)

to avoid the singular behaviour in the integrals 111..-Z2 close to 9 = 7.

Thus u du dB _(3.196) ( 1 u2 1 _k

=

= 0 0 = u = 1 ' (3.197) (3.189) (3.190) = (2±1)2)2

(i