Time-domain for predicting ship motions

rEDOMAIN

ANALYSIS FOR PREDICTING SHIP MOTIONS Robert F. BECK and Allan R. MAGEE

Department of Naval Architecture and Manne Engineenng, The University of Michigan, Ann Arbor, Michigan

The use of time-domain analysis for predicting ship motions is investigated. In the method, the hydrodynamic problem is solved directly in the time domain as an initial value problem starting from rest rather than the more conventional frequency-domain approach. For linearized problems the time-domain and frequency-domain

results are Fourier transforms of one

another and are therefore complementary. For fully nonlinear

simulations the time-domain approach is preferred.

In this paper, both linear and body-nonlinear problems will be discussed. The body-nonlinear problem requires the body boundary condition to be satisfied on the instantaneous position

of the body while

maintaining the linearized free surface boundary condition. In the linear problem, both the free surface arid body boundary conditions are tineanzed. The body boundary condition is linearized about the mean position of the body. Because the free surface condition is linearized about the calm water level, a time-domain Green function approach is used to solve both problems.

Results of linear time-domain calculations are

presented for the Wigley hull form and compared with

experiments. Body-nonlinear computations are shown for a submerged ellipsoid. In both cases, the influence on the time-domain results of the singularity in the frequency domain at t = 1/4 is discussed.

1. INTRODUCTION

Time-domain analysis has proven to be a useful and enlightening tool to analyze ship motions. In

the method, the hydrodynamic problem is solved

directly in the time domain as an initial-value

problem. Its simplicity and applicability to a wide

variety of problems are major advantages.

The computational algorithms remain substantially the same regardless of the amplitude, direction or time

history of the motion. The major disadvantage is

that the solution must be time-stepped with memory in

the system.

This can lead to numerical

instabilities for a fully nonlinear approach, and

requires the evaluation of convolution integrals in the Green function approach.

The

conventional

approach

to

solving

seakeepirig problems is to develop a

frequency-domain solution.

For linearized problems at

constant or zero forward speed the time-domain and frequency-domain solutions are Fourier transforms of one another and are, therefore, complementary. One method or the other might have advantages for a particular problem.

At the present time, fully nonlinear time-domain solutions for arbitrary three-dimensional bodies are

under development but are not yet practical.

Therefore, an intermediate approximation, the so-called Neumann-Kelvin approach, has been used

by many researchers.

In the Neumann-Kelvin

'(

T.

_{LEborabrium 'icor}

Scheepshydompj

kctilet

Mekelweg 2,2828 CD Deft

1&

approach the fluid is considered incompressible

and inviscid so that the Laplace equation governs the flow. The body boundary condition is applied on the exact body surface, but a linearized free surface boundary condition is used. These assumptions

allow the development of numerical solution

techniques using a Green function approach. The

so called panel methods have been used on a

variety of problems.

At zero forward speed the

frequency-domain panel method enjoys widespread popularity in the offshore industry. The confidence level in the results is quite high, (see for example Korsmeyer et al., 1988). At steady forward speed, the freque ncy-domai n panel method encounters

difficulties

because the Green function

is complicated and difficult to compute. Nevertheless, results have been obtained by several researchers

including Chang (1977), Inglis and Price (1981),

and Guevel and Bougis (1982).

The use of time-domain methods is not new. The

solution for the fundamental hr singularity

is credited to Finklestein (1957). Discussions of direct

time-domain solutions are presented by various

authors such as: Stoker (1957), Cummins (1962),

Ogilvie

(1964), and Wehausen (1967). As

computational power has increased, it has become practical to study actual solutions and investigate

the computational advantages of time-domain

methods.

Adachi and Omatsu (1979), Yeung

(1982), Newman (1985), Beck and Liapis (1987),

Korsmeyer (1988), Korsmeyer etaL (1988), King et aI. (1988), and Ferrant (1989) are among those who have successfully obtained results.

For linear problems at zero forward speed, the time-domain computations are not as fast as the

conventional frequency-domain approach because

many time steps are needed (rather than a few

frequencies) to obtain an adequate representation

of the results.

However, at forward speed the frequency-domain Green function becomes very difficult to compute and the time-domain method appears to be significantly faster. For problems where the body boundary condition is applied on

the instantaneous exact body surface the

time-domain method is the only alternative; frequency-domain solutions are limited to a few simple cases.

In this paper the theoretical development for the time-domain approach given arbitrary body shape and motion will be presented. The reduction of the formulation to the case of linear motions at constant forward speed will also be outlined. Comparisons with strip theory and experiments will be given. The

influence of the singularity around r UOú)e

i

g 4

where Uo = the body velocity, 0e is the encounter frequency and g is the acceleration of gravity, will also be discussed.

2. THEORETICAL FORMULATION 2.1. Radiation Forces

The origin of the axis system is fixed on the free surface of an infinitely deep, incompressible, ideal fluid, initially at rest. The z-axis is positive upwards and the

x-y

plane corresponds to the calm water level.

The direction of motion of the body is

generally ¡n the positive x direction. The governing equation in the fluid for the velocity perturbation potential is the Laplace equation:

V2Ø(x,y,z,r)=O (1)

where

u=V

(2)

The boundary condition on the free surface

is linearized to yield:

onz=0

(3)

The body boundary condition is satisfied on the instantaneous position of the body surface:

on Sh(t) '(4)

dn

VØ-+0

VØ-0

The initial conditions are that: and

0-+0

dØ_90

dr

In the usual manner, an integral equation that must be solved to determine is found by applying Green's theorem to the fluid domain yielding:

$

N(0V2G_GV20)=$dS(Ø_G

dn

_dnJ

(8)

V s

with the volume V bounded by S

, where

S Sh+Sf+S00 and Sh= body surface, S1= free

surface, and S00 = surface at infinity.

The Green function for the time dependent

problem is given by G(P,Q,t, r)=

-

!)8(r

-

'r) + H(t - 'r)Ô(P, Q,r, 'r) 00 Ô(P,Q,t, 'r) =$

