Natural frequencies of free vertical ship vibrations

MATHEMATICAL INSTITUTE TECHNOLOGICAL UNIVERSITY DELFT

JULIANALAAN 132 DELFT

STUDIECENTRUM T.N.O. VOOR SCHEEPSBOUW EN NAVIGATIE

Netherlands' Research Centre T.N.O. for Shipbuilding and Navigation

SHIPBUILDING DEPARTMENT MEKELWEG 2, DELFT

*

NATURAL FREQUENCIES

OF FREE VERTICAL SHIP VIBRATIONS

(EIGENFREQUENTIES VAN VRUE VERTICALE SCHEEPSTRILLINGEN)

by

Ir. C. B. VREUGDENHIL

Issued by the Council This report is not to be published

unless verbatim and unabridged

This report is a joint publication of the Mathematical Institute

of the Technological University Deift and the Netherlands' Research Centre T.N.O. for Shipbuilding and Navigation.

page

List of notations 4

Summary 5

i Introduction 5

2 Basic equations 5

3 Some methods of solution ₆

3.1 Introduction 6

3.2 Systems of orthogonal functions 6

3.3 Energy-method 7

3.4 Method of direct integration 7

4 Hoizer-Myklestad method 8

5 Method of integral equations 8

5. 1 Integral equations ₈

5.2 Influence functions 9

5.3 Iteration-process 11

6 Applications 12

7 _{Discussion and conclusion 13}

8 References 14

9 Appendix 14

9.1 Hoizer-Myklestad method for a beam, divided into

prismatical parts 14

9.2 Theoretical investigation of the integral equation . 15

9.3 Expansion of a function into a series of eigenfunctions 19 9.4 Choice of the computation-method 20 9.5 Algol-program for the Hoizer-Mykiestad method . 21

9.6 Algol-program for the method of integral equations 23

List of notations

a(x) unit generalized coordinate

A(x) shearing stiffness B(x) bending stiffness G (x,u) influence function Ck transfer matrix E modulus of elasticity

j(x) eigenfunction

g1(x,u) kernel of integral equation G shear modulus

h step in difference equation H(x) Unit-step function

(1 x>Q

x<O

1(x) area-moment of inertia J (x) mass-moment of inertia k'(x) factor in shearing stiffness K(x.u) kernel of integral equation

length

m(x,t) bending moment

M(x) amplitude of bending moment q(x,t) shearing force

Q(x) amplitude of shearing force

s(t) generalized coordinate S(x) area of cross-section S(x,v) part ofK(x,u) t time T kinetic energy vk transfer vector V potential energy x space-coordinate y(x,t) deflection

Y(x) amplitude of deflection

= fx(x)dx = JJ()dX Dirac-function +

t5(x) =0 forx0

fò(x)dx= 1

-

w f i for k = I â Kronecker-symbol = _for k LI determinant p(x,t) angle of deflection

'(x) amplitude of angle of deflection

íe(x) for O <x < i

(x) = J(xl) _for I <x ₂₁

eigenvalue, j = w2

r number of elements

Q(x) mass per unit length w frequency

w. eigenfrequency

[] reference to list of literature

differentiation with respect to x differentiation with respect to t

NATURAL FREQUENCIES OF FREE VERTICAL SHIP VIBRATIONS

by

Ir. C. B. VREUGDENHIL

Summar)

The present investigation deals with some problems in the calculation of ship-hull vibrations. For practical purposes only the determination of natural frequencies of free vertical vibrations is considered.

Especially one assumption, that the influence of the mass-moment of inertia is negligible, is investigated. Therefore a number of computation-methods is considered on their possibility of taking that influence into account.

Two methods appear to be useful:

(I) the Holzer-Myklestad method (Prohl's method);

(2) a method using influence functions, that has not yet been applied to ship vibration.

Computations with these two methods were carried Out on an electronic computer for three cases: a uniform beam, the SS. "Gopher Mariner" and the MS. "Naess Falcon". This made clear that the Holzer-Myklestad method requires a division of the ship into about 60 elements to obtain a reasonable accuracy. The method of influence functions gives satisfactory results when using only 20 elements.

For the MS. "Naess Falcon" computations were carried Out, both taking the influence of the mass-moment of inertia into account and neglecting it. There was no significant difference between the results of these cases, so that the neglection is permitted.

Besides, as could be expected, it was found that beyond the fourth or fifth mode the approximation of the ship-hull by a ,,Timoshenko-beam" is no longer justified.

i

Introduction

q(xt)

In the literature on ship-hull vibration several Orn

U m(xt) x x+x m+yax I x-axis

methods of computation have been proposed.

Usually they are based ori the assumption (among others) that the influence of the mass-moment of

Ox

inertia is negligible. It appears useful to investigate

the validity of this assumption and the merits of y-axis some methods of computation. Because of the

many difficulties in this problem, among which

the uncertainty about the exerted forces; It is subject to

the uncertainty about the damping factors; _{a. shearing forces -q(x,t) on the section x and} the coupling between horizontal and torsional

Lx on the section x+x.

vibrations, q(x, t) +

a

this investigation is confined to the computation

of natural frequencies of free vertical hull vibra- b. bending moments -m(x, t) on the section x tions. The problem of the added mass caused by am

the entrained water is not treated here. and m(x,t) +axAx on the section x+Ax. Our procedure will be to find methods of corn- _{Then the equations for translation and rotation}

putation that take the influence of the mass-

_{of the element are:}

moment of inertia into account and to investigate

that influence experimentally, i.e. by means of a _translation:

a(xt)

(X)AXa2(xt)

. (a)

computation. ax at

where (x) = mass per unit length

2

Basic equations

y(x,t) = deflection of the beam

t = time

As a mathematical model of the ship we take an elastic beam with the same properties in bending

[1], [2]. Consider the part of the beam between

x and x-f-Ax.

Figure 1. Forces and moments in the beam

rotation:

am (x t) a (x, t)

Axq(x,t)Ax = J(x)Ax

. . (b)

6

where

J (x) = mass-moment of inertia of the section

relative to the neutral axis; q(x,t) = angle of deflection.

From the elementary theory of bending [3] we

know: (x, t)

m(x,t) = B(x)

(c) ax

J(x,t)

q(x,t) = A(x) (x,t)} . (d) I a where

B(x) = bending stiffness E1(x);

A(x) = shearing stiffness /c'(x)S(x)G;

E = modulus of elasticity;

1(x) = area-moment of inertia

of the section

relative to the neutral axis;

k'(x) = factor, denoting the influence of the shape of the section on the shearing stiffness; S(x) = area of the cross- section;

G = shear modulus.

We consider simple harmonic vibrations with

fre-quency O) and introduce amplitudes of the

varia-bles in the following way:

m(x,t) = M(x) cos wt etc.

One of the purposes-- of this investigation is the

determination of the possible values of the

fre-quency, i.e.

the natural frequencies or

eigen-frequencies of the system.

The equations (a)(d) become: dQ = _w2(x)Y(x) (2 1) = Q(x)+w2J(x)0(x) dx dO M(x)

dx -

B(x) dY _O(X)+Q(x) (24) dx A(x)

which we will call

the basic equations. The

corresponding boundary conditions are those of a beam with free ends:

M(0) = M(l) = O

Q(0) = Q(l) = o

This system of equations is only valid for small

amplitudes of the vibrations because of the

assump-tion on which (c) and (d) are based: that a plane

section of the beam remains plane in bending.

3 Some methods of solution

3.1 Introduction

One of the main difficulties in solving the

equa-tions (2.l)(2.5) is the fact that the frequency (J

appears in a complicated way. This will be clear when we derive a differential equation for 0(x):

d 1 d21

dol

-

B(x)

í+

dx o(x)dx2t dx) d

ldI

I

jJ(x)0(x)j +

dx (x)dxt

o2 dl

_dol

+

B(x)--1 +

A(x)dxk dxj w4J(x) 0(x)

Ax

= o

. (3.1.1)

The term with w makes it very difficult to

deter-mine the eigenfrequencies from this equation.

Differential equations for Y(x), M(x) or Q(x) show

the same feature and it is also impossible to derive a simple integral equation for one of the variables. We will now investigate some other methods of solution, but we omit methods that confine them-selves

to special problems. Examples of such

methods are those, used by KUMAI in his papers

[4] and [5]. He imposes a certain shape to the

distributions of mass and mass-moment of inertia.

In the following pages we will show that this is

not necessary for our purpose.

3.2 Systems of orthogonal functions

The system of differential equations (2.l)-2.4)

can be converted into a set of algebraic equations

by expanding the variables in a series of orthogonal

functions (a system of functions Xk(x) is called

orthogonal if

/Xk(x)Xm(x)dx = Ckòkm (3 2.1)

where Ck is a constant and ôjtrn is the

Kronecker-symbol).

