FINITE ELEMENT ANALYSIS OF THE FORCED

OSCILLATIONS OF SHIP HULL FORMS.

David Arthur Smith

X LIBRARY
iRAOUATE SCHOOL
CALIFORNIA 93940
 :

NAVAL POSTGRADUATE SG L
Monterey, California

T P $

FINITE ELEMENT ANALYSIS OF THE FORCED

OSCILLATION OF SHIP HULL FORMS

by

David Arthur Smith, Jr.

J une 197^

Th esis Advisor R .E. Newton

Approved for public release; distribution unlimited.

T160131

* \
SECURITY CLASSIFICATION OF THIS PAGE (When Data Entered)

REPORT DOCUMENTATION PAGE READ INSTRUCTIONS

BEFORE COMPLETING FORM
1. REPORT NUMBER 2. GOVT ACCESSION NO 3. RECIPIENT'S CATALOG NUMBER

4. TITLE (and Subtitle) 5. TYPE OF REPORT 4 PERIOD COVERED

Finite Element Analysis of the Forced Engineer's Thesis;
Oscillations of Ship Hull Forms June 197^
6. PERFORMING ORG. REPORT NUMBER

7. AUTHORf*) B. CONTRACT OR GRANT NUMBERf*)

David Arthur Smith, Jr.

9- PERFORMING ORGANIZATION NAME AND ADDRESS 10. PROGRAM ELEMENT, PROJECT, TASK
AREA 4 WORK UNIT NUMBERS
Naval Postgraduate School
Monterey, California 939^0

11. CONTROLLING OFFICE NAME AND ADDRESS 12. REPORT DATE

Naval Postgraduate School June 197^
Monterey, California 939^0 13. NUMBER OF PAGES
102
14. MONITORING AGENCY NAME 4 ADDRESS) // -
different from Controlling Office) 15. SECURITY CLASS, (of thla report)

Naval Postgraduate School Unclassified

Monterey, California 939^0
15a. DECLASSIFI CATION/ DOWN GRADING
SCHEDULE

16. DISTRIBUTION STATEMENT (of thla Report)

Approved for public release; distribution unlimited

17. DISTRIBUTION STATEMENT (of the abatract entered In Block 20, If different from Report)

IB. SUPPLEMENTARY NOTES

19. KEY WORDS (Continue on reverae aide If neceaaasy end Identify by block number)

Added mass
Finite element
Nonlinear

