Below are the first 10 and last 10 pages of uncorrected machine-read 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, chapter-representative 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 three-dimen-
sional body beneath the free surface is solved in the
time domain by a semi-Lagrangian 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 bi-cubic B-spline 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 Runge-Kutta 4-th 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
Longuet-Higgins and Cokelet (1976) presented
a mixed Eulerian-Lagrangian method for following
the time-history of space-periodic 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 time-step 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
free-surface 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.
Two-dimensional 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 428
and computing time. Applications to breaking water
waves over finite-bottom topography and interacting
triads of surface and interracial waves were given.
The semi-Lagrangian 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
three-dimensional free surface flows. The free surface
shape and the forces acting on a three-dimensional
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
three-dimensional problems. The three-dimensional
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 4-th 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(Longuet-Higgins
& 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) =-fix-Y) (6)
2. MATHEMATICAL FORMULATION
Consider an ideal fluid below the surface given
by F(x,t) = 0, where x(z,y,z) is a right-handed
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 fix-y) 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 429
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 (Longuet-Higgins &
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 Sommer-feld condition gave reasonable
results for an axisyrnmetric cylinder heaving in the
still water. Also it was extended to 2-D 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 B-spline 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 B-spline
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 B-spline approximation. There are sev-
eral methods to impose end conditions according to
the geometrical characteristics (Barsky(1982)). The
derivative of B-spline 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 B-spline
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 430
OCR for page 432
OCR for page 433
OCR for page 434
OCR for page 435
OCR for page 436
OCR for page 437
OCR for page 438
Representative terms from entire chapter:
time derivative
(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 u-v 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
non-singular 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
jEG-F2dudv
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 u-axis is not perpendicular to the
vector to (go) along v-axis. 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 three-dimensional 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' d-t
= 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 r-w 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 · Vie-2 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) = llpodS-milk
431
-pun(d)-V. V++ 2V¢e V++gZ)dS
s
-mgk
M = |ip_ x ndS.
s
For three-dimensional bodies the force and mo-
ment are rearranged as follows:
F = F. + F2 + (pgV-mg~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 . Ad-2 Vet . V¢) dS, (33)
SB
Ml = -p If r x n dads, and
SB
M2 = -p || _ x n(V Vet-2V' · 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. 5-6. 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
Saw-tooth 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
7. CONCLUSION
In this paper the nonlinear hydrodynamics of
a three-dimensional 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 4-th 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 Free-Surface Flow Pro-
blems," Journal of Fluid Mechanics, 123, pp477-
501, 1982 .
[3] Barsky, B.A. & Greenberg, D.P., Determining
a Set of B-Spline Control Vertices to Generate an
Interpolating Surface," Computer Graphics and
Image Process 14, pp203-226, 1980.
.
[4~ Barsky, B.A., End Conditions and Boundary
Conditions for Uniform B-Spline Curve and Sur-
face Representations," Computers in Industry 3,
pp. 17-29, 1982.
[51 Dommermuth, D.G. & Yue, D.K.P. 1986, "Study
of Nonlinear Axisymmetric Body-Wave Interac-
tions," In Proc. 16th Symp.on Naval Hydrodyna-
mics, Berkeley.
[6] Dommermuth, D.G., Numerical Methods for
Solving Nonlinear Water-Wave Problems in the
Time Domain," Ph.D. Thesis, MIT, 1987.
[7) Faltinsen, O.M., "Numerical Solution of Tran-
sient Nonlinear Free-Surface 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 3-Dimensional
Free-Surface 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 Surface-Body 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 , Addison-Wesley 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] Longuet-Higgins, 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, ppl-26, 1976
[18] Newman, J.N., Transient Axisymmetric Mo-
tion of a Floating Cylinder," J. Fluid Mech., 157,
ppl7-33, 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
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~
I-TH COS SIN
o
1
2
3
-0. 2850707E-O1 0. OOOOOOOE+OO
0. 1843049E+01 0. 2665849E+OO
-0. 1033978E+00 0.6141333E-O1
0. 6449880E-01 0. 8355246E-02
-0. 511481 lE-02 -0. 1019468E-O1
0. 2418429E-01 0. 2957444E-02
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
I-TH COS SIN
o
2
3
4
5
-0. 1308621E-02 0. OOOOOOOE+OO
0.1927655E+01 0. 1237885E+OO
0.2138726E-03 -0. 1002954E-02
0. 3005543E-01 0.1583521E-02
-0. 1381397E-04 -0. 4408678E-03
0. 2590462E-01 -0. 4266882E-03
Heave Force F/pgKaR3
I-TH COS SIN
-0. 1274143E-01 0.0000000E+OO
-0. 4001656E-03 0. 4211906E-04
0. 2037149E-01 0. 2325045E-O1
-0. 3133378E-03 0. 8355940E-04
-0. 3846344E-03 -0. 2575103E-03
-0. 2204935E-03 0. 9557988E-04
t
SF
y
( )
Fig. 1 Coordinate System
1~ '
1
Fig. 2
Comparison of Heave Force Acting on
the Sphere by 3-D 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
~=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
~rat --------- r2
Appendix 1. The 4th order Runge-Kutta
method [9]
When y' = fishy) is nonlinear, this can be
solved by Runge-Kutta methods. The most com-
monly used Runge-Kutta methods are fourth order
accurate and there are a number of these. The best
known such method (sometimes called the fourth or-
der Runge-Kutta 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)
Z-Force
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 B-spline
Surface Representation [3]
436
A B-spline 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 B-spline 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 B-spline 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 three-dimensional x-y-z space, they are or-
ganized as a two-dimensional 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, ·, m-1; j = 0, · . ., n3
U{(V`;, Vi+~,j) Pi = 0, · , m-1; j = 0, · . . , nil
The interior vertices are the vertices Vij, where
1 < i < m-1 and 1 < j < n-1 And the boundary
vertices are Voj, j = 0, ...,n-1, Vi",i = 0, ...,m-1,
Vmj, j = 1, ..., n, and V`O, i = 1, ..., m . To emphasize
this graph-theoretic interpretation, the term control
graph to describe the set of control vertices is chosen
(Fig. 10). A major advantage of the B-spline formu-
lation is that it is a local representation. A bicubic
B-spline 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 B-Spline 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 B-spline 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 B-Spline Surface is a Mosaic
of Surface Patches.
Thus, the formulation for the patch Qij(u,v)
can be rewritten as
The univariate uniform cubic B-spline basis func-
tions are graphed in Fig. 11, and can be written in
matrix form as
437
[ b-2(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 B-Spline 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)]~/EG-F2dudv
(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 body-wave 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