It is possible to use sines and cosines for %k(x)

and construct Fourier-series for Y(x), 0(x), A-1(x)

and Q (x). This has several disadvantages:

it is difficult to determine the number of terms required for a certain accuracy;

it is difficult to construct Fourier-series for

terms like (x) Y(x), using the expansion for

generally the number of equations, resulting

from this procedure will be very large, viz. 8N

when N terms are used in the Fourier-series

(there are four variables with 2N coefficients

} (25) dx dM (22) (2 3)

each). It requires much work to determine the

eigenvalues of such a system.

Another method has been proposed by RICHARDS

[8]. His system of orthogonal functions consists of the natural modes of vibration of a homogeneous

beam. This method not only has the same

dis-advantages as the preceding one, but also

the expansion of terms like dQ/dx is difficult;

it is necessary to derive criteria for the

con-vergence of a series of these functions.

Because of these problems we will not use this

method.

3.3 Energy-method

Because of the influence of the mass-moment of inertia the well-known energy-method results in

a very complicated system of equations. This is

shown as follows.

The kinetic energy of a vibrating beam is [7]:

T =

1/2/TÇ(x) ()2+J(x)

(a)2}

dx (a)

and the potential energy (strain energy):

1m2(x j

_{2x t1}

V = /2 /

'

+

' ' dx. . . . (b)

i, B(x) A(x)

o

As usual in this method we introduce generalized

coordinates s1(t), i.e.

y(x,t) =

a1(x)sj(t) (c)

i=O

where the functions ai(x) satisfy the boundary conditions for y(x,t) and form an orthonormal system (i.e. satisfy

(3.2.1) with all Ck=l). By

using equation (2.4) and the fact that we

con-sider harmonic vibrations, it is possible to express

(x, t) in terms of si(t) X ay

lfa2y

(x,t) =

_{- Açx)Jat2 =}

X 1'(x)sl(t) -

=

X { w2 1 aj'(x)+ f(x)ai(x)dxJts1(t) . (d) A(x) i=0 0

where ' means differentiation with respect to x

andO with respect to t. In a similar way m(x,t) and

q(x,t) are treated and so (a) and (b) are

trans-formed into

-{h15 + w2kj j + w11s}J (t)

{h15 + (02k15 + (,)4l11}c02S5 (t)

j0

because s5(t) is supposed to be a harmonic func-tion. Similarly

av

as1 15 + w2n5 + w4rjj}sj (t) 5=0

so that (e) can be written as

{m15+o2(n15-h15) +w4(rij-kij) _w615}sj(t) = O

J0

_{(i = 0,1,2,..)}

It will be clear that it is very difficult to determine the eigenvalues w2 of the system. For that reason we will not use this method.

3.4 Method of direct integration

In the simple case that there is no influence of

shearing forces and rotational inertia (so 5(x) 0

and A(x)oo) it is possible to derive:

X U V W

Y(x) = G, xC2 +

w2fdu/)/dwf

(z) Y(z) dz (a)

with boundary conditions

/(x)Y(x)dx = O

(b)

fx(x) Y(x)dx = O

For this system there exsts a simple

iteration-method, known s the Stodola-method or method of direct integration [2], [8]. It is possible to take

into account the influence of either shearing forces

or rotational inertia, but not of both. For this

reason we will not consider this method anymore,

but we will derive a generalization of it in section 5

of this report.

T= /2

h15 + oi2kj +D4l1i}3'1 (t).tj (t) 1=01=0 co co V = '/2 mj5 + w2njj + w r iJ}s1 (t) s5 (t) i0 j=0

where h, k, i, m, n, r are constants, determined by

a(x) and the properties of the system. The

com-putation of eigenvalues is based on the Lagrangian equations dt as1

) + =

i as1 0