20. ABSTRACT (Continue on reverae aide It neceaaary and Identify by block number)
dynamic, ideal fluid
A study of the first and second-order,
effects on the forced, harmonic oscillations of a rigid cylinder
in a free surface is presented. The solutions are obtained with
a FORTRAN IV computer program based on the finite element method
using isoparametric elements.
First-order added mass and damping in heave for a semi-circular
and bulb hull form including finite depth effects are presented.
Second-order solutions for the semi-circular hull form in heave
RM
DD ,^ N 73 1473 EDITION OF NOV 65
1 IS OBSOLETE
(Page S/N 0102-014- 6601
1) SECURITY CLASSIFICATION OF THIS PAGE (Whan Data Bnta
|

1
i'liLUHITY CLASSIFICATION OF THIS P AGEC^on Data Ent«fd)

(20. ABSTRACT continued)

also including finite depth effects are obtained. The solutions

are verified by comparing with existing theory and experiment.

DD Form 1473
, Jan 73
1
S/N 0102-014-6601 SECURITY CLASSIFICATION OF THIS P AGE(Whon Dmtm Enffd)
 Finite Element Analysis of the Forced
Oscillations of Ship Hull Forms

by

David Arthur .Smith, Jr.

Lieutenant, United States Navy
B.S., University of Missouri, 1966

Submitted in partial fulfillment of the

requirements for the degrees of

MASTER OF SCIENCE IN MECHANICAL ENGINEERING

and

MECHANICAL ENGINEER

from the

NAVAL POSTGRADUATE SCHOOL

June 197 1*
 ABSTRACT

A study of the first and second-order, dynamic, ideal

fluid effects on the forced, harmonic oscillations of a

rigid cylinder in a free surface is presented. The solutions

are obtained with a FORTRAN IV computer program based on the

finite element method using isoparametric elements.

First-order added mass and damping in heave for a semi-

circular and bulb hull form including finite depth effects

are presented. Second-order solutions for the semi-circular

hull form in heave also including finite depth effects are

obtained. The solutions are verified by comparing with

existing theory and experiment.

i
 TABLE OF CONTENTS

I. INTRODUCTION 16

A. HISTORICAL BACKGROUND 16

1. First-Order Solutions 17

a. General Potential Theory 17

b. Finite Element Solutions 19

2. Second-Order Solutions 20

a. General Potential Theory 20

b. Finite Element Solutions 22

B. OBJECTIVES OF PRESENT RESEARCH 23

II. FORMULATION OF THE MATHEMATICAL MODEL 24

A. THE NONLINEAR BOUNDARY VALUE PROBLEM 24

B. PERTURBATION ANALYSIS 29

1. The Linearized Free Surface Condition — 30

2. The Ship Interface Condition 31

3. Perturbation Formulas for Pressure,

Force, and Wave Amplitude 35
(1
C. THE BOUNDARY VALUE FOR <J)
^ 39
2
D. THE BOUNDARY VALUE PROBLEM FOR <^
^
39

E. THE SECOND-ORDER BOUNDARY VALUE

PROBLEM WITH FINITE DEPTH Hi

1. The Second-Order Boundary Value

Problem for ^ (2) 42

F. FORMULAS FOR PRESSURE, FORCE, AND

WAVE AMPLITUDE 4 3

III. FINITE ELEMENT DISCRETIZATION 45

A. THE DISCRETIZED FIRST-ORDER PROBLEM 47

i \
 B. THE DISCRETIZED SECOND-ORDER PROBLEM 48

1. Calculation of Nodal Values of the

Field Derivatives of <fr 50

IV. THE COMPUTER PROGRAM AND NUMERICAL SOLUTION 52

A. MAIN FEATURES OF THE COMPUTER PROGRAM 52

B. NUMERICAL SOLUTIONS 54

1. First-Order Numerical Solutions 5^

a. Mesh Requirements 5^

b. Mesh Data Preparation 56

c. Convergence of First-Order
Solutions 57

2. Second-Order Numerical Solutions 57

a. Mesh Requirements 57

b. The Second-Order Boundary Condition

on S 58

c. The Ship-Fluid Interface Second-

Order Boundary Condition 60

d. Convergence of Second-Order
Solutions 61

V. NUMERICAL RESULTS 62

A. DEFINITIONS :
62

1. Added Mass and Damping Coefficients 62

2. Dimensionless Force Coefficients 64

3. Wave Amplitude at Infinity 65

B. HULL FORMS 65

C. FIRST-ORDER NUMERICAL SOLUTIONS 67

1. Semi-Circular Hull 67

2. Bulb Hull 73

3. Wave Generator Studies 73

» \
 D. SECOND-ORDER NUMERICAL RESULTS 77

1. Semi-Circular Hull - Infinite Depth 77

2. Semi-Circular Hull - Finite Depth 79

VI. CONCLUSIONS 88

APPENDIX A. Complex Algebra 90

APPENDIX B. Isoparametric Elements 91

APPENDIX C. Computer Program Economization 95

BIBLIOGRAPHY 98

INITIAL DISTRIBUTION LIST 102

I
 LIST OF TABLES

TABLE I. DIMENSIONLESS FORCE COEFFICIENTS 64

TABLE II. CONSTANTS FOR HULL FORMS 66

\
 LIST OF FIGURES

Figure Title

1. Hull at Free Surface 24

2. Reduced Region 25

3. Hull Coordinate System 32

4. Finite Element Mesh 45

5. Hull Forms 66

6. Added Mass and Damping Coefficients in

Heave, Semi-Circular Hull, Infinite Depth 68

7. Variation of Added Mass Coefficient with

Depth, Semi-Circular Hull in Heave 69

8. Variation of Damping Coefficient with

Depth, Semi-Circular Hull in Heave 70

9. Added Mass and Damping Coefficients in Sway,

Semi-Circular Hull, Infinite Depth 72

10. Exciting Force Amplitude in Heave, Bulb

Hull, Infinite Depth 74

11. Wave Amplitude for Heave, Bulb Hull,

Infinite Depth 74

12. Variation of Added Mass and Damping

Coefficients with Depth, Bulb Hull in Heave — 75

13. Reduced Region for Wave Generator Problem 76

14. Wave Amplitude Generated by a Plunger Type

Wave Generator 78

15. Wave Amplitude Generated by a Flapper Type

Wave Generator 78

16. Second-Order Exciting Force Amplitude in

Heave, Semi-Circular Hull, Infinite Depth 80

17. Second-Order Exciting Force Phase Angle

in Heave, Semi-Circular Hull, Infinite Depth - 8l

\
Figure Title

18. Variation of Second-Order Exciting Force

Amplitudes with Depth, Semi-Circular
Hull in Heave 82

19. Variation of Second-Order Exciting Force

Phase Angles with Depth, Semi-Circular
Hull in Heave 83

20. Ratio of Second-Order Exciting Force

Amplitude to First-Order Exciting Force
Amplitude, Semi-Circular Hull in Heave 84

21. Total Exciting Force vs. Time 86

22. Second-Order Wave Form Effects 87

Bl. Element Mapping from £,ri Plane to x,y Plane — 91

CI. Schematic of Coefficient Matrix 95

10 ;
 . .

LIST OF SYMBOLS

Symbol Definition

A One-half of ship's cross-sectional area

a Amplitude of ship's motion (sway or heave)

b Ship's half-beam measured at y =

C C, Added mass and damping coefficients,

m, respectively

c. Celerity of first-order gravity wave and

Stokes second-order wave

c Celerity of second-order gravity wave

-(2)
from velocity potential <{r

d Ship's draft

e, ,e_,e,- Real coefficients in Lewis hull form

^
equation

F(t) Total nonlinear ship's exciting force

P ,P First and second-order components of F,

respectively

F^ >Fg »
ri
Nondimensional first and second-order
force coefficients

F^ ,F^ Resultant first-order hydrodynamic force

y acting on the ship (sway or heave,
respectively)

f ,f First and second-order force amplitudes,

respectively
(2) (2)
*s **fl
Amplitude of second-order static and
dynamic components, respectively

g Acceleration of gravity
G, ,Gp, H,,Hp Frequency dependent coefficients

h Fluid region depth

i,j Unit vectors in cartesian coordinate

11
Symbol Definition

i The imaginary unit J^l

k, ,k ? First and second-order wave numbers,

respectively (k. = a/c.)

I One-half of the total arc-length of a

symmetric ship hull form

m Total number of nodes in fluid region

finite element mesh

ma Added mass
-»-

n Unit outward normal vector from region R

p Hydrodynamic overpressure (referred to

atmospheric)

P >P jP Zero, first, and second-order components

of p

Q(x,0) ,P(x,0) Asymptotic functional forms

q(y) Function describing mode shape of wave

generator

q (y) Function describing mode shape of plunger

p type wave generator

q f (y) Function describing mode shape of flapper

type wave generator

r-,,r
?
,6 ,6 Moduli and arguments of complex numbers

R(x,0) Function representing second-order Stokes

free surface boundary condition

R The domain of the boundary value problem

S Boundary of the Region R

S. Denotes the ith region boundary

s Arc-length along ship's hull measured from

ship's centerline

t Time

u,v x,y components of fluid velocity vector,

respectively

12
 (

Symbol Definition

V Ship's velocity vector

w Reduced region width

w . Reduced region width parameter

min
x,y Two-dimensional inertial cartesian
coordinates

x(s),y(s) Parametric hull coordinates

x,,y h Functions of time such that ex^ey^

represent ship's harmonic motion in sway
and heave, respectively

y * Complex form of y.

z, ,Zp,w, ,Wp Complex numbers

a Reduced region width parameter

Y >Y First and second-order exciting force

phase angles measured with respect to
ship's "top dead center" position

6 Dimensionless frequency parameter

2
(6 = a b/g)

C Complex mapping parameter in Lewis hull

form equation

e Perturbation parameter (e = a/b)

n(x,t) Wave amplitude measured from y =

n (x,t),iy (x,t) First and second-order components of n(x,t)

rP
'00
Dimensionless first-order wave amplitude
at infinity

n ,n First and second-order dimensionless wave

amplitude coefficients
(^(i) = ^(i)
X ,t))

y Maximum amplitude of wave generator

£,n Cartesian coordinates in C,n plane

P Fluid density

13 (
Symbol Definition

a Circular frequency of ship's motion

$(x,y,t) Real-valued nonlinear velocity potential

Real-valued nonlinear velocity potential

i> (x,y,t) of a traveling gravity wave

$ ,$ First and second-order components of

$, respectively

$*^ ,$* First and second-order components of

$*, respectively

c}r ,<J>
First and second-order complex, time-
independent velocity potentials,
respectively
v
4> Modified complex, second-order velocity
potential

<f>
A known test velocity potential
(2)
v
ip Second-order, complex Stokes wave
velocity potential (time-independent)

tyi^ ,?^ m x 1 column vector of nodal values of

9 or 9

H,D,Q m x m real, symmetric coefficient matrices

J Jacobian matrix (2x2)

K^ Total coefficient matrix (m x m) of the
~
discretized first-order system

N 1 x m row vector of cubic isoparametric

shape functions

N 1 x 12 row vector of cubic isoparametric

shape functions at the element level

r ,r ,a,B,y m x 1 column vectors

V The vector gradient operator

2
V Laplacian operator

m '
 ACKNOWLEDGEMENTS

The author wishes to express his deep appreciation to

Dr. Robert E. Newton, Professor of Mechanical Engineering,

for his advice and perceptive guidance throughout the

development of this research. His patience and insight were

key factors in the successful conclusion of this work.

The author is indebted to Professor C. J. Garrison and

Professor T. Sarpkaya for their helpful discussions during

the course of this study.

The W. R. Church Computer Center of the Naval Postgraduate

School provided the excellent computer service and technical

advice for the numerical solutions.

Finally, the author wishes to thank his wife, Wendy, for

her patience, understanding, and encouragement during the time

of this research.

15
 I. INTRODUCTION

Naval architects and marine engineers have long been

greatly interested in the study of the dynamic fluid response

to the forced oscillations of a two-dimensional floating body

The applications of the studies include studying the rigid-

body response of the ship due to an incoming wave train, and

determining the natural frequency and vibration mode shapes

of a ship's structure in connection with propulsion systems

design [3^].

Marine engineers, in modern times, are concerned with,

among other things, the natural frequencies and stability

characteristics of floating platforms near a coastline.

Also, they are concerned with the related problem of finding

the dynamic loads on submerged, moored, or fixed bodies [8],

The modern techniques for solving these types of problems

involve seeking the solution of a steady-state, periodic,

potential flow problem with a free surface in a fluid of

either infinite or finite depth. This problem is made

tractable by linearizing the free surface boundary condition.

A. HISTORICAL BACKGROUND

The historical development of potential flow problems

involving a free surface with wave motion, and the related

problem of the dynamic response of fluids due to the motion

of rigid bodies on and under a free surface is lengthy; but

of interest. Several authors have reviewed this history in

16
some detail [3,17,24,31]. The history is repeated in the

following in the interest of providing a foundation upon

which this author endeavors to add additional information in

an attempt to provide a more complete history.

1 . First Order Solutions

a. General Potential Theory

The foundation of the modern day theory of fluid

motions begins, of course, with the classical works of Newton,

Lagrange, Euler, Stokes, and others.

However, in 1879, Lamb [15] provided one of the

first complete studies of the theory of wave motions; which

he updated five times until 1932. It appears that the next

major contribution was made by Havelock [10] in 1929; wherein,

he studied the solution of the wave forms created by a wave

generator.

The modern beginnings of the study of the inter-

action of a fluid and a harmonically oscillating rigid body

at the free surface start with the work of Lewis [19] in 1929.

Lewis studied the determination of added mass with a boundary

condition requiring zero pressure at the mean position of the

free surface, thus eliminating the effects of frequency of

body oscillation and of damping. Lewis considered the high

frequency case of the more general problem.

The first work which included the effects of

frequency and damping was done by Ursell [38] in 19^9. This

work solved the problem of a right-circular, semi-immersed

cylinder heaving on the free surface of a fluid infinite in

17
extent. Ursell presented added mass and damping coefficients

as a function of frequency. His method of solution involved

placing an "infinite" number of singularities of all orders

at the intersection of the free surface at rest and the

cylinder centerline . The strengths of the sources were

determined from the kinematic boundary condition applied on

the immersed surface of the cylinder. Grim [9] arrived at a

similar solution in 1953.

A period of about ten years passed before any

additional work was reported on the subject. There is a

rather obvious reason for this time lag. Ursell' s solution

involved laborious numerical calculations which had to be

done by hand. The advent of the digital computer made it

feasible to pursue further solutions.

In I960, Tasai [36] extended the work of Ursell

to include roll and sway motions of Lewis forms. At the same

time, Porter [30] compared linearized pressure calculations

with experimental results, and extended the range of solutions

to include other hull forms.

It appears that Yu and Ursell [42] provided the

first theoretical study of the effect of finite depth on

two-dimensional solutions in 1961. This was followed by

additional work by Kim [14] in 1969, who also included finite

depth. In the first work, added mass coefficients and the

properties of the resultant linear wave were reported. The

second work, by Kim, gave added mass and damping coefficients

as a function of frequency and depth. Work by Paulling and

18
 .

Richardson [28] and Paulling and Porter [29] provided some

experimental verification of the theory in 1962. Vugts [40],

at the Netherlands Research Center in 1968, reported additional

extensive, experimental research on seven different ship hull

forms

b. Finite Element Solutions

Zienkiewicz [46], in 1965, provided the first

successful application of the finite element method to Poisson's

equation. Although the finite element method was originally

developed to solve problems in the theory of elastic mediums,

Zienkiewicz ' s work "opened the door" to possible solutions to

boundary value problems of many different types, including the

solution of fluid flow problems; this follows since Laplace's

equation is contained in Poisson's equation.

In the same year (1965) 9

Zienkiewicz, Irons, and

Nath [44] first used the finite element method to find added

nfass . Just as in the case of Lewis and others before, surface

waves were excluded from this solution, as well as in subse-

quent works by Rftren [32], Holand [12], Matsuura and Kawakami

[22], and Matsumoto [21].

Zienkiewicz and Newton [43] presented the first

general theoretical treatment of the fluid-structure inter-

action problem using the finite element theory in 1969. Their

work provided a means of including surface waves and deter-

mination of energy loss due to wave propagation. The latter

was accomplished by providing a radiation boundary condition.

This radiation boundary condition was essential to the

19
successful modeling of infinite regions under steady-state

conditions. In addition, compressibility effects could be

included in the analysis.

Chenault [4] 3 in 1970, appears to be the first to

apply the finite element analysis developed by Zienkiewicz

and Newton to find added mass and damping coefficients as a

function of frequency and hull form. The fluid regions he

studied were chosen to simulate infinite depth. Some of his

results are published in Ref. [23]. Bai [3] also included

surface waves in his 1972 analysis of a wide range of problems.

His work demonstrated the flexibility of the method by studying,

in addition to the problems of Chenault, the additions of sway

and roll motions, finite depth, irregular fluid bottom geometry,

axisymmetric geometry, and fluid stratification. Bai also

provided criteria for properly placing the radiation boundary

and developed a new form of the radiation condition for the

axi-symmetric case.

The author has extended the studies of Chenault

and Bai in the solution of the first-order problems with an

improved version of Chenault 's computer program. The results

are presented later in this work along with results separately

reported by Newton, Chenault, and Smith [23].

2. Second Order Solutions

a. General Potential Theory

The origins of the second-order theory of gravity

waves begin again with the work of Stokes and his classic

solution of nonlinear traveling waves. Others have extended

20
 .

the work of Stokes, but it appears that a modern version

of Stokes' theory was first compiled by Stoker [35] in 1957.

The generally accepted terminology associated with Stoker's

work is "perturbation theory" which of course has applications

in many fields

Fontannet [5], in 1961, appears to be the first

to attempt to use the basic nonlinear perturbation theory as

described by Stoker to solve the problem of finding the

characteristics of waves generated by a plane wave-maker using

Lagrangian coordinates. Ogilvie [25], in 1963, found the

first and second-order forces on a fixed or oscillating

circular cylinder submerged in a fluid with a free surface.

However, his second-order forces were only those due to the

first-order solution, and only the time-independent part was

reported. Tuck [37] found the second-order forces on a

submerged cylinder in a uniform stream in 1965. Salvesen [33]

extended Tuck's work to include wii^g-shaped bodies in his

Ph.D. dissertation in 1966.

Lee [17] and Parissis [27] appear to be the first

to apply second-order perturbation theory to attempt to solve

the fluid-structure interaction potential flow problem with

a free surface in 1967. Lee confined his studies to the

heaving of floating circular and U-shaped sections, while

Parissis studied only the circular shape in heave; both

studies were in fluids of infinite extent. These investiga-

tions, as well as those previously mentioned, required the

development of a fluid-structure second-order boundary

21 \
 .

condition (discussed in more detail later) in addition to

the free-surface equations compiled by Stoker.

The work of Lee and Parissis was followed by that

of Potash [31] in 1970. Potash extended the second-order

studies to include the additional degrees of freedom of roll

and sway and the attendant coupling between modes of motion.

Recently, Garrison [7] completed a second-order solution of

the problem of determining the dynamic loads experienced by

a "fixed" three-dimensional body due to an incoming wave

train. His theory includes submerged and surface-piercing

bodies

All of the second-order solutions previously

mentioned used some form of singularity distribution in the

flow field. The functional forms of the singular functions

used were usually taken from the work of Wehausen and Laitone

[41].

b. Finite Element Solutions

The author has been able to find only one

solution to the second-order, of the class of problems being

discussed, by the finite element method. Allouard and Coudert

[2] presented a paper at the International Symposium on Finite

Element Methods in Flow Problems in January of 1974. They

presented a method of solving problems related to floating

bodies in water waves under unsteady-linear or stationary-

nonlinear (second-order) conditions.

22 I
 .

B. OBJECTIVES OF PRESENT RESEARCH

The experimental studies of Paulling and Porter [29],

Paulling and Richardson [28], and Vugts [40], previously

mentioned, indicated good agreement with the developed linear

theory if the amplitude of motion was small. Potash [31]

observed that small variations in amplitude caused large

discrepancies between Vugt ' s experimental results and those

of linear theory. Potash [31] has shown that second-order

effects become significant at higher frequencies and are more

significant when roll and/or sway modes are present. However,

he did not examine the effect of finite depth on second-order

solutions

The present work has two primary goals. The first is to

determine if the finite element method can be successfully

applied to finding solutions to second-order boundary value

problems of the type previously described. And, if success

is achieved in the first objective, the second major goal is

to study the significance of nonlinear effects in heave in

water of finite depth.

23
 ;

II. FORMULATION OF THE MATHEMATICAL MODEL

The physical problem is to determine the force required

to sustain steady, vertical harmonic motion of a floating

cylindric ship hull form. The ship will be considered

oscillating in an "ideal" fluid of infinite horizontal extent

and arbitrary constant depth.

^ ee Surface

Bottom
T/y-rr/ ;;//;///;;;;;;;;//
jT

FIGURE 1. Hull at Free Surface

A. THE NONLINEAR BOUNDARY VALUE PROBLEM

The domain of the boundary value problem is shown in

figure 1. In the figure, b is the ship's half beam, h is the

mean fluid depth, d is the ship's draft, and 2A is the total

submerged cross-sectional area of the ship.

The fluid is assumed to be incompressible, inviscid,

homogeneous, and without surface tension. It is well known

that under the assumption of irrotationality , a velocity

potential exists. Furthermore, the potential will be a

2k
\ \
 S .

single-valued function in any simply-connected region. The

convention for the potential used herein is

$$ = ut + vj , (1)

where u is the x-component of fluid velocity, v is the

y-component of fluid velocity, and 1 and j are the unit

vectors in the positive right-handed cartesian coordinate

system directions x and y respectively. The symbol V is the

standard "gradient" used in vector calculus

The nature of the problem dictates that at a large

horizontal distance from the ship's hull (cross-hatched area

fig. 1) only an outgoing one-dimensional traveling wave will

be present. Therefore, it is advantageous to consider the

reduced region shown in figure 2. In figure 2, w is the

region semi-width. The meanings of the other symbols shown

in figure 2 will become apparent in the following discussion

W )
Free Surface,S i

