Below are the first 10 and last 10 pages of uncorrected machineread text (when available) of this chapter, followed by the top 30 algorithmically extracted key phrases from the chapter as a whole.
Intended to provide our own search engines and external engines with highly rich, chapterrepresentative searchable text on the opening pages of each chapter.
Because it is UNCORRECTED material, please consider the following text as a useful but insufficient proxy for the authoritative book pages.
Do not use for reproduction, copying, pasting, or reading; exclusively for search engines.
OCR for page 427
A Numerical Solution Method for
Three Dimensional Nonlinear Free Surface Problems
C.G. Kang, I.Y. Gong (Ship Research Station, KIMM, Korea)
ABSTRACT
The nonlinear hydrodynamics of a threedimen
sional body beneath the free surface is solved in the
time domain by a semiLagrangian method. The
boundary value problem is solved by using the bound
ary integral method. The geometries of the body and
the free surface are represented by the curved panels.
The surfaces are discretized into the small surface
elements using a bicubic Bspline algorithm. The
boundary values of ~ and ~ are assumed to be bi
linear on the subdivided surface. The singular part
proportional to R are subtracted off and are inte
grated analytically in the calculation of the induced
potential by singularities.
The far field flow away from the body is rep
resented by a dipole at the origin of the coordinate
system. The RungeKutta 4th order algorithm is
employed in the time stepping scheme. The three
dimensional form of the integral equation and the
boundary conditions for the time derivative of the
potential is derived. By using these formulas, the
free surface shape and the equations of motion are
calculated simultaneously. The free surface shape
and the forces acting on a body oscillating sinu
soidally with large amplitude are calculated and com
pared with published results. Nonlinear effects on a
body near the free surface are investigated.
1. INTRODUCTION
The free surface effects are considered in the
design of submerged bodies operating near free sur
face. The linearized theories were developed during
many decades. Recently the nonlinear free surface
problem is solved in the time domain by the semi
Lagrangian method
LonguetHiggins and Cokelet (1976) presented
a mixed EulerianLagrangian method for following
the timehistory of spaceperiodic irrotational sur
C.G. Kang and I.Y. Gong, Korea Research Institute of Ships and Ocean Engineering
171 Jangdong Yuseonggu Daejeon, Korea
427
face waves. The only independent variables at the
beginning of each time step were the coordinates
and velocity potential of marked particles on the
free surface. At each timestep an integral equation
was solved for the new normal component of veloc
ity. This method was applied to a free, steady wave
of finite amplitude, and was found to give excellent
agreement with calculations based on Stokes's series.
It was then extended to unsteady waves, produced
by initially applying an asymmetric distribution of
pressure to a symmetric, progressive wave. The re
sults showed the freely running wave then steepened
and overturned.
Using a technique similar to that of Longuet
Higgins and Cokelet (1976), Faltinsen (1977) solved
a nonlinear two dimensional free surface problem in
cluding a harmonically oscillating body. The body
intersected the free surface and was constrained to
move in the vertical direction. The numerical cal
culations were reduced by representing the flow far
away from the body as a dipole located at the center
of the body. A formula to calculate the exact force
on the body was presented. It was only necessary to
know the velocity potential on the positions of the
free surface and the wetted body surface.
A numerical method for the time simulation of
the nonlinear motions of two dimensional surface
piercing bodies of arbitrary shapes in water of finite
depth was presented by Vinje & Brevig (1981~. Pe
riodicity in space was assumed. At each time step,
Cauchy's integral theorem was applied to calculate
the complex potential and its time derivative along
the boundary. The solution was stepped forward in
time by integrating the exact kinematic and dynamic
freesurface boundary conditions as well as the equa
tion of motion for the body. They solved the problem
of capsizing in beam seas, caused by extreme waves.
Twodimensional nonlinear free surface prob
lems by a dipole (vortex and source) distribution
method were solved by Baker, et al. (1982) . The
resulting Predholm integral equation of the second
kind was solved by iteration which reduced storage
OCR for page 427
and computing time. Applications to breaking water
waves over finitebottom topography and interacting
triads of surface and interracial waves were given.
The semiLagrangian method was extended to
vertically a~cisym~netric free surface flows by Dom
mermuth & Yue (1986) . Since they solved the fi
nite depth problem, a far field closure was imple
mented by matching the linearized solution outside
a radiation boundary. The intersection line between
the body and free surface was treated by extending
Lin's(1984) method.
The nonlinear hydrodynamics of an axisym
metric body beneath the free surface in the time do
main were solved by Kang & Troesch (1988~. The
free surface shape and the forces acting on a sphere
oscillating sinusoidally with large amplitude are cal
culated and compared with published results. The
far field flow away from the body is represented by a
dipole at the origin of the coordinate system similar
to Faltinsen (1977~. This is only valid until waves
arrive. Waves generated by the numerical error at
the truncation boundary are not observed.
In this paper, the method used for axisymmet
ric flows by Kang and Troesh(1988) is extended to
threedimensional free surface flows. The free surface
shape and the forces acting on a threedimensional
body moving forward and oscillating sinusoidally with
large amplitude are calculated and compared with
published results. When the body motion is un
known, the time derivative of the potential on the
body is needed for the time simulation. In two di
mensions, Vinje & Brevig (1982) derived the inte
gral equation and the boundary conditions for the
time derivative of the potential and stream func
tion. However their formulas may not be extended to
threedimensional problems. The threedimensional
form is derived in this work. By using these for
mulas, the free surface shape and the equations of
motion are calculated simultaneously. The Runge
Kutta 4th order algorithm is employed in the time
stepping scheme (See Appendix 1~.
Numerical calculations are performed for the
following cases:
(a) A body oscillating vertically near the free sur
face
origin located at the mean free surface. The fluid is
assumed to be inviscid and incompressible and the
Dow is assumed to be irrotational. The fluid domain
is bounded with the following surfaces, the free sur
face, SF, the body, SB, and the surfaces at infinity,
SOO (Fig. 1~. The surfaces, taken as a whole, will
be denoted as S. The governing equation and the
boundary conditions are as follows(LonguetHiggins
& Cokelet (1976) and Dommermuth & Yue (1986~:
Laplace equation:
V2¢ = 0 in the fluid domain (1)
Kinematic free surface condition:
D1; = VI on F(x,t) = 0 (2)
Dynamic free surface condition
D} =  9Z + 2 V'· Vie on F(x, t) = 0 (3)
Body boundary condition:
V¢.n~x,t) = V.n on B(x,t) = 0 (4)
Radiation condition:
~ ~ O as Axe ~ oo,t < oo (5)
where Din= a&t + V¢. V) iS the substantial deriva
tive, F(x, t) = 0 is the function representing the free
surface geometry at time t, V includes both trans
lational and rotational velocities, and B(x,t) = 0 is
the function representing the body surface geometry
at time t.
The Green function, G~xjy), satisfies the fol
lowing equation.
(b) A body oscillating horizontally near the free V2GtX;Y) =fixY) (6)
2. MATHEMATICAL FORMULATION
Consider an ideal fluid below the surface given
by F(x,t) = 0, where x(z,y,z) is a righthanded
coordinate system with z positive upwards and the
where x is the vector to the field point, y is the vector
to the source point, and fixy) is the Dirac delta
function. Through the application of Green's second
identity in the fluid domain, the potential is given as
a~x,t)~(x,t) = i/t ~¢¢~ Aids (7)
where a is an included solid angle at x. In this prob
lem, a is 27r on the surface.
The Green function that satisfies Eq. (6) is
428
OCR for page 427
G(x,y) = R = ~(8)
where x is the position vector of a field point and y
is that of a source point.
The solution of Eq. (3) gives the potential ~ on
the free surface F(x,t) = 0. Also in on the body is
known from the body boundary condition, Eq. (4).
The normal and tangential velocities on the free sur
face are needed to solve Eq. (3). The normal velocity
on the free surface is a solution of Eq. (7). Details to
calculate the tangential velocity is given in the sec
tion 3. Consequently, a Fredholm integral equation
of the second kind on the body and of the first kind
at the free surface may be solved.
Far Field Condition
The far field condition is important in the non
linear free surface problem. It can be resolved by
using periodic boundary conditions if the physical
problem has spatial periodicity (LonguetHiggins &
Cokelet, 1976~. Faltinsen (1977) and Kang (1988)
assumed that the behavior of the potential is like
that of a dipole at the origin of the coordinate sys
tem.
A far field closure by matching the nonlinear
computational solution to a general linear solution
of transient outgoing radiated waves was used by
Dommermuth & Yue (1986~. This method is math
ematically complete.
A numerical radiation condition was posed so
as no waves reflected from the truncated surface
(Yang & Liu (198933. They found that the usual one
dimensional Sommerfeld condition gave reasonable
results for an axisyrnmetric cylinder heaving in the
still water. Also it was extended to 2D case for the
cylinder swaying in the still water.
In this work, the far field closure similar to that
used by Faltinsen (1977) and Kang (1988) is posed.
It is simple and it works well until waves arrive at
the truncation boundary.
At the far field, the velocity potential, ¢, and
the wave elevation, a, are small from the radiation
condition, Eq. (5~. For example,
¢(z = q) = +(x = 0) + ~ .~¢ (z = 0) + H.O.T. (9)
Assuming the behavior of the potential, ¢, is
like that of a dipole at the origin of the coordinate
system, it follows that as r · oo
¢(z = 0) = 0
a~ (z = 0) ~ 3
a=/ ozfz=0)dt~~ 3dt (10)
¢(Z = Q) ~ FIZZ = 0) ~ 6
where r is `~. This is only valid until waves
arrive at the truncation boundary.
If we take a large value of r, the potential
on the free surface must be relatively small to the
vertical velocity ~ . Therefore the effect of the po
tential at the far field can be neglected. The far field
condition is approximately satisfied by including the
effect of the vertical velocity ~ at the far field.
3. NUMERICAL IMPLEMENTATION
The body surface and the free surface are dis
cretized into the small surface elements AS`j using
a bicubic Bspline algorithm ( Barsky & Greenberg
(1980) ) . The surfaces I\Sij(x,y,z) can be repre
sented by the parameters, u and v. Thus
1 1
ij(U, V) = ~ ~ bJ(U)Vi+J,j+tat(V)
a=  2 t=  2
1 1
ydi(uiv)= ~ ~ bJ(U)Vi+s,j+tbt(V) (11)
s=2 t=2
1 1
z,}(U,v)= ~ ~ bJ(U)ViX+s,j+tbt(V)
a=  2 t=2
for 0 ~ u < 1 and 0 < v < 1
where bJ(u) and lot(u) are the uniform cubic Bspline
basis functions and V`j are vertices (See Appendix
2). This allows the curved panels.
The end condition should be imposed to get
a complete Bspline approximation. There are sev
eral methods to impose end conditions according to
the geometrical characteristics (Barsky(1982)). The
derivative of Bspline interpolation at the end is set
to get the tangent of the given geometry if the tan
gent is known. If the tangent is not known, the
derivative at the end is set to be the slope between
two vertices at the end obtained by using Bspline
algorithm.
The boundary values of ~ and ~ are assumed
to be bilinear on the subdivided surface AS`j as
shown below .
= aO + alu + a2v + a3uv
~3¢ = ho + be + b2v + b3uv (12)
for 0 ~ u < 1 and 0 < v < 1
To evaluate the integrals over the segments the
two point Gaussian Quadrature formula (Perxiger
429
OCR for page 427
(1981), Abrnmowitz ~ Stegun (19643) is used when
the field point is not a corner of the pannel. In Eq.
(7), Gn is not singular but G has (R) type singu
larity in the transformed uv domain as the field and
point approaches the source point. The singular
ity is integrable and can be integrated by numeri
cal quadrature. But since an accurate integration
of the singularity requires a higher order quadrature
formula, the method following Ferziger (1981) and
Dommermuth & Yue (1986) can be used. The in
tegral can be factored into the sum of the (R~ type
singular part which is integrable analytically and the where
nonsingular part which requires numerical quadra
ture (Ferziger(1981~.
Removal of (R~ type singularity
In Eq. (7), G has (R~ type singularity as the
field point approaches the source point. The (R)
type singularity is integrable in the surface integra
tion.
/l^sij R (13)
First, consider the induced potential at (/Oo, 900, hoo)
,which is one of the corners of panel, by the source
panel AS`j . Eq. (11) can alternatively be repre
sented by the following equations:
~ 3
x' = ~ ~ fijuiv]
i=0 j=0
~ 5
Y' = ~ ~ 9iiuivi (14)
`=o j=o
z' = ~ ~ hijUiV]
`=o j=o
By using Eqs.~14) and (40), dS and R can be
transformed and expanded into Taylor's series about
u = 0, v = 0 as follows :
dS = IJ~dudv
jEGF2dudv
J ~ (~ ) + (39 ) + (68 ) ]
'~(~) +(BY) +(aZ) 1
where
, , ~
{3X' ~X' + ~y ~y + ~X ~Z ~ 2 ~ ~ d
v respectively. Generally, the tangential vector to
(fU,gu~hu) along uaxis is not perpendicular to the
vector to (go) along vaxis. Therefore Gram
Schmidt orthogonalization is needed to get the or
thogonal tangential vectors on the surface.
4. CALCULATION OF THE TIME DERIVA
TIVE OF POTENTLY
For the time simulation, ~ should be known to
calculate the forces and moments acting on the body.
In two dimensions Vinje & Brevig(1982) derived an
integral equation and boundary condition for ~ by
using the ~ and ~ formulation. However their results
can not be extended to the threedimensional case.
Since `~&n(~) can not be calculated by using the given
motion, a boundary value problem for A, the time
derivative of the potential in body fixed coordinates,
is derived as follows:
dt On
= n dtV' + V' dt
= n ldtV'+ (y. V)V¢] + V) ~ (w x n)
= n ~ [V dt + V(V V¢) + w x V+] + V} (w x n)
an dt
which can be expressed as
.g ~ dt ~ = n ~ ~ mitt + w x rw x VT) (23)
Following the nomenclature of Vinje & Brevig (1982),
the operator `~` is ~ 98' + V V), V = V T + w x _, VT is
the translational velocity of the center of mass of the
body, _ is the position vector to the boundary from
the center of mass of the body, and w is the angular
velocity vector of the body. Eq. (23) is useful in that
most quantities of interest are expressed in the body
coordinate system rather than a fixed, inertial one.
Since V Vie satisfies Laplace equation, the time
derivative of the potential, A, can be calculated by
using Green's theorem. The limiting behavior of V.
V' at red oo can also be checked, or
V V' = 0(r¢) = 0((~)) as r ~ oo (24)
Applying Green's theorem for ~ instead of
in Eqs. (1), (6), and (7), the following equation can
be obtained:
R.H.S. of Eq. (23) may be represented as follows:
n . ( ,9fT + W X _~ X VT) = ~ anti (26)
where n = (n~,n2,n~,) and _ x n = (n`,n5,n6) .
The time derivative of the potential, A, can be
decomposed as follows:
die = Thai deli + d¢7 (27)
The auxiliary terms, A; and "7, are solutions
of the following boundary value problems:
deli) = ni and
and
6¢i = 0 and
dt
~n( dt ) = 0
on B(x,t) = 0 (28)
d¢7 = V · Vie2 V' · V'gz
on F(x,t) = 0 (29)
The time derivative of the potential on the free
surface,"7, is calculated by using solutions of the
integral equation, Eq. (7).
5. THE PRESSURES AND FORCES
Once the time derivative of the potential is
known, the pressures are found by applying Bernoulli's
equation. Bernoulli's equation is derived for the vari
ables relative to an inertial coordinate system. How
ever, it is convenient for the purpose of solving the
boundary value problem to use body fixed coordi
nates. Under these circumstances, spatial differen
tiation is invariant with coordinate transformation,
but temporal differentiation is not. Bernoulli's equa
tion can be expressed as
 =  a'  1 V' V'  go
dt + V V' 2V¢) V.d)go. (30)
The term ~ in the above equation is calculated
by Eq. (27). With the pressure known, the force and
moment become
F = mV
_ _
at lli'~n(~t)  (at) ]GdS (25) = llpodSmilk
431
OCR for page 427
pun(d)V. V++ 2V¢e V++gZ)dS
s
mgk
M = ip_ x ndS.
s
For threedimensional bodies the force and mo
ment are rearranged as follows:
F = F. + F2 + (pgVmg~k (32)
M = M~+M2
where V in Eq. (32) is the displaced volume of the
sphere,
F. = p  n d,¢~;7dS
SB
F2 = p AL n(v . Ad2 Vet . V¢) dS, (33)
SB
Ml = p If r x n dads, and
SB
M2 = p  _ x n(V Vet2V' · V+) d S.
SB
6. NUMERICAL CALCULATION
IIeave motion
To demonstrate the usefulness of the technique
shown in the previous section, the force acting on
a sphere oscillating beneath the free surface is de
termined. The motion of a sphere is given by z =
h + a cos at for t greater than zero. Initially the
potential and wave elevation at the free surface are
zero.
The number of elements on the body is 200 and
that on the free surface is 40 x 40. The truncation
boundary is the position from the origin of the cm
ordinate system where waves reaches in four periods
of the body motion ~ 16 < x/R < 16,  16 <
y/R ~ 16~. So, it depends on the group velocity of
wave. Even spacing is used on the body and free
surface. The typical time interval is approximately
0.05 period of motion for the time simulation of the
sphere.
The mean depth of immersion for the center of
the sphere, h, is h/R=2.0. The time history of the
force acting on the oscillating sphere with a large ra
tio of motion amplitude, a, to radius, R. (a/R=0.5)
was calculated and is compared with the results of
the axisymmetric free surface problem (Kang &
Troesch (19883) in Fig. 2. They show good agree
ment. In the reference~l3], the calculation results
for the sphere oscillating vertically with small am
plitude showed good agreement with those given by
Ferrant(19873. This means the AD algorithm in this
(31) paper works well.
In Fig. 3, the time history of force components
which consist of F. and F2 is shown. The forces are
nondimensionalized by pgKaR3, where K is a wave
number, w2/g, and 9 is the gravitational constant
and the time t is nondimensionalized by ,/i~. The
harmonic distributions of the total force are shown
in Table 1. The second order amplitude of the force
is 6.5% of the first order one. And the mean force
is 1.5% of the first order. Fig. 4 shows the three di
mensional wave profiles at four different times. All
the wave profiles are exaggerated by factor 50 in z
direction. In the figures T is a period of the motion.
Surge motion
The surge motion of the sphere is given by x =
a cos wt for t greater than zero. The mean depth
of immersion for the center of the sphere is h/R =
2.0. The amplitudes of the motion is a/R = 0.5. The
nondimensionalized wave number, KR, is equal to
1.0.
The time histories of the forces acting on the
sphere are shown in Fig. 56. The harmonic distri
butions of the horizontal and vertical forces are given
in Table 2. The three dimensional wave profiles at 4
different times are shown in Fig. 7.
In case of surge motion, the first order surge
force is dominant unlike the heave motion. However
nonlinear effects appear only in the vertical force.
The first order vertical force is negligible, but the
mean and second order vertical forces are not. The
mean vertical force is important for a submerged
body to keep a constant depth.
Advancing motion
Sawtooth instability is not observed in the com
putation of the oscillatory motions. But it seems to
be inevitable and break down the numerical time
stepping in case of advancing sphere. It may be due
to short waves generated by the body. The length of
short waves is less than the mesh size in this compu
tation. The numerical error does not die out but was
accumulated continuously. Eventually the numerical
scheme breaks down. Thus a simple numerical filter
ing scheme are tried to avoid break down, but still
does not work well. Fig. 8 shows the wave profile
before breakdown.
All the calculations were carried out on CRAY2S
and each solution time was approximately 50000 sec
onds.
432
OCR for page 427
7. CONCLUSION
In this paper the nonlinear hydrodynamics of
a threedimensional body beneath the free surface is
solved in the time domain. The free surface shape
and forces acting on a sphere oscillating sinusoidally
with large amplitude are calculated and compared
with published results. The far field flow away from
the body is represented by a three dimensional dipole
at the origin of the coordinate system. This is only
valid until waves arrive at the truncation bound
ary. Any numerical instability is not observed in the
computation of the oscillatory motions. The inter
gral equation and boundary conditions to calculate
the time derivative of the potential on the body are
derived. By using these formulas, the free surface
shape and forces are calculated simultaneously. A
Rung  Kutta 4th order scheme is employed in the
solution method. Nonlinear effects on the oscillating
body submerged in infinite water depth are studied.
ACKNOWLEDGMENT
This work was supported by the Basic Research
Program, contract ED469 of the Korea Research In
stitute of Ships and Ocean Engineering (KRISO).
Acknowledgement is also given to the Ministry of
Science and Technology of Korea (MOST). The au
thors should appreciate Dr. Choung Mook Lee for
his encouragements.
REFERENCES
[1] Abramowitz, M. & Stegun, I.A., Handbook of
Mathematical Functions, Government Printing
Office, Washington, 1964.
[2) Baker, G.R., Meiron, D.I. & Orsag, S.A., "Gen
eralized Vortex Methods for FreeSurface Flow Pro
blems," Journal of Fluid Mechanics, 123, pp477
501, 1982 .
[3] Barsky, B.A. & Greenberg, D.P., Determining
a Set of BSpline Control Vertices to Generate an
Interpolating Surface," Computer Graphics and
Image Process 14, pp203226, 1980.
.
[4~ Barsky, B.A., End Conditions and Boundary
Conditions for Uniform BSpline Curve and Sur
face Representations," Computers in Industry 3,
pp. 1729, 1982.
[51 Dommermuth, D.G. & Yue, D.K.P. 1986, "Study
of Nonlinear Axisymmetric BodyWave Interac
tions," In Proc. 16th Symp.on Naval Hydrodyna
mics, Berkeley.
[6] Dommermuth, D.G., Numerical Methods for
Solving Nonlinear WaterWave Problems in the
Time Domain," Ph.D. Thesis, MIT, 1987.
[7) Faltinsen, O.M., "Numerical Solution of Tran
sient Nonlinear FreeSurface Motion Outside or
Inside Moving Bodies," Proc. 2nd Intl. Conf.
on Num. Ship Hydra., U.C. Berkeley, 1977.
[8) Ferrant, P., "Sphere immergee en mouvement de
pilonnement de grande amplitude," Premiers
Journes De L'hydrodynamique, Nantes, 1987.
[9I Ferziger, J.H., Numerical Methods for Engineer
ing Application, John Wiley and Sons,Inc., 1981.
[10] Forbes, L.K., "An Algorithm for 3Dimensional
FreeSurface Problems in Hydrodynamics, "
Journal of Computational Physics 82, pp.33() 347,
1989.
[11] Gradsteyn, I.S. and Ryzhik, I.M., Table of Inter
grails, Series, and Products, Academic Press, 1980.
[12) Kang, C.G., Bow Flare Slarruning and Non
linear Free Surfac  Body Interaction in the Time
Domain," Ph.D. Thesis, Univ. of Michigan, 1988.
[13] Kang, C.G.& Troesch, A.W., "Nonlinear In
teration betwwen Axisyrnmetric Bodies and The
EYee SurfaceBody in Water of Infinite Depth,"
Proc. the Seminar on Ship Hydrodynamics to

honor Prof. J.H.Hwang, Seoul, Korea, 1988.
.
[14] Kaplan,W., Advanced Mathematics for Engine
e , AddisonWesley Publishing Co.
[15] Lin, W.M., "Nonlinear Motion of the Free Sur
face near a Moving Body," Ph.D. Thesis, M.I.T.,
Dept. of Ocean Engineering, 1984.
[16] Lin, W.M., Newman, J.N. & Yue, D.K.P., "Non
linear Forced Motions of Floating Bodies," Proc.
15th Symp. on Naval Hydra., Hamburg, 1984.
_ . .
[17] LonguetHiggins, M.S. & Cokelet, E.D., UThe
Deformation of Steep Surface Waves on Water, I.
A Numerical Method of Computation,"
Proc. R. Sac. Land. A., 350, ppl26, 1976
[18] Newman, J.N., Transient Axisymmetric Mo
tion of a Floating Cylinder," J. Fluid Mech., 157,
ppl733, 1985.
[19] Vinje, T. & Brevig, P., Nonlinear Ship Mo
tions," Proc. 3rd Intl. Conf. on Numerical
Ship Hydrodynamics, Paris, 1981.
[20] Vinje, T. & Brevig, P., ~Nonlinear, Two Di
mensional Ship Motions," Seminar on the Norwei
gan Ships in Rough Seas ISIS) Project, 1982.
21] Yang,C. & Liu,Y.Z., ~ Tim  Domain Calcu
lation of the Nonlinear Hydrodynamics of Wavy
Body Interaction," 5 th Intl. Conf. on Num. Ship
433
OCR for page 427
Hydro., Hiroshima, 1989.
Table 1 Harmonic Distributions of the Total
Force for Heave Motion
(a/R=0.5,h/R=2.0,KR=l.o)
Heave Force F/pgKaR~
ITH COS SIN
o
1
2
3
0. 2850707EO1 0. OOOOOOOE+OO
0. 1843049E+01 0. 2665849E+OO
0. 1033978E+00 0.6141333EO1
0. 6449880E01 0. 8355246E02
0. 511481 lE02 0. 1019468EO1
0. 2418429E01 0. 2957444E02
Table 2 Harmonic Distribution of the Total
Force for Surge Motion
(a/R=O.5,h/R=2.0,KR= 1.O)
Surge Force F/pgKaR3
ITH COS SIN
o
2
3
4
5
0. 1308621E02 0. OOOOOOOE+OO
0.1927655E+01 0. 1237885E+OO
0.2138726E03 0. 1002954E02
0. 3005543E01 0.1583521E02
0. 1381397E04 0. 4408678E03
0. 2590462E01 0. 4266882E03
Heave Force F/pgKaR3
ITH COS SIN
0. 1274143E01 0.0000000E+OO
0. 4001656E03 0. 4211906E04
0. 2037149E01 0. 2325045EO1
0. 3133378E03 0. 8355940E04
0. 3846344E03 0. 2575103E03
0. 2204935E03 0. 9557988E04
t
SF
y
( )
Fig. 1 Coordinate System
1~ '
1
Fig. 2
Comparison of Heave Force Acting on
the Sphere by 3D and Axisymmetric
Solutions (a/R=0.5, h/R=2.0, KR=1.O)
Fig. 3 Time History of the Heave Force Compo
nents Acting on the Heaving Sphere
(a/R=0.5, h/R=2.0, KR=1.O)
434
OCR for page 427
~=1.99T
t=1 . 99T
Fig. 4 Wave Profiles for Heave Motion
(a/R=0.5, h/R=2.0, KR=1.0)
t=3 2 4T
Fig. 7 Wave Profiles for Surge Motion
(a/R=0.5, h/R=2.0, KR=1.0)
435
OCR for page 427
~rat  r2
Appendix 1. The 4th order RungeKutta
method [9]
When y' = fishy) is nonlinear, this can be
solved by RungeKutta methods. The most com
monly used RungeKutta methods are fourth order
accurate and there are a number of these. The best
known such method (sometimes called the fourth or
der RungeKutta method) is
Fig. 5 Time History of the Surge Force Compm Yn+l/2 = Yn + 2J(Zn,Yn)
Dents Acting on the Surging Sphere (Euler predictor  half step)
(a/R=0.5, h/R=2.0, KR=1.0) * h
Yn;~/2 = Yn + 2 /(Xn+~/2, Yn+~/2)
ZForce
Fig. 6 Time History of the Heave Force Compo
nents Acting on the Surging Sphere
(a/R=0.5, h/R=2.0, KR=1.0)
Fig. 8 Wave Profile Generated by an Advancing
Sphere (U = 2.557m/sec for 2.557 < t,
U = t for O < t < 2.557,
h/R = 2.0, t = 3.6sec)
(Backward Euler corrector  half step)
Yn;1 = Yn + h f (fin+ 1/2 ~ An+ 1/2 ~ (34)
(Midpoint rule predictor  full step)
Yn+1 = Yn + 6 [fern, Yn) + 2f~xn+l/2,yn+ll2)
+2f~=n+ll2,yn+ll2) + f~xn+l~yn+l)]
(Simpson's rule corrector full step)
Looking at this method one can see that derivation
of such methods is not an easy task. An analysis of it
for the general case is also difficult. It is not too dif
ficult to analyze, however, when applied to y' = ay.
We find that
Yn+1 = (1 + ah + 2 + 6 24 )
(35)
so that the method is indeed fourth order accurate
and the error is of order (ah)5/120 . It is interesting
to note that the steps that comprise this method
are of order one,one,two, and four, respectively, and
the method has inherited the accuracy of the final
corrector.
Appendix 2. Parametric Uniform Bspline
Surface Representation [3]
436
OCR for page 427
A Bspline surface is defined in a piecewise man
ner, where each piece is a segment of the surface
called a surface patch. The entire surface is a mo
saic of these patches sewn together with appropriate
continuity (Fig. 93. A bicubic Bspline surface con
sists of patches which are cubic in each of the two
parametric directions and it is everywhere continu
ous along with its first and second derivative vec
tors, in both directions. This continuity constraint
reduces to requiring continuity of the first and sec
ond parametric derivative vectors across the borders
of adjacent patches. The Bspline surface is defined
by, but does not interpolate, a set of points called
control vertices. These control vertices form a two
dimensional array. Although the vertices actually
exist in threedimensional xyz space, they are or
ganized as a twodimensional graph. Each vertex is
either an interior vertex or a boundary vertex. This
notion can be formalized quite elegantly by drawing
on graph theory. The set of control vertices can be
considered as a graph V, E whose vertices form the
set
V={V`j ~ i=O,.~.,m;j=O,.~.,n)
and with the set of edges
E = {(V`;, V`,j+~) I = 0, ·, m1; j = 0, · . ., n3
U{(V`;, Vi+~,j) Pi = 0, · , m1; j = 0, · . . , nil
The interior vertices are the vertices Vij, where
1 < i < m1 and 1 < j < n1 And the boundary
vertices are Voj, j = 0, ...,n1, Vi",i = 0, ...,m1,
Vmj, j = 1, ..., n, and V`O, i = 1, ..., m . To emphasize
this graphtheoretic interpretation, the term control
graph to describe the set of control vertices is chosen
(Fig. 10). A major advantage of the Bspline formu
lation is that it is a local representation. A bicubic
Bspline surface patch is controlled by 16 control ver
tices and is unaffected by all other control vertices.
I\ An interior
conrro! vertex
/\ / \ / \/ ~ A boundary
 ~1 corporal ve rcex
Fig. 10 A BSpline Control Graph
Conversely, a given control vertex exerts influence
over only 16 surface patches and has no effect on the
remaining patches. This means that the effects of
moving a control vertex are limited to 16 patches. A
point on the (i,j) th uniform bicubic Bspline sur
face patch is a weighted average of the 16 vertices
V`+rj+,, ~ ~ =  2,  1,0,1 and s =  2,  1,0,1 . The
mathematical formulation for the patch Qij(u,v) is
then
1 1
Qij(U,V) = ~ ~ bbre(U,V)Vi+r,j+e
r=  2 A=  2
for O < a, v < 1 (36)
The set of bivariate uniform basis functions is
the tensor product of the set of univariate uniform
basis functions. That is,
bbr,,,(u,v) = b,(U)be(V)
for r =  2,  1, O. 1, and s =  2,  1, O. 1 (37)
('`, I = ~ ~ b,(.)V`+,,,+,b,(~)
for O < u,v < 1 (38)
Fig. 9 A BSpline Surface is a Mosaic
of Surface Patches.
Thus, the formulation for the patch Qij(u,v)
can be rewritten as
The univariate uniform cubic Bspline basis func
tions are graphed in Fig. 11, and can be written in
matrix form as
437
OCR for page 427
[ b2(U) b_l(U) loo(U) bl(u) 1
1 3 3 1
= [u3u2ul]1/6 3 6 3 0 (39)
1 4 1 0
1, ~
b1 (u) loo(u) b ~ (u) b (us
Fig. 11 Graphs of the Univariate Uniform
Cubic BSpline Basis Functions
Appendix 3 Calculation of the Surface
Integral and the Normal Vector
The surface integral can be calculated as fol
lows (Kaplan, 1981~:
[; Hair, y, x) dir
s
sat..
  Of (u, v) ~ g(u, v), h(u, v)]~/EGF2dudv
(40)
where
x = f(u,v), y=g(u,v), z=h(u,v)
P1 = xu~+y~+zuk
P2 = xv~+yuj+zvk

E  IP112 = (0~)2 + (~Y)2 + (~Z)2
XX + Byway + ~Z~Z
MU TV MU TV MU TV
= ( 3~)2 + ( ~Y)2 + ( ~Z)2
F = P1 P2
G = lP2 12 =
The normal vector on the surface S is calcu
lated by using the following formula.
n= Apt pa (4~)
DISCUSSION
Tor Vinje
Norwegian Contractors, Norway
It seems to me that you have treated the computation of the forces in
a way that is much more complicated than necessary. The splitting
up of the B¢/3t values, according to Vinje & Breirg, does not make
sense for forced motion. In this case it can be computed by means
of a central difference scheme at a later stage. If you, on the other
hand, are dealing with problems where the dynamic equations of
motion of the body have to be integrated in time, you have to split
the forces (and B¢/3t) in the way you did. This is due to the
occurrence of numerical instability if you apply, for instance,
backward differences in time to compute B¢/8t, and thus the forces
acting on the body.
AUTHORS' REPLY
Even if the calculation examples in this paper are restricted to the
forced motion, this motion can be directly applied to the problems
where the body motion is unknown. The calculation method of the
forces in this paper is not much complicated because the integral
equation for the potential (7) has the same influence coefficients that
the integral equation for the time derivative of the potential (25).
Also, this method gives more accurate results for the calculation of
the forces than a central difference scheme.
DISCUSSION
Pierre Ferrant
Sirehna, S.A., France
Your computations were run only for w = 1.0. According to my
computations on the submerged heaving sphere, w = 1.0 is a local
minimum for the higher harmonic terms in the force. Furthermore,
this frequency is too high to make higher harmonics appear in the
wave field, and your result for the wave field is purely harmonic and
certainly very close to the result of a fully linearized models as
demonstrated in my paper presented in this conference. For this
geometry and for heaving motion, major nonlinear effects appear
especially low frequency, say about w = 0.4. Future runs of your
code should be situated in this frequency range in order to obtain
significant nonlinear bodywave interaction.
AUTHORS' REPLY
I agree with you. I have some experiences for the computation on a
submerged heaving sphere by using the axisymmetic code (13). The
nonlinear effects appear strongly at low frequencies and small
submerged depth (21). Reference: (21) Kang, C.<, Nonlinear
Free Surface Flows for an Axisymmetric Submerged Body,
submitted to J. of the Society of Naval Architects of Korea.
438