(i=0,l,2...(e)

d jaj

av

aT

Now

d faT

50

i

dt ia1 I

8 Vk =

(0

_\

O v±i Yv+i!

=

0\

O Yo! (4 1.7)

This is a system of four equations, the first two of which are:

Ci3*cIo+C14*Yo = O C2i*Jio±C24*Yo = O

This can be solved only

if the determinant

= Ci3*C2i_Ci4*C23* is equal to zero. The

eigenvalues are those values of w2 for which this condition is satisfied. They are determined by an

iteration-process.

It is not very sensible however to execute the

process as described here, because the matrix C

depends on w and consequently has to be

cal-culated in each iteration-step again. This requires theoretically 64v multiplications per iterationstep.

This can be reduced to 32v by applying the follow-ing procedure.

First take Vo = (0,0,0,1) and calculate

using (4.1.5) repeatedly. According to (4.1.6)

= (C14*, C24", C*, C44*). In the same way

we determine C13* and C23* by repeating the calculation with Vo = (0,0,1,0). Then z can be

calculated. The zero's of the thus defined function (w) are determined by means of a simple linear interpolation method, known as the "regula falsi".

This method has been programmed in Algol for

a computer. The program

is reproduced in

appendix 9.5.

5

Method of integral equations

As far as the author knows, the method of integral equations has never been used in

ship-hull-vibra-tion calculaship-hull-vibra-tions, at least not in the comprehensive

form that is presented here. One simple method

has already been mentioned in section 3.4, another

is suggested in KocH's book [10], but without real

applications. The method can be formulated in

two different ways, that will be shown to be

equivalent.

Ji

Integral equations

The first way is the construction of a system of two

simultaneous integral equations for Y(x) and From (2.1) and (2.2) we derive:

dM

-w2/(u) Y(u)du+ w2J(x) (x)

Integration with respect to x gives (because

M(0) = 0): dx M(x) = - w2fdvf (u) Y(u) du + w2[J(v) (v) dv = X X = _U)2/ (xv)(v) Y(v)dv+w2/J(v)(v)dv ô d

Divide by B(x) and integrate again, now using

(2.3):

4

Hoizer-Myklestad method

The Hoizer-Mykiestad method (sometimes called Gümbel- or Prohi-method) is based on division of the mathematical model into small parts. This can be done in several ways. Among them are:

a method in which the system of differential

equations (2.l)(2.4) is replaced by a system of difference equations [2], [8]. In this section we will explain this method.

a method in which the mathematical model is

replaced by a beam that consists of a finite

number of prismatical parts [8], [9]. In

appen-dix 9.1 we will show that this method is in-ferior to the preceding one, because it

intro-duces a relatively large systematic error.

In method (a) we divide the interval

[0,1] in v+l equal parts with lengths h = l/(v+l). When

we introduce the notation Y(kh) = Yk etc. and use the approximation , the equations

dx h

df

_fk+ifk

(2.1)(2.4) can be written as:

(41.1 fl 91. V k+1 = kW1L(k1 k (4.1.2) Mk+1 = Mk+hQk+w2hjkbk.

k+i = Ii1k_B

(4 1.3)

= Yk+hk+

(4 1.4)

For convenience we introduce a vectornotation as

follows:

¡1 0 0

hw2k

h i hw2Jk O

O h/&1

O

\h/AkO h i

Then Vk+j = Ckvk (Ic = 0, 1,2, .. .,i') . (4.1.5) By applying this rule repeatedly we obtain

= C'v0

(4 1.6)

where the matrix C* depends (among others) on w. Inserting the boundary conditions (2.5) gives:

X W dw

f(wv)(v)Y(v)dv+

(x) = B(w) o b X W

(02

fJ(v)(v)dv

(a)

jB(w)

o

This is the first integral equation. The second one

is derived from (2.1) and (2.4) with Q(0) = O:

X IdY } = 2J;(u)Y(u)du A(x) dx o or dv

Y(x) = Y(0)

+fv)dv_wz/

/(u) Y(u)du

A(v)ì

o o ò

Substitution of (a) and partial integration with

respect to r gives: Y(x) = + w2

dvJ)

J (w - u) (u) Y(u) du + ' dw w2 /dv /

fJ(u)(u)du+

.j jB(w) o o o

rdv

/(u)Y(u)du - Wz I .1 A(v).o

=

V

r(xv)dv f

+ w2J

B(v)

J{ (vu)(u) Y(u) J(u)(u)}du+

o o

X V

rdv

"

(02 f f(u) Y(u)du (b)

ô o

This is the second integral equation.

The equations (a) and (b) can be written in a

more intelligible way by reversing the order of the integrations:

Y(x) = Y(0)+x(0)+

+gi(x,u)(u)Y(u)+gz(x,u)J(u)(u)}du

(5.1.1) (x) =

X

+2/{g3(x,u)(u) Y(u) +g4(x,u)J(u)(u))du (5.1.2)

where gi(x,u) = g2(x,u) g3(x,u) g4(x, u) and X

¡Ç(x_u)(vu)

I J 1 B(v) A(v)

f(vx)dv

J

B(v) [(v_u)dv B(v) ç. dv

J

B(v) A = (02 dr (5.1.3)

The functions g(x,u) are not defined for u > x.

The two remaining boundary conditions are

Q(1) = O and M(1) = O or

They can be identified with the conditions for

equilibrium for the system.

The present formulation of the problem has two important advantages:

The functions gi(x,u) depend only on the stiff-nesses A(x) and B(x) and not on the mass and mass-moment of inertia. This means that they

have to be computed only once for a certain

ship and that they can be used then regardless

of the loading-conditions.

the main unknown A appears in such a way

that it is possible to apply a simple

iteration-process for the determination of the natural

frequencies. This will be shown in one of the

following sections. 5.2 Influence functions

The equations (2.l)(2.4) have the same shape as those describing the static bending of the beam when subjected to an external force A(x) Y(x)

per unit

of length and a

bending moment

AJ(x)(x) per unit of length. This will be described

now by means of influence functions. The boundary

conditions have to be adjusted in a somewhat artificial but not unusual way: suppose that the

beam is clamped in x = O and give the point

where the beam is clamped a translation and

rotation such that the boundary conditions for the whole system are satisfied (figure 2).

/(x)Y(x)dx = O

(5.1.4)

lo

Y

o

This means:

Y(x) = Dix+D2+

+J{Ci(x,u)(u) Y(u) +C2(x,u)j(u)(u) }du

Di+

-DEFLECTION OF CLAMPED BEAM

Figure 2. Influence functions

ò(x-u) is the Diracfunction (cf. list of notations). The case (r = 1, s = O) results in a Y(x) and (x) (5.2.1) representing Ci(x,u) and Ca(x,u); the case (r = O, s = 1)

results in a Y(x) and

(x) representing

C2(x,u) and C4(x,u). The integration of the system (a), (b), (c) is a simple process, resulting in

(a)

x Y(0) = (0) = O (b)

M(l) = Q(1) = O

(e)

describing the clamped beam with a force r arid a bending moment s at the point u (figure 3).

>x Figure 3 u Ci(x,u)

I(xu)(uv)dv

=

H(xu)+

B(v)

[(xv)

(uv)dv

B(v)

H(ux)+

o ±J)H(xu)

±f)H(u

x)

C2(x,u) = u X

[(xv)dv

'(xv)dv

=1 B(v) H(x_u)+] B(v) ô o

C(x,u) =

u X

,[(uv)dv

(5 24)

/'(u_v)dvH()

B(v) H(ux) B(v o C4(x,u) = X

[dv

[dv

J

I

B(v) o O

H(xu) represents the unit-step function (cf. list

of notations).

Apparently there is some relation between gi (x, u) and C1 (x, u); thisrelation is found to be:

Ci(x,u)

_{gl(x,u)H(xu)gl(x,O)ug2(x,O) i}

C3(x,u) = ga(x,u)H(xu) g3(x,O) ug4(x,

C2(x,u) = g2(X,u)H(Xu)g2(X,O)

(d)

C4(x,u) = g4(x,u)H(xu) g4(x,O)

We substitute this in-to (5.2.1) and (5.2.2) remem-bering the boundary conditions (5.2.3). Then the terms containing gl(x,O) vanish and so we obtain

(5.2.5)

+Af C3(x,u)(u) Y(u) +C4(x,u)J(u)(u) }du (5.2.2)

o

where D1 and D2 are the rotation and translation

of the point where the beam is clamped. They

are determined from [(x)Y(x)dx = O

(5.2.3)(5.l.4)

/{J(x)(x)+x(x)Y(x)}dx = O

The influence functions C1(x,u) are defined as follows:

C1 (x, u) = deflection at x caused by a

unit force at u;

C2(x,u) = deflection at x caused by a unit moment at u;

C3(x,u) = angle of deflection at x caused

by a unit force at u;

C4(x,u) = angle of deflection at x caused

by a unit moment at u.

Of course all deflections and angles of deflection

are relative to a clamped beam. Considering these definitions, the influence functions can be found

from the equations dQ dx

= râ(xu)

dM dx

= Q+sò(xu)

d M(x) dx B(x) dY Q(x) dx

=

<u<i

ZuZ21

21, 0 u < i

< 21, i < u 21

(5.2.8) the two equations (5.2.6) and (5.2.7) are

com-bined to one: f(x)

=

[K(x,u)f(u)du (5.2.9) Ki(x,u) for fK2(x,u_1) for K(x,u)

=

ÌK3(x-1,u) Cor

1K4(xi,ui)

for 5.3 Iteration-process

Unfortunately the kernel K(x,u) is not symmetric. According to the theory of integral equations it is then generally impossible to prove the existence

of eigenvalues. However in appendix 9.2 the

following theorem is proved:

The kernel K(x,u) has an infinite number of

single positive eigenvalues Ai with

correspond-ing eigenfunctionsj(x) that satisfy the relation

(9.2.8):

(x)

₌

(53.1)

and the orthogonality-relation

21

frn(x)f1(x)fi(x)dx = (5 3.2)

_Í(x)

forOx<l

where (x)

- tJ(xi)

for i < x 21

Using these properties we can construct and justify

the following iteration-process [12].

Assume that the eigenvalues are ordered

accord-ing to

Ai<A2<A3

. (53.3)

Take an

initial approximation f(°)(x)

of the

desired eigenfunction; f (1 (x) must satisfy the

boundary conditions. Then according to section 9.3

f(0)(x)

=

b1fi(x)

=1

If y eigenfunctions are known already, compute

21

b1 =/i,(x)f(0)(x)JTj(x)dx and subtract their

con-tributions bjfj(x) fromf(°)(x), so that

f(0) (x) =2Jifi(x)

In the real computation this elimination has to be

repeated after each iteration-step, because the

approximations that are used, can introduce small contributions of fi(x) . . .j,(x) to f(°)(x).

Now starting withf(°)(x) we compute a sequence

of functions f(n)(x) from

21

f(z)(x)

=

fK(x,u)fn_1(u)du Then from (5.3.1) it will be clear that

f(n)(x)

=

j=v+l

and because of (5.3.3) after a sufficiently large

number of steps:

f(n)(x)

bif()

so by norma1izingf(z)(x) we find the new

eigen-function f-1-i(x). Moreover from the sequence

f ()(x) the eigenvalue A±i can be computed. One

method is: b1 jv+1 b1 f (n) (Xi)

-fi(xi)

i= V+ i i

b1 (A+i-1 fi(xi)

f (n-1) (xi) 1+ b+1\A1

=

'v+1 2+\ b+2 /'2l,1\n_hfv+2(xi)

1tl+(1_A

)b

f1(x1)J

b1 (av+i)n A1

(5.1.1) and (5.1.2) again. This means that the two methods are equivalent.

For our theoretical investigations the system (5.2.l)(5.2.3)

has some advantages over the

system (5.l.l)(5.1.4) of the preceding section.

As we mentioned already, the constants D1 and

D2 are determined by substituting (5.2.1) and

(5.2.2) into the boundary conditions (5.2.3). This

results in a system of two equations for D1 and D2.

The solution is substituted back into (5.2.1) and

(5.2.2), with the result:

Y(x)

=

2/{Ki(x,u)Y(u)+K2(x,u)t1i(/)Idu (5.2.6)

(x)

=

A/{Ka(x,u) Y(u) +K4(x,u)(u) }du (5.2.7)

O

For the derivation of these formulae we refer to

appendix 9.2.

By the substitutions

JY(x)

forOx<l

f(x)

=

12

so the quotient off (')(xi) and f(n)(xi) gives )±i

with a relative accuracy of approximately

\,a2)

A better method is to form the so-called

Ray-leigh-quotient f(x)f (n-1) (x)f (n) (x) dx i=v+l f(x){f(n) (x) }2dx i=v+1 1+

1+

b12 (2+i)2n_1 b+12\ 1 v-l-1 b2

/

i ')1 'v+1

\ /i

Al2n-1 b12 2n j=v+1 ° 7.a \2n-fl v+1\

so in this way we find .ar+1 with a relative accuracy (A+i\2fl_l

of approximatelyt I

From experiments we know that very roughly

i2Ai. This means that using the same number

of iterations n, the two methods give an accuracy

of

(v+

1' 2n-2 ' 4n-2

+

2) and + 2 respectively, the last

of which is

preferable.

After this procedure the whole process can be repeated to compute v±2 and so on.

An Algol-program, based on this method (but

differing slightly from it; cf. appendix 9.4) is given

in appendix 9.6. The number of iteration-steps is

found from j2, _{as follows:}

an accuracy of 1% is obtained when ¡A._1\2n-1

/i- 1'\'

6 Applications

The computing-methods, described in sections 4

and 5.3 have been used to compute the natural

frequencies in three cases: a homogeneous beam;

the SS. "Gopher Mariner" with data taken

from the paper [13] by MCGOLDRICK & Russo; the MS. "Naess Falcon", the data of which have been provided by the Netherlands'

Re-search Centre TNO for

Shipbuilding and

Navigation.

The data are reproduced in appendix 9.7.

The computations have been carried out on the

Electrologica-Xl computer of the Netherlands'

Ship Model Basin, Wageningen, The Netherlands.

TABLE I. Frequencies for the homogeneous beam, computed

by the Hoizer-Mykiestad method

TABLE II. Frequencies in Hz. computed by

Holzer-Mykiestad method, using 20 elements, compared with

measured frequencies

)

<0.01

For the first five eigenvalues n

=

i+ I is sufficient.

TABLE ¡ II. F eclucncies in Hz computed by method of integral equations using 20 elements, compared with experimental values

Remark: for future computations the computing time can be reduced considerably, because in thc present ones some super

fluous work was included.

of mode

number of elements used

theor. [6] 20 30 40 0.36 2 1.09 1.03 1.01 0.98 3 2.39 2.11 .... 1.93 4 4.64 3.72 . .. . 3.18 no. of mode

"Gopher Mariner" "Naess Falcon"

comp. meas. comp. meas.

1.39 1.33 1.03 ? 2 2.90 2.57 2.31 2.08 3 4.91 3.47 4.25 3.28 4

....

4.25 6.05 4.22 5

....

4.92

....

4.92 no. of mode

homogeneous beam "Gopher Mariner"

"Naess Falcon"

heavy condition light condition

-- computed

comp. theor. comp. meas. meas. comp. meas.

J0

J=0

0.35 0.36 1.48 1.33 1.00 1.01 1.08 0.90 2 0.96 0.98 1.90 2.57 1.82 1.86 2.08 1.99 2.30 3 1.85 1.93 3.88 3.47 3.27 3.30 3.28 3.72 3.68 4 .. .. 3.18 5.34 4.25 4.50 4.52 4.22 5.10 4.62 5 .. .. 4.74 . . . . 4.92 5.76 4.92 6.31 5.13

For the MS. "Naess Falcon" it was possible to

use data for the mass-moment of inertia so that its influence could be determined.

Fig. 4. Theoretical and computed frequncies for the uniform beam (table III)

7

Discussion and conclusion

As the tables I and II show, the results of the

com-putations with the Hoizer-Mykiestad method are

far from satisfactory when using 20 elements.

Errors up to 50% are obtained in this way. Only

for the homogeneous beam it was possible to

carry out a computation with more than 20

ele-ments. The results (table I) suggest that a number of about 60 elements might be satisfactory. This

is in agreement with the views, CSUPOR expresses

in his paper [8]. It has not been possible to show

no.of mode

2 3 4 5

no.of mode

Fig. 6. Measured and computed frequencies for the "Naess

Falcon" (heavy condition) (table III)

The results of the computation are reproduced in the tables no. I, II and III and in the figures

4, 5, 6 and 7. Hz 6 f req. 5 4 3 2 .5 4 no.of mode

that the inaccuracy of the method is caused by

the inaccuracy of the difference-formula that was

used to derive (4.l.1)(4.l.4). An attempt to

im-prove the process by using the more accurate

formula

df

_{fk+1fk_1}

+ 0(112

dx 2h

failed; in this case the method did not even

converge. Although this point was not investigated

more completely, it is possible that it is caused by instability of the method; a similar experience has

- -..-measured computed

Pl

1A

-.measured -o--comp uted o--.- measured computed 2 3 4 5 6 no of mode

Fig. 7. Measured and computed frequencies for the "Naess Falcon" (light condition) (table III)

Fig. 5. Measured and computed frequencies for the "Gopher Mariner" (table III)

Hz freq. 5 4 3 2 Hz 6 f req. A 5 4 3 Hz 6 f req. A 5 4 3 2

14

been obtained in some techniques for solving

partial differential equations.

Although consuming more computing-time, the

method of integral equations gave much better

results (table III). The division into 20 elements is sufficient here. For the MS. ,,Naess Falcon" an accuracy of 15% or better is obtained. Probably it

is not possible to obtain a better accuracy using

the present mathematical model.

An important fact is shown by the computation

for the heavy loading-condition of the ,,Naess

Falcon" with and without influence of the mass-moment of inertia. The results do not differ signi-ficantly so that we may conclude: the influence of the mass-moment of inertia is negligible.

This means that it is not necessary to use such

elaborated methods as the two that have been

used here. Some methods mentioned in section 3

are possible now.

It

will be clear anyway that a

frequency-computation on an electronic computer is worth wile. With one of the described methods it will be

possible to compute five frequencies in a few

minutes.

Computation of more than five frequencies is

possible but the obtained accuracy will be poor. This is suggested by the figures 6 and 7 and is in agreement with the ideas, expressed in the litera-ture. For the computation of higher modes other

methods must be found. An attempt in

this direction is given by GREENSPON [14].

Acknowledgements

The author wishes to express his gratitude to

Prof. Dr. E. VAN SPIEGEL for his direction and

important suggestions during the investigation,

and to the Netherlands' Ship Model Basin for the possibility to use the computer.

8 References

TIMOSHENKO, S.: Vibration problems in engineering,

New York 1955.

9 _{MCGOLDRICK, R. T.: Ship vibration, David Taylor}

Model Basin (Department of the Navy), Report 1451 (Dec. 1960).

TIMOsHENKO, S.: Elements of strength of materials, London 1935.

KUMAI, T.: Vertical vibration of ships, Rep. Research

Inst. Appl. Mech. Kyushu Univ. (Japan) Iii, 9 (April 1954).

KUMAT, T.: Shearing vibration of ships, Rep. Research Inst. Appi. Mech. Kyushu Univ. (Japan) IV, 5

(Jan. 1956).

RICHARDS, J. E.: An analysis of ship vibration using

basic functions, Trans. N.E. Coast Inst. Eng. &

Shipb. 68 (1951) pp. 51-92.

LANGHAAR, H. L.: Energy methods in applied

me-chanics, New York 1962.

CSUPOR,D.: Methoden zur Berechnung der freien

Schwingungen des Schiffskörpers, Jahrb. Schiffbaut. Ges. 50 (1956).

BIGRET, R.: Quelques examples d'application du calcul électronique en méchanique industrielle, Ing. de

l'automobile 11(1961) pp. 584-595.

KocH. J. J.: Eenige toepassingen van de leer der eigen-functies op vraagstukken uit de toegepaste mechani-ca, Thesis Delft, 1929.

li. SCHMEIDLER, W.: Integralgleichungen mit

Anwen-dungen in Physik und Technik, Leipzig 1955. V0OREN, A. I. VAN DE: Theory and practice of flutter

calculations for systems with many degrees of free-dom, Thesis Delft 1952.

MCGOLDRICK, R. T. & Russo, V. L.: Hull vibratíon

investigation on SS. Gopher Mariner, Trans.Soc.

Naval Arch. Eng. 63 (1955) pp. 436-494.

GREENSPON. J. E.: Theoretical developments in the vibration of hulls, J. Ship Res. 6,4 (1963) pp. 26-47.

Mir /

I -1H ELEMENt

Q,, Figure 8

To facilitate the computation it is assumed that the inertia-termw2eiY of a part does not work on the

center of gravity, but on the left end section (figure 8).

Then the transfer from the right end of the (i l)-th part to the left end of the i-th part is described by

Q.z = Q1_1.r+co21Y1.1 (a)

M1,1 = M1_î,r+w2ji1.j (b)

while 1'i1,r (c)

i-1,r (d)

9 Appendix

9.1 Holzer-Myklestad methodfor a beam, divided into prismatical parts

Divide the beam into prismatical parts [9] with

length h1

mass-moment of inertia Jt

mass

bending stiffness B1

[e(x)Y(x)dx = O /{j(x)(x) +x(x) Y(x) }dx = O (x) = D2 _{H-)/{C3(x,u)e(u)Y(u) +C4(x,u)J(u)(u)}du} o h12_Qi,z B1 1 '2 B1 (g) (52.2)

Temporarily we will call f 13(X) = fGi,3(x,u)(u)Y(u)du F2,4(x) = AfC24(x,u)J(u)(u)du

o

Substitute (5.2.1) and (5.2.2) into (5.2.3). This gives

ß0D1

+ßiD2 =

-/' (y){F1_{(o) +F2}(o) }dv

O

ßiD1+{yo+ß2JD2 = _[v(v){Fi (y) ±F2(o) }dv _/(v){F3 (o) +F4(v) }do

where

ß =fxi(x)dx and

y

=/'J(x)dx

With = ßo(yo+ß2)ß12 we find forD1 and D2:

where the subscripts I and r denote left and right respectively. The transfer from the left end to the right end of the i-th part is now a matter of a prismatical beam that is loaded at the ends only.

the shearing force is

Q = Qi,l, so _{Qir = Qi,i} (e)

the bending moment is

M =

M11+xQ1z so Mir = Mi i+hiQ1 z (f)

according to (2.3) the angle of deflection is

XJlljj Qi.t

=

-

1/2 X- SO i,r =

Dj according to (2.4) the deflection is

AI1 z x3Q1 j xQj i

Y = Yí,l+Xl,l1/2X

B1 6B1

+

A1

M

1h h.31

SO Yi,r = Yi,z+híi,i_h/2hj2

+

-

Qij (h)

B A1 6B1;

From the equations (a)(h) it is possible to construct a vector-equation similar to (4.1.5). However the shifting of the inertia-forces to the left end of the parts introduces a serious error. This error results in a

moment of approximately

1/2 hiw2iYig

(where Y1 is the deflection at the centre of gravity) and this is an error of order h1 whereas the error

in (4.1.5) is of the order h2. This means that the present method is inferior.

9.2 Theoretical investigation of the integral equation

First we will show the procedure that leads to the equations (5.2.6) and (5.2.7). We have the equations

Y(x) = Di+xD2+A/{Ci(x,u)(u)Y(u)+C2(x,u)J(u)(u)}du (52.1)

o

j

16

D1

₌

-f{ (o+ ß2) vßj}{Fj (u) +Fz (u) } (u) dv + ißiJ(v){F3 (u) +F4(v) }dv

D2

=

+F2(v) }dvo/j(v){F3(v) +F4(v) }dv

Substitute these expressions back into (5.2.1) and (5.2.2):

Y(x) = Fi(x)+F2(x)-f{yo+ß2-(v+x)ßi+xvßo}e(v){Fi(v)+F2(v)}dv+

±k/(ßl_xßo)J(v){F3(v) +F4(v) }dv

(x) = Fa(x)+F4(x)-f(ßov-ßi)e(v){Fi(v)+F2(v)}dv+

_/'oJ(v){Fa(v)

+F4(v) }dv

This can be rearranged in the desired form

Y(x) = ¿/{Ki(x,u)Y(u)+K2(x,u)(u)}du (92.1)

(x) = )/{K3(x,u)Y(u)+K4(x,u)(u)}du (92.2)

o

where ÀfK1(x,u)Y(u)du

_{= Fi(x)/yo+ß2(x+v)ßißoxv}e(v)Fi(v)dv+}

+kf(ßlxßo)J(v)F3(v)dv = ÂiCi(x,u)(u) Y(u)du+

/{yo+ß2 (x+v)ßißoxv}o(v)/Ci(v,u)(u) Y(u)dudv+

+f(ß - ßox)J(v)fC3 (v,u) (u) Y(u) dudv =

(reverse the order of the integrations)

=

À[Ci (x,u) (u) Y(u)du

4[Y(u)

_{du/Ì(ßi xßo)J(v)e(u) Cs(v,u) +}

{yo+ß2 (x+v)ßi +ßoxv}(v)(u) C1 (v,u)]dv

so K1 (x,u) = Ci(x,u) (u)

_+k

/Ì(ßi xfio)J(v) (u) C3(v,u) +

{yo+ß2 (x+v)ßi+ßoxv}o(u)(v)Ci(v,u)]dv (9 2.3)

Similarly

Kz(x,u) = C2(x,u)J(u)

+

{yo+ß2 (x+v)ßì+ßoxv}(v)J(u)Cz(v,u)]dv (9 2.4)

Ku(x,u) = C3(x,u)(u)

_+ki

[(ßi-vßo)0(v)0(u)Ci(v,u) -ßoJ(v)(u)C3(v,u)]dv . . . (9.2.5)

Unfortunately the kernel K(x,u) as described in (9.2.3)(9.2.6) and (5.2.8) is not symmetric. In such

cases it is usually impossible to prove the existence of eigenvalues. However the kernel has a favorable shape viz.:

The validity of this expression is easily shown by working back from (9.2.7).

When we investigate the kernel C(x,u) we find that C(x,u) = C(u,x) so C(x,u) is symmetric. With this

knowledge it is possible to apply a theorem from SCHMErDLER'S book [11] (p. 317) that says:

If a kernel K(x,u) is ,,left-sided symmetrisable" (,,Iinksseitig symmetrisierbar") i.e. there exists a symmetric positive definite kernel D(x,u) for which

21

H(x,u) = /D(x,v)K(v,u)dv

is a nonzero symmetric kernel, then K(x,u) has (generally an infinite number of) single

eigen-values A1 with corresponding eigenfunctionsj)(x) and gi(x) that satisfy the relations

J(x) = A1/K(x,u)jì(u)du (9 2.8)

gi(x) = AifK(u,x)g1(u)du (9 2.9)

g(x) =fD(x,u)fj(u)du (9 2.10)

21

K(x,u) = C(x,u) r (u) + /S(x,v) (y) C(v,u) ì (u) dv (92.7)

o

where C(x,u) is composed from Cj(x,u) in the same way as K(x,u) from K1(x,u) (cf. (5.2.8));

0<u<1

(u)

= IJ(ul)

i < u < 21

{yo+ß2(x+v)ßi+xvßo}

0< x <i,

0< r <i

(ßixßo}

O<x<l,

1<v<2i

S(x,v)

=

-(ßivßo}

1<x<21, O<v<l

ßo

l<x<21, l<v<21

and the relation of bi-orthogonality

21

/J(x)g5(x)dx = (92.11)

If the number of eigenvalues is n, the kernel can be written as

K(x,u) = fi(x)g(u) In our case we put

D(x,v) = (x)C(x,v)i(v)

21

Then H(x,u) = /D(x,v)K(v,u)dv = o

21 2121

=f? (x) C(x,v) (r) C(v,u)i1 (u) dv+f/ (x) C(x,v) (v)S(v,w) (w) C(w,u) (u) dvdw

and it is easy to show that this is symmetric.

gi(x) =/(x)C(x,u)j(u)J(u)du O

where we put

Yj(u)

0<u<1

Pj(u-1)

I < u <

21 according to section5.2. With

jp(x)

0<x<I

g(x) =

Iql(xl)

l<x<21

we find pj(x) = (x)fiCi(x,u)(u)Yj(u) +C2(x,u)J(u)i(u)}du

21

qj(x) = J(x)/C3(x,u)(u)Yj(u)+C4(x,u)J(u)(u)}du

o

and combined with (5.2.1) and (5.2.2) this gives

p(x) = )(x){Yl(x)_Di_xD2}

qi(x) =

Now we will show that for the kernel K(x,u) the eigenvalues are positive. This is accomplished in three

steps:

1. C(x,u) is positive definite, i.e.

21 2/

//C(x,u)(x)q(u)dxdu > 0 for all (x) of integrable square. For

[/Ì1x_1

C3(x1,u)(u) Y(u) +J(x C4(x-1,u)J(u)j(u) }du 1<x

<21

(92.12) j(x)

0<x<1

with q(x) = we find Iq2(X-1)

1<

X < 21 18 21 I{e(x)Ci(x,u)o(u)Yí(u)+e(x)C2(x,u)J(u)i(u)}du 0< x <1 21 21

ì[c(x,u) (x) (u) dxdu =f/{Gi(x,u)i(x)i(u) +C2(x,u)i(x)2(u) +

+C3(x,u)l(u)2(x)+C4(x,u)q2(x)z(u)} dudx = (using (5.2.5.))

!!!

1

B

od5 (o)

+2(x)2(u) }+

_A(o)

i(x)i(u) {H(xu)H(uv) +H(ux)H(xv) }dudxdv

Now H(x-u)H(u-v)+H(u-x)H(x-v) = H(x-v)H(u-v) so we obtain:

1 1

/ / _B {(x_0)1(x)+2(x)}{(11_0)1(11)+2(u)}+A 1(X)q111) .H(xv)H(uv)dudxdv =

(o) _(V)

!

i

!

B(v) (xv)i(x) + p2(X) }H(x_v)dx]2+A{/1(x)H(x_v)dx}2)dv

and this is non-negative, q.e.d.

2. AiJ(x)gi(x)dx O (a)

For Ai/(x)gi(x)dx = (according to (9.2.12)) = /(x)Yi2(x)+J(x)i2(x)}dx> O

3. fri(x)g(x)dx> O (b)

1i

For /J(x)gí(x)dx

_{= /}

/j(x)1(x)C(x,u)j(u)fi(u)dudx > O

ô

because C(x,u) is positive definite.

From(a)and(b)wefind:Aj> O

(i= 1,2,3,...).

Finally the relation (9.2.11) can be rewritten

Aif(x)gj(x)

=

A[{ Yi(x)pj(x) +i(x)q5(x)}dx

=

Or f(x)fi(x)fj(x)dx

=

(9 2.14)

where we used the boundary conditions (5.2.3) to eliminate the terms with D1 and D2 from (9.2.13).

The normalization of the functionsf1(x) is here slightly different from (9.2.11). We introduced a factor A. 9.3 Expansion of a function into a series of eigenfunctions

In this section we will prove that it is possible to expand a function that satisfies the boundary conditions (5.2.3) into a series of eigenfunctionsfk(x).

First we prove that the integral equation

(x) = 2/(x)S(v,x)(v)dv

has only one eigenvalue A = - I.

For by inserting the expression for S(v,x) we find

ppi(x) = Ae(x){Ci+xC2} O < x < i

?p2(x-1) = AJ(x-1)C2 I x 21 (-o-ß2+vß1) '1(v) +ß12(v) }dv

C2

ßo2(v)}dv

We substitute (a) into this and find:

C1 = {(_yo_ßz)ßoGj+(_yo_ßz)ßiC2+ßi2Ci+ßiß2Cz+ßiyoG2} = AC1 G2 = {ßißoCi+ßi2G2_ßofliCi_ßoß2G2_ßoyoG2} = AC2

p(x) =

where C1

so A must be - I. The corresponding eigenfunction is

í(x){Ci+xCz}

O

x < I

I'(x) =

tJ(xl)Cz

i Z x 21 _J

According to FREDHOLM'S theorem ([11] p. 269) the integral equation 21

=

(x)

has a solution if q(x) is orthogonal to the eigenfunctions of S(v,x)(x), i.e. if

[(x)(x)dx

= O

}

(a)

20

21

=/K(x,u) (u)du

-This condition is satisfied if q(x) satisfies the boundary conditions (5.2.3).

Next, with the function (x) that we found, we consider the integral equation

21

[C(x,u)(u)(u)du

= (x)

O

This has a solution for every (x) if

/'(x){f(x,u)(u)h(u)du}2dx> /1> O

O O

for every h(x) that satisfies /(u)h2(u)du = 1 ([11] P. 47). (In accordance with (9.2.14) we introduced

a weight-function (x) in the integrals). Now according to SCHWARZ' inequality ([11] p. 18),

21 21 /;ì(x)h(x) /(x,u)(u)h(u)dudx /(x){[C(x,u)1(u)h(u)du}2dx> 0 0 ô ° _{f21(x)h2(x)dx} o 2121

=/[C(x,u)

r (x) h(x) (u) h(u) dudx (c)

In section 9.2 we found that this expression is non-negative and that it is zero if and only if i1(x)h(x) O. But this is not the case here, so the expression (c) is positive. In other words: the expression (c) represents a set of positive numbers, which consequently has a largest lower bound j,, such that

2121

[fC(x,u)(x)h(x)ìj(u)h(u)dudx

> t>

O (q.e.d.) Ô O

Finally by combining (9.3.1) and (9.3.2) we find

92(x) =/k(x,u)(u)du

if 92(x) satisfies the boundary conditions.

Now we calculate the "Fourier-coefficients"

92k =/92(x)gk(x)dx and remember that K(x,u) f,(x)gi(u)

Then 92(x)

92kfk(x) =/(x,u) (u)du-

fk(x)fcp(v)gk(v)dv =

k=1 ¡ç=1 O

21

= /K(x,u)(u)du-

((x)gk(v)/(v,u)z(u)dudv =

o _k1

and to prove this was the purpose of this section. 9.4 Choice of the computation method

Two methods for computing the eigenvalues of (5.2.9) suggest themselves. They are theoretically

equivalent so that the proof of convergence, given in section 5.3 applies to both. But there are some

practical differences.

First method. In this methd we use the integral equation as given by (5.2.6)(5.2.7) or by (5.2.9).

We use the iteration-process as described in section 5.3.

Second method. We use the integral equation as given in (5.1.1) and (5.1.2) together with the boundary conditions (5.1.4). Here too we use the iteration-process, but now we have to compute the con-stants Y(0) and (0) after each step.

21 gk(u) Ifk(x)x(u) du

_{= f}

1K(x,u Ak _O fk(x)gk(u) (u)du O (93.2)

In section 5.3 one method is described for the elimination of previously computed eigenfunctions. This

is necessary because the process converges to the smallest eigenvalue. Another possibility is the so-called

deflation method [12]. Here the elimination takes place by means of the kernel. After the computation

of n eigenvalues, the kernel K(x,u) is replaced by K(x,u) j(x)g(u). This kernel has the eigenvalues

An+1, 2n+2, . ., the smallest of which can be determined with the process of section 5.3.

Because the kernels g(x,u) do not depend on mass and mass-moment of inertia, it is advantageous to

keep them unmodified in the memory of the computer. Then the corrections of the kernel in the deflation

method have to be stored separately, which requires much space, or have to be added by including

them in the formulae (5.1.1) and (5.1.2) which has the same effect as the elimination method as to the work that has to be executed.

For a small computer we come now to the following comparison.

We conclude that for the available small computer, the second method is preferable. Because of item (4)

above we have to use the elimination method. The computation method that is obtained in this way, is described in section 9.6 in the form of an Algol-program.

9.5 Algol-program for the Holzer-Myklestad method For a flow-chart we refer to fig. 9.

begin integer m, n; real length, freq, epsilon, pi, norm;

m: = read; n: = read; length: = read; freq: = read;

epsilon: = read; pi: = 4arctan (1);

comment m is the number of elements into which the beam is divided, n is the number of loading

conditions, freq is the step in the process of determining the approximate location of a zero, epsilon is the required relative accuracy;

begin integer k,l; array a,b[1:m],c,d[l:n,l:n];

comment a, h, c, anddare 1/A(x), l/B(x), J(x) and e(x) respectively;

for I: = 1 step 1 until m do

begin a[kJ: = read; a[k]: = l/a[k];

= read; b[k]: = length 2/b[kI;

for 1: = 1 step 1 until n do

begin c[l, k]: = read; d[l, k]: = read; d{1 k]: length ' 2 x d[l, k] end

end input;

for 1::-- 1 step 1 until n do

begin integer i; real h; array mu[l :3], delta[l :3];

First method Second method

contra:

complicated kernels

2. contra:

kernels depend on mass and mass-moment oF

inertia

pro:

kernels do not depend on mass and mass-moment of inertia

3. pro:

he integrals in (5.1.1) and (5.1.2) can be

con-structed from a number of functions of a single variable, which saves memory-space (cf. section

9.6)

4. pro:

deflation method useful; so all eigenvalues

re-quire an equal amount of work per

iteration-step

contra:

as stated above, the deflation method is not

use-ful,

so higher modes require more work per

iteration-step, which is the disadvantage of the

22 BY ,u3+1O N TIMES t N T REPLACE1L/, NO READ DATA MAKE z1O MAKE _!'2 /11+FREQ.-STEP T FIND REGION OF THE ZERO I-PRINT RESULT NO REPLACE ,uBY ,LI3AND 2BY,

Fig. 9. Flow-chart for the Holzer-Myklestad method

REPLACE/A, BY /A3AND REPLACE p.BY1L11 AND , BY COMPUTE 2Z

()

comment mu is o2, delta is the determinant of section 4;

real procedure determ (mu); value mu; real mu;

comment determ computes the determinant as a function of o2;

begin array i[l:2, 1:4], p[1 :4];

procedure transf (mu); value mu; real mu;

comment transf determines Vm for given vo;

begin integer r,s;

for r: = I step 1 until în do

begin v{2,1]: = v[1,l]mu>hd[l,r] ><v[1,4J;

v[2,2]: =v[l,2]+h><(v[l,1]+muxc[l,r]xv[1,3]);

v[2,3]: =v[l,3]h<b[rJ<v[1,2];

v[2,4]: =v[1,4]+hx(v[l,3]+a[r]xv[l,1]);

norm: =sqrt(v[2,l]t2+v[2,2fl'2+v[2,3] t2+v[2,4] t2);

for s: = I step I until 4 do v[l,s] : = v[2,s]/norm; end end transf;

v[I,l]: = v[l,2]: = v[1,3j: = 0; v[l,4]: = 1;

transf (mu); p[lJ : = v[l,1] ; p[2] : = v{l,2]

v[l,l]: = v[l,2]: = v[1,4]: = 0; v[1,3]: = 1;

transf (mu); p[3] : = v[l,1] ; p[4] : = v{1,2] determ: = p[3] xp[2]p[1] xp[4] end determ;

Ii: = 1/rn; mu{l]: = 10; delta[1]: = determ (mu[1]);

NLCR; print (delta[1]) ; i: = 0;

scan: mu[2]: = mu[1J-1-freq; delta[2]: = determ (mu[2]);

NLCR;print (delta[2]);

if cign(delta[2]) = sign(delta[I]) then

begin mu[1]: = mu[2]; delta[l] : = delta[2J; goto scan end;

bc: comment this part of the program localizes the zero that has been passed somewhat

more accurate. This is necessary because the convergence of the following process is

rather slow;

delta[2]: = determ (mu[l]+10);

if s(gn (delta[2]) = sign (delta[lJ) then

begin mu[l]: = mu[1]+ 10; delta[l]: = delta[2] ; goto loe end

else begin mu{2]: = mu[1]-Fl0;goto inter end;

inter: comment this part of the program determines the zero that has been localized, with a

relative accuracy ;

mu{3] : = (mu[2]xdelta[1] mu[1]x deita[2])/(delta[1]

delta[2])

NLCR;print(mu{3]);

if abs((mu[3] mu[2])/mu[2]) epsilon then goto exit; if abs((mu[3] mu[1])/mu[l]) < epsilon then goto exit;

delta[3] : = determ (mu[3]) ; print(delta[3])

if sign(delta[3]) = sign(delta[2]) then goto equal; mu[l]: = mu{3]; delta[l]: = delta[3]; goto inter; equal: mu[2]: = mu{3]; delta[2]: = delta{3]; goto inter;

exit: NLCR; PRINTTEXT( frequency );

print(sqrt(mu{3])/(2xpi)); NLCR;

i: = i+1; ¡fi <5 then

begin mu[1]: = mu{3] + 10; delta[l]: = determ(mu[1]); goto scan end

end end

24

9.6 Algol-p rogram for the method of integral equations

All integrations are executed by means of the trapezoidal rule. To simplify the process of computation,

the following functions are introduced:

az(x) o A(a) du X ui2 b(x)

du (i = 2, 3,4)

ai(x) rfirst I/A (x)

I later on a2(x) +b4(x) -xb3(x) j first 1/B(x)

later on xb2(x) b3(x)

= /(u)Y(u)du

o

V2(X) =/u(u) Y(a) +J(u)(u) }du

v3(x) = fai(u)e(a)Y(u)du

o

X

n4(x) =/{bi(u)(u)Y(u) +b2(u)J(a)(u)}du

V5(X) =/b3(u)J(u)P(u)da By means of these functions we can write:

/{gi(x,u)(u) Y(a) +g2(x,a)J(a)(a) }du =

o

= aj(x)vi(x) bi(x)vz(x) +V3(x) +x04(x) v5(x) /{g3(x,u) (u) Y(u) +g4 (x,u)J'(u) P(u) }du =

= b3(x)vi(x)b2(x)v2(x)+v4(x)

We will call

J(x) = ci(x), e(x) = di(x), c2 = fj(u)du

o

and di = fi_2o(u)du (i = 2, 3, 4)

C)

The flow-chart is given in figure lO.

A

K1(1)N

COMPUTE A() AND BJ

COMPUTE C2 02 03 AND 04 =1(1)5 INITIAL APPROXIMATION Y(X) :COS(I+1) TTX/L Øcx- (I+1)SIN(I+1)i J: 1(1)1+1 ELIMINATE PREVIOUS EIGENFU NC TI ONS COMPUTE NEW (X) ANO (X) ADJUST BOUNDARY CONDITIONS ENO

Figure 10. Flow-chart for the method of integral equations

begin integer m,n; m: = read; n: = read;

begin integer i,j, 1, t,p; real length, h,pi, delta, c2, d2, d3, d4;

array a,.y,fi[l :2,1 :m+ 1], b, v[l :4, 1 :m± l],c,d{l :n, i :m± l],f,g[l :m+ 1], z,psi[l :5,1 :m+1], mu[l :6], e{1 :5];

comment f and g are in- and output for the integration subroutine, z and psi are eigenfunctions, mu is w2;

procedure integrate (alfa);

boolean alfa; 'I COMPUTE ANO EIGENF UNCTION PRINT

D

comment integrate determines the integral g(x) off(x) with g(0) O

when alfa = true, and the definite integral g(l) off(x) when alfa =

false;

begin if alfa then

begin integer k; g[l]: = O;

fork: = 2 step i until m+l dog[k]: =g[k-1]+hx(f[k--1]+f[k])/2

end else

begin integer k; g[in+l]: = O;

fork: = 2 step i until m dog[m-f-I]: =g[rn+l]+f[k];

g{m--l]: =hx(g[m-]-lJ+(f[1]+f[m-]-i])/2)

end

end integrate;

length: = read; h: = length/rn; pi: = 4arctan (1);

for i: = i step I until m

l do

begin a[l,iI : = read; a[1,i] : = l/a[l,i]

b[l,i]: = read; b[I,i]: = l/b[l,i];

forj: = i step I until n do begin c[j,i]: = read; d[j,i]: = read end

end input;

for j: = i step 1 until 5 do

begin for i: = i step i until ln+ I do

begin ifj= 1 thenf[i]: = a[l,iI else

ifj = 2 thenf[i]: = b[l,i] else

ifj < 5 thenf[i]: = ((iI) xh) t (j-2) xb[I,i];

ifj = 2 then a[2, i]: = g[i] else

ifj> 3 then b[j-1,i]: =g[i]

end;

if j < 5 then integrate (true) end;

for i: = i step i until in + I do

begin a[1,i]: = a[2,i]+b{4,i](il) xhxb[3,i];

b[1,i]: = (i-1) xhxb[2,i]b[3,iJ

end;

for i: = i step i until n do

begin for j: = i step 1 until 5 do

begin for i:

step i until in+ i do

begin if j = I thenf[i]: = c[l,i] else

ifj = 2 thenf[i]: = d[l,i] else

ifj < 5 thenf{i]: = ((ii) xh) t (j-2) xd[l,i]

end;

ifj = 2 then c2: = g[m+ 1] else

ifj = 3 then d2: = g[m+1] else

ifj = 4 then d3: = g[m+ 1] else d4: = g[rn+ 1];

if j < 5 then integrate (false)

end;

delta: = (c2+d4) xd2d3 t 2;

for i: = i step i until 5 do

begin for j: = 1 step i until m+i do

begin5[1,jJ: = cos((i+l) xpix(ji)/m);

fi[l,jJ:

(i+i) xpixsin((i+l) xpix (ji)/m)/length

end;

for j: = I step I until 1+1 do

begin fori: = i step i until i i do

comment now follows the elimination of preceding eigenfunctions;

begin for p: = i step I until rn+ i do

26

end

integrate (false);

eit]: =g{m+1] xmu[t];

for p: = i step i until m + i do

beginy[ l,p]: =y[ l,p] e[t] X z{t,p] fi[1,p] : =fi[1,p] e[t] xpsi[t,p] end

end elimination;

for t: = i step i until 5 do

begin forp: = i step i until m+l do

begin if t> i then u[t-1,p]:

= g{p];

= if i = I then d[1,pJ xy{i,p]else

if t = 2 then (pl) x h xf[p] +c[l,p] xfi[i,pJ else

if t = 3 then a{l,p] x d[l,p] xy[1,p] else

if t = 4 then b[1,p] x d[l,p] xy[1,p] --b[2,p] x c[l,p] >jl[l,p]

else b[3,p] xc[1,p] xfi{1,p]

end;

integrate (true)

end:

for p:

I step i until m+ I do

beginy[2,p]:

_{= a[l,p]xv[l,p]b[i,p]xv[2,p]+v[3,p]_g[p]+(p_flxhxv[4,p];}

Jl[2,p]: = b[3,p] xv{l,pJb[2,p] xv[2,p]+v[4,p]; = d[1,p] xy[2,p] end; integrate (false);

e[l]: =g{m+1];

for p: = i step i until m + i do

= c[1,p] <fi[2,p]+(p-1) xhxf{p];

integrate (false);

e[2: =g[m+l];

e[3]: = ((c2+d4) xe[l]d3xe[2J)/detta;

= (d3xe[1I+d2xe[2])/Jelta;

for p: = i step i until m+i do

comment here we find a new approximation of the desired eigenfunction;

beginy[2,p]: =y[2,p] e[3] e[4] x (p-1) x h;

fi[2,p]: =fi[2,pJ e[4];

ifj < i+1 then begin y[i,p]: ==,[2,pi;fi[l,pJ: =fi[2,p] end

end end iteration:

for p: = i step i until m +1 do

= d[1,p] xy[l,p] xy[2,p] +c[1,p] xfi[1,p] xfi[2,p];

integrate(false);

= g[m+1];

for p: = i step i until m+1 do

f[PI: = d[í,p] >Çy[2,p] t 2+c[1,p] xfi[2,p] t 2;

integrate(false);

e[2]: =g[m+1];

mu[i]: = e[l]/e[2];

NLCR; PRINTTEXT ( frequency ) ;print(sqrt(mu[i])/(2 xpi)); e[1] : = 1/sqrt(eíl])

for p: = i step I until m +1 do

begin z[i,p]: =y[2,p] xe[1];

psi{i,p] : =fi[2,p] x e[l];

end

end integral equations

9.7 Data for the computations

TABLE IV. Data of the homogeneous TABLE VI. Data of the MS. "Naess Falcon" beam

/- 100 A(x) 00 9(x) 100

B(x) = 108 3(x) 0

TABLE V. Data of the SS. "Gopher Mariner" {l3

1= 170.7m E=2.l2xlO7ton/m2 G=83x10ton/m (lton-l000kg

The added mass of the entrained water has been computed by the usual Prohaska-method.

3(x)O

I = 160 m shearing stiffness x 10-s bending stiffness x 10-v

heavy condition light condition

mass-moment of inertia mass per unit length mass-moment of inertia mass per unit length ton tonrn tonsec2 tonsec2/m2 tonsec2 tonsec2/m

1 1.0 1.0 2 0.79 2 0.79 2 18.85 25.44 50 7.55 50 6.08 3 25.56 75.26 188 11.75 188 9.56 4 34.28 147.34 252 19.44 252 17.24 5 38.76 152.64 276 20.29 276 20.03 6 50.38 125.08 256 35.71 150 21.10 7 47.56 110.25 277 45.41 90 22.42 8 47.73 114.48 207 34.37 201 32.59 9 47.89 115.54 287 41.35 204 34.23 10 48.06 116.60 263 40.05 207 35.04 11 48.31 116.60 224 37.15 215 36.07 12 48.97 122.96 438 42.32 392 36.20 13 55.20 159.00 517 50.90 440 38.53 14 47.73 116.60 328 48.06 201 34.91 15 47.97 116.60 319 49.36 181 34.20 16 48.22 114.48 200 37.54 67 22.74 17 45.48 102.82 192 33.86 162 26.09 18 44.99 93.20 175 27.12 176 25.92 19 33.86 78.44 130 15.75 130 14.43 20 28.88 57.24 77 8.14 77 7.42 21 1.0 1.0 34 1.09 34 1.09 I0°A ton 10 8B tonm tonsec2/m o 2.32 7.1) 3.1 1 3.10 4.8 5.8 2 3.71 6.9 9.3 3 3.94 8.6 16.3 - 4.65 10.1 23.6 5 4.62 11.4 30.4 6 4.38 12.4 28.8 7 3.97 13.2 31.6 8 3.68 13.0 35.6 9 3.49 14.1 37.9 10 3.44 14.2 39.7 11 3.48 14.2 33.5 12 3.61 13.9 31.6 13 3.77 13.1 28.0 14 3.95 12.2 21.2 15 4.20 10.7 17.2 16 4.44 9.1 10.7 17 4.82 7.6 11.1 18 5.25 6.5 9.9 19 5.36 4.7 6.2 20 0.1 7.6 2.3

No. 1 S No. 3S No. 4S No. 5S No. 6S No. 7M No. 8M No. 9S No. 10 S No. 11 M No. 12 M No. 13 M No. 14 M No. 15 S No. 16 M No. 17 M No. 18 M No. 19 M No. 20 5 No. 21 S No. 22 S No. 23 S No. 24 M No. 25 S No. 26 M No. 27 S No.28 M No. 29 M No. 30 5 No. 31 M No. 32 s No. 33 M No. 34 S No. 35 5 No. 36 S

PUBLICATIONS OF THE NETHERLANDS' RESEARCH CENTRE T.N.O. FOR SHIPBUILDING AND NAVIGATION

Reports

The determination ofthe natural frequencies ofship vibrations (Dutch).

By prof. ir H. E. Jaeger. May 1950.

Practical possibilities of constructional applications of aluminium alloys to ship construction.

By prof. ir H. E. Jaeger. March 1951.

Corrugation of bottom shell plating in ships with all-welded or partially welded bottoms (Dutch).

By prof. ir H. E. Jaeger ana' ir H. A. Verbeek. November 1951.

Standard-recommendations for measured mile and endurance trials of sea-going ships (Dutch).

By prof. ir J. W. Bonebakker, dr ir W. J. Muller and ir E. j. Diehi. February 1952.

Some tests on stayed and unstayed masts and a comparison of experimental results and calculated stresses (Dutch.

By ir A. Verduin and ir B. Burghgraef. June 1952.

Cylinder wear in marine diesel engines (Dutch).

By ir H. Visser. December 1952.

Analysis and testing of lubricating oils (Dutch).

Bv ir R. N. M. A. Malolaux and irJ. G. SmiLJuly 1953.

Stability experiments on models of Dutch and French standardized lifeboats.

By prof. ir H. E. Jaeger, prof. ir J. W. Bonebakker and J. Pereboom, in collaboration with A. Audigé. October 1952.

On collecting ship service performance data and their analysis.

By prof ir J. W. Bonebakker. January 1953.

The use of three-phase current for auxiliary purposes (Dutch). irj. C. G. van Wijk. May 1953.

Noise and noise abatement in marine engine rooms (Dutch).

By " Technisch-Physische Dienst T.N.O.- T.H." April 1953.

Investigation of cylinder wear in diesel engines by means of laboratory machines (Dutch).

By ir H. Visser. December 1954.

The purification ofheavy fuel oil for diesel engines (Dutch).

By A. Bremer. August 1953.

Investigation of the stress distribution in corrugated bulkheads with vertical troughs.

By prof. ir H. E. Jaeger, ir B. Burg/zgraefand I. van der Ham. September 1954.

Analysis and testing of lubricating oils II (Dutch).

By ir R. N. M. A. Malotaux and drs J. B. Zabel. March 1956.

The application o.uew physical methods in the examination oflubricating oils.

By ir R. N. M. A. Malotaux and dr F. van Zeggeren. March 1957.

Considerations on the application of three phase current on board ships for auxiliary purposes especially with

regard to fault protection, with a survey of winch drives recently applied on board of these ships and their

in-fluence on the generating capacity (Dutch).

By ir J. C. G. van Wijk. February 1957.

Crankcase explosions (Dutch). By ir J. H. Minkhorsl. April 1957.

An analysis of the application of aluminium alloys in ships' structures.

Suggestions about the riveting between steel and aluminium alloy ships' structures. By prof. ir H. E. Jaeger. January 1955.

On stress calculations in helicoidal shells and propeller blades.

By dr ir J. W. Cohen. July 1955.

Some notes on the calculation of pitching and heaving in longitudinal waves. By ir J. Gerriisma. December 1955.

Second series of stability experiments on models of lifeboats. By ir B. Burghgraef. September 1956.

Outside corrosion of and slagformation on tubes in oil-fired boilers (Dutch). By dr W. J. Taat. April 1957.

Experimental determination of damping, added mass and added mass moment of inertia of a shipmodel.

By ir J. Gerritsrna. October 1957.

Noise measurements and noise reduction in ships.

By ir G. J. van Os and B. van Steenbrugge. May 1957.

Initial metacentric height of small seagoing ships and the inaccuracy and unreliability of calculated curves of

righting levers.

By trof. ir J. W. Bonebakker. December 1957.

Influence of piston temperature on piston fouling and piston-ring wear in diesel engines using residual fuels.

By ir H. Visser. June 1959.

The influence of hysteresis on the value of the modulus of rigidity of steel.

By ir A. Hoppe and ir A. M. Hens. December 1959.

An experimental analysis of shipmotions in longitudinal regular waves.

By ir J. Gerritsma. December 1958.

Model tests concerning damping coefficients and the increase in the moments of inertia due to entrained water on ship's propellers.

By N. J. Visser. October 1959.

The effect of a keel on the rolling characteristics of a ship. By ir J. Gerriisma. July 1959.

The application of new physical methods in the examination of lubricating oils. (Continuation of report No. 17 M.)

By ir R. N. M. A. Malotaux and dr F. van Zeggeren. November 1959.

Acoustical principles in ship design. By ir J. H. Janssen. October 1959. Shipmotions in longitudinal waves.

By ir J. Gerritsma. February 1960.

Experimental determination of bending moments for three models of different fullness in regular waves.