3
p-Radiation
Ship Hull,S 2 Boundary,S 4

—
\a Plane of
Symmetry.S, /-Bottom, 5

FIGURE 2. Reduced Region

25
 Conservation of mass requires:

2
V $ = , in R . (2)

p
R is the fluid domain and V is the Laplacian operator.

There are two nonlinear free surface conditions, one dynamic

and one kinematic. The dynamic condition, from the Bernoulli

equation, is

1
• + \ [*£ + S>p + E + yg = . (3)
t

Equation 3 must hold on the moving free surface (p = 0) . In

the equation, g is the acceleration of gravity, p is the

overpressure (referred to atmospheric), p is the fluid

density, and t is real time. The subscripts denote partial

differentiation with respect to the indicated variable. The

kinematic condition comes from a basic assumption from

continuum mechanics simply stated by Stoker: "A particle

once on the free surface remains on it." This assumption

translates into the following expression

D[n(x,t) - y]
Dt
Q
U
*
m
v ;

In equation 3, the possible pure function of time,

which may appear, is assumed to be identically zero without
loss of generality (See Stoker [35], p. 10).

26
 .

Again, equation *! must be satisfied on the moving free

surface. In equation 4, D/Dt denotes the co-moving derivative

(sometimes called the substantial or total derivative) and the

function n(x,t) is defined by the following expression

y = n(x,t) . (5)

Therefore, n is defined to be the wave height above the line

y = at any point x and any time t.

Equation 4 can be rewritten using the definition of $ in

p
equation 1 and the definition of the co-moving derivative

on the moving free surface

On and "near" the bounding surface Sk (radiation boundary)

<Kx,y,t) -* $*(x,y,t) , as x+oo . (7)

In equation 7, $* denotes the velocity potential of a

traveling wave moving away from the ship. It must be mentioned

here that inherent in the statement of equation 7 is the

assumption that the surface Sk is "far enough" away from the

ship's hull. This point will be discussed in more detail

later.

Turning to the interface Sp between the fluid and the

ship's hull, the kinematic condition is

2
The co-moving derivative Is defined as
D( )/Dt = u( ) + v( ) + ( ) .
x t

27
 ;
$-n = v-n . (8)

The vector n is the unit outward normal vector (from the

region R) to the surface at the point of interest . The dot

is the standard scalar product and V represents the vector

velocity of the ship at any instant of time.

Since the motion of the ship is assumed to be restricted

to the vertical direction (heave), the bounding surface S-

becomes a plane of symmetry. The surface S,-, represents the

bottom of the fluid region and is rigid. Therefore, both

surfaces have the same kinematic boundary condition

$$.n = . (9)

An examination of equations 3 and 6 indicates the nonlineari-

ties inherent in the boundary value problem. An additional

problem, not as obvious, is the fact that the free surface

and the hull are moving boundaries. In its present form the

problem has not yielded solutions. Therefore, some approxi-

mate theory is applied.

In the development which follows, a perturbation analysis

will be applied which when carried to the second-order will

produce two linear boundary value problems from the one

presently posed. Further, the boundary conditions of both

linear problems will be referred to fixed boundaries . In

addition, appropriate equations for determining force, and

wave profiles will be obtained.

28 \
B. PERTURBATION ANALYSIS

First it is assumed that $ and n each possess a perturba-

tion series representation as follows

4> = e$ (1) + eV 2) + 0(g 3 ) , (10)

n = en (1) + eV 2) + 0( e
3
) , (11)

where z is a small non-dimensional parameter which is a

measure of the size of the ship's motion and is defined as

e = a/b ("a" is the amplitude of the ship's motion and "b"

is the ship's half-beam). Also in equation 10, $ and $

are defined to be the first and second-order velocity

potentials, respectively, and similarly for r\ and r\ .

Note that as e-*0 there is no disturbance and therefore there

are no zero-order functions. Substitution of equation 10

into equation 2 yields

eV $
2 (1)
+ e
2
vV 2) + 0(e 3 ) = . (12)

Equation 12 implies that

2 (1)
V $ = , and (13)

2 (2)
V <D = , in R . (14)

A similar result holds for equation 9.

29
 1 . The Linearized Free Surface Condition

The following development of the linearized free

surface conditions follows in general that of Stoker [35].

The final linearized equations for the free surface boundary

condition are obtained in the following steps. Equations 5,

10, and 11 are substituted into equation 3, noting that p =

on the free surface. Also, it is observed that $/: , $ ,

and $ may be represented by the following form of the

Taylor series

$ (x,y,t) = $ (x,0,t) + n(x,t)$ (x,0,t)

v v
2
+ \ <D
vyy
(x,o,t) + o(n 3
) . (15)
3

After collecting coefficients of like powers of e in the

resulting equation, the following set of equations results

(1) 1}
gn + <^ = o , (16)

gn
(2) + #
(2)
+n (D # (l) + | [{$
(D } 2 + {f
(l)
}2] = f

on S . (17)
3

Applying a very similar procedure to equation H yields the

following result to the second order

The subscript v denotes partial differentiation with

respect to x, y, or t.

30
 ( 1) X)
n - - and, (18)
t ^ ,

<2) _ # (2) + n
(l),(l) . (D,(l) = o on 5, (19)
n't y x x
n
yy ,
" 3
.

Differentiating equation 16 with respect to time, and

eliminating n. from equations 16 and 18 yields the well

known linearized free surface condition

1) + =
g$ ° on y = ° (20)
y *tt* ' •

In equation 20, because of the use of the Taylor series

expansions in n, the functions are evaluated for y = 0,

rather than on the actual free surface.

Following the same procedure, as before, for equations

17 and 19 yields

g
S
*(2) + ,(2) = 1 v
#
(D (# (1) + s
g$
(D);
y tt g t ^tty yy

- 2(^ (1)
x
^,1)
xt
+ t^U^h
y yt
. (21)

Note the nonlinear function of $ on the right-hand side of

equation 21.

2. The Ship Interface Condition

Since the ship boundary is also in "motion", a

similar procedure must be applied there to "fix" the boundary

Sp. The development follows to some extent that of Potash

[31]. Consider figure 3 in which the general ship hull form

31
 8-1 s °«
s s>

y(s) A
x(s)

Ship Hull.S,

V
S-0

FIGURE 3. Hull Coordinate System

is described in two-dimensional cartesian coordinates in

terms of the parameter s (arc length) in the following way

x(s) y = y(s) . (22)

The inertial coordinate system is the same as shown in

figures 1 and 2. The parametrically represented body

coordinates (equation 22) correspond to the inertial coordi-

nates of the body when the ship is at rest (neutrally buoyant

position) . The motion of the ship in this reference frame

can be defined in the following way

x(s,t) = x(s) , (23)

y(s,t) = y(s) + e y h (t) (24)

32
1
\
where the function y, is defined to be
h

lat
yh = Re{be } . (25)

In equation 25, a is the circular frequency of the harmonic

motion and 1 is the imaginary unit.

The velocity vector for the ship referred to the

inertial reference frame is

V = e yh j . (26)

The dot denotes ordinary time differentiation. The unit

outward normal vector (from the region) in terms of the

parametric coordinates becomes

n = -y'l + x'j , (27)

where the prime indicates differentiation with respect to

the variable s. Substituting equations 26 and 27 into

equation 7 yields the ship-fluid interface boundary condition

?*«n = ey
h
x' . (28)

Now expanding $ in a Taylor series in e (similar to free

surface expansion) yields

2
*[x(s,t),y(s,t),t] = $(x,y,t) + ey $ (x,y,t) + 0(e ) . (29)
h y

Using equations 10 and 29, the normal derivative of $ on S

2

may be expressed as

33
 (1 (1)
$$«n = (-y'<l> ' + x'<i> )e
x y

2) 2) 5
+ [-y f

^ + ^ ,$ - y'y^xy
+ 7 '
yh$ )]6:2 + 0(e3) »
y yy

on S . (30)
2

Substituting equation 30 into equation 28 and equating like

powers of e yields the following set of equations valid on

the fixed boundary S ? :

(1)
V<D -n = y x"'
h , (3D

v*
(2
U = -x'y
h «<J>
+ y'y
h ^\ (32)

Equation 32 may also be written

(1)
V* (2) -n = -y (h.V$ )y . (33)
h

It remains to consider the boundary S^ (radiation

boundary) and the function $*. For consistency, $* must

also be expanded in a perturbation series

<i>* = e $* (1) + e
2
$* (2) + 0(e 3 ) . (3*0

Potash [31] and others have shown that $* and $* , for

a region of infinite depth, each represent a simple harmonic

traveling wave of appropriate frequency and wave length.

Zienkiewicz and Newton [43] have provided the appropriate

homogeneous boundary condition for the surface Sh to success-

fully pass a wave form without reflection. The equations are

34
 $M<D.S = - J- .(!> , (35)
C
l

^« (2) .S = - i
C
•(25 . (36)
2

In the above equations c-, and c (the wave celerities) are

?

defined by the following implicit relationships

c, = (g/a) tanh (2&) , (37)

la {s -.

c
p
= (g/2a)tanh(^-) . (38)

3. Perturbation Formulas for Pressure,

Force and Wave Amplitude

The determination of dynamic loads on the ship

requires the following equations which follow from the pertur-

bation series assumed for $ and n

(0)
(x,y) + ep (1) (x,y,t) +
2 2)
p(x,y,t) = p e p
(
(x,y ,t ) + 0(e 3 ) ,

(39)
(2)
F(t) = eF (1) (t) +
2
e F (t) + 0(e 3 ) . (40)

In equation 40, F(t) is defined as the total force per unit

length required to sustain simple harmonic vertical trans-

lation (see figure 3).

Following a similar procedure as before, evaluating

the Bernoulli equation (equation 3) on the ship's hull using

equations 23, 24, and 29 yields

35
 1 1
p[x(s),y(s),t] = gy(s) + e[gy (t) + •f *(x,y,t)]
P h

+ e
2
[y h
^ }
+ *£
2)
+ | ({^
1}
}
2
+ {^
1}
}
2 3
)] + 0(e ) . (41)

By comparing equation 39 with equation 41 we find

(0)
P = - Pgy(s) , (42)

(1) 1}
P = - P[gy h + *£ ] , (43)

P
(2) -- P [y
h *<J>
+ |({^} 2 + C^} 2
) + *< 2 >] .

(44)

At this point observe that the second-order pressure is made

up of contributions from both the first and second-order

potential solutions.

Potash [31] and others have provided relations for

the force per unit length acting on the cylinder. However,

in this treatment, it is desired to obtain relations for the

exciting force as previously defined. The motive is that the

theoretical result obtained should correspond to that

measured by an appropriate physical experiment. The relation

between the exciting force F, harmonic displacement y, and

hydrodynamic pressure p, is

F(t) = ep2Ay, - / pn-j ds ,

(45)
-I

36
where I is one-half of the total symmetric arc-length of the

ship's hull (see figure 3). It should be noted here that

Potash [31] has examined the question of the time-dependent

total arc-length of the actual wetted hull and its effect on

the limits of integration in equation 45. However, It will

be stated without proof that, provided the sides of the ship

are vertical at y = 0, the second-order effect is zero. The

only hull form examined to the second order satisfies this

stated requirement. Therefore, equation 45 is exact for the

case considered.

In the calculation of the force F, the constant

hydrostatic lifting force which comes from p will not be

considered since it is not in the definition of F. Accordingly,

the following expression for F(t) results using equations 27,

39, 43, 44, and 45

1)
F(t) = e[p2Ay - -p(gy + ds
h ~
/
At
h ^ )x« ]

2 1)}2 + 1)}
+ e [-P
~*
/ -(y
n ^y
}
+
l [U x {
*y ^ + $
t
2))I ' ds] *

(46)

Comparing equation 46 with equation 40 yields the following

set of equations

F (1) (t) = 2Ay + (gy h +

1}
ds (47)
P
h p /
^ )x' ,

F
(2)
( t) - P Avg }
+ 1«4 1} >
2
+ <*< 1} >
2
) +
^
2)
]x« ds .

(48)

37
 The wave amplitude may be written by inspection of

equations 16 and 17.

