A three dimensional method for the calculation of the un-steady ship wave pattern using a Neumann-Kelvin approach

b

ATkee Dimensional Method for

the

_{lc_afion of the Unsteady}

Ship Way Pattern

using

a

Neumann

-

Kelvin Approach

Carlo van der Stoep

A Three Dimensional Method for

the Calculation of the Unsteady

Ship Wave Pattern using

a

Neumann

-

Kelvin Approach

TECHNISCHE UNIVERSITEIT Laboratorium voor I _{ScheepshydromechanlCa} Archiof Mekelweg 2, 2628 CD Deift Tel.: O15-786873-Fax 015-78183G

A Three Dimensional Method for

the Calculation of the Unsteady

Ship Wave Pattern using

a

Neumann

-

_{Kelvin Approach}

P roefsc h rift

ter verkrijging van de graad van doctor aan de Technische Universiteit Delft, op gezag van de Rector Magnificus, prof. drs. P.A. Schenck,

in het open baar te verdedigen ten overstaan van een commissie

aangewezen door het College van Dekanen

op 6 februari 1992 te 14.00 uur

door

Carlo van der Stoep

geboren te Rotterdam, Wiskundig Ingenieur.

Promotiecommissie:

Prof. dr. ir. A.W. Heemink Prof. dr. ir. A.J. Hermans Prof. dr. ir. J.J. Kaker Prof. dr. ir. G. Kuiper Prof. dr. H.G. Meijer Prof. dr. ir. J.A. Pinkster

Contents

Nomenclature 3

i

General Introduction 5

2

Mathematkal Formulation

11

2.1 Introduction 11

2.2 Linearization of the boundary conditions 14 2.2.1 The free surface condition 14 2.2.2 The body boundary condition 18 2.2.3 The radiation condition 19

2.3 The leading equations 20

2.3.1 The singularity distributions 20 2.3.2 Expansion of source strength and potential function. 26 2.3.3 Expansion of body boundary condition. 28

2.4 (Added) Resistance 30

3 Computational Aspects 31

3.1 Solution of the integral equation 37

3.2 Generating the grid 40

3.2.1 Hullform definition 40

3.2.2 Panel division 41

3.3 The Green function 41

3.3.1 Behaviour of the Green function 41 3.3.2 Calculation of the Green function. 44

3.3.3 Double integral 45

3.3.4 Single integral 47

3.4 Padé approximations. 50

3.4.2 Padé approximations in the single integral 54

4

Solution of the Matrix Equation

57

4.1 Introduction 57

4.2 The matrix equation 57

4.3 Convergence of the solution algorithm 59

4.4 Conclusions 64

5 Computational Results 65

5.1 Introduction 65

5.2 Steady motion 65

5.2.1 Wave profiles of Wigley hull 65 5.2.2 Comparison with measurements I 67

5.2.3 Convergence aspects 68

5.2.4 Comparison with measurements II 71

5.2.5 Resistance 73

5.3 Unsteady motion 78

5.3.1 Wave profile of Wigley hull 78 5.3.2 Comparison with measurements. 78

5.3.3 Dynamic swell-up 78

5.3.4 Added resistance 80

6 Concluding Remarks 85

A Expansion of the Free Surface Condition

87

B The Differential Operator

93

C Computer Program Outline

95

D Measurements: Wigley HuH

99

References 103

Acknowledgements 107

Summary 109

Samenvatting 113

Nomenclature

Coordinates:

x = (x,y,z)

Coordinate system moving with velocity U x Direction of the forward velocity

y Lateral direction

z Vertically upward

n

Unit vector normal to in outward direction

a

Motion vector

Physical properties:

B

Beam at midship

Total wetted surface of the body

L Length at water line

H

Draft at midship

V

Ship volume

g Gravitational acceleration

Non dimensional parameters:

E Principal small parameter

CB Ship block coefficient(Vs/BLH)

Fn