dk.Jjsin(iJj(r

-o P=(x(r),y(r),z(t)) Q=

(er,iier,cer)

1/2

=[(x_ )2

_{+(_}

17)2

r

+(z_)]

00

H(t) = unit step function

= O

r<0

=

i t>O

(9) I ..2 1/2

r'={(x_ç)

+(y_T7)2(z+Ç)]

R=[(x )2_{+ (} oc

6(t)= deltafunction where _{f ô(r)f(:)dr}=

f(0)

asR-co

(5)

as z*oo

(6)

ast*oo

(7)

where n

= inward unit normal to the body

surface, out of the fluid.

V = instantaneous velocity of a point on

the body surface including angular velocity effects.

At infinity the fluid velocities must all go to zero such that:

It is easily shown that the Green function has the following properties: V2G

= 4irö(P - Q)S(t -

r) d2G G

G0

_9t

VGO

Ofl z=O

As shown in Magee (1990), substituting the Green function into (8), integrating both sides with respect to r and then reducing the integrals over the free

surface using Stokes' Theorem results in the

following integral equation for Ø(P,t):

Ø(P,)+_L Jf

dSQ

Ø(Q,)_L. !

1' Sh(t)

dnQ(r

-2ir dSQ(!

_{-i:.) -Ø(Q,z)}

=JJ

Sh(t)

r

r

dnQ

_j

5cir f5cisQ{øQr)

f ò(°:r)

-00 _{Sh(t) (10)}

-Ô(P,Q,t,r)

Ø(Q,r) dflQ _.J:_. _{Ç dr} dJQ _{ø(Q

_{f Ô(P,Q,t,t)}

-00

-Ô(P,Q,t,t)

±ø(Q.r)}

VN(Q,r)

where F(t) is the curve defined by the instantaneous intersection of the mean hull position and the z=O plane and VN

is the two-dimensional normal

velocity in the z=O plane of a point on

r.

lt should

be noted that in the linearized problem the line

integrals are zero at zero forward speed and reduce

to tnose given by King et aL (1988) for constant

forward speed.

For the maneuvering case of

unsteady, large amplitude motion in the z=O plane, (10) reduces to the equation given by Liapis (1986), Appendix A. For the body-nonlinear problem the

line integrals have nonzero values even at zero

forward speed unless the body is wall sided for all points on

r

for all time.

In many situations a source formulation is more convenient because it leads to easier computations for the tangential velocities on the body surface. In the usual manner of potential theory, it is possible to derive the following integral equation for the source strength on the body surface:

a(P,t)

JJ dSQ

c(Q,r)--(!

1_

_c9Ø(P,t) 2 4x Sh(t)

dnpr

r')

dnp

....L f

dr 1f

dS

4r

np -00 Sh(r)

__!_ 1dr f

dt

a(Q,v)V(Q,r)

4irgJ

s _OO

VN(Q, r)--Ô(P,Q,r,

r) dnp

where V

is the three-dimensional normal velocity

of a point on r(t), o(P,t)

is the unknown source strength, and the potential on the body surface is given by Ø(P,r) =

IS ''

_a(Q,r)1' 4,r

\r

r')

Sh(t) 4ir J

1dt J,' dSQ a(Q,r)Ô(P,Q,t,r)

00

S(r)

₍₁₂₎ +-.__

_"di

_{dtQ a(Q,r)V(Q,r)}

4rg J -00

r(r)

VN(Q, r)Ô(P,Q,t, r)

The hydrodynamic forces acting on the body due to a prescribed body motion are found by integrating the pressure over the body surface. Neglecting the hydrostatic pressure, the unsteady hydrodynamic pressure is given by Bernoulli's equation:

(13)

The generalized force on the body in the jth

direction is then given by:

F=

5J'dSPnj

(14)

Sh(t)

where n ,

representing the generalized

unit

normal, is defined as

(n1,n2,n3) = n (n.j,n5,n6) = rx n

r =(z,y,')

(i,y,)=

body axis system

and j=1,2,3 corresponds to the directions of the axes, respectively.

t- r cz O

r-LIIT(x,y,z,t) = U0x + 0(x,y,z)+ Ø0(x,y,z,t)+ i1(x,y,z,t)

+ O

where as

r)--oo k=1,2,...,7

(19)

U0x + cI =potential due to steady translation = incident wave potential

7

c1(x,y,z,t) =

k=1

k = 7 is the diffracted wave

k=1,2,...,6

are the potentials due to the body

motions surge, sway, heave, roll, pitch and yaw, respectively.

To meet the appropriate body

boundary condition,

on the mean body

surface S0 , the following boundary conditions are

specified for the various potentials:

= U0n1

dn

dØ7 _dØ0

an dn

-=nkÇk+mkÇk

_dn

k1,2,...,6

where Çk(t) is the displacement in thekth mode of motion, and the overdot represents the derivative with respect to time, n is the generalized unit

(16)

normal given in (14) and _mk, resulting from the steady forward motion, is given by

(mi,m2,m3)=(n.V)W

(m4,m5,m6) = (n. V)(r x W) ₍₁₇₎

W = V(U0x+0)

W is the fluid velocity due to the steady forward motion of the vessel in the ship fixed coordinate frame.

The linearized free surface condition is

written as:

(f_uo -I;)2øk+føk

=0

On z=0 (18) The initial conditions for the unsteady potentials are

Since the disturbances generating the unsteady

potentials originate in the neighborhood of the

origin

V(x,y,z,r)-0

on Sc,, k=1,2,...,7 (20)

lt can be shown that the Green function given by (9) also satisfies these conditions if the appropriate transformation is made to the moving coordinate

system.

In developing the linear boundary value problem implied by equations (15) through (20), there is an implicit assumption that the body geometry is such that the steady disturbance potential P0 is small. This is a consequence of the free surface boundary condition that has been linearized around the free stream velocity U0. In addition, the amplitudes cl motion must be small because the body boundary

condition has been expanded about the mean

position of the body surface.

As with the body-nonlinear problem, an integral equation to determine the unknown linear potentials is found by applying Green's theorem and using the Green function (9). The final result may be found in King etaL (1988) and is

The above formulation is for the body-nonlinear

satisfied on the instantaneous position of the body problem in which the body boundary condition is surface. Because of this, the linear system theory

normally invoked in seakeeping analysis cannot be used. The body-nonlinear approach is primarily useful for nonlinear simulations.

A linear time-domain analysis may be developed for the constant forward speed case and can be compared directly with frequency-domain analysis. Either impulsive or nonimpulsive input can be used. Liapis and Beck (1985) developed a theory for the impulsive radiation problem. King et aI. (1988) developed a nonimpulsive approach to both the radiation

and exciting-force

problems. By

appropriate choice of the nonimpulsive input,

numerical errors in the computation can be reduced. In

linear time-domain analysis

it

is more

convenient to work in a coordinate system fixed to the moving vessel. In this case, the total velocity potential is defined as:

Ø(P,r)+.Lffa5QØ(Q,

d (1 1'\

2jrJj

')!-

-so

=__JIdSQ (!

-

!

_--Ø(Q,t)

2ir

r

r')dnQ

__L

$

d$ItLS

Ø(Q,)--Ô(P,Q,r)

2r

J -00

a

G(P,Q,r

r)-dflQ

__L

_{5d4d {u(Ø(Q)}

_{.ô(P,Q,t_}

22rg

Ô(P,Q,t

-

r).Ø(Q,

r))

_uo(Ø(Qr)._ Ô(P,Q,t

₎

Ô(P,Q,t

-

.r)fØ(Q. .r))}

where F is the curve defined by the intersection of the mean hull position and the z= O plane.

In the linear case (13) and (14) may be

linearized to yield:

Fjk(t)=

_P$$iS[À+W.VØk]fl.

(22)

The gradient of _4k may be eliminated from (20) by employing a theorem derived by Tuck (cf. Ogilvie,

1977)

$fiS[mjØk +n(W.

_VØk)]=

_4dt

økzf(lx n

s0

r

where ¡ is the tangent vector to the waterline curve T. Substituting (23) into (22) results in

4k(t)

=

_$f dSqLfl +$f dS

økmJ

so so

(24)

pd1

øk'2J(lx n). W

(21)

lt should be noted that since W is tangential to the

ship hull it is almost perpendicular to (ixn) and the line integral term is of higher order. If W is parallel

to the waterline this is exactly true and the line

integral is zero. For the results shown in this paper

the contribution

of the

line integral to the

hydrodynamic forces acting on the body wilt be

neglected.

For the linear problem, the hydrodynamic forces

acting on the vessel can be related to the more

traditional frequency-domain added mass and

damping. The radiation force in mode j due to an imposed motion in mode kmay be written in general as (cf. Cummins, 1962):

1jk(t)=_/1Jk:k(t) - bfk4k(r) -cfktÇk(r)

(25)

fdKJk(:-r) k()

where Kik represents the memory effect due to the free surface and the hydrostatic restoring forces

have been neglected.

lt is shown in King (1987) that bfk is zero. Cjk is

a hydrodynamic force that is proportional to the

body displacement and is given by

Cjk =

_pj5ds

tPk00mJ pdt øk'J(1 x n).W

(26)

s0

r

where øk represents the large-time limit of the potential Øk øk has a non-zero large-time limit

because

_t9

1im k_=mk which is nonzero. The dn

integral equation that must be solved to determine øk is given in King et aL (1988).

The factor represents the infinite frequency added mass and is given by

where 'Yk is the solution to the integral equation

d(1

1)

1f

V'k+$fdSQ'Yk\;

-

=fdSQnk(

-For sinusoidal motions, the radiation forces are

usually given in terms of the added mass and

damping coefficients. The equivalence between the time-domain force formulation and the frequency

(23)

jk

=p$$ds

lVknj (27)

domain is found by substituting Çk(r) = e03f into (25)

and equating

it with

the frequency-domain

representation:

[w2PJk

-

Cjk -

iwJdi

KJk()e0]ei

oit

=[ W2AJ (w) - iwB (w)]e4»t

Equating real and imaginary parts yields

00 C_k AJk(w)=JLJk

_!j'dwKjk(r)sincot__4

o 00 BJk(w) =

Jdw Kk(t)coswr

The term c is frequency independent. Since it multiplies the motion amplitude, it could be added to the hydrostatic restoring force coefficients in the usual frequency-domain equations of motion. In the typical strip theory (cf. Salvesen, Tuck, and Faltisen, 1970 ) this term does not appear. When comparing

time-domain predictions with

strip theory

it, therefore, appears reasonable to retain c with the added mass. In experiments, the hydrodynamic

part of the spring constant term is sometimes

subtracted off and at other times it is not. Thus,

when comparing time-domain predictions with

experiments one must be careful. The influence of

this term on the time-domain predictions will be

shown in section 4. 2.2. Exciting Forces

Consistent with the radiation force problem, the exciting forces for the body-nonlinear Neumann-Kelvin problem are computed assuming the incident wave system is linear. The boundary condition on the free surface is still given by (3) but the normal

velocity on the body is modified to include the

induced velocity of the incident wave system. As shown in King et aI. (1988) the induced velocity on the body is given by

V0(P,t)= JdrK(P,r_ 'r)0(r)

(30)

K(P,t) =

!Re{

where

Y=xcosß + ysinß

Ç0(t) = arbitrary incident wave amplitude at origin ß = wave propagation angle (ir = head seas) The function K may be identified as the impulse

response function for the velocity field in the fluid resulting from an impulsive wave elevation at the origin. The waves in (30) are long crested; short crested seas could be simulated by adding waves propagating in different directions. The required

body boundary condition for the body-nonlinear

problem is then given by

= V. n - n VØ0 (31)

The remaining boundary conditions on the velocity potential (5) - (7) are unaltered. Therefore, the computer code that solves the radiation problem can be easily modified to include the exciting force problem. A simple addition of the second term in (31) to the body boundary condition is all that is necessary.

King et aL (1988) give the solution of the exciting force problem in the linear case at both zero and

non-zero forward speed.

3. NUMERICAL METHODS

The principal numerical task in the time-domain

method is to solve the integral equation for the

perturbation potential ((10) or (21)) or the source strength (1 1). The hydrodynamic forces acting on

the body are then found from the integration of

pressure over the body surface at each time step so

that in the future, when full simulations are

attempted, these forces will be available to solve the equations of motion. To solve the integral equation a panel method is used in which the instantaneous

wetted surface of the body

is divided into quadnlateral panels. On each panel the unknown potential or source strength is assumed constant. A

time marching scheme is used to solve the

governing integral equation in the time domain. At

each time step a new value of the unknown is

determined on each panel.

The convolution

integrals over the time history are evaluated using a simple trapezoidal wie. The details of the numerical method may be found in Magee (1990) for the

body-nonlinear problem or King (1987) for the linear

problem.

For the linear problem, the calculations of the m1 terms (17) are done using two methods. The first (and most common) is to use simplified m1's in which they are computed assuming W=-U0i. In this case, the simplified mçs are given by:

i cos ß

j sin ß

$ do.)w

where P= (x,y,z) =kr' /.1 =(z+ Ç)/r' /3=

/(ii(

-Vr')

211'2

r=[(x_)2(y_1)2(z;)

_j

m

=(O,O,O,O,U0n3,U0n2) (32)

The second technique is to calculate the m's using a finite difference scheme. The required terms are

derived from the normal derivatives of the fluid

velocities due to the steady flow on the body surface (cf. Ogilvie, 1977). The value of the fluid velocities on the body and in the fluid a small distance, n, away from the body along the normal vector can easily be calculated from the known solution to the steady Neumann-Kelvin problem (cf. Doctors and Beck, 1987). Once these velocities are known, a simple finite difference can be used to determine the

required normal derivatives, and these can be

extrapolated to the limit as n-4O.

The procedure works quite well for the

semi-infinite fluid problem of the double-body flow.

However, when the wave terms are included in the computations for the gradients of the velocities, the results are highly oscillatory, and do not converge for this finite difference scheme. As expected, it is

the panels near the free surface that cause the

difficulties.

Therefore, in section 4 results are

presented only for simplified and double body mj's.

The major part of the computer time

is consumed evaluating the Green function.

In a

typical run of the body-nonlinear program order NP2

x KT2 Green function evaluations are required,

where NP is the number of panels and KT is the

number of time steps.

To evaluate the Green

function efficiently, a vectonzed technique has been developed by Magee and Beck (1989). Integrals involving the i Ir terms are performed using methods similar to Hess and Smith (1964). To evaluate the wave terms they are first written in a two-parameter,

nondimensional form: Ô(P,Q,t, r) = 00 Cu'ß) =

d'Jsin(ßV)L)eJ0

_(1Vf(1

p2))

The parameter j relates the vertical to horizontal distance between source and field points, and /3 is time-like and relates to the phase of the generated waves. The function is oscillatory for large /3 and is sharply peaked (though not singular) for i near O. A plot of is given in figure 1.

A bicubic interpolation technique has been

developed to compute quickly using a vector processor. To reduce memory requirements and increase the accuracy of the interpolation, a two-step approach is taken. First, (33) is written as

=exp'

':

](Oß)+

G1(z,ß)

where (34)

(O,ß) =

{J;/4(Ç].

J_i,[Ç]

+Jy_J.

3i4(ÇJ]

is given in Wehausen and Laitone (1960) and is the value of the Green function when both the source and field points lie on the free surface (i.e., i=O). The function

(0,ß) may be precomputed and

stored for simple one-dimensional interpolation. The

second step is an interpolation of the G(bL,ß)

function in z and /3 space. Figure 1 also shows a plot of G(p,ß). Note that G(jz,ß) is smooth and a small percentage of (j,ß) , thus allowing a much

coarser grid spacing for equivalent accuracy.

G1(jt,ß) is interpolated using bicubic interpolation on a nonuniform grid spacing. On a CRAY Y-MP the Green function routine runs at approximately 140

Million Floating Point Operations per Second

(MFLOPS) and

takes

approximately

2.7

microseconds per call for Ô and its derivatives

including the return to physical space for a vector (33) _{length of 200.}

In order to understand the difficulties associated

with time-domain analysis,

it is

instructive to

examine the form of the potential due to an isolated source traveling below the free surface. Figures 2a and 2b show the time history of the potential and its time derivative for a field point travelling on the free surface located a horizontal distance R = i behind the source point. The source is travelling at a depth

of z = -2 below the free surface with a Froude

number, F = u0/.J equal to 0.4.

In this case, the critical wave number where t = 1/4 is kR = .3906. The source has unit strength and is created at t = 0; thereafter,

its strength remains constant.

The potential has the form of a typical impulse response

Figure 1.

The time-domain Green function

G(p,B) and the remainder, Gi(p,B) which must

be interpolated.

Figure 2.

The potential and Its time derivative

for a unit source.

Undamped

(a=O.0); Damped (a=O.05).

0.004

-0 .004

100.

TI(G/R)

160.

Figure 3.

Greatly expanded view of the

large-time tail of (a) the potential and; (b) its

time derivative for the source of figure 2.

a=O.O; cz=O.05.

Figure 4.

Fourier transform of the time

derivative of the potentials of figure 2b.

a=O.O;

a=O.05; Frequency domain

Green function by Wu and Eatock Taylor, 1987

i . 5

0.0 0.

function: there is an initial large response followed a decaying oscillation about the large-time limit.

The solid lines

in figures 3a and 3b are

expanded views of the oscillatory tail at large time. The oscillations are at a frequency equivalent to

t = 1/4 and lead to the singular behavior in the

frequency domain. Recall that the Fourier transform of a sine wave of frequency 0)0 in the time domain is a delta function at w in the frequency domain. Thus, the t= 1/4 singularity appears as a decaying sine wave at the appropriate frequency in the large-time tail of the large-time-domain results.

Figure 4 shows the real and imaginary parts of

the Fourier transform of the derivative of the

potential (figure 2b). The Fourier transform of the time derivative is equivalent to the values in the

frequency domain of the Green function for a

translating, pulsating source. The small squares in figure 4 are the values of the Green function taken from Wu and Eatock Taylor (1987). As can be seen,

the agreement between the frequency-domain

computations and the Fourier transform of the

time-domain response is

excellent.

The singular

behavior of the frequency-domain Green function at t = 1/4 is clearly visible.

To eliminate the anomalous behavior in time-domain simulations due to the t = 1/4 singularity, artificial damping can be introduced into the

large-time tail of the large-time-domain Green function.

In particular, the large time asympotic expansion (cf.

King, 1987) that is used in the Green function

subroutines is multiplied by , where a is

an arbitrary constant which determines the strength of the damping. Typical values of a are in the range 0.0 to 0.1 The artificial damping has the effect of forcing the large-time tail of the time-domain Green function to decay exponentially to zero. The dashed

curves in figures 2, 3, and 4 are the results

computed using the damping.

In figure 2 the

dashed curve cannot be seen because on the scale of the figure the damped and undamped curves are identical. The effects of the damping are clearly seen, however, in the expanded scale of figure 3. In this figure, the oscillations in the damped curve are substantially lower than in the undamped case and they decay much faster.

The effects of the damping in the frequency

domain can be seen in figure 4. The damped and undamped curves are effectively identical except around t= 1/4. In this region, the damping smooths out the singularity. As the damping coefficient a is

reduced, the results become more and more

singular. Ir, section 4, the effects of the damping on motion simulations will be shown.

4. RESULTS

An example of computed results using linear time-domain analysis is shown in figures 5 and 6.

These figures show the heave and pitch added

mass and damping as a function of frequency for a modified Wigley mathematical hull form. The Wigley form has a length-to-beam ratio of 10 and a beam-to-draft ratio 1.6. The half beam is given by the equation:

Z

₌

b

where b = the half beam of the model

L = model length

T = model draft

For the calculations, 240 panels (30 lengthwise x 8 girthwise) were used on the half-body. The non-dimensional time step, defined as &'=&

i,

was

VL 0.088 and the total number of time steps was 256. Figure 5 presents results for zero forward speed and figure 6 is for a Froude number of 0.3. Also shown in the figure are strip theory results computed using

a Salvesen, Tuck, Faltinsen (1970) based strip

theory (cf. Beck, 1989). The experimental results

are due to Gerritsma (1988).

No experimental results are available for F=0.0.

The Wigley hull

is fore-and-aft symmetric.

Consequently, for F=0.0 (figure 5) the cross

coupling coefficients are zero. As can be seen in the figure 5, the strip theory and time-domain results

agree very well for high frequencies.

For low frequencies,

however, there

are significant differences.

This is presumably due to

three-dimensional effects that are neglected in strip

theory.

Note the oscillations in the damping

coefficients in the frequency range 6. - 8. which are a result of irregular frequencies. Inglis and Price (1981) used a box barge with the same dimensions as the ship to estimate the irregular frequencies. By this method, the first irregular frequency occurs at 5.7 and the second at 7.9. The singular nature of the coefficients around the irregular frequencies is not properly captured because the record length and time step size are insufficient to give adequate resolution to the Fourier transform. This is both a weakness and strength of the time-domain method. As with the t = 1/4 singularity, in the time domain the irregular frequencies appear as oscillations in the

large-time tail and do not cause any particular

problems until the Fourier transform is taken.

Because of the finite record length and time-step size, the irregular frequencies are effectively filtered out. On the other hand, in the frequency domain the results are singular and calculations in the region of the irregular frequency are meaningless.

(

\.T)

A

\.T)(

.L) )

()8

"y

'...T)).

+0.2(

_\L)

-KL))

(35)

Figure 5.

Added mass and damping for a

Wigley hull, Fn=O.O. Time

domain; - - Strip theory.

Figure 6 presents the results at forward speed. It should be noted that the added mass coefficients include the Cjk/CO term as presented in (29). This ¡s consistent with strip theory and Gerritsma's (1988)

experimental results since he used only the

hydrostatic coefficients in reducing the measured

data to obtain his added mass and damping

coefficients. The solid line was computed using the double body mfs and the long dash curve results from using the simplified version of the mjs (32). As

can be seen, for this smooth slender body the

different m's do not change the results very much. The most prominent feature of figure 6 is the singularity at t = 1/4. The strip theory results are not affected by t = 1/4 and the experiments were not conducted at such low frequencies. As previously discussed, the t = 1/4 singularity will affect all linear three-dimensional theories. lt

will be shown in

figures 9 - 13 that it also affects the body-nonlinear problem and random sea simulations.

In figure 6 the agreement between theory and experiment is mixed.

For pure heave, the strip

theory appears better while for the cross-couplings

and for pure pitch the three-dimensional

time-domain predictions give better agreement. For this ship, using the double body m's is not worth the extra effort. Perhaps the complete m1's might improve the agreement for all modes of motion. lt

should be noted that

in the body-nonlinear

computations the mj terms are automatically

accounted for. The agreement between strip theory,

linear time-domain analysis, and experiments

shown in figures 5 and 6 is typical of the results that have been calculated for many different ship types (see Magee and Beck, 1988).

As discussed previously, the character of the curves shown in figures 5 and 6 around the irregular frequencies and at t = 1/4 is the result of the large-time tail in the large-time domain. This is demonstrated by figure 7 which shows an expanded view of the large-time tail of the pitch moment time-domain record for the Wigley hull calculations. At a Froude number of 0.3 the tail is dominated by the t = 1/4 frequency and at zero forward speed the irregular frequencies appear. Note that the zero speed tail has been multiplied by 1 OOx in order to plot it on the

same curve as the forward speed result.

The apparent beat frequency in the zero speed tail is probably the result of interference between the two

lowest irregular frequencies. We have seen no

evidence of irregular frequencies in calculations at forward speed. This might be numerical because of the dominance of the t = 1/4 frequency or it might be the result of forward speed effects; at this point we do not know.

Examples of body-nonlinear calculations are shown in figures 8-13.

The results are all for a

submerged ellipsoid of length-to-diameter ratio of 5. At the present time, the computer code for

body-0.7

In t

e' 0.0 i i I I I i i i

0.0

4(L/G)

9.0

Figure 6.

Added mass and damping for a Wigley hull, Fn=O.3.

double-body mj's;

- Tkme domain,

approximate mi's;

+ Experiments, Gerritsma, 1988.

0.3 0.0 i i I r i r

'I(L/G)

9.0 Time domain,

Strip theory;

2.5

'N

0.0 i r ¡ i i I I 00

I(L/G)

90

Figure 7.

Large-time tail of the time-domain

record of the pitch moment for the Wigley hull.

Fn=O.3; - -

Fn=O.O.

Fn0.O solution has been expanded 100X.

0.020

-0.004 -j j J I j

0. 70.

TI(G/L)

Figure 8.

Surge and heave force for a

submerged ellipsoid started from rest at two

different depths of submergence.

Time domain, H/L=O.16;

-

Doctors and Beck, 1987 H/L=O.16; Time domain, H/L=O.245;

- Doctors and

Beck, H/L=O.245.

Figure 9. Greatly expanded view of the

approach to steady state for H/L=O.245.

Time domain;

- Doctors and

Beck.

nonlinear calculations in the case of floating bodies is not complete. In order to save computer time, the ellipsoid has been approximated by 36 panels on one-half of the body (12 lengthwise x 3 girthwise). This is the same configuration used by Doctors and Beck (1987) for some of their steady Neumann-Kelvin calculations. Using more panels will alter the absolute value of the forces, but the character of the

curves remains the same.

For numerical

investigations of a submerged body, the 36 panel ellipsoid seems sufficient. The calculations were all made starting from rest, using a smooth start-up in both forward speed and heave motion. The time histories of the forward speed and forced heave motion are given by:

=

uo(i

-

e_0121) (36)

Ç3(t)

= 1e01

(

2\N

Iz

Ajsin(wjr) )i=1

Note that this motion has zero initial velocity and acceleration, but contains a jerk at t = 0.0.

Figure

8 shows the surge force

(wave resistance) and heave force acting on the ellipsoid as it starts from rest and approaches a constant forward speed. The ultimate Froude number is 0.35. The two depth to length ratios are 0.16 and 0.245. The straight lines are from the steady Neumann-Kelvin code used by Doctors and Beck (1987). The same body geometry was used for both codes, but the time-domain results were extrapolated using three values of the time step size. As can be seen, the time-domain results quickly approach the steady state values. The large-time values for both heave and surge agree with the steady Neumann-Kelvin results to within 0.5% even for the smaller depth of submergence. This is a good verification for two entirely different computer codes.

If figure 8 is examined very closely, it can be

seen that there are oscillations in the time-domain results as they approach steady state. Figure 9 shows a greatly exaggerated scale of the approach to steady state for the surge and heave forces at the deeper depth. The surge force is expanded 125x and the heave force 400x. The oscillations are at a

frequency equivalent

to

t = 1/4 and decay

approximately as lit. These results are in complete

agreement with Wehausen (1964) in which he

investigated the effects of the starting transient of a thin ship started from rest. Figure 9 shows that even a smooth start-up leads to some oscillations. The more abrupt the start-up, the bigger the oscillations.

To investigate the effects of t = 1/4 on

simulations of sinusoidal motions, figure 10 shows

the heave force on the ellipsoid as a function of

time.

The body has forward speed (ultimate

(37)

F=0.35) and is also heaving sinusoidally at a

single frequency corresponding to t = 1/4.

The amplitude of motion is AIL = 0.085 and the mean depth of submergence is H0/L = 0.245. In figure 1 Oa, the total heave force with and without artificial damping and the contribution to the force from just the (hr -1k') terms is shown. The curves for the total force with and without artificial damping are almost coincident. The (hr - 1k') terms do not vary from cycle to cycle and are responsible for a large part of the force, particularly on the bottom of the stroke. As

expected, the wave terms make the biggest

contribution near the free surface.

The difference between the total force and the (1fr -1fr') component is shown in figure lOb. As can be seen, the peaks in the force curves are growing slowly for no artificial damping and reach constant amplitude for the case with damping. A line is drawn through the peaks of the curves to emphasize the growth without artificial damping. Since there is a singularity in the frequency domain response at t = 1/4, it is expected that the results without artificial damping will continue to grow and there will be no

steady state solution.

The artificial damping

apparently eliminates this growth. The growth rate for the case without artificial damping is extremely slow. Dagan and Miloh (1980) show that in three-dimensions the singularity in the frequency domain behaves in the vicinity of t = 1/4 as Inico

-

_{wI where}

o. is the t = 1/4 critical frequency. While we have

not yet worked out the asymptotic form for the

growth, it is probably on the order of tn(t).

Figure 11 presents the surge force for the 'r = .20 case with the same parameters as used for figure 10, except that the frequency was lowered. No

artificial damping was used in the computation.

Several interesting features of this figure should be noted. First, the peaks in the unsteady force do not grow as opposed to the 'r = 1/4 case. There is an initial starting transient and then the peaks exhibit a small oscillatory behavior as they approach steady state. For other combinations of frequency and forward speed the oscillatory behavior can be much more pronounced. lt is easily calculated that this oscillation is the beat frequency between t = 1/4

and t = .20.

The starting transient must induce some t = 1/4 response and this forms a beat with the t = .20 force. Also shown in the figure is the approach to steady forward speed if there is no forced heave (a replot of the curve in figure 8 as long dashes at small times) and the mean of the total surge force near the end of the record (small

dashes at large times).

The vertical difference between these two lines is the mean added drag due to heave. This is another advantage of the body nonlinear time-domain method: the mean shift

and slowly varying forces are automaticaIy

computed, no special calculations are necessary.

0.04

-0 . 02_0.

T

t--t

-t-

s-TI(G/L)

Figure 10.

(a) Heave force on a heaving

ellipsoid at Fn=O.35, t=l/4, amplitude

A/L=O.085, mean depth H0IL=0.245.

«=0.05;

(hr - l/) force.

(b) Expanded view of total force minus

(hr - h/r') force.

0.010 -0 .01 6_0.

T(G/L)

i

80.

Figure 11. Surge force on a heaving ellipsoid

at Fn0.35,

0.20, A/L=0.085, H0/L=0.245;

Total;

-Forward speed only; Mean.

Figure 13. Greatly expanded view of the

diflerence forces from Figure 12.

With t=h/4;

-

Without r=h/4. o. -0.03 vf1 _II o. O. 02 -0.03_o.

ÌJ

i.

! t I t t I

'I

,

T'I(GiL)

¶ 60. 180.

T'I(G/L)

I

180. -t

r

r

waves in heave.

-between these two.

Figure 12.

(a) Surge force on a heaving

ellipsoid at Fn=0.35, with a sum of five sine

waves in heave.

(b) Surge force due to 4 sine

The effects of t =

1/4 on more general

simulations are shown in figures 12 and 13. The random sea simulations were made using a sum of

five sine waves in heave as in (37) with forward

speed. Figure 12a has one component placed at

t = 1/4 and in 12b this component has been

eliminated. Each of the simulations were run with and without the artificial damping. The dashed line near the axis is the difference between the two runs.

Without the t = 1/4 component (figure 12b) the

difference is so small it cannot be seen in the figure.

For the t = 1/4 case, the difference continues to

grow.

Initially, the difference is exactly zero

because the artificial damping affects only the large

time asympotic form of the time-domain Green

function. A greatly expanded view of the differences is shown in figure 13. As can be seen, there are some differences between the artificially damped

and undamped calculations even for the case

without the t = 1/4 component; however, they

remain very small.

With the t = 1/4 component

present, the differences are much larger and more importantly, they continue to grow. lt appears that any long-time simulation which uses the linear free

surface boundary condition and includes t = 1/4

forcing is doomed to failure because this component will eventually dominate the solution. Our proposed

artificial damping fix-up is simple but effective.

Undoubtedly, improved methods can be developed.

5. CONCLUSIONS

The primary conclusion from the work to date on time-domain analysis is that it is a viable alternative to the traditional frequency-domain analysis. For

linearized

problems

with

forward

speed, computations in the time domain appear to be much

easier than the equivalent frequency-domain

calculations.

The Wigley model six degrees of

freedom calculations for linear radiation and

exciting force coefficients using 240 panels on a

half-body and 256 time steps required 25 CPU

minutes for F1, = 0.3 on a CRAY Y-MP. For radiation forces only, at zero forward speed the required time

was 9 CPU minutes.

In comparison with

experiments and strip theory, linearized

time-domain analysis gives mixed results. At times the

time-domain calculations compare better with

expenments and at other times the predictions are worse. For the mathematical Wigley hull form, the use of doLble body or simplified m's does not seem to make much of a difference in the added mass and damping predictions.

Time-domain analysis is easily extended to the

body-nonlinear problem

in

which the body

boundary condition is always satisfied on the exact instantaneous wetted surface of body. Research into the body-nonlinear problem is continuing at The University of Michigan. To date computations have been performed only on submerged bodies.

Calculations with a submerged ellipsoid have shown

excellent agreement with a steady

Neumann-Kelvin code. In addition, the calculations have indicated the importance of the frequency domain singularity at t = 1/4. Because of the linear

free surface boundary condition, singularities

appear in the frequency domain results at 'r = 1/4 and the results grow slowly without bound in time-domain simulations. A simple fix-up using artificial damping on the large-time tail of the time-domain

Green function is proposed.

In time-domain simulations, the fix-up prevented the unbounded

growth and otherwise did not alter the solution

significantly. Much more research needs to be done on the large-time asymptotics and rates of decay ¡n order to obtain the exact solutions in both the time and frequency domains.

ACKNOWLEDGEMENTS

This research was funded by The Office of Naval

Research, Contract

No. N00014-88-K-0628. Computations were made in part using a CRAY

Grant, University Research and Development

Program at the San Diego Supercomputer Center.

REFERENCES

(1) Adachi, H. and Ohmatsu, S., "On the Influence

of Irregular Frequencies in

the Integral Equation Solutions of the Time-Dependent Free Surface Problems," Journal of engineering Mathematics, Vol. 16, No. 2, 1979, pp. 97-119.

(2) Beck, R.F., "Motions in Waves," Section 3 of ?rinciples of Naval Architecture,, Chapter 8, edited by E. Lewis, Society of Naval Architects and Marine Engineers, Jersey City, N.J., Vol. III, 1989.

Beck, R.F. and Liapis, S.J., "Transient Motions of Floating Bodies at Zero Forward Speed," Joumal of Ship Research, Vol. 31, No.3, 1987, pp. 164-176.

Chang, M.-S., "Computation of

Three-Dimensional Ship Motions with Forward

Speed,"

Proceedings of

the Second

International Conference on Numerical Ship

Hydrodynamics, University

of California, Berkeley, CA, 1977, pp. 124-1 35.

Cummins, W.E., "The Impulse Response

Function and Ship Motions," Schiffstechnik, Vol. 9, 1962, pp. 101 -1 09.

Dagan, G. and Miloh, T., "Flow Past Oscillating Bodies at Resonant Frequency," Proceedings 13th Symposium on Naval Hydrodynamics, Tokyo, Japan, October, 1980, pp. 355-373.

Doctors, L.J. and Beck, R.F., "Convergence Properties of the Neumann-Kelvin Problem for

a Submerged Body,"

Journal of Shio, Research. Volume 31, Number 4, 1987, pp. 227-234.

(3)

Ferrant, P., "Radiation d'Ondes de Gravite par les Mouvements de Grande Amplitude d'Un Corps Immerge," Ph.D. Thesis, Universite de Nantes, 1989.

Finklestein, A., "The Initial Value Problem for Transient Water Waves," Communications or

Pure and Applied Mathematics.. No. 10, 1957, pp. 511-522.

Gerritsma, J., "Motions, Wave Loads and

Added Resistance in Waves of Two Wigley Hull Forms," Report No. 804, Technical University of Deift, Netherlands, November, 1988.

Guevel and Bougis, "Ship Motions with

Forward Speed in Infinite Depth," International Shipbuilding Progress. No. 29, 1982, pp. 103-117.

Hess, J. L., and Smith, A.M.O., "Calculation of Nonlifting Potential Flow About Arbitrary Three-Dimensional Bodies,"

Journal of Ship

Research, Vol. 8, No. 2, 1964, pp. 22-44. Inglis,

R.B., and Price, W.G., "A

Three-Dimensional

Ship Motion

Theory

--Comparison Between Theoretical Prediction and Experimental Data of the Hydrodynamic

Coefficients

with

Forward

Speed,"

Transactions of the Royal Institution of Naval Architects, Vol. 124, 1981, pp. 141-157.

Inglis,

R.B., and Price, W.G.. "Irregular

Frequencies in Three Dimensional Source

Distribution

Techniques,"

International

shipbuilding Progress, Vol. 28, No. 319, 1981, pp. 57-62.

King, B.K., Beck, R.F., and Magee, A.R.,

"Seakeeping Calculations with Forward Speed Using Time-Domain Analysis," Proceedings of the Eighteenth Symposium on Naval Hydrodynamics, Delft, Netherlands, 1988. King, B.K., "Time-Domain Analysis of Wave Exciting Forces on Ships and Bodies," Report No. 306, Department of Naval Architecture and Marine

Engineering, The

University of Michigan, 1987.

Korsmeyer, F.T., "The First and Second Order

Transient Free-Surface Wave Radiation

Problems," Ph.D. Thesis, Department of Ocean

Engineering, Massachusetts

Institute of Technology, January, 1988.

Korsmeyer, F.T., Lee, C.-H., Newman, J.N., and

Sclavounos, P.D., "The Analysis of Wave

Effects

on Tension Leg

Platforms,"

Proceedings of the International Conference on Offshore Mechanics and Arctic Engineering, '88, Houston, 1988.

Liapis, S.J., "Time-Domain Analysis of Ship Motions," Report No. 302, Department of Naval

Architecture

and Marine Engineering, University of Michigan, April, 1986.

Liapis,

S.J., and Beck, R.F., "Seakeeping

Computations Using Time-Domai n Analysis,"

Proceedings of the Fourth International

Symposium on Numerical Hydrodynamics, National Academy of Sciences, Washington,

D.C., 1985, pp. 34-54.

Magee, A.R., "Nonlinear Time-Domain Ship Motions," Ph.D Thesis, University of Michigan, Department of Naval Architecture and Marine Engineering, 1990.

Magee, A.R., and Beck, R.F.,"Vectorized

Computation of the Time-Domain Green

Function," The Fourth International Workshop on Water Waves and Floating Bodies, Oystese, Norway, 1989, pp. 139-144.

Magee, A.R. and Beck, R.F., "Compendium of Ship Motion Calculations Using Linear Time-Domain Analysis," The University of Michigan, Department of Naval Architecture and Manne Engineering, Report Number 310, 1988.

Newman, J.N., "Transient Axisymmetric Motion

of a Floating Cylinder," Journal of Fluid

Mechanics, Vol. 157, 1985, pp. 17-33.

Ogilvie, T. F., "Singular-Perturbation Problems in Ship Hydrodynamics," Advances in Aoplied Mechanics, Vol. 17, 1977, pp. 91-1 88.

Ogilvie, T.F., "Recent Progress Toward the Understanding and Prediction of Ship Motions," Proceedings of the Fifth Symposium

on Naval Hydrodynamics, Office of Naval

Research, Washington DC, 1964, pp. 3-128. Salvesen, N., Tuck, E.O., and Faltinsen, O., "Ship Motions and Sea Loads," Transactions of the Society of Naval Architects & Marine Engineers, Vol. 78, 1970, pp. 250-287.

Wehausen, J.V., "Effect of the Initial Acceleration Upon the Wave Resistance of Ship Models," JournaI of Ship Research, Vol. 7, No. 3, January, 1964, pp. 38-50.

Wehausen, J.V., "Initial Value Problem for the Motion in an Undulating Sea of a Body with Fixed Equilibrium Position," Journal of Engineering Mathematics, Vol. 1, 1967, pp. 1-19.

Wehausen, J.V. and Laitone, E.V., "Surface Waves," in J-landbuch der Physik, edited by S. Flugge, Springer-Verlag, New York, Vol. 9, 1960, pp. 446-778.

Wu, G.X. and Eatock Taylor, R., "A Green's Function Form for Ship Motions at Forward Speed," Intern -tional Shipbuilding Progress, Vol. 34, Number 398, 1987, pp. 189-196. Yeung, R.W., "The Transient Heaving Motion cl

Floating Cylinders," .Journal of Engineering Mathematics, Vol. 16, 1982, pp. 97-119.