(1) I 1}
(x,t) = -
n
^ , (H9)

2 2 2
+ n (1) 1) 1} 2
n< >(*,t) = - I[*< > )
+ |({»< ) +
(^ > )]
^J .

(50)

The boundary value problems for $ and $ may now

be formally stated. However, at this point, a standard

separation of variables technique is applied to eliminate

time from the solution (since the solution is assumed

periodic) . The appropriate definitions are

lat
*
(1)
(x,y,t) = ReU (1) (x,y) e } , (51)

iat
*
(2)
(x,y,t) = ReU (2) (x,y) e } . (52)

Note the newly defined spatial potential functions <}> and

(2)
<f>
are in general complex. Also because of the definitions

in equations 51 and 52 it is necessary to replace y, by the

complex form

iat
y* = be . (53)

38
C. THE BOUNDARY VALUE PROBLEM FOR (1)
<t>

The boundary value for <fr '

is formally stated using

equations 9, 13, 20, 25, 31, and 35.

2 (1)
V <f>
= , in R ,
(5l|)

g4»y
1}
- a V 1}
= , on S
3
(55)

(1)
V"<J> .n = iabx' , on S (56)
2

^(l).S = i2|^l , onS ,H (57)

C
l

(1)
V<f> .n = , onS 1 and S
5
. (58)

2)
D. THE BOUNDARY VALUE PROBLEM FOR 4>^

(2)
Similarly, the boundary value problem for 4> may be

formally stated using equations 9, 1^, 21, 25, 32, and 36.

2 (2)
V (j> = , in R , (59)

- ia[{0>^ 1) }
2
+ {^
1)
}
2
] , on S (60)"
2

V 2
>-n = ^-^ , o„S 2 (61)

The factor H on the right-hand side of equations 60 and

6l is due to complex algebra involved and is explained in
Appendix A.

39
 V 2
>-S = -2i|iHi , onS, (62)

2
?V -n = , onS 1 and S
5
(63)

Note the striking similarities between the first and second-

order problems. Note also, that the nonlinearities in the

original problem have been approximated by the nonlinear

terms in <J>
in equation 60, and the nonlinear terms in y,

and cj) in equation 6l (y, does not explicitly appear) .

2
The additional nonlinear amplitude dependence comes from e

in the original perturbation expansions (equations 10 and

11) . The structure of the problem requires the solution of

<J>
first in order to define the solution for <J>

The conversion to complex algebra allows another observa-

tion to be made. Since only the real part of <j> is

meaningful, Appendix A demonstrates the required complex

manipulation. Close examination of that result indicates

(2)
that the nonhomogeneous -boundary conditions on <p contain
(2)
time-independent terms. This implies <t> has a time-

dependent and a time-independent part. However, since the

physical quantities of interest in the solution all involve

(2)
derivatives of f(2) , the time-independent solution of <J>

5
will not be pursued further. The above comments are also

true for the second-order functions of pressure, force, and

Lee [18] comments on the time-independent solution in

terms of mass transport.

HO
wave amplitude. However, the time-independent quantities

are of interest and have physical meaning as will be seen

later.

The boundary value problem for <p

K
is valid only in
i

water of infinite depth. The following discussion will

explain the problem and develop a boundary value problem

(2)
which will allow for the solution of <Jr in finite depth.

E. THE SECOND-ORDER BOUNDARY VALUE PROBLEM

WITH FINITE DEPTH

Potash [31] has shown that the asymptotic form of the

velocity potential $ as x -* °° becomes (in water of infinite

depth)

f(-.fl.t) = Re{ E G e i C° t
1
- (k x
l ^ (1) " + ^^^t-i^^y^hl >

(64)

where k.. denotes the wave number of the first order wave

(k 1 = a/c-, ) and G-, and G 2 are coefficients which depend on

frequency a. However in the more general case of finite

depth, an additional term due to the Stokes' second-order

potential appears as follows

*/ a *.\ u u i[ot-(kix+Y )]
,J 2 ru -i(kpX+Y )
<K°°,.0,t) = f
Re{eH-.e 1
- '
+L e [Hpe £

2
- 3 _3£^ H 1 cosh(2k 1 h) 2i(k 1 x+Y (1) )n 2iat, . (65)
e e
8c^ sinh 2 (k h) J '
1

The Stokes potential to the second order is well known

to be zero in infinite depth.

41 I
\
In equation 65, k~ is defined as the second-order wave number

(kp = 2a/c
2
). H, and Hp are constants similar to G-, and G ? .

Close examination of equation 65 reveals that the radiation

boundary condition (equation 62) would fail since the

"Stokes wave" which is present is traveling at the same

celerity as the first-order solution. This means that the

radiation boundary would "reflect" the Stokes* wave since

the boundary condition defined by equation 62 will only pass

a wave traveling at a celerity of Cp, defined by equation 38.

A successful resolution of the above described difficulty

requires the formulation of a different boundary value

(2)
problem. We define the complex velocity potential i> as

«£_ H
2
c osh[2k (y +h)] ^(j^U))
t
(2) . .
2
1
2
1
^ ( ^ &)
Qc
±
sinh (k h)
1

which corresponds to the third term in equation 65. Then,

following previous conventions, the complex velocity potential

_(2)
<|) is defined as

*<*> - <*> 2
- *< > . (66)

1. The Second-Order Boundary Value __

Problem for <j> ;
J

(2)
The definition of "$"
in equation 66 allows the

formal statement of the following boundary value problem,

(2)
which follows from equation 66, the definition of ^ and
(2)
the boundary value problem for <}>

42
 2 (2)
V ,f>
=0 , in R , (67)

bj<
2
> - toV 2) = i ^ 2g
[«\!
Y
y
1)
* &y
8 il!»]
yy

- iali^h 2 +
{J
1 '} 2
] - g*( 2 ) 4aV 2) , on 83,

(68)

btlV b*<!>
V 2)
.S-f (^S- +
,„
^)) _x'(^ ^ + 2)
) , onS 2>
(69)

*<2>.J - _ 2i£i^ , on s , ( 7 0)
4

+_(2); > ± (2) •>

Vcj)^ -n = - Vij/^-n , onS. , (71)

(2)
^cf) -n = , on S^ . (72)

Equation 70 will now properly behave because the asymptotic

—(2)
v
form of (j) using equations 65 and 66 is

2i0t - -2i(k 2 x +Y (2) ) 2iat

T<2>(., y) e H e e (?3)
2

Note the celerity of the complex wave form ~§ e is c~

and the radiation boundary condition will now work.

F. FORMULAS FOR PRESSURE, FORCE, AND WAVE AMPLITUDE

The development which follows requires additional defini-

tions, necessary because of the separation of variables scheme

previously introduced. They are:

43
 (1) lot 2iot
Re{f (1) e
(2
p (x ( y,t) = } , p > = Re{p< 2) e } , (7*)

lat 2iat
F (1) (t) = Re{f (1) e } , F
(2)
(t) = Re{f (2) e } , (75)

(1) (1 lot (2) 2) 21ot

n (x,t) = Re<n >(x)e > , n = Re{n< (x)e }. (76)

The functions p^ , p , rf , n" are complex. This is

also true for the constants f and f

Using the above equations and equations 43, 44, 47, 48,

49, and 50 the equations for pressure, force, and wave

height for both first and second-order solutions follow.

(1) (1)
p = -p[bg + ia<{> ] ,
. (77)

p(2) . _ p[ ^4^ + ^(< 1) >

2
+ {^
1}
}
2
) + 2icV
2
>] , (78)

f
(1) = -2pbAa
2
+ p / (gb + lo^ 1 ')** ds , (79)
-I

l0
*y 1)b
f (?) = _p /( ~2 +
'
2
l
Llv x h
kui 2
• + i^h
l
*y
2
i + 2io^ 2 ')x- d S>

(80)

(1)
-(1) = _ lof (81)
g
1)}2 1)}2
(1) (1)
-(2) = - i1 (2ia<f>
f2l
(2) +
ian
iqn fy +
d)
k—^
1
{
*x
+ —*v5L
{

]) .

(82)

Again, equations 78, 80, and 82 include the special considera-

tion of the complex algebra involved which is explained in

Appendix A.

44 \ \
 III. FINITE ELEMENT DISCRETIZATION

The Finite Element Method and its application to scalar

field problems is described in various texts, including one

by Zienkiewicz [45]. The essential ideas involved follow.

The region of figure 2 is replaced by a grid or mesh of

finite elements as shown in figure k.

FIGURE 4. Finite Element Mesh

Then it is assumed that the field function sought can be

represented by the following equation

4>(x,y) = N , in R . (83)
J

The total number of nodes in the region will be defined to

be m. Then N is a lxm row vector of shape (interpolating)

functions and <J>

is an m x 1 column vector of nodal values of

the complex potential $ (in this problem) . Equation 83

implies that knowledge of the vector <J>

defines <j>(x,y)

everywhere in R.

45 1 \
 There are several choices of shape functions N. which may

be used in equation 83. This work selected the cubic shape

functions described by Zienkiewicz [*J5] as the "serendipity

family". The elements used are isoparametric. A brief

description of the 12-noded isoparametric elements is given

in Appendix B for the unfamiliar reader. The key advantage

of the isoparametric element lies in the fact that regions

with curved boundaries may be readily represented.

The development that follows defines the set of linear

equations which determine the vector <J>

of equation 83.

One general approach to discretization of several types

of field problems in finite element analysis is to define a

functional whose integral over the domain of the problem is

to be minimized. The Calculus of Variations shows that if the

Euler equation of the integral to be minimized is the governing

equation (i.e., Laplace's Equation) then the two problems are

7
equivalent.

However Zienkiewicz and Newton [43] have shown the

discretization may be readily accomplished by applying the

Galerkin process. The weight functions are chosen to be the

nodal shape functions. Mathematically stated

T
/ N U XX + d> ) dR = . (84)
R

This statement assumes the same boundary conditions are

used.

46
The superscript T denotes transposition and the function <J>

is considered to be any of the complex potential functions

whose solution is sought. Application of the divergence

theorem transforms equation 8*1 into

' [ fyl dR = / N V<J)-n dS , (85)

A
!?x y ~
R S

where S denotes the boundary of the region R.

A. THE DISCRETIZED FIRST-ORDER PROBLEM

Application of equations 84 and 85 to the first-order

boundary value problem (<j) ) stated in Chapter II yields

the following matrix equation

(1)
iabr (1)
2 i£ D
(_ Q + + H) * = . (86)
1 =

The definitions of the matrices in equation 86 follow. The

matrix H is m x m, real, symmetric, and defined as

T T
H = / [N N ] N dR (87)
X ~y ~x
R ~
N„

The matrix H comes from the integral on the left-hand side

of equation 85. The matrix D is mxm, real, symmetric, and

defined by

T
D = / N N dS . (88)
z ~
S„ ~

47
 .

D is contributed by the boundary integral on the right-hand

side of equation 85 and represents the discretization of the

radiation boundary condition (equation 57)

The matrix Q Q is real, -m x m, symmetric and given by

T
QQ = / N N dS . (89)
z I s
3

Q parallels D in its origin and represents the contribution

~0 as

from the homogeneous free surface boundary condition (equation

55).

The vector r is real and defined by

KX) mT-,
v (1) = fr N x» dS . (90)
S.
2

The above vector is the discretized ship-fluid interface

kinematic boundary condition (equation 56). The surfaces

S, , S~, S|, and S,- have homogeneous boundary conditions and

therefore make no contribution to r . It follows that r

is non-zero only at nodes along the ship's hull (surface Sp) .

B. THE DISCRETIZED SECOND-ORDER PROBLEM

-(2)
The problem for <j>
v
only will be presented since it
K 2)
(
contains the solution for ' in the limit. The boundary
<p

—(2)
value problem for <fr yields a linear system of complex

equations similar to equation 85.

2 2ia -(2) = (2)

(-Ho Q„ +
^O Cq~D + H) r
~ 2~ ^ (91)

H8 (
 )

The right-hand side vector r (2) in equation 80 receives

contributions from the surfaces S, , S~ and S^. Accordingly

a further decomposition is necessary:

^(2) . ^(2) + r U) + ^(2) _ (92)

(2)
The vectors r^ are defined below. First to be consistent,
(2) —(2)
the function ty must be represented in the same way as <$>

(2) {2)
* = N i> (93)

( 2)
Then v^. becomes, from equations 71 and 85,

2) =
T (2)
r{ / N N ij, dS . (94)
^1

(2)
The nodal values of ip are determined by partial differen-
(2)
tiation of the function ty and then representing the
(2)
function \p in similar manner as equation 93.
(2
The vector v\ is defined by equations 69 and 85

(1) (1)
r<
2) = / N
T
Ny"'[^—
cj> b
+ ij/
2)
]dS - /
T
N Nx'
<f>

l=M — b
+ ^ 2) ]dS .

2 2
(95)

Similarly r~ ' follows from equations 68 and 85

j(2) = I /

3
T
n N[- ^ a + * g - lay -
g^
2)
+ 4c
2
^
2
)] dS .

(96)

49
1
\
The i component of a is

i Yi y
i,y (97)

where <J>. and <f>. are the i components of <Jr ^

and ~y ,
l l ,y ~ <J>

respectively. Similarly

H
i y
i
y
i,yy
and (98)

Tl <>> a
<>> 2
(99)

The nodal values of the various partial derivatives of <f>

are determined by differentiating equation 83. The method

of calculation is discussed in the next section. The reader,

unfamiliar with isoparametric elements, should refer to

Appendix B prior to reading the next section.

1. Calculation of Nodal Values of the

Field Derivatives of ft
OT
1
The vectors ft^ \ K >
*xv >
and
^vv
re Q uired in

equations 95 and 96 are determined by applying standard

concepts of finite element analysis. The appropriate relation

for the i node of a given element is

(l)e
* i,x
= J
-1 rf (De (100)
4>
(De N
e

50
«
J is the 2x2, real Jacobian matrix of the transformation.
The vector N is a 1 x 12 row vector of element level shape

functions. The superscript "e" denotes all quantities in

brackets are at the element level. In equation 100, the

-1 e
elements of J and N are evaluated at node i.

Application of equation 100 to all twelve nodes of a

(l)e (l)e
given element determines the vectors $ and <}>

~x ~y
Then, replacing <$>
(De by
(De in equation 100 yields the
<J>^

components of the vectors

(De and (l)e
. Mathematically
<f>
v
<f>Jt

stated

<$>
(De N
e
i,xy
r -i (De (101)

4>
(De
i,yy

51
 .

IV. THE COMPUTER PROGRAM AND NUMERICAL SOLUTION

A. MAIN FEATURES OF THE COMPUTER PROGRAM

The formation and solution of the linear system of equations

defined by equation 86 ( <jr system) or equation 91 (<J>

system) was accomplished on an IBM 360/67 computer at the Naval

Postgraduate School. The computer has a maximum word length

of sixteen bytes.

A comparison of equations 86 and 91 demonstrates the close

similarity of the coefficient matrices for the first and

second-order solutions. Note that only the coefficients of

Q and D change. An additional observation follows from

examination of equations 87, 88, and 89 defining H, D, and

Q respectively. Note that all three matrices are purely a

function of mesh geometry. Therefore, once a mesh is chosen

to define the region of interest, the three matrices previously

mentioned need only be calculated once regardless of the forced

motion

A well known property of the finite element method is that

it usually produces a banded, symmetric, coefficient matrix.

That is the case here and this has obvious computational

advantages. However, in this problem, the final system of

equations is a function of the frequency of motion as for

example in the first-order problem

K^V 1 ' = iobr (1) . (102)

52
The matrix K is an m x m complex matrix defined by equation

86 to be

K (1) = -a Q^ +
z z°
2
—D
c
i ~
+ H
z
. (103)

Although K^ is complex and frequency dependent, it is

possible to perform most of the triangular decomposition of

K ' (by Gauss elimination) only once, using real arithmetic,

for a given problem geometry. This economization was first

employed by Chenault [4] and is described in more detail in

Appendix C. In addition, this author employed a vector

storage of the upper triangular portion of K , which greatly

reduced the storage requirements of the program. The vector

storage scheme made it possible to store only the portions of

K ' which are non-zero.

~ o
The net advantage of employing Chenault ' s scheme and the

compact vector storage was that the system of 733 complex-

equations which is generated by the mesh shown in figure 4

could be solved for both the first and second-order problem

twenty times in twenty minutes. Further, only ten percent of

the total matrix was actually stored and the computer solution

was accomplished in core.

o
The Chenault scheme described in Appendix C was developed
by Professor R. E. Newton at the Naval Postgraduate School.

53
 .

B. NUMERICAL SOLUTIONS

The solution of the first-order problem for a given mesh

geometry is required before the second-order solution can be

obtained. Further, the first-order solution must be highly

accurate if a good second-order solution is to be expected.

This fact is made more obvious by re-examining equations 68

and 69 in part II

1. First Order Numerical Solutions

Several first-order problems were solved and closely

examined for accuracy. The numerical solutions are presented

in part V.

a. Mesh Requirements

One of the primary considerations in choosing a

proper mesh concerns the positioning of the radiation

boundary (Sj,). It may be recalled that the surface Sh must

be far enough away from the ship hull (location of disturbance)

so that the local effect of the disturbance has decayed and

only an outgoing wave is present. Bai [3] has addressed this

question analytically. He has shown that the eigenfunctions

associated with the non-propagating velocity field decay

exponentially with x. Bai recommends the criterion

w
"
= ah + b ,
(104)
min *

where a is a monotonically increasing function of frequency.

The range of a is from 1.5 (low frequency) to 3. The values

of a were determined based on the decay rate of the most

persistent eigenfunction and an attenuation factor of 0.01.

5^
 .

This investigation has found the above criterion somewhat

conservative. A value of a = 2/3 was found to be adequate

for the velocity potential to numerically approach its

asymptotic limit for all frequencies studied for the semi-

circular hull (h/b = 6). Similarly, a value of a = 9/10 was

found to be adequate for the bulb hull form described in

part V (h/b = 10) for all frequencies studied. The values

of a specified above were determined for the deepest regions

studied for each hull form. The region widths used in this

case (w/b = 6, w/b =10) were found to be more than adequate

for shallower depths

Another requirement on w which tends to restrict

its maximum possible value comes from the practical limits on

mesh fineness and the resulting computer storage requirements.

The mesh must be sufficiently fine to represent properly a

traveling surface wave. The finer the mesh, the larger the

computer storage requirements and the longer the computer

solution time. Therefore, this latter restriction tends to

oppose the former one.

Bai [31 used quadratic elements and recommended

at least 10 surface nodes per wave length. Visser and Van

der Wilt [39] also used quadratic elements and recommended

that an element could not span more than one quarter of a

wave length. This investigation examined this question in

detail. As previously mentioned, cubic isoparametric elements

were used. It was found that if the gravity wave did not have

to be propagated more than two wave lengths, the elements

55
could be allowed to span six-tenths of a wave length.

However, as in the case Visser and Van der Wilt [39] studied,

the successful propagation of a gravity wave over many wave

lengths requires a finer mesh spacing. The one quarter wave

length Visser and Van der Wilt recommend was found to be a

practical upper limit.

A parallel problem to positioning the radiation

boundary concerns the location of the bottom (surface S._)

to simulate infinite depth. Comparison of the finite element

solutions obtained in this study with analytic solutions

based on infinite depth led to the conclusion that a keel to

bottom clearance (h - d) of 5b will adequately simulate

infinite depth. However, satisfying the classical hydro-

dynamic requirement that the depth must be at least one-half

of a wave length becomes impractical at low frequencies.

Ignoring this requirement did not seem to affect solution

accuracy.

b. Mesh Data Preparation

The meshes for a given boundary geometry were

prepared using a very versatile mesh generator computer

program developed by Adamek [1] at the Naval Postgraduate

School. This tool greatly expedited mesh preparation. For

example, the entire data for the mesh shown in figure H

(733 nodes, 132 elements) could be generated in twenty

seconds. A modification was used to supply accurate values

of the nodal coordinates along the ship-fluid interface.

56
I
\
 c. Convergence of First-Order Solutions

Convergence of the first-order solutions could

not be determined by standard numerical techniques (i.e.,

mesh refinement) due to computer size limitations. However,

the excellent agreement between the first-order solutions

obtained in this study and both theory and experiment left no

doubt that the solutions were quite good.

2. Second-Order Numerical Solutions

The successful solution of the first-order problem

made study of the second-order problem possible.

a. Mesh Requirements

All of the requirements on mesh geometry previously

mentioned for the first-order solutions of course apply to the

second-order solution. However, the second-order problem

demands even more stringent requirements. One: reason for this

follows from the fact that the wave lengths of the propagating

portion of the second-order velocity potential are .much

shorter (one-fourth of the first-order wave lengths in infinite

depth solutions). Therefore, since it is impractical to use

a different mesh for each solution, the mesh must have a free

surface element spacing approximately one-fourth of that

required for the solution of <j> for a given frequency range

to be studied. The mesh shown in figure k was used to solve

the infinite depth second-order problem. It contains 132

elements and 733 nodes.

57
 b. The Second-Order Boundary Condition on S_

Additional constraints on mesh geometry were

found to be required. These constraints arise from the

requirements to calculate second-order partial derivatives of

<J)
which appear in equations 95 and 98. It was essential

to have rectangular elements with equal node spacing along

the free surface for success in representing the non-

homogeneous boundary condition given by equation 96.

The further discussion of this point will be

aided by the definitions of some complex functions. Define

Q(x ,o) = isgl r-c

2
^ + g $>] - icru^) 2 + u^) 2
: .

(105)

R(x,0) = g^ 2) - 4oV 2) .
-
(106)

In terms of these definitions the right-hand side of equation

68 is

P(x,0) = Q(x,0) - R(x,0) . (107)

The reader should recall that the successful solution of the

(2)
boundary value problem for <J>
using the method of finite

elements required the definition of a new boundary value

—(2)
problem for the function <$>
K
. The reason was that in the

finite depth problem a second-order Stokes' wave is produced

with celerity c, as
,S well
W' as a second-order wave with celerity
(2)
Cp (contained in IJr ). The presence of the two waves with

different celerities meant the radiation boundary condition

58
would fail. The problem was mathematically eliminated by

subtracting the well-known analytic solution for the Stokes

(2)
wave portion of the second-order solution of <p ' on every

boundary of the region R. The mathematical steps are straight-

forward as shown in Chapter II. The asymptotic value of the

function P(x,0) is zero as x gets large. But, the function

R(x,0) which represents the classic Stokes second-order wave

free surface boundary condition has the form

310* (l) 2 -2i(k 1X+ Y (1)

R(X) 0) = |
t (Wj0 )| e )
> (108)
g sinh k-,d

everywhere on S~. Therefore, if the numerical representation

of the function P(x,0) is to approach zero near the radiation

boundary; the numerical representation of the function Q(x,0)

must approach R(x,0) or

Q(x,0) -*-
R(x,0) as x * °°
(109)

Note that the amplitude and phase of <f>

must be known to
(2)
completely define R(x,0) and therefore ijr .

The above discussion points out the necessity of

a very accurate numerical determination of the function Q(x,0).

— (2)
The second-order solutions for <f>
were closely examined near

the radiation boundary and it was found that on the average

the function Q(x,0) was within one-percent of the analytic

function R(x,0). This accuracy was achieved even though

formation of Q(x,0) requires calculation of second derivatives

59
of <f>
using element shape functions which guarantee only

C° continuity across inter-element boundaries.

c. The Ship-Fluid Interface Second-Order

Boundary Condition

The accurate representation of the ship-fluid

interface (S
?
) boundary condition for the second-order

boundary value problem (equation 69) is of course essential.

Potash [31] encountered difficulties in making a numerical

calculation on this boundary because the potential singulari-

ties which define his solution lie along S ? . This results

in the required functions 4> and <J>

being only piecewise

continuous. This problem required an additional approximation

to be made.

The method used herein also has the problem that

T
<f> , Y
<P > Y
<p j and are not continuous across element
T
<J>
x ' y » xy '
yy
boundaries. However, no singularities are present. Extensive

investigations in this work have shown that it is essential

that the elements along the hull be as square as possible,

and further that the node spacing on the elements along the

hull should be uniform. It was found that if the hull

elements were distorted (not square or rectangular) the

numerical representation of the functions <J>

and <b was

highly unstable. However, the mesh requirements for

numerically calculating <j> and <K along the hull were

much less stringent.

The adequacy of the hull elements used in this

investigation (see figure 4) was tested by defining a known

60
1
 . .

potential function on the elements along S

?

* =
| y
2
+ i xy (110)

and then applying the computer program's numerical scheme

t t
to calculate the functions <J>
and $ . Observe that

*L
xy
= i > (in)

and

4=1. (112)

It was found that the cubic, isoparametric elements along S

?

produced an average error of about one percent in calculating

the functions and T

<\> . The maximum error observed was
<J>
xy yy
five percent at isolated nodes.

d. Convergence of Second-Order Solutions

Again, computer capacity limitations prevented

the application of systematic mesh refinement to determine

convergence of the second-order solutions. However, it was

possible to refine the elements next to the ship's hull and

no significant change in the solutions was observed (the

elements along the free surface were performing properly)

Further, it is shown in Chapter V that second-order results

for infinite depth are in satisfactory agreement with those

of other researchers. No comparison is available for finite

depth solutions, but a smooth transition is observed as the

depth becomes effectively infinite (See Chapter V)

61
 .

V. NUMERICAL RESULTS

A. DEFINITIONS

The definitions which follow are useful for presentation

of results

1. Added Mass and Damping Coefficients

The complex amplitude F^ of the first-order

resultant hydrodynamic force acting vertically on the ship's

hull is

pd) = / p^x' ds . (113)