Froude number(U//TE)

Re Reynolds number(UL/v)

CD

Total resistance coefficient(D/pU2)

CF

Frictional resistance coefficient(F/pU2)

Cw Wave resistance coefficient(W/pU2E) Flow data:

U Ship speed

p Pressure

p Density

li

Viscosity

((x,y)

Elevation of the free surface

o(x,y)

Steady wave elevation (1(x, y) Unsteady wave elevation

Za Oscillator amplitude

SA relative motion composed of pitch, heave and

undisturbed wave amplitude

SA swell-up

F.S. Free surface

Steady wave potential Unsteady wave potential Total potential w Frequency of motion We Frequency of encounter Dipole strength a Source strength Mscellaneous:

C(x,)

Green function

T(x)

Chebychev polynomial of order n

Chapter 1

General Introduction

Fluid motions are basically governed by three kinds of forces: Inertial, Viscous and Gravitational. In order to estimate the actual full-scale flow from the scale model, the actual force have to be corrected by a suitable multiplicative factor.

The problem is that all three forces act simultaneously and they each have their own multiplicative factor. For instance, using the length scale L, the

three forces are of order of magnitude L2, L and L3 respectively. To predict

the forces of the full-scale model from tests with the small-scale model, it is

helpful to form the ratios of the three forces. Because any two of these ratios are sufficient to define the third, only the well-known Froude number, the ratio of the inertial to the gravitational force and the Reynolds number, the ratio of the inertial to the viscous force will be used (Newman [28]).

At the end of the last century, Froude stated his hypothesis that the ship resistance coefficient can be expressed as the sum of two components: the

residual drag coefficient, depending on the Froude number and the frictional drag coefficient, depending on the Reynolds number:

CD(Re, Fn) = CF(Re) + CR(Fn)

(1.1)

For almost every hull form, provided the ship velocity is not too low, the dom-inant part of the residual drag is associated with wave resistance. The viscous forces in the fluid are confined to a thin boundary layer close to the ship hull

and can be treated separately.

The neglect of viscosity in the Navier-Stokes equations results in the Euler

equations. Under the additional assumptions of constant density p and irrota-tional flow this leads to the class of Potential flows. Under the assumption of a potential flow the stationary problem of calculating the ship wave resistance is described by the Laplace equation, subject to conditions on the ship hull and the free surface. Many theorems and computer programs have been developed to solve this problem. For instance computer algorithms based on finite differ-ence methods and algorithms based on boundary element methods. The latter can be classified in two categories, based on the choice of the singularity:

Kelvin wake sources. These are originally due to Brard(1972). One of the

major advantages of this approach is the elimination of the free surface

from the domain over which the resulting integral equation is defined.

Also the so-called radiation condition is automatically satisfied.

Rankine sources. Dawson(1977) was one of the first to use the Rankine

source as the elementary singularity. One of the advantages of this method is the flexible use of different free surface conditions. One of the disadvantages is the introduction of a distortion of the wave sys-tem due to the free surface discretization. Also errors are introduced by imposing the radiation condition numerically. See also Sciavounos & Nakos(1988), Piers(1983) and Raven(1988).

The estimation of a ship's speed and power was usually based on still water

performance. However, in order to be able to predict ship performance n seaway, it is also desirable to be able to solve the instationary problem of a

sailing ship. For example a ship sailing in waves or an oscillating vessel. In preliminary design studies the use of a fast computer algorithm could help to assist the traditional model testing. With the considerable increase in computer speed and resources, the first numerical solutions of 3-D seakeeping problems appeared. Recently Nakos(1990) and Bertram(1990) reported work with the

Rankine source approach.

All of the other existing algorithms are based on the Neumann-Kelvin

ap-proach. According to Nakos(1990):

"Further developments of the unsteady scheme, stumble, however, on the addi-tional and numerical difficulties ¡n the evaluation of the unsteady Kelvin source potential, which serves as the Green function."

One of the drawbacks of the Rankine panel method was not being able to model the propagation of the complicated unsteady ship wave pattern by a

discretized free surface. Recently Nakos and Sclavounos(1990) proposed a free surface discretization scheme that is free of numerical damping and introduces

small numerical dispersion, while providing flexibility in the choice of the

lin-earization of the free surface condition.

In this thesis a different approach will be used. For the solution of the unsteady problem, the solution method of the steady state problem will be used. Partic-ularly the additional numerical difficulties are reduced so they can be handled more easily. The unsteady wave pattern will be calculated, especially the part near the vessel itself. This is important for the relative motion between the ves-sel and the free surface, shipping of green water, bow slamming and propeller clearance. Shipping of green water, for instance, is assumed to occur whenever the relative motion amplitude is greater then the freeboard (see Blok [6]). For the prediction of the wave resistance use can be made of far-field approxi-mations for the evaluation of the (added) resistance. Also a method based on direct integration of pressure can be used. This direct pressure method is the approach which will be followed here and is similar to Pinkster [30].

One of the main problems in this thesis will be the calculation of the dynamic swell-up (and as a result the added resistance) of a ship in seaway. Dynamic swell-up is the effect of water being pushed up around a bow higher than can be accounted for by considering heaving, pitching and incident wave alone (see

also figure (1.1)). The swell-up s given by the quotient of the amplitude of

relative motion ((i) and the osciUator amplitude (Za).

Figure 1.1: The dynamic swell-up.

was the first to introduce the term swell-up. On the basis of experiments

he obtained an empirical formula for the calculation of the dynamic swell-up, provided that the wave frequency, ship speed and ship block coefficient CB are between some given limits.

Overview:

This thesis deals with the hydrodynamical problem of a ship sailing at forward speed. Of practical interest will be features such as shipping of green water,

(added) resistance, sinkage, trim and water height. The mathematical

for-mulation based on smallness of the unsteady movements will be given in the

first chapter. A new free surface condition will be given, that connects the

steady flow to the unsteady flow. A linearization round the 'known' steady solution will lead to a free surface condition. The body boundary condition

will be expanded round a mean position with respect to the small displacement

parameter , to be defined appropriately in Chapter 2.

For the solution of the two coupled integral equations as a result from a

Neumann-Kelvin approach an expansion in terms of the small oscillation

fre-quency w has been chosen. The body boundary condition is also expanded

in this small parameter w. A modified Neumann-Kelvin formulation has been found, using only the steady state Green function for the calculation and

eval-uation of the ship wave resistance problem.

In Chapter 3 the computational aspects are dealt with. The integral equa-tion will be solved using a boundary element method. The wetted body will

be divided into triangular panels to solve the equations. One of the most time

consuming aspects of the problem is the calculation of the Green function.

This is one of the major problems when using a Neumann-Kelvin approach (as has been mentioned by Nakos in his Ph.D. Thesis).

In Chapter 4 the stability and convergence of the resulting matrix equation will be analysed. Some of the aspects of these equations, like singular values

and convergence of the solution algorithm will be dealt with.

In Chapter 5 some computational results will be presented. The forward speed

problem has been solved for a specific hull form: the Wigley hull.

At the

performed. The results from these measurements will be compared with the theoretical results presented here.

In the last chapter conclusions and recommendations for future developments will be given.

The research has been carried out in cooperation with MARIN and the Deift

Chapter 2

Mathematical Formulation

21

Introduction

We consider a ship moving horizontally in still water of infinite depth and a constant velocity U. Viewed from an inertial frame (x, y, z) attached to the ship, there is an incident flow of velocity U in the direction of the negative

x-axis. lt is assuming that the fluid is inviscid, incompressible, irrotational and

free from surface tension effects, wave breaking and cavitation. By virtue of

these assumptions, the velocity vector V of the fluid can be represented as the

gradient of a potential funtion :

V=V

(2.1)

The problem will be formulated in a right handed Cartesian coordinate system

Oxyz moving with the ship. The x-axis is pointing astern and the z-axis

vertically upward.

Figure 2.1: The moving coordinate system.

The free surface elevation is defined by a function z = ((x,y,t).

For practical applications, the vessel is performing in seaway. This means that there are incoming waves present or some sort of forced oscillation exists.

This means that there is not only a steady (time independent) but also an

unsteady (time dependent) part of the velocity potential present.

The total velocity potential 4ttal can be split ¡n a steady (the sum of the

uniform flow and the stationary wave part) and an unsteady part:

= Ux + (x)+

(x,t)

(2.2)

steady tinsteady 4tota1 has to satisfy the following conditions: alt has to satisfy the Laplace equation.

A1totai = O

¡n the fluid region ,z < ((x,y)

(2.3)

The physical nature of the free surface requires two boundary conditions: firstly the kinematic boundary condition, which states that the normal veloc-ity of the fluid equals the normal velocveloc-ity of this free surface and secondly the dynamic condition which requires that at the free surface the pressure

should be atmospheric and independent of the position on the free surface(see

figure 2.2).

undisturbed water surface

Figure 2.2: The free surface conditions. The free surface S is given by the equation:

T(x,y,z; t) = 0

(2.4)

In our case the free surface is given by the equation z = C(x,y;t), so

T =

z - ((x,y;t).

This free surface should have the property that any particle which is once on

Remark

dT ÔT 0T ÔT aT

onS

y ôz at

Thus equation (2.5) results in the so-called kinematic boundary condition:

ay

---=O onS

The dynamic boundary condition can be derived from Bernoulli's law:

¿N)

i

P = Pam - PYZ - P

-

. V + C(t)

when Vn is the ncwmal velocity of the body.

(2.5)

The solution also has to satisfy a so-called radiation condition. Fo.- this, the reader is referred to section (2.2.3).

n

1One has to assume that the fluid is ¡deal and that the surface tension is negligible, which

indeed is the case. The surface tension for air-water is about 0.07 N/m (see Newman [28]). (patin

So at the

conditions

.Finally

may be taken equal to zero and C(t) = pU2 free surface ((x,y,t) the exact1 dynamic

read:

V4

U2

without loss of generality) and kinematic boundary

at z = (

(2.8)

at z = (

(2.9)

+

.

+

gz =

a

a

a

a

äatax

ôx _Oy ôy the boundary condition on the ship hull:

at

vn,

E

=

Un (2.6) (2.7)

2.2

Linearization of the boundary conditions

As can be observed from equation (2.8) the main difficulty in this problem is to find a solution of the Laplace equation satisfying boundary conditions at a free surface which is still unknown. One way of dealing with this problem could be in using a direct discretization scheme. For the complete time dependent problem this would involve some sort of iterative solver: at every time step a new solution must be computed andbefore this is possible, the free surface has to be calculated iteratively as well. Obviously this would lead to an algorithm which is very time consuming.

In 1898, Michell was the first to obtain an approximating analytical solution

for this problem. His 'thin' ship theory has become worldwide famous. Michell and others after him obtained his perturbation scheme using the slenderness2 of the ship. Another small parameter could be thevelocity of the ship. If the Froude number, based on the length of the ship, is a small parameter, an ex-pansion with respect to this parameter can be obtained. Brandsma ([7J,[8]) has developed a computational method to obtain the wave resistance of full ship forms. Further work can be found in Baba [21, Sakamoto [331, Hermans [17] and many others.

To derive an approximating solution of the problem we continue by linearizing the free surface condition (2.8) and the body boundary condition (2.9).

2.2.1

The free surface condition

The free surface condition (2.8) consists of two parts: The dynamic boundary condition:

+VV4 +gz =

U2

and the kinematic boundary condition:

o'1

ô(

o o a

Oz ot Ox Ox _Oy Oy

at z =((x,y,t)

(2.10)

2Note the difference between thin and slender; Thin means that the beam is small

com-pared to all other characteristic lengths of the problem, by slender we mean that the transverse

dimensions of the body are small compared to its length.

The complicated non-linear nature of the two free surface conditions (2.10)

and (2.11) prohibits the development of an exact solution. Therefore a method of approximation is required.

First we give a simple example of the linearization procedure. For instance in the case of small unsteady potentials and the forward speed equal to zero one

way of dealing with this problem could be neglecting the second and higher order terms in 1 and (.

Then the kinematic boundary condition (2.11) leads to:

212

azat

(.)

This approximations means that the vertical velocity of fluid equals the vertical velocity of the free surface.

Substituting the (for z in equation (2.10) leads for the linearized free surface height to:

i ô4

at

(2.13)

combining equations (2.12) and (2.13) leads to the well known free surface

condition:

at2

+

= O

at z = 0

(2.14)

Of course this free surface condition should be imposed on the actual free sur-face z = (, but in this linearized theory it makes no difference.

D

We now proceed with the more general case. In order to be able to solve

the steady and its perturbed wave problem, we introduce a more complicated

linearization by making further assumptions. We assume that the unsteady wave part is due to a small displacement of the ship. So when the body is displaced from its equilibrium position, the deflection of any point of the hull

is assumed to be small (see also section 2.2.2).

In many cases a small parameter is involved. For instance can be the

slenderness of the ship or the magnitude of the instationary disturbance. In this thesis we consider the case that the instationary potential is small com-pared to the stationary potential. Therefore we replace

{x,t) by cç(x,t) in

16 CHAPTER 2. MATHEMATICAL FORMULATION

We will now obtain an asymptotic solution for the free surface condition. Let

« i and expand the free surface elevation around a 'known' solution ( =

((x,y,t)

(o(x,y) + c(1(x,y,t) + ...

(2.15) Denote the total velocity potential by total, and split this total potential in a

steady and an unsteady part as follows:

4)totai(a,t) = Ux +

(x) +

EÇb(X,t) (2.16) All the terms in the free surface condition (2.10) and (2.11) have to be

ex-panded (see Appendix A for details). This leads to the following free surface condition:

Foi the zeroth order non-linear problem (the steady state solution):

(0(x,y) =

U-

ç (x, y, 0)

g

and for the first order linear problem (the unsteady state solution):

-

-

(oy) -

zz(t + V(Ux +

V) +

2V(Ux +

)

Vçf. + V(Ux + ç) V[V(Ux + ç) V4] +

(2.19)

ctt

= O

atz=(o

with the first order wave elevation given by (see equation A.14):

(1 =

_{+V(_Ux+YV}

(2.20)

(2.18)

+

V(Ux +

) V[V(Ux + Y V(-Ux +

)] = O

_(2.17)

at z = (o

Regrouping of equation (2.19), by leaving the terms which are independent of the steady state solution to the left, gives us:

g& - 2U

+

tt +

U2rr =

at z = Co (2.21)

Where £(U;) denotes a linear differential operator acting on as defined in

equation (2.19). This operator can be found in appendix B.

The linearization procedure (of the F.S.C.) can be performed with respect to ei-ther the undisturbed uniform flow or the double-body flow (see also Raven [32]). The former approach has been used in this thesis, the latter has been used by Sakamoto and Baba [33]. After some rather strong demands (the unsteady

wave part must be of order O(U') with 3 < n < 7) Sakamoto and Baba

obtained two independent linearized equations for the unsteady steady free

surface conditions: 1 O

02

+

<arO + VrQ O =

D(x,y)

z g Ox OyJ

1(0

0 0 2

+UrØ+VrO

1+0

g 1.Ot Ox OyJ

where (see Sakamoto and Baba [33]):

1r velocity potential representing double body flow in calm water. This

is regarded as the base flow on which the following perturbation poten-tials are superimposed.

o: velocity potential representing steady wave making. velocity potential representing unsteady wave making. ttro,Vro: velocity components of double body flow at z=O.

D(x,y)

= Ox

+

_ay

Cr: wave height of double body flow.

In this thesis we are especially interested_in the formulation of the influence of the steady part ç on the unsteady part ç, so we will use equations (2.17) and

(2.19).

at z =CT (2.22)

2.2.2

The body boundary condition

The body boundary condition for the stationary potential reads as: (see also

equation (2.9))

= V

at the body E (2.24)

On

or equivalently for the steadilysailing ship:

n- V(Ux +

) = O

at E (2.25)

Next we consider the calculation of the instationary part (see also forinstance

Timman [391 et al.). The displacement of the ship hull is assumed to besmall

compared to the length of the ship. We write:

z - z' =

&(z,t)

(2.26)

where is a small parameter, z denotes the coordinate system moving with velocity U in the direction of the positive x-axis, and z' denotes the

coordi-nate system fixed to the ship. The body boundary condition states that: at

the instantaneous position of the hull, the norma! velocityof the fluid is equal to the normal velocity of the hull.

Expansion of 4 around the mean position with respect to a small parameter

yields:

4moncntary = toia1 +

''totai

So for the body boundary condition this leads to:

n V(40k +

4totai) c n

Since on E we have that n- V(Ux +

) = O equation (2.28) reduces to:

= n(

-

V(à-V(Ux+)))

Now using the assumption of a(z,t) and (z,t) to be oscillatory:

&(x,t) =

{a(x)

-(z,t)

=

{(x) e}

(2.30)

This finally leads to the following body boundary condition:

(2.29) (2.27)

= n.(iw +V(aV(Ux+)))

and from the free surface condition (2.21): _w2

_{+ 2iwU}

_{+ U2}

_{+ g}

₌

at z = Co

where r(U;) denotes the linear differential operator acting on çb.

2.2.3 The radiation condition

In this thesis we are dealing with the situation of a body which is oscillating

with a frequency ' and moving steadily with velocity U. Besides the boundary conditions given in section (2.2) an additional condition has to be imposed: a condition at infinity to ensure a unique solution.

For instance a body moving steadily has no waves (fluid motion) far ahead and far below of the body. An oscillating body has no waves far below of the body,

but outgoing waves at infinity in all directions (radiation condition).

In our case, the body is both oscillating and moving steadily, so the situation is more complicated. For a small velocity U the situation will have the tendency to behave like in the case of an oscillating body. (i.e. a radiation condition). However for larger velocities U, the situation will have the tendency to behave like a steadily moving body and the condition for the steady state case has to

be imposed.

This situation will be taken care of by formulation of an initial value prob-lem.

The solution obtained as we let t -

is the desired solution (see

Wehausen [43] and others).

Remark:

A simplification of the steady free surface condition( 2.17) can be obtained if is small compared to Ux. This is true for thin or slender ships. In this case the steady state problem is called a_Neumann-Kelvin problem. The linearized free surface condition (U2ç5 _{+ g} = O) is that used by Lord Kelvin. The

exact body boundary condition on the ship hull is of Neumann's type.

(2.31)

(2.32)

In the next section we will solve the instationary problem by using Green's theorem.

2.3

The leading equations

2.3.1

The singularity distributions

A Green function G(x,; U)) is introduced as a solution of Laplace's equation, representing a source of oscillating strength in uniform motion. In this section

we will construct a solution of the Neumann-Kelvin problem in terms of an

integral distribution of Kelvin wave sources. We will apply Green's theorem to the domain as can be seen in figure (2.3).

Figure 2.3: The mathematical domain.

Let denote the wetted surface of the hull. Cf s the closed contour

(water-line) intersected by in the free surface. C1 is the closed contour intersected by the free surface in a certain vertical prism. Let p, be the part of the free surface located inside the hull and the free surface outside , with

the part of

Fe inside the prism. Let D denote the interior of the closed

surface E + F, and De the interior of the closed surface E +

+

+

the free surface and an infinite prism.

We choose the function G such that it fulfills the following adjoint

homoge-neous free surface condition:

w2G - 2iwUG + U2G

+ gG(

= O at =0 (2.33)

with the function G written as:

G = (

- -- +

W(x,;U))

(2.34)

where r denotes the distance between the source and the field point and r1 denotes the distance between the field point and the image of the source above the free surface. In the following we will apply Green's theorem to a problem in

D inside E and to the problem in De outside E. Combining the formulation inside and outside the ship we obtain a description of the potential function defined outside E by means of a source and vortex distribution (see equa-tion (2.47)_).

In the same way as Brard [9J we consider the following integrals:

and

1EJf

F, +

Ee+IT

(e

- G Oe)ds

0e

At EF and EFe equations (2.32) and (2.33) are valid:

U0

=

{r(ç) +

- 2iwU

-g

0G

_'{w2G

0G

-u2 02G

+2iwU-0e2}

(2.35) h

=

(q5-

-

G)dS

(2.36) +

then the following is valid:

I 4irç&e

X E D

xED

(2.37)

Ij -

{ 4ir

X E D

-

O

XEDe

(2.38) (2.39) (2.40)

Combnng equations (2.39)and (2.40)leads for equation (2.35) to:

8G

G

=

ÇeG + 2i 0G U2 02G

(&

Ofle ôfle g

-

8C2

+

2iw ÔeG

U202çe

--

eG+U

g g g

+

G (2.41)

Regrouping the expressions at the right hand side of equation (2.41) leads to

the following formulation:

0G

2iw

_0(&G)

U2O(G -

_(e)G

_(2.42)

Ô4

¿k

g

and almost the same for equation (2.36). The EF, part of equation (2.35)

leads to: 8G 8Çe_)dS

=

I

--U

e G d71 +

(,- - G

äfle

äfle

J 9 coo - Cj 8G U2 ôç

J (-VG -

q5e--)di

--JJ

£()Gds

(2.43)

Eç,

and for the EF1 part:

JI

EF,

[f

8G

Ji

(çti EF, + f 2iWu(, (2.44)

cf

-

G)dS

=

J

i8

_-

0G

Adding equations (2.43) and (2.44) results in:

a

_-

a

a

- afl

- ¿371

'E

+ Ij

=

-

f_.)GdS

-

Uf(e

- qf)Gdi7 +

E C1

-

1f

_{r(e)Gds -}

_{fJ £()GdS + JJ(ci - ç)dS +}

gli

EFe EF E (with 7

=

eÇbi

aj

a

Dn

-a =

u2 ¿3G U2

a

e

+

_{- J}

( )Gdi1 + J(&

-g Ox CI Cf

The source and dipole strength are defined as:

equation (2.45) will now transform into:

4(x)

_JJ 7()DG(xE)ds

_-

_{ffa(e)G(xe)ds}

₊

E E 2zw U2

0G

-

gf

U

7()G(x,)d + J7(e)

'd+

g Cf Cf U2 ' _{D-y ¿37}

+

J (at-- - cr + ana())G(x,e)dui +

t

Dr

Cf

f[)G(x,)dS

gJi

F.S. (2.45) (2.46) (2.47)

r = tXn

cos(Ox, t)

= cos(Ox,i-) =

cos(Ox,n)

where t is the tangent to the waterline.

87 O-y

- od- - ch.- per

unit arc length. Or

(2.48)

So the potential is generated by the following five singularity distributions(see

also Brard [9]):

A distribution of simple sources over with strength a per unit area. A distribution of double sources over with strength y per unit area.

U2 A distribution of double sources over the waterline with strength --y

per unit arc length.

A distribution of simple sources over the waterline with strength

g

A distribution of simple sources over the entire free surface with strength

£( ).

There is still a free choice in -y and a. We can make

çj ((e) = O) at

the boundaries. This means that the tangential velocities are now continuous

at the boundaries, but the normal velocities are discontinuous (see also the remark on the next page).

With a choice of y()

O, the following expression will be obtained:

4(x)

= -

ffa()G(x)dS

+

+

fi £()G(x,)dS

9JJ

F.S. (2.49) with:

The potential (x) is expressed in simple sources only!

Using the body boundary condition (2.31) we may obtain a description of the potential function by means of a source distribution of the following form:

ôçf

[f

aG(,e)

4ir = 2ira(x)

_-

_jj

o(e)

_{dS +}

E

u2

+

fn)dii_ !ff £()8'dSe

_(2.50)

gi

an

C1 F.S.

with (see equation (2.31)) given by:

an

n (iwa + V(

V(Ux + j)))

Re ma rk

We consider the following potential 4):

4) =

ffdS

The following is valid (Kellogg [21]): if the densityo of the distribution on

is continuous at x, the normal derivative of the potential 4) approaches limits as X approaches x along the normal to at a from either side.

These limits are:

()+ = 2rcr(x) +

ffa()"'dSe

(2.53)

() = +2ira(x) +

ff)TdS

(2.54)

D

So when using the sources the resulting equations looks like:

rr

_ôG(x,)

+2'rcT(x)

-

_ji

I/a

_{dSe +}

_{... =}

an,x E (2.51) (2.52) (2.55)

and when the sources are used (as in Brard [9]) the resulting equation

transforms into:

2.3.2

Expansion of source strength and potential function.

We will now consider small values of . From model test experiments one

generally concludes that a large part of the dynamic swell-up of a ship originates at the lower frequencies. So the vessel is in a nearly steady state motion, but still swell-up occurs.

On the basis of experiments Tasaki [38] obtained an emperical formula to

evaluate the height of the dynamical swell-up (see also Blok [6]):

ir

8G(x,)

27ra(x)

- I/a

dSe +

...

=

JJ

onx

E

SA g

with:

C =

(CB -

0.45)

this formula should only be valid for the following: 0.16 <

Fn

< 0.29

2_L

1.6 < s-- < 2.6

(2.59)

g

0.6 CB

<0.8

Now the following approach is justified. We can solve the unsteady state prob-lem by expansion of all the relevant physical parameters in this small parameter. We will use an expansion of a, and G in the small parameter :

a(x)

_{= ao(r) + wcy(x) +}

=

o(x) + w1(x) +...

G(x,) = Go(x,) + wGi(x,) +...

(2.60) (2.56)

(2.57)

Substitution of these expansions in equations (2.49) and (250) leads for the first order problem to:

-

4n

V(Ux + )) = 2o(x)

_ff

()oco(xe)dS

+

onz

E

+ - I

_gJ

!JJ

0(ç5)

°'dS

(2.61)

C1 F.S.

where is given by:

U2

- 4o(x) =

_ffao(e)Go(x)dS

+

_gi

f0(e)G0(x,)diì+

E C1

+

_g

_JJ

ro()Go(x,)dSe

F.S.

and for order w:

ôGo(x,)

_{dS + - [ai(e)}

U2 ÔG0(x

'd =

2iri(x) - ffcTi(e)

ôn

g

8n

E C1

-

_{47rn.ic+fJao()°"'dSe}

U2

onz

--

f

oao()

ôflv

+

E Cf

+

-

i

(r

1+ro

OG0(x,e)

+

Ow1(x,e)dS

g

fi

onz

Oflz

F.S.

where 'i is given by:

= -

_{ff(ao(e)Wi(x) + i()Go())dSe +}

+

f a(ao()W1(x,) + ai(E)Go(x,))d

+

Cf

-

_{H (r +}

L1 )Go(x,

) +

(2.62) (2.63) (2.64)

2.3.3 Expansion of body boundary conciftion.

The right hand side of equation (2.31)contains a vector

k(x t)

f Ok(t) tk

k = 1,2,3 (2 65)

'S

Uk(t) (ik Xx)

k = 4,5,6

Here k(t) is the deflection in translational motion for k = 1,2,3 and for k =

4,5,6, ck(t) represents rotation angles about the Xk_3 -axis. See figure (2.4).

heave

yasv

roll

surge

-Figure 2.4: Body motions.

Combining equation (2.65) with equation (2.31) leads to the body boundary

condition.

So for k=1 (surge, translation parallel to the longitudinal axis) this looks like:

= (2.66)

k2 (sway, translation in the lateral direction):

k=3 (heave, translation in the vertical direction):

ôç - - .

-=

fli&jzfl2yzfl3(W.'+&z)

k=4 (roH, rotational motion about xi-axis)

On

+ n3(yiw +

+ z-J

- yxz)

k=5 (pitch, rotational motion about x2-axis):

On

= ni(ziw - zçb

+ & + x&) + n2(zqxy + yÇyz) +

-

y) + n2(ziw + zyy -

- Yyz) +

ri

0G0 0G1

2irai- //(ai

ji

.

+ a0 )+=ni

On '----' On

E =0

When 4 = 0(1) then the expansion will start with leading order zero (a a0 + wa1 + ...) and equation (2.61) will look like:

fr

0Go -

-

-27ra0

-

J]

a0----

= 21Ç5xz

-n

E

lt should also be noted that for the calculation of the steady wave potential

the same matrix equation (261) as for the unsteady potential can be used!

But of course the right hand side for now looks like:

= nU

(2.74) (2.68) (2.69) (2.72) (2.73)

+ n3(xiL) + U -

- zçL + Xqzz) (2.70)

k=6 (yaw, rotational motion about x3-axis)

=

On

+ Yxy -

+ fl3(j

-

(2.71)

When ç5 «

i then the expansion (2.60) will start with leading order w

(a = wa + w2cr2 + ..) and then equation (2.61) will transform into (for

2.4

(Added) Resistance.

If the potential '1 is known, the pressure in a point in the fluid can be calculated using Bernoulli's equation:

8cI

i

pV.V4+C(L)

In order to determine the resistance of a vessel, two different lines of approach can be followed (see Pinkster [30]):

Using momentum considerations. The change of momentum of the fluid

is equated to the mean force. Placing the momentum control volumes

at infinity use can be made of the knowledge of the far field behaviour of the potential flow (see for instance Maruo [24] and Huijsmans [18]). Direct integration of pressure. This approach will be followed in this thesis. The underlying approach is similar to Pinkster [30].

The fluid forces can be determined by direct integration of the pressure:

F

= -

JJ ndS

(2.76) and for the moment:

M = -

JJp(a

x

n)dS

(2.77) in which:

p = fluid pressure

n

= outward normal S = wetted surface

x = coordinate of surface element

(2.75)

The oscillating movement of the ship is defined by the sum of the translational and the angular motions (see also equation (2.65)):

X = Xg(t) +

(x,t)

(2.78) Where x9(t) denotes the motion vector of the centre of gravity of the body. Now the first order (steady state) force is given by:

F= _JJ_pgz

+pU

+ p(+ +)ndS

(2.79)

s

If the vessel is in an oscillatory motion, the pressure in each point can be found using a Taylor expansion of this point round the mean position:

= Po +

(x(')

-

x(°)) Vp +

... (2.80) Using this, Bernoulli's equation (2.75) leads to:

PIxrx(1) = pgz° - pgz'

5)o)

5) (0)

'J ₍₁₎

-

₍₀₎ 'J

pot

us

-

- at

pV°

.

v(°) -

p(x(1)

-

x°)) V(V(°)

.

Using Pinkster [30] equation (Iii-68) and further, the following equation for

the second order force can be found:

F

=

-

fJ p ndS

(2.82)

s'n$t

The instantaneous wetted surface is split in two parts: a constant part

So up to the static waterline on the hull and an oscillating part s between the

static waterline on the hull and the wave profile. The force due to the part s is given by: F5

=

_JfpndS

Z=o +(] (2.81)

=-J J

(2.83) WL z=(o--z

At the waterline WL the following is valid:

a

i

pg((o +

_{(i) =}

p

-so equation (2.83) transforms into:

F3

=

-

J J (pgz + pg((o + () n dzdl

WL z=t0+z3

The time average of this force is given by:

=

çg

- pg(z)ndl

(2.86)

Now the second order (unsteady state) force is given by the sum of the following six components (see Pinkster [30]):

Wave elevation:

pgf (n dl + pgf ClZafl dl

(2.87)

The first order velocity:

- :;(

_pIVI2ndS

(2.88)

Product of first order motion ((1)) and gradient of first order potential: (2.84)

(2.85)

Product of first order motion and gradient of squared velocity:

-

Jf

1)

V(V

V))ndS

Product of first order angular motions and inertia forces:

aX (Mg)

Unsteady potential:

The first order wave elevation is given by (see equation A.14):

=_{+V_Ux+.V}

with the use of the summation convention

3

=

where the subscript i denotes differentiation with respect to x.

This can be written as:

cl =

_

{

-

+

(2.90)

(2.91)

-

Jf

p'ftndS

(2.92)

In our case only the time average is of importance, so for the first contribution this is: i i I.

(-pg

_{J (?ndl)}

(2.93) WL (2.94) (2.95) (2.96) Now consider heave motion (mode = 3). The motion and potential are given

by:

a(x,t)

₌

{c(x) e}

so

(x,t)

= {cos(t) - isin(wt)}

(2.98) For heave mode the quantity is pure imaginary, so this leads to

(We replace ç by i 5):

=

sin(t)

(2.99)

c(x) = (0,0,1)TZaCOSWt

so the first order wave elevation (equation 2.96) is now given by:

(1 =

2{wcoswt

- Usint +

(2.100)

Most of the unsteady forces (2.87) - (2.92) can be calculated easily. For

instance equation (2.90) produces:

o

_'\

((-U +

+

+

sinw) \

0

_1.21

1(2.101)

Za

cos.t J

k

(7 + & sin wt)(

+ zz sin.'t)

J

In this equation only products like sinwtcoswt occur.

In order to calculate the time average, the following three relations should be

remembered:

Jsintdt=

ifcostdt=0

(2.102)

-

_{f Cos2tt =}

_ f sin2 tdt =

(2.103) 2ir

---

[sintcostdt=0

2ir j

o (2.104)

((1)

.

V)) = O

So contributions in the mean force can only originate from the (sin2 t)-terms and the (cos2 1)-terms.

For pitch motion (mode = 5) is given by:

= çcos(wt)

(2.105)

c(x) = (z,0, _X)TOa cos.'t

So equation (2.90) leads to:

(

z

\

((-u +

₊

₊

OaCOSWtI 0 1.21

-x )

\,

Now the time average is given by:

8a{Z((U +

+

-

(2.106) In the case of equation (2.87) - (2.92) this leads to:

HEAVE:

1{1,22

₊

+

1()2

-i

(iza =

)ZÇb 2g 2

22)

2x

((1) v)

=

ZaW5z

PITCH:

=

{tw22

_{+ U2ç +}

- U(11)}

=

±WXOa 2g

222)

2x

Y

(X(') V) = O

((1) V(V4. V))

=

Oa{Z((U +

+ xx) -

zz)z}

The dependence on the different parameters w and Za can be found in the

following table:

Table 2.1: The dependence of added resistance on amplitude and frequency.

heave pitch

frequency amplitude frequency amplitude Force-1

(,2

_{+ ...)} -.- (w2 + ...) Force-1h w2 - -Force-Ip - - w Force-Il w2 -. ' w2 Force-Ill w2

-

- -Force-IV - - (U + w) Force-V - - - -Force-VI - - -

-Chapter 3

Computational Aspects

3.1

Solution of the integral equation.

The solutions obtained by the singularity distributions are the result of solv-ing two coupled integral equations. Equations (2.61) and (2.62) can be

solved using an iterative scheme. In this scheme use has been made of the numerical evaluation of the wave resistance Green function as done by

New-man ([26]+[27]).

ek

Figure 3.1: Panel influence.

The numerical solution of equation (2.61) and (2.62) is carried out using a

finite element method. The wetted body is divided into N triangular panels (see figure (3.1)).

n( x)

X

Integration is done using a piecewise constant variation of the source strength

a(e).

JI

(JJ

a(C)dS

onx

E j=1 e

In this way a set of N linear equations for the N source strengths a is obtained. So the main structure of the integral equation (2.50) looks like:

= 2(x')

(E(l))8G(xU,())s(l)

₊ (1)

'9f'

i) rr (2)

_{S(1)+ ...}

Un1 = 2iw(x2)

-Un1(N)

2((N))

((l))'9G(X())

(N) S(1)+ (N))

_an,

(N) S(N)

with the normal vector n at element number (i) denoted by:

S(i) denotes the area of element number (i) and G(x,C) the Green function.

The matrix equation now reads as:

Aa

(k)

_(3.3)

where

ÇI')

depends on the (k)-th oscillation mode.

(3.1) <N) 7( )

_an,

(1) (N)

i(

)

_an,

(2)

=1

I

/

(i) I i (i) 2 I (i)

\ fl3

\

I I I J

(3.2)

So for each element we need to evaluate the following integrals: For the evaluation of x) when a(x) is known:

fdSe

fG(xe)dS

JG(x)diì

(3.4)

For the calculation of the integral equation:

J1dS,

ônx ek

fDG(x,)

[ôG(x,)

i

dS,

i dij J Ufla; J Ofla; 1k (3.5)

with e an area element and 1k a waterline element. Numerical errors occur at different places:

Ship hull panelling.

The numerical integration of the Green function.

Numerical integration of the waterline integral; errors occur here as a

result of the singular behaviour of the Green function at the free surface.

The Green function G contains a term. The integration of these

1/r-terms is carried out by a subroutine developed at MARIN [121 using analytical

expressions of those integrals in order to avoid large errors from numerical

3.2

Generating the grid

3.2.1

Hullform definition

In this thesis, most of the computations have been performed on two hull forms: a Wigley hull and a modified Wigley hull. In this section the generation of the

grid for these two hull forms are discussed. Of course, other ship hull forms

can be implemented without much difficulties.

-1.50 -1.25 -1.00 -0.75 -0.50 -0.25 0.00 0.25 0.50 0.75 1.00 1.25 1.50

-1.50 -1.25 -1.00 -0.75 -0.50 -0.25 0.00 0.25 0.50 0.75 1.00 1.25 1.50

Figure 3.2: Wigley hull contour plot. A right-handed coordinate system is defined by:

L

amidships in the waterplane longitudinal, positive forwards lateral, positive to port side vertical, positive upwards

Then the non-dimensional hull form of the model (a two-sided parabolic hull form) is defined by (see figure (3.2)):

=

(1-(2).(1-C2)

(3.7)

I

with: O

e:

q:

H (3.6)

The principal dimensions of this specific Wigley hull are given by: = 10 and

L H

In modern analysis also a modified Wigley hull definition is used:

= (1(2).(1-2).(1+a22+a4.4)+

for all models for all models

, for models with CM = 0.9090

for models with CM = 0.6667

The mid-body of this modified Wigley hull is fatter then the mid-body of the standard Wigley hull. The shape of the modified Wigley hull looks more like

the shape of commercial vessels then the standard Wigley hull.

3.2.2

Panel division

In our calculation the wetted body in divided into triangular panels . This

has been done in such a way that a symmetrical grid with respect to the y O plane has been obtained. The panel arrangement can be observed in figure (3.3).

3.3

The Green function

3.3.1

Behaviour of the Green function.

In Wehausen and Laitone [43] the Green function of an oscillating source is

given:

=

I

_{+ 4'(x,e;U)}

(3.10)

with the function 11 given by:

=

1ffF(O,k)dOdk

O L1

+

]JF(9,k)dodk

(3.11) 1, L2 with: a2 = 0.2 , a4 = 0.0

A =

1.0

=

0.0 (3.8) (3.9)

3 x=O 4 5 6 _{z = O}

Figure 3.3: Panel arrangements.

where

F(9, k)

is given by:

kek(Z+(+ i(x-)cosO)co5(k(y_)sjflg)

F(O,k)

-gk

- (

+ kUcos9)2

the paths L1 and L2 are given in figure (3.4).

With k1 - k4the poles of F(O, k):

Jgk4

=

+ 'J1-4rcos9

U)

OE(,ir)

2r cos G

z = -H

(3.12) (3.13) (3.14) (3.15) (3.16)

With

r =

(so r « i

in this thesis) this leads to:

Jgk1 = \/gk3

w + 0(w2)

(3.17) 1 3

45

9

442l

22

iì°v

i -

/1-4rcosO

=

w

OE(O,)

i +

2r cos O

v"1-4rcosû

=

U)

9E(O,)

i

-2i- cos O

./1-4rcosO

=

U

9E(,7r)

2r

cos O

zLy

Figure 3.4: Integration contours L1 and L2.

g

+w+Q(w2)

U cos O (3.18)

When for instance U -* O then k2 and k4 willgo to a and the paths L1 and L2

will coincide with their poles located atk1 = k3 w2/g, see also figure (3.5).

Figure 3.5: The Iocition of the poles when U - O.

When w - O the poles k1 and k2 will move to the origin and the paths L1

and L2 can be seen in figure (3.6). However this is not correct. When w -+ O a factor k can be removed from the function

F(9, k)

and looks like this:

ek((+(X_)CO5O)cos(k(y -ij)sinO)

Fo(9,k) =

_g

-kU2cos2O

with only one pole, located at k2 = k4 _{J2 cos2}g O. (3.19) k4 k3 k1 k2 o k1 = k3 L1

=

L2

k3 k4

o k2

A

Figure 3.6: The location of the poles when w -* O.

3.3.2

Calculation of the Green function.

Newman ([261+ [27j) has published two papers in the Journal of ship research

concerning the evaluation of the wave resistance Green function: one for the

calculation of the double integral and one for the calculation of the single integral on the centreplane. The Green function is written in the following form (see Wehausen and Laitone [43]):

i

i

4

¡.

Je'cos(kxcos9)cos(kysin9)dkdO

R0 R

7rJ J

kcos20-1

00

J-Z

sec2 0 sin(x secO) cos(y sec2 O sin O) sec2 OdO

o

(3.20) The quantities are defined as can be seen in figure (3.7).

with R

=

Jx2 + y2 + z2.

As has been done by Newman, this integral can be split in a Double and a Single integral as follows:

+,r

2

j

jkz+ikIxIsectL'+kytaniI,

Double = i

dkd,I (3.21)

image

/

s field source z X

Figure 3.7: Location of field-, source- and image point.

+ ir

Single = 4iH(x)

_f

sec2 6e sec2 O + ix sec O + ijyl sec2 O sin °io

7T

(3.22)

3.3.3

Double integral

The double integral as given by equation (3.21) will be approximated by

Cheby-chev polynomials as done by Newman. In order to approximate a function f

of one variable x in the normalized range [-1,+1], the Chebychev polynomial of order n. is defined by:

T(x) = cos(n arccosx)

-1 < x < +1

(3.23) so the function 1(x) is approximated by:

N

f(x) =

cmT(x)

m=O

The coefficients c1 can be found according to:

N m

Cm =

_{La "f@mn)Tm()}

1.0 f(

LO

ri=O

with t

= 1,Em = 2, the double prime

indicates that the first and the last

terms in this summation are multiplied by . The coordinates x are given by:

=

cos(i)

(3.26)

(see Fox and Parker [14]).

For the 3D-case the situation is completely equivalent. In equation (3.21) logarithmic singularities are present when R = O. These singularities must be subtracted and approximated first in order to improve the convergence of the

approximation. The final approximation is given by:

16 16 8

D

S +

CIkT[f(r)]T(1 +

_ir

ir

i=O j=O k=o

The function f(r) is defined so that the interval (O,x) transforms into

the interval (-1,+1), see figure (3.8).

10 15 20 25 30

Figure 3.8: The transformation function f(r). and O are defined as:

(3.25)

(3.27)

x =

Rsin6

T, T3 and T2k are Chebychev polynomials. S is the logarithmic part of the

double integral. The Chebychev coefficients Ck are calculated and tabu-lated by Newman. Also the differentiated Green function has to be evaluated,

which contains terms like ¿iG(x,)/Dn (see equation (2.61)). Each part of

the expansion (3.27) has to be differentiated and evaluated analytically. The following terms have to be evaluated:

¿ix ôc (3.29)

with the transformation matrix A given by:

¿ix ¿Ix ¿ix

¿19 ¿ic ¿IR

A (3.30)

06 ¿Ia ¿IR

¿IZ ¿IZ ¿iz J

And also for the singular part S. For example the next figure (3.9) is obtained. The singular character is well shown here.

3.3.4

Single integral

After changing the variable of integration to s = sec 6 and i = tan 9, the

following expression equivalent to (3.22) is obtained.

f(x,y,z) = 4H(x) J sin{(x + yt)I1 +

t2}_Z(l + t2)dj

(3.31)

Numerical integration of this equation is rather difficult because of the rapidly oscillatory integrand as can be observed from figure (3.10).

In order to avoid time-consuming numerical integration Bessho [5] has derived two complementary Neumann series (a summation of different kind of Bessel functions) to evaluate these integrals.

The single integral as given by equation (3.22) is evaluated by Newman at the centreplane, i.e. a special case where the source and field point are in the same longitudinal plane. This is especially important when analyzing thin ships. The case y O will be dealt with later. For instance with the use of

ctv'rP

i*XL4Z

Figure 3.9: Function plot of the Double integral.

F(x) =

e2

Jt2

O < x < (3.34)

S(x, y, z)1=o = 8H(x)

_J

sec2 0e_Z sec2 O

sin(x sec0)dG (3.32)

o

In each of the different x-z regions (see figure (3.11)), the integral will be ap-proximated differently.

Region A: (small z) an expansion involving differentiated Bessel functions of the second kind:

S

=F() -

[Yi(x)+

-J

(3.33)

where Yi(x) denotes the Bessel function of the second kind, equals and F() Dawson's integral:

4.0

0.4

0.4

-Figure 3.10: Integrand of single integral.

Region B: an expansion in Neumann series, products of Bessel functions of the first kind and modified Bessel functions of the second kind:

i

s

(_l)?2J2+1(x) [K(z)

₊

_Kn+i(Z)]

4.0

(3.35)

where .J(x) denotes the Bessel function of the first kind and K(z) the

mod-ified Bessel function of the second kind.

Region C: large distances form the origin, steepest-descent expansion. The final expansion looks like this:

Figure 3.11: Domains for the approximation of the single integral.

where only the coefficient B71 has to be evaluated.

For all the expansions in the different domains (A)-(C), the differentiated Green function has to be evaluated. An example for the Green function can be seen in figure (3.12). The wave character is well observed here.

3.4

Padé approximations.

3.4.1

Introduction

The basis for the Padé approximation technique is the formal Taylor series

expansion. From this basis a Padé approximation can be found.

lt

is also

possible for a Taylor expansion to be divergent and the Padé expansion to be convergent and vice versa.

As an example we will use the Taylor expansion of the exponential function and the Euler function.

eX

=

Definition

Figure 3.12: The Single integral at y O.

XII i

= i + X +

x2 +

+

n! 2

et

E(x)

=

f

i + xt

dt = i - x + 2x2 - 6x3 +

r1 = 0

(3.38) o

The idea of Padé expansion is to approximate the function by a rational

func-tion of the following form:

(see Baker [3],[4J)

We denote the [L, M] Padé approximant of 1(x) by:

[L/M}

PL(x)

QM(x)

rc =

(3.37)

where PL(X) is a polynomial of degree at most L and QM(x) is a polynomial of degree at most M. The formal power series of f(x) reads as:

f(x) =

(3.40)

When we require:

f(x) - [L/M] = O(xM+1)

(3.41)

Then the coefficients of PL and QM can be found according to the following

scheme (Baker): a0

= Po

a1

₊

a0q1

= Pi

a2

+

a1q1

_{+ aoq2}

=

P2 a1

₊

a1_1q1

₊

aoq

=

Pl (3.42)

a1+

+

ajq1

+

a1_miqm =

O

ai+m +

aim_q1 +

alqm

=

O

(3.43) inevitable a pole occurs at x = 2. So the Padé approximants seems worse here.

If equation (3.43) is evaluated, the meaning of equation (3.41) is better un-derstood:

a E

O ii<O

qi

=1

q3

Oj>M

for instance the [1/1] approximant for eX reads as:

eX 2 + x

2+x

1 ₂ 1

2_=l+x+;x +

t...

The difference with the Taylor expansion begins at the O(x3)-term.

The Euler function evaluated at x = i is wildly oscillating ( when using Taylor

expansion). The answer however for E(i) is known: E(i)

0.5963. The

Taylor expansion never reaches this value. The Padé approximants however do. The [2/21 approximant reads as:

i + 5x + 2x2

[2/2](x) =

i + 6x + 6x2

So E(i)

0.6154 , only five Taylor terms have been used to get this accuracy. The Taylor expansion leads to the result:

Taylor5(1)

= 20. (3.46)

The [6/61(x) approximant is even better: E(i) 0.5968. Padé approximations

can also be used when calculating continued fractions, using the [M/M] and the EM/Al + 1] Padé approximants (see Baker [3]). For instance the Taylor approximation of the following function reads as:

(x+i+1)1+1

X X

i

2+

_r=1

4 16

and the continued fraction approximation is given by:

(3.44)

(3.45)

(3.47)

The advantage in mainly due to the fact that the series (3.47) has a range of convergence

r = 1,

whereas the continued fractions (3.48) have a

r

of the whole complex plane. This is true as long as x is not on the real axis for

x E (oo,-1] see figure (3.13).

So Padé approximant could lead to valuable results when Taylor expansion fails (see also Hermans and van Gemert [16]). In the next paragraph use will be

made of this when evaluating the single integral at y 0.

Figure 3.13: Range of convergence.

3.4.2

Padé approximations in the single integral

The single integral to be approximated is given by equation (3.22) repeated

here:

Single = 4iH(x)

_J

sec2ûe_Z5&20+ ixsecO+ ilylsec2OsinOdû

(3.49) Only the real part will be used:

_8JI(_x)Jsec2Oe_Z529sin(xsecO)cos(ysec2OsinO)d9

(3.50)

o

which will be denoted by:

changing of variable to s = secO leads to:

zs2

t se

f(x,y,z)

=

J

'/s

-

sin(sx)cos(ys/s2

-

1)ds

Taylor expansion of f(x, y, z) leads to:

af(x, y, z)

f(x,y,z) = f(x,O.z)1=o+y

a

1y°

+

ay2 Iy=°

applying this to (3.52) leads to:

00 2

se5

sin(sx)ds

.J2

_{- i}

12

f

3/2

-

1e2

sin(sx)ds + 0(y4)

(3.52)

(3.53)

(3.54)

Newman [27] uses only the first term of the Taylor approximations (the thin ship term) in his article. So expanding further we might be able to get more

accurate results.

Equation (3.54) is the formal Taylor series expansion. Numerical investigation reveals the fact that this series is not converging very well. Perhaps Padé

ap-proximation could lead to valuable results. Numerical experiments give raise to the following table (see table (3.1)):

Examine for instance the first row in this table. The Taylor(Newman) expansion of the function leads to a value of the calculated integral of 0.5179. Using Padé approximation leads to surprisingly good results in some cases.

Table 3J: Comparison between different approximation techniques.

(x,y,z)

f(x,y,z)

Taylor

Pade[2/2]

Pade[O/4]

(.1,.1,.1)

.30124 .5179 .3019 .3031 (.1,.5,.1) .0053 2.7 102 .0073 -.0018 (.1,.9,.1) .0029 1.1 10 .006 -.00012 (.1,.9,.9) .00252 .0044 .00254 .00255

(.5.1,.9)

.2132 .2132 .2132 .2132 (.5,.9,.5) .1608 2.11 .168 .177

Chapter 4

Solutipn of the Matrix

Equation

4.1

Introduction

In order to obtain accurate results for the forward speed problem a sufficient

number of panels for the discretization of the ship hull has to be used. In

most of our calculations we used 48 panels in the x-direction and 8 panels in the y-direction. So the total number of element will be 768. This is also

the matrix size. If we increase the number of panels, a direct solver (like

LU-decomposition) for the matrix equation cannot be used anymore. So some

sort of iterative method has to be used. In this chapter some suitable iterative methods will be discussed.

4.2

The matrix eqUation

In order to solve the ship wave problem the following two equations are impor-tant (see also equation (2.49) and equation (2.50)):

4ir(x)

= -

ff

a(e)G(x,)dS +

(4.1)

2xa(x) -

+ ... =

(4.2)

Once the potential çt(x) is known, the water height can be found according to equation (2.18):

K=

So at (A) we have:

N. N

K(i,j) = K(i +

₊

and at (B) we have:

K(i,j)=K(i-Fo the symmetrical problem the following block-structure can be observed:

(see also figure (4.1)).

(4.7)

(4.8)

(4.9)

(4.10)

ab e

f

b1\

g h b2

f

a b X3

-

b3 c d

b41

(o(z)

g

(4.3)

After discretization of these equations the following matrix equations are

ob-tained:

The system influence matrix for the determination of the source strength a

(equation (4.2)) is given by:

Ka=b

(4.4)

with a the source strength at panel number i. The potential ç can be written as (equation (4.1)):

(4.5)

with çj the potential at panel number i.

The wave height ((equation (4.3)) becomes:

(4.6)

( is the waterheight at the waterline panel. The matrix structure for the matrix

K looks like (with K the influence of a source at panel j uponpanel i.):

(

Figure 4.1: Symmetrical structure in the influence matrix.

with a right hand side vector of (b1, b2, b1,b2)T, the solution will be the vector

(x1,x2,x1,x2)T _{for the matrix now reduces to:}

( a+e

c+g

b+f

(x1

(b1

d+h )

_{X2 ) -}

b2

4.3

Convergence of the solution algorithm.

In order to investigate the convergence behaviour of the matrix the singular

values have been plotted (see figure 4.2). The singular values of a matrix K are the square roots of the eigenvalues of the matrix

KTK (a2(K)

=

)(KTK)).

From this figure a few things can be noticed:

(4.11)

The condition number (amax/crmjn) of this matrix is very small (so is very good), it is about 3.2.

25.0

20.0 15.0 -10.0 *

5.0

-t I I I I I I I 0 100 200 300 400 500 600 700 800

Figure 4.2: Singular values of Matrix K.

There is an enormous range in singular values which are almost the same. For instance the 150th singular value is 12.68 and the 170th is 12.66. We

think that this behaviour is typical for integral operators of this type.

A sudden change in singular values near the begin and the end of the

spectrum and slowly varying values in the middle. For a differentia

operator just the opposite can be observed: slowly changing values at the begin and at the end of the spectrum and a relatively rapid change in the middle.

Most of the problems we have been dealing with in this thesis have been solved using only a relatively small number of panels (about 800). If we increase the number of panels, a direct solver (like LU-decomposition) for the matrix

equa-tion K

= cannot be used anymore. So some sort of iterative method has to be used. Not only does the error increase rapidly using LU-decomposition,

but also the amount of work. LU-decomposition needs 0(n2) and iterative methods like CG, ART or SIRT only Q(3/2).

Because of the special structure of the singular values a method like CG could be most successful. The rate of convergence of CG is given by:

with C given by:

C=

+ i

and e1: Amax

cl =

Amin

But each time a singular value is approximated by CG, this value can be 're-moved' from the spectrum and the value of c1 can be adjusted. If there exists a large cluster of singular values which are almost the same, the convergence will increase rapidly. Some startup delay may occur, but this can be fixed

by preconditioning the matrix K (see also van der Vorst and van der Sluis

[36},[41]).

Due to the block structure of the matrix K (see figure (4.1)) the LU decom-position of K will almost have the same structure as the matrix K itself (see figure (4.3)). Of course this specific block structure of the matrix K only

occurs when the panels are generated as indicated in figure (3.3).

So the matrix equation

K .

= b can be preconditioned using only the lower

part of the matrix K. What you are doing is to approximate the matrix K

using the following:

K = LU LU

(4.15)

with L the actual lower triangular part of the matrix K and with U the identity matrix I. The matrix L is easy to invert because of its lower triangular form.

So instead of solving the original equation:

Ko=b

(4.16)

we now solve the preconditioned equation:

LK

=

(4.17)

This preconditioned matrix .Î1K should be much easier to handle using a CG

algorithm. A first approximation for the solution a is L1. In the case of a

symmetrical problem, the same preconditioning could be applied. The matrix (4.13)