y -I

The real part of F is in phase with ship's acceleration

and the ship's added mass is defined by

Re{F^ 1} }
ma (114)
=
c
/a ,

2
where the ship's acceleration amplitude is given by a a.

The added mass coefficient is defined to be the ratio of

added mass to ship's displaced mass

m Re{F (1) }

m 2pA 2
2pAaa

The imaginary part of F^ is in phase with the

Q
velocity of the ship and is the "non-conservative"^

The potential theory is, of course, conservative. However,

the waves which propagate to "infinity" represent work done by
the ship which is effectively lost.

62
 .

component. Accordingly, the damping coefficient is defined

to be

1 *)
-ImCF^
C = 1 . (116)
d
2pAacT

C and C, are frequency dependent for a given hull form and

region depth.

A dimensionless frequency parameter is defined as

follows
2.
6 - 2g2. . (117)

The parameter 6 is commonly used by other researchers and is

used herein for ease of comparison of solutions . The

parameter 6 is equivalent to —j— in infinite depth (A is the

wave length)

The added mass and damping coefficients previously

defined are for heave. The same coefficients for sway motion

are similarly defined. The only difference in the definitions

is F is replaced by the horizontal force amplitude F

y x
defined by

(1) y~'
F (1) = - / p ds , (118)
x
-I

where p^ '
is the hydrodynamic pressure due to the swaying

motion. There is a coupling between sway and rolling motion

for a general ship cross-section. However, the only hull

form studied in the sway mode was semi-circular in

63 i \
cross-section and the coupling force for this case is zero.

Therefore, the appropriate coefficients typically used are

not defined herein.

One further comment is warranted here. The second-

order time-dependent hydrodynamic forces are harmonic with

frequency 2a. Therefore the definition of second-order added

mass and damping coefficients is meaningless.

2. Dimensionless Force Coefficients

The force coefficients f and f are also

presented in dimensionless form. The non-dimensionalizing

p
factor for the various force amplitudes is 1/pgb . Table I

presents the correspondence between the dimensional and non-

dimensional force quantities. The second-order quantities

(2) (2)
f and fj, which appear in table I are amplitudes of the
s d
time-independent (static) and time-dependent (dynamic)
(2) (2)
portions of f respectively. The existence of f^ is

readily verified by examining equations 75 and 80 in light of

the discussion in Appendix A.

TABLE I

DIMENSIONLESS FORCE COEFFICIENTS

(1) 2) 2)
Force Amplitude f f£ f^

F (1)
2) (d:
Dimensionless F^ F.
d
Force Amplitude

64
 3 . Wave Amplitude at Infinity

Another suitable measure of the damping of ship's

motion (to the first-order) is the wave amplitude at

"infinity" (i.e., near the radiation boundary). Further,

experimentally it has been found to be an easier quantity

to measure (particularly for high frequency motion since the

damping force is very small then). Accordingly, the dimen-

sionless first-order wave amplitude at infinity is defined

by

-(l>_]^_!™ ( ll 9)

B. HULL FORMS

Other investigators have studied a variety of hull forms

both by experiment and analytically. The semi-circular hull

was' chosen because of the wealth of data available for

"comparison of both first and second-order solutions. The

bulb-shaped hull studied by Paulling and Porter [29] was

also examined for first-order solutions. The shape of the

bulb hull in comparison to the semi-circular hull is shown

in figure 5.

The coordinates of the bulb hull studied are defined by

the following mapping

e e e
x + iy"=^ + -^+-4+-| • (120)

65
 1

Centerline
Waterline

Semi-Circle, d/b -

FIGURE 5. Hull Forms

The complex number r, represents the points on a circle in

the £-plane. The choice of radius in the c; -plane defines

the scale of the hull cross-section. The values of the e.

are given in table II.

TABLE II
CONSTANTS FOR HULL FORMS

Hull e e e A/bd
l 3 5

Semi-Circular 0.0 0.0 0.0 0.785^

Bulb -0.7132 -0.02096 0.0 0.6950

66
C. FIRST-ORDER NUMERICAL SOLUTIONS
1. Semi-Circular Hull

The semi-circular hull form has been studied by many

researchers mentioned previously in part I. Specifically,

the original solution by Ursell [38] for heaving motion has

been confirmed analytically by Porter [30], Vugts [40], and

others. Ursell's solution has also been experimentally

confirmed by Paulling and Porter [29] and Vugts. The

analytic
J results of Ursell for C and C-, are shown by
J the
m d
continuous curve in figure 6. The numerical results obtained

by the finite element method (PEM) in this work are repre-

sented by the circled points in figure 6. The FEM points

shown were obtained using a mesh similar to that shown in

figure 4, but using only 62 cubic elements and 359 nodes.

The extension of the solutions for C and C, to

m d
finite depth is easily accomplished by simply redefining the

mesh geometry (decreasing the vertical dimension) . Kim [14]

and Yu and Ursell [42] reported analytic results. The

results obtained by the FEM in this work for h/d = 4 agree

with those reported by Kim for C and n*^ within three

percent over the range 0.25 < 6 < 5. The results of Kim are

in disagreement with those obtained by Yu and Ursell for

h/d = 2. The results obtained herein fall between the two.

Bai [3] reports inability to verify Kim's result, but gives

no details. Bai used the finite element method to obtain

his solutions as mentioned in part I.

67
 FIGURE 6. Added Mass and Damping Coefficients in
Heave, Semi-Circular Hull, Infinite Depth

Figure 7 shows the effect of finite depth on C

obtained in this work. The figure shows that significant

changes do not occur until very shallow water is encountered

(h/d <_ 2) . Further, the effect of shallow water on C is

more pronounced at very low and very high values of 6 . Note

the unusual variation with depth for 6=0.5. C decreases

with decreasing depth (opposite the trend at higher values

68
\
 1.25

s
4.5-
1.00
3.0 iL
2.510^.
\^^^4^_
20 rC^^r? 1
0.75 l- 5 T^^^_^
m l.oi _, J

0.5<i
0.50

0.25

4 8
h/d

FIGURE 7. Variation of Added Mass Coefficient with

Depth, Semi-Circular Hull in Heave

69 •
 .

of 6). Although it does not appear true in figure 7, C

(for <5 = 0.5) has reached its asymptotic value at h/d = 6.

Figure 8 shows the corresponding effect of finite depth on

C,. The figure indicates that C, is only sensitive to

shallow water effects at very low values of 6

0.75

Cd

0.50

0.25

FIGURE 8. Variation of Damping Coefficient with

Depth, Semi-Circular Hull in Heave

70 \
 .

The semi-circular hull was also studied in sway

motion to the first-order. The only changes in the boundary

value problem for <jr defined in part II are in equation 56

which defines the surface S, . The changes are

(1)
V<J) .n = iabx' , on S , (121)
2

for equation 56, and

(1)
<J>
= , on S , (122)
1

for equation 58. The boundary condition on S,- defined by

equation 58 remains unchanged. The ship's motion is defined

in a manner analogous to equations 23 and 24

x(s,t) = x(s) + e x (t) (123)

h ,

y(s,t) = y(s) . (124)

The function x, (t) is defined the same as y (t)

h
a
i
1 t
z
x
h
= Reibe ] , (125)

and e is defined in the case of sway motion to be

= a
e (126)
b

The symbol "a" is the amplitude of the sway motion in this

case

Vugts [40] presents both experimental and analytic

solutions for C and C, for the semi-circular hull in sway.

m d
Figure 9 shows a comparison of Vugts' theoretical solution

71
 1.0
Theory(40)
FEM
-C
0.8 in
-a h- C
d
h/d = 6 w/b=5
0.6

c m> c d

0.4

0.2

FIGURE 9. Added Mass and Damping Coefficients in Sway,

Semi-Circular Hull, Infinite Depth

72 v \
 .

with the FEM results obtained herein. The FEM calculated

points are shown by the solid circles and squares in the

figure. The solid curve (from Vugts' theory) was scaled

from a very small graph and may be somewhat in error.

2. Bulb Hull

Paulling and Porter [29] reported both analytic and

experimental results for the bulb hull form. They show

excellent agreement by comparing first-order exciting force

amplitude coefficients F^ and wave amplitude coefficients

rj"^ respectively. Figures 10 and 11 compare the values of

F^ ' and rp obtained by the FEM herein (circled points)

and Paulling and Porter's theoretical calculations (solid

line). The FEM mesh parameters were h/d = 2 and w/b = 25.

This value of w/b was found later to be quite conservative

(w/b = 10 would have been adequate)

Figure 12 shows added mass and damping coefficients

for the bulb hull vs h/d for various values of 6. The curves

have a similar trend to those of the semi-circular hull form.

-(1)
The corresponding curves of force F and wave amplitude
v

n^ are not presented. However, the maximum change in

exciting force amplitude does not exceed five percent for

h/d = 1.2 as compared to infinite depth (h/d = 2). Similarly,

the maximum change in rf^ was only ten percent over the same

range of h/d.

3. Wave Generator Studies

MacCamy [20], in 1957, reported an analytic theory

using a Green's function approach for solving the problems

i \
73
 . .
Thrnru (29)

— FEM
O h/d = 2, w/b= 25

p(.)
—

— Ll rf

1 1 1 1

0.5 i.o

8
FIGURE 10. Exciting Force Amplitude in Heave,
Bulb Hull, Infinite Depth

0.5
Theory (29)
0.4 FEM
O h/d= 2, w/b= 25
-in °- 3
r"°~° O O ( r ,

0.2

0.1

1 1 1 1 1 1
i i i i
: 1 I

0.5 1.0 1.5

FIGURE 11. Wave Amplitude for Heave,

Bulb Hull, Infinite Depth

71
 FIGURE 12. Variation of Added Mass and Damping
Coefficients with Depth, Bulb Hull in Heave

of gravity wave production by a moving partition and gravity

wave reflection from a horizontal strip. The problem of

gravity wave production by a moving partition was studied in

this work to further help determine proper mesh characteris-

tics for gravity wave propagation.

The wave generation problem is readily placed in the

framework of the class of problems being studied herein with

75
 S S

only minor modifications. ° All of the boundary conditions

posed for <J>

in part II remain unchanged except those on

S and S . The surfaces S- and S become a vertical line.

1 2 [ 2

The region therefore is rectangular. The left-hand boundary

will be designated as S' as shown in figure 13.

Ay
W -H
x~¥
/ Free Surface,
I
p-Radiation
~^V-Wave Generators! Boundary, 4
/ (flapper type) I

I I

I
{~ B0 " Qm S I
>
5\
-U.
777777777777777777

FIGURE 13. Reduced Region for Wave Generator Problem

The boundary condition on Si becomes

v> (1) -n = q(y) , on Sj_ , (127)

The function q(y) describes the mode of displacement of the

wave generator. The plunger type and flapper type (shown in

figure 13) wave generators were studied and compared with

This observation is not true in the potential theory

singularity distribution approach. The source singularity
must be replaced by doublets and the problem completely
reformulated.

76
i \
MacCamy's results. The function q (y) for the plunger type

generator is simply

q (y) = U , (128)
p

where y is the maximum amplitude of the wave generator

motion. Similarly, the function q (y) for the flapper type

f
generator is

q f (y) = ^j-21 y . (129)

Q

The comparison of MacCamy's results with those obtained

herein is shown in figure 1H for the plunger. The circled

points were obtained by the FEM and the solid line is MacCamy's

theoretical curve. Similarly, the comparison of MacCamy's

result for the flapper type wave generator vs the FEM results

obtained in this work is shown in figure 15. The excellent

^agreement allows two observations to be made. First, the FEM

approach to solution of problems of the type being presented

is highly versatile in its application. Second, the radiation

boundary placement was further verified.

D. SECOND-ORDER NUMERICAL RESULTS

1. Semi-Circular Hull — Infinite Depth

The first-order solutions just presented established

confidence that second-order solutions could be attempted.

The only hull form studied, to the second-order, was semi-

circular. The accuracy of the second-order solutions was

77
 2.0

1.5

HH 1.0

Theory
0.5
o FEM
0,
0.1 0.2 0.3 0.4 0.5
cr 2 h
2irg

FIGURE lH. Wave Amplitude Generated by a

Plunger Type Wave Generator

0.4

( .

0.1 0.2 0.3 0.4 0.5

2
cr h

FIGURE 15. Wave Amplitude Generated by a

Flapper Type Wave Generator

78 I
\
 )

established by comparing with the work of other researchers.

Lee [17,18], Parissis [27], and Potash [31] have reported

numerical solutions for second-order exciting force coeffi-

_(2) —(2) (2)
cients F and phase angle y
, Fg Figure 16 compares the .
d
—(2) —(2)
values of F^ and F^ ' reported by Lee, Parissis, and

Potash with the numerical results obtained by the FEM in this

work. The points presented from the other researchers' works

were all obtained from Potash [31]. Potash indicates Lee's

points are corrected in some manner from his original work

[17] and otherwise unpublished.

Figure 16 shows excellent agreement for F^ and very

5
—(2)
good agreement for F^ . Only Potash's points show signi-

ficant deviation and then only for 6 near 1.5. Potash

reported some numerical difficulties for a range of 6 near 1

and therefore Potash's results are not presented for those

values of 6 which are suspected to be in error. Similarly,

( 2
figure 17 shows a comparison of y (phase angle) results.

The figure shows the agreement is not as good. However, none

of the other works shown (Lee, Parissis, Potash) agree in any

consistent fashion over a large range of 6. The reason for

the discrepancy has not been resolved.

2. Semi-Circular Hull — Finite Depth

—(2)
Figure 18 presents the effect of depth on F^ and
~(v 2)'
F as obtained by the FEM in this work. The author knows
s

of no other solutions. Solutions were obtained for values

of h/d from 6.0 (infinite depth) to 1.5. The figure shows

79
 Theory FEM
Potash(3l)
e(2)
Lee (17,18)

* Parissis(27) —V- F
"V (2)
fi

FIGURE 16. Second-Order Exciting Force Amplitude in

Heave, Semi-Circular Hull, Infinite Depth

80 I
\
 Theory FEM
• PotashOl)

-f Lee(l7,l8) +- r (Z)
* Parissis(27)

-50

-100

-150

-200
0.5 1.0 1.5

FIGURE 17. Second-Order Exciting Force Phase Angle

in Heave, Semi-Circular Hull, Infinite Depth

81 * \
2.0

FIGURE 18. Variation of Second-Order Exciting

Force Amplitudes with Depth,
Semi-Circular Hull in Heave

82
 — (2) —(2)
that Fl and F v
are most sensitive to finite depth at high
J

d s

values of 6 or in very shallow water (h/d < 1.5).

(2)
Figure 19 shows the corresponding effect on y of

finite depth for the same range of h/d. Note that at

(2)
h/d = l.M 6 has negligible effect on y .

Also the effect of decreasing depth causes a monotonic

(2)
variation of y at fixed 6.

FIGURE 19. Variation of Second-Order Exciting Force Phase

Angles with Depth, Semi-Circular Hull in Heave

83
 Figure 20 presents the ratio of F^ 2 '
to f'
1 '

plotted vs 6 for various values of h/d. The curves peak

close to 6 = 1 for all values of h/d studied. The curves

peak because the first-order force F^ '

passes through a

minimum very near 6=1 for all values of h/d studied.

FIGURE 20. Ratio of Second-Order Exciting Force Amplitude

to First-Order Exciting Force Amplitude,
Semi-Circular Hull in Heave

84
»
\
 A quantitative illustration of the physical signif-

icance of the nonlinear effects (second-order) is readily

made using figure 20. Assume a ship with an approximately

circular hull form is being excited such that the maximum

acceleration of the motion is 0.2 g's. The maximum time-

dependent amplitude of the force is given by

(2)
e|f (1)
2
F(t)| = + e |f (130)
i

max i
|
i i |
i ,
>

which may be written

F(t)l py
'max m P.
w(2)
-77m-- iTb£ 1 +
JflT
'
^ (131)

Assume h/d = 1.5 and 6 = 1.0, the maximum ship acceleration

(in g's) may be written

HiLtl)] = e 6 . (132)
g

Equation 132 implies e = 0.2 under the assumed conditions.

Using figure 20 and equations 131 and 132

P(t) I

jj2££ = 1 + (0.2X1.75) = 1.35 (133)

e f

Therefore, neglecting the second-order contribution could

result in a thirty-five percent error in l

F (t)l
max •

Figure 21 further demonstrates the nonlinear (second-

(1)
order) effects on F(t). The function F(t)/ef is plotted

85
 .

vs t. The ship's displacement as a function of time is also

shown on the plot for reference (cosine function of unit

amplitude). The parameters chosen for the plot are h/d = 1.5,

6 = 1.0 and e = 0.2 (corresponds to an equivalent acceleration

of 0.2 g's)

2.0

1.0

(i)
ef

-1.0

Ship's Displacement

-2.0
1.0 3.0 4.0
time

FIGURE 21. Total Exciting Force vs. Time

86
\
 .

Figure 22 illustrates, for the same values of the

parameters, the nonlinear effects on the gravity wave

generated by the ship's disturbance. The figure presents

the sum of the dimensionless first-order wave amplitude and

the second-order Stokes' wave (both have celerity c, )

Separately, the second-order wave from "<p '

(celerity c ) is
p

plotted. Both curves are plotted vs x/b for fixed time.

X2 t O»I.(

h/d=l.5

{2) {
77 (f rom ')
<f>

... (Sfokos second

^7 order correction)
-2

-4
2 4 6 8
x/b

FIGURE 22. Second-Order V/ave Form Effects

87
 . .

VI. CONCLUSIONS

It is concluded that this work has demonstrated the

feasibility of using the finite element method to solve the

linear and nonlinear problems of the type presented. Further,

the finite element method is shown to be quite versatile in

that various boundary conditions may be applied to the

boundary value problems studied, with only minor changes in

the computational scheme. In addition, the isoparametric

elements chosen are capable of representing arbitrary hull

shapes

It is further concluded that a first-order study of the

effect of finite depth on the ship motions studied demon-

strates that little change in the solutions occur. However,

it is asserted that the second-order results obtained

demonstrate that nonlinear (second-order) effects on heaving

motions are most significant in very shallow water; particu-

larly for motions at a frequency corresponding very closely

to the natural frequency of free oscillations of the floating

body

The present state of digital computer capacity appears to

preclude the extension of the present problem to three-

dimensions. However, the extension of this study to other

modes of motion (sway and roll) is possible as demonstrated

for the linear sway problem. Further, the possibility of

studying non-symmetric, two-dimensional hull forms is

suggested.

88
 An adequate finite element mesh to represent the region

for the nonlinear problems studied tends to tax the core

capacity of the computer used. Nevertheless, meaningful

results were obtained through careful design of the computa-

tional scheme.

89
 —

APPENDIX A

Complex Algebra

The development which follows determines the correct

form of the complex numbers w, and w ? in the following

expression

iat iat 2iat

Re{z
1
e }Re{z
2
e } = Re^ + w e
2
} (Al)

z-. and Zp are complex numbers with moduli r-, and r ? and

arguments 8-, and 0p. It follows using appropriate trigono-

metric identies that

1 iat 2 iat
Re{r e e }Re{r e e } =
1 2

r r ? [cos6 1 cose 2 ^(l+cos2at) + sine,sin6p ^(l-cos2at)

p -p

- sin(0 +6 )^sin2at] =-^ -2- [cos(e -0 ) + cos (9- +e +2ct ) ] (A2)

2 1 ? 1 2 L 2

Equation A2 implies that

Z Z Z Z
,.-.
w
1
= — 1p 2
, w2 = —1 —2 2 . (A3)

The symbol (~) implies conjugation.

Equation A3 shows the origin of the factor h in various

places in the text.

90
 APPENDIX B

Isoparametric Elements

H,l) (£.1) (i.l) (1,1)

H.t)

H,1f)<M

(-1,-1)^

£,r) plane x f y plane

FIGURE Bl . Element Mapping from £,n Plane to x,y Plane

The essential idea in the isoparametric element lies

in the following equations which define a map from the £,n

plane to the x,y plane.

e e
x = *NT x , (Bl)

e e
y = N y (B2)

The vector N e is a row vector of element shape functions

.
r
1 x 12

e e
one for each node in figure Bl. The vectors x and y are

91 \
12 x 1 column vectors representing the respective x and y

coordinates of the nodes in the x,y plane (figure Bl).

The shape functions are conveniently defined using the

conventions of Zienkiewicz [45]. Let

= =
«o «i > % ™± •
(B3)

where E, and n are the coordinates of any point of the square

element in the £,r) plane. £. and n . are the coordinates of

the i node. The shape functions are given by the following

set of equations.

Corner Nodes:

N
i
=
J2
(1 +
^o
)(1 + % )[ " 10 + 9( ^ + n2)]
'
(B2])

Edge Nodes

for q = ±1, n
±
= ±| ,

2
N® = -^ (1 + C )(l " n )(l + 9n ) ,
(B5)

for K = ±| , n± = ±1 ,
±

2
N| = ^ (l + n )(l - C )(l + 9? ) (B6)

92 \
 The element is defined to be isoparametric when the same

shape functions are used to achieve the element mapping as

are used to define the field function in the x,y plane.

There are certain consistency requirements on the element

shape functions as discussed by Zienkiewicz [45]. The shape

functions chosen satisfy them all.

The mapping has the advantage of being able to represent

a region with curvilinear sides. The area and line integrals

given in part III are calculated by the expressions given

below.

From the calculus, the following relation is known to hold

ft ~x (B7)
e
N r

J is the real 2x2 Jacobian matrix of the transformation

,,e e
N x
r hi
(B8)

r?e
N x
e
V
Me e

It follows, provided J is not singular, that

Vx" .-1 £ (B9)

e
N
_~y

93 «
\
Then, for example, equation 86 in part III may be written at

the element level as

,e
N
T ~x
T
e
H = / [N® N® ] dRe (BIO)
z ~x ~y
Re
~y

Using equations B8, B9, and

dRe = dxdy = |J| d£dn , (Bll)

equation B9 may be written

e
+1+1 T T
1 7 1
»t
H = / / [N? N® ] [J" ] [J" ] Jl d£dn . (B12)
~*> ~n z ~
~ -1 -1

Equation B12 demonstrates the method of the evaluation

e
of the matrix H or in general the matrix H. The integrations

are accomplished using six-point Gauss quadrature.

91 \
 APPENDIX C

Computer Program Economization

FIGURE CI. Schematic of Coefficient Matrix J iCH

The economizations in the computer solution of the linear

system of equations defined by equation 102 in part IV are

achieved by observing certain properties of the matrix K

First all interior nodes of a given mesh contribute only real

numbers to K . This is also true of nodes on S.,, Sp and S~

Only the nodes on Sj. contribute complex entries in K .

Further only nodes on S~ and Su contribute entries which are

frequency dependent.

The nodes on S„ and S^ typically comprise no more than

fifteen percent of K^ . Further, the nodes on S^ alone

95
comprise nor more than five percent of the total rows in K .

Therefore if K is rearranged with equations for nodes on

Sjj at the bottom, equations for S^ nodes directly above the

Sn equations, and all other rows of the matrix filling the

remaining part of K , the resulting matrix would appear as

schematically represented in figure CI. The region designated

with a Roman numeral I (uncrosshatched) indicates the portion

of the matrix completely real and frequency independent.

Region II indicates the portion of the matrix which has

elements that are frequency dependent. And, Region III

indicates the portion of the matrix which is complex and

frequency dependent. If the matrix K is formed as shown

in figure CI, eighty-five percent (approximately) of it may

be eliminated only once using real arithmetic. This is

2
possible because the matrix -a Q Q does not have to be added

to Region II until after Region I has been reduced by Gauss

elimination, and the same fact applies to the matrix ,4-CtJ

(ia/c )D in Region III. Then, after Region I is eliminated,

n

p
the matrix -a Q Q is added to Region II (after saving the

previous results) and K is reduced to upper triangular

form through Region II using real arithmetic. Now, if this

result is saved and then the matrix (ia/c-.)D is added to

Region III only a small set of complex equations (usually

about 25) is left to solve.

After the final elimination, the resulting matrix is in

complete upper-triangular form and ready for a back-

substitution process. The right-hand side vector is complex

96 (
 .

and frequency dependent in general. However, the required

information to perform the forward substitution process on

it is contained in the final form of K (upper triangular)

Therefore, the forward substitution on ibar may be done

after K is processed. Then a back-substitution of the

resulting system using mixed-mode computer arithmetic yields

the final solution vector <jr '

.

A close examination of the process just described reveals

that for a whole family of frequencies only a small portion

of K has to be reassembled and eliminated again. These

matrix properties contribute toward a sizeable saving in

computer time.

i
\
97
 BIBLIOGRAPHY

1. Adamek, J. R. An Automatic Mesh Generator Using Two

,

and Three-Dimensional Isoparametric Elements M. S.

,

Thesis, Naval Postgraduate School, Monterey, California,

1973.

2. Allouard, Y., and Coubert J. P., "Numerical Study of

,

Transitory or Non-Linear Waves," Proceedings of the

International Symposium on Finite Element Methods in
Flow Problems Swansea, January, 1974, PP. 283-288.
,

3. Bai, K. J., "A Variational Method in Potential Flows

with a Free Surface," University of California, Berkeley,
College of Engineering, Report No. NA 72-2, 1972.

4. Chenault, D. W. II, Motion of a Ship at the Free Surface

, ,

M. S. Thesis, Naval Postgraduate School, Monterey,

California, 1970.

5. Fontannet, P., "Theorie de la Generation de la Houle

Cylindrique par un Batteur Plan (2 e Order d' Approxima-
tion)," La Houille Blanche Vol. 16, 1961, pp. 179-197.
,

6. Frank, W., "Oscillations of Cylinders in or Below the

Free Surface," Naval Ship Research and Development
Center, Report No. 2375, October, 1967, pp. 1-48.

7. Garrison, C. J., The Consistent Second-Order Theory

of Wave/Structure Interaction paper to be presented
,

at the 14th International Conference on Coastal Engi-

neering, Copenhagen, 24-28 June, 1974; To be published
in the Proceedings of the above: September, 1974.

8. Garrison, C. J., "Hydrodynamics of Large Objects in the

Sea Part I - Hydrodynamic Analysis," Journal of Hydro-
nautics, Vol. 8, No. 2, January, 1974, pp. 5-12.

9. Grim, 0., "Berechnung der durch Schwingungen Eines

Schiffskorpers Erzeugten Hydrodynamischen Krafte,"
der Schiffbautech. Ges. Vol. 47, 1953, PP. 277-296.
,

10. Havelock, T., Philosophical Magazine, Vol. 8, No. 51,

October 1929, pp. 569-576.

11. Havelock, T., "Waves Due to a Floating Sphere Making

Periodic Heaving Oscillations," Proceedings Royal
Society of London, Vol. 231, 1955, pp. 1-7.

98
 .

12. Holand, I., "Finite Elements for the Computation of

Hydrodynamic Mass," Proceedings of the Symposium on
Finite Element Techniques Stuttgart, 1969, pp. 509-531.
,

13. John, F., "On the Motion of Floating Bodies,"

Communications on Pure and Applied Mathematics , Vol. 3,
1950, pp. 45-101.

14. Kim, C. H., "Hydrodynamic Forces and Moments for Heaving,

Swaying, and Rolling Cylinders on Water of Finite Depth,"
Journal of Ship Research Vol. 13, No. 2, 1969,
,

pp. 137-154.

15. Lamb, H., Hydrodynamics , Cambridge University Press,

6th edition, 1932.

16. Landweber, L., and Macagno, M. "Added Mass of a Three

,

Parameter Family of Two Dimensional Forms," Journal

of Ship Research Vol. 2
, No. 4, 1959, pp. 36-4 8
,

17. Lee, CM., "The Second-Order Theory of Cylinders

Oscillating in a Free Surface," University of California,
Berkeley, College of Engineering, Report No. NA-66-7,
September, 1966.

18. Lee, C. M., "The Second-Order Theory of Heaving Cylinders

in a Free Surface," Journal of Ship Research Vol. 12, ,

1968, pp. 313-327.

19. Lewis, F. M., "The Inertia of Water Surrounding a

Vibrating Ship," Transactions of the Society of Naval
Architects and Marine Engineers, Vol. 37, 1929,
pp. 1-20.

20. MacCamy, R. C, "A Source Solution for Short-Crested

Waves," La Houille Blanche No. 3, July-August, 1957,
,

PP. 373, 379-389.

21. Matsumoto, K., "Application of Finite Element Method

to Added Virtual Mass of Ship Hull Vibration," Journal
of the Society of Naval Architects of Japan Vol. 127. ,

1970, pp. 83-90.

22. Matsuura, Y., and Kawakami , M., "Calculation of Added

Virtual Mass Moment of Inertia of Ship Hull Vibration
by the Finite Element Method," Journal of the Society
of Naval Architects of Japan Vol. 124, 1968, pp. 281-
,

291.

23. Newton, R. E., Chenault D. W., II, and Smith, D. A.,

,

Jr., "Finite Element Solution for Added Mass and Damping,"

Proceedings of the International Symposium on Finite
Element Methods in Flow Problems Swansea, January,,

1974, pp. 159-170.

99
24. Newton, R. E., Finite Element Analysis of Two-Dirnensional
Added Mass and Damping, to be published in Collection
of Papers presented at the International Symposium
on Finite Element Methods in Flow Problems, Swansea,
J. Wiley and Sons, London.

25. Ogilvie, T. F., "First and Second-Order Forces on a

Cylinder Submerged Under a Free Surface," Journal of
Fluid Mechanics Vol. 16, 1963, pp. ^51-^72.
,

26. Ogilvie, T. P., and Tuck, E. 0., "A Rational Strip

Theory of Ship Motions," University of Michigan,
Department of Naval Architecture, Report No. 013, 1969.

27. Parissis, G., "Second-Order Potentials and Forces for

Oscillating Cylinders on a Free Surface," Massachusetts
Institute of Technology, Report No. 66-10, 1966.

28. Paulling, J. R., and Richardson, R., "Measurement of

Pressures, Forces, and Radiating Waves for Cylinders
Oscillating in a Free Surface," University of California,
Berkeley, Report No. 82-28, 1962.

29. Paulling, J. R., and Porter, W. R., "Analysis and

Measurement of Pressure and Force on Heaving Cylinders
in a Free Surface," Proceedings of the Fourth National
Congress of Applied Mechanics 1962, pp. 1369-1382.
,

30. Porter, W. R. , "Pressure Distributions, Added-Mass,

and Damping Coefficients for Cylinders Oscillating in a
Free Surface," University of California, Berkeley, •

Report No. 82-16, i960.

31. Potash, R. L., "Second-Order Theory of Oscillating

Cylinders," University of California, Berkeley, Report
No. NA70-3, 1970.

32. Roren, E. M. Q., "Impact of Finite Element Techniques

on Practical Design of Ships Structures," Proceedings
of the Symposium on Finite Element Techniques ,Stuttgardt,
1969, PP. 509-531.

33. Salveson, N., On Second-Order Wave Theory for Submerged

Two-Dimensional Bodies Ph. D. Dissertation, University
,

of Michigan, Ann Arbor, 1966.

3^. Saunders, H. E., Hydrodynamics in Ship Design Vol. 2,

,

The Society of Naval Architects and Marine Engineers,

New York, 1957.

35. Stoker, J. J., Water Waves Interscience Publishers,

,

Inc., New York, 1957, pp. 19-5*1.

(
\
100
 . . , ,

36. Tasai, P., "Damping Forces and Added Mass of Ships

Heaving and Pitching," Kyushu University, Research
Institute of Applied Mechanics, Report, Vol. 8,
pp. 131-152, 1959; Vol. 8, pp. 36-39, I960.

37. Tuck, E. 0., "The Effect of Non-Linearity at the Free

Surface on Flow Past a Submerged Cylinder," Journal
of Fluid Mechanics Vol. 22, 1965, pp. 401-414
,

38. Ursell, F. , "On the Heaving Motion of a Circular Cylinder

on the Surface of a Fluid," Quarterly Journal of Applied
Mechanics and Applied Mathematics Vol. 2, No. 2, 1949.
,

pp. 218-231.

39. Visser, W., and van der Wilt, M., "A Numerical Approach
to the Study of Irregular Ship Motions," Proceedings
of the International Symposium on Finite Element
Methods in Flow Problems Swansea, January, 1974
,

pp. 277-282.

40. Vugts, J., "The Hydrodynamic Coefficients for Swaying,

Heaving, and Rolling Cylinders in a Free Surface,"
Netherlands Ship Research Center, Report No. 1125,
1968.

41. Wehausen ,J.V. , and Laitone E. V. ,Surface Waves

,

Encyclopedia of Physics Vol. 9, edited by S. Flugge

,

Springer-Verlag, Berlin, i960.

42. Yu, Y.S., and Ursell, F., "Surface Waves Generated

by an Oscillating Cylinder on Water of Finite Depth:
Theory and Experiment," Journal of Fluid Mechanics ,
Vol. 11, No. 4, 1961, pp. 529-551. *

43. Zienkiewicz, 0. C, and Newton, R. E., "Coupled Vibra-

tions of a Structure Submerged in a Compressible Fluid,"
Proceedings of the Symposium on Finite Element Techniques ,

Stuttgart, 1969, pp. 359-379.

44. Zienkiewicz, 0. C, Irons, B. M., and Nath, B., "Natural

Frequencies of Complex Free or Submerged Structures by
the Finite Element Method," Proceedings of the Conference
on Vibrations in Civil Engineering Institution of
,

Civil Engineers, London, 1965

45. Zienkiewicz, 0. C, The Finite Element Method in

Engineering Science McGraw-Hill, London, 1971.
,

46. Zienkiewicz, 0. C, and Cheung, Y. K., "Finite Elements

in the Solution of Field Problems," The Engineer ,
September, 1965, PP 507-510.
.

101 i
\
 INITIAL DISTRIBUTION LIST

No. Copies

1. Defense Documentation Center 2

Cameron Station
Alexandria, Virginia 22314

2. Library, Code 0212 2

Naval Postgraduate School
Monterey, California 939^0

3. Professor R. E. Newton, Code 59Ne 2

Department of Mechanical Engineering
Naval Postgraduate School
Monterey, California 939^0

4. Commander, Naval Shipyard. 1

Att: LT David A. Smith, Jr., Code 338.22
Long Beach, California 91701

5. Professor J. R. Paulling 1
Department of Naval Architecture
University of California
Berkeley, California 9^720

6. Naval Ship Research and Development Center 1

Engineering Research Project in Hydromechanics
Bethesda, Maryland 2003^

7. Naval Ship Systems Command Headquarters 1

Code 03
Navy Department
Washington, D. C. 20360

102
 Thesis
S5759 Smith
Finite element anal
1

c.l os
ysis of the forced
ci nations of
ship
hull forms.

1 C{
JU W U \J ^3
*J

Thesis
S5759 Smith
c.l Finite element anal-
ysis of the forced os-
ci llations of ship
hull forms.
thesS5759
F e 6lement analysis
of th e forced OS
',n,n

3 2768 002 00898 9

DUDLEY KNOX LIBRARY