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 403
Numerical Computations
for a Nonlinear Free Surface Flow Problem
K. J. Bai and J. W. Kim
Seoul National University
Seoul, Korea
Y. H. Kim
Dae Woo Shipbulding and Heavy Machinery Co.
Jangseungpo, Korea
ABSTRACT
A finiteelement method is presented for solving a nonlin
ear free surface flow problem for a ship moving in a towing
tank. This problem is formulated as an initial/boundary
value problem governed by Laplace's equation. The varia
tional functional used in the present paper is basically sim
ilar to the variational functional given by Luke (1967~.
In the numerical procedure, the variational problem is
reduced to a set of nonlinear ordinary differential equation
where the unknown functions appear as the first deriva
tives with respect to time. In this numerical procedure, the
governing equation, i.e., Laplace's equation, is treated as a
constraint to the ordinary differential equation. In solving
this problem, we employed the lumping and the unwinding
schemes to maintain the numerical stability in the integra
tion with respect to time. To illustrate this method, a simple
wedgeshape ship stretched to the bottom is treated. Com
putations are made for the range of depth Froude number,
Fh = 0.7 ~ 1.1. The computed results show a good agree
ment with the previous results obtained by other methods.
Even though the present computations are made mainly in
the neighborhood of Fh = 1, this method can be applied to
general Froude number.
1 Introiluction
A free surface flow of an inviscid, incompressible fluid
past a ship moving with a constant velocity in a towing tank
is described by an initial/ boundary value problem governed
by Laplace's equation with a free surface as a part of solu
tion.
III the past, the problems of this type were generally
treated after the boundary condition on the free surface was
linearized. Recently, however, there are growing interests in
treating the nonlinear free surface flow problems by various
403
numerical schemes. For example, twodimensional problem
was treated by Washizu et al. (1977), Hess (1977), Korving
& Hermans (1977), Salvesen & van Kerczek (1978), Yen,
Lee & Akai (1977) and Mori & Shin (1988~. Recently, a three
dimensional nonlinear problem is treated by Dommermuth
& Yue (1988~.
There is another line of investigations on the nonlinear
free surface flow problem based on the shallow water approx
imation which results in the Kortewegde Vries, Boussinesq
and KadomtsevPetviashivili equations and GreenNaghdi
formulation. Many references along this approach can be
found in Choi & Mei (1989~. To name few, Choi & Mei
(1989), Ertekin, Webster & Wehausen (1984,1986) treated
threedimensional problems.
In the present paper, we initially planned to investigate
two different computational methods, i.e., the finiteelement
method and the fundamental singularitydistribution method
for a three dimensional nonlinear problem. However, the
progress along the second approach is not quite satisfactory
at this time. Thus we describe only the application of a
finite element method in this paper. One may find the ap
plication of fundamental source distribution method for a
sloshing problem in three dimensional rectangular tank in
Kim (1989~. To illustrate the present numerical method,
computations are made for a wedgeshape ship extended
from the bottom to the free surface. Main computational ef
forts are made for the critical Froude number, i.e., Fh = 1.0,
just because of less computation time compared to the small
depth Froude numbers. To treat the case of small Froude
number, one simply has to take more finite elements. The
computed results show two dimensional solitons periodically
generated in the upstream and complicated three dimen
sional waves in the downstream. The comparisons show a
good agreement with the previous results.
OCR for page 403
OCR for page 403
and then with respect to by,
~ J{o ~ I/SF
—/~ID V47 V§g} dV ~
=  dt ~ J~J. (O— ~ And) 8+ dS (19)
+ fff v2¢ b~ ~ ~ .
Here, [J = [Js + [J'
Equation (18) shows that the stationary condition on
J for the variation of ~ recovers the dynamic free surface
boundary condition in each time and that the wave elevation
at t = 0, t* should be specified as the constraints.
Equation (19) shows that the stationary condition on J
for the variation of ~ recovers the kinematic condition on
SF and the governing equation.
The variational form written above is previously given
by Miles (1977) and is slightly different from that given
by Luke(1967~. In the present variational formulation the
wave elevation ~ is assumed to be known at t = 0, t*, whereas
Luke assumed the potential ~ to be known at both initial
and final times additionally. Specifically, the present varia
tional form is obtained from Luke's form by subtracting the
volume integral of the potential resulted in the process of the
integration by parts with respect to time. The present vari
ational functional has more advantage over the original Luke
variational functional in treating the nonlinear free surface
boundary conditions.
3.2 Finit  Element Discretization.
The original initial/boundary value problem is well de
fined for fib of which the admissible solution should be twice
continuously differentiable in space. However, in the above
variational method, it suffices that the admissible trial func
tions have the square integrable properties of the function
~ and only their first derivatives in space. This enables us
to look for an approximate solution in a wider class in the
variational method.
As in usual finite element analysis we discretize the fluid
domain into finite number of finite elements. In the present
computations, the finite elements are generated such that
the projections on the horizontal plane, i.e. z=0, is un
changed while the other coordinate, i.e. the zaxis, is allowed
to change in time. This simplifies the mesh generation and
the computations considerably. However, this restriction is
not necessary for the present method in applying to a gen
eral problem. Then we approximate ~ in N dimensional
function space whose basis is continuous in D and has con
tinuous derivatives in each element. However, it is allowed
to have a finite discontinuity in the normal derivative across
the common boundaries of the adjacent elements. We de
note the basis of this space by {Ni~i=i,...,N and approximate
by the span of the restrictions of {Ni~i=~,...,N on SF which
is also continuous and piecewise differentiable on SF;
¢(x, y, z, t) = Ditty Ni(z, y, 3; S) (20)
<(x,y,t) = Skits Mk~x,y) (21)
where
Mi~x,y) = ~i`~,Y,8;~z=: k= 1,.,NF (22)
and NF is the number of nodal points on SF and i' is the
nodal number of the basis function Ni of which the node
coincides with that of the free surface node k. Summa
tion conventions for the repeated indices are used here. It
should be noted that the basis function {Ni~i=~,...,N is de
pendent on the free surface shape z = §(x,y,t) but its re
stricition on SF is the function of (x,y) and independent
of (. This special property of {Mk~k=~, ,NF is due to the
simplicity of the fluid domain D which has no variation in z
direction such that we can change the position of the nodal
points only in the zdirection in each time step.
Once the trial function is approximated by using the
above basis function, we obtain the Lagrangian L for this
trial solution as
L = ~i,~Til~ (231)
 2¢iKij~j 2fkPHn
where
Tkl
=  MkM' dS
SF
Kij =  VNi VNj dV (24)
Pal = ,/ MkMI dS.
SF
The tensors Kij, Pal are the kinetic and potential energy
tensor and Tk' is a tensors obtained from the free surface
integral, which can be interpreted as a tensor related to
the transfer rate between these two energies. It is of in
terest to note that in Eq. (24), Tk' = Phi. However Tkl will
be defined differently from this in the computation through
lumping.
405
OCR for page 403
The stationary condition on J = f Ldt is equivalent to
the following EulerLagrange equation
To ~ = Kim ~i, (25)
T,c, ~i' =  2¢i phi JO—Pit 0, (26)
fork=l, ,NF
Kij Hi = 0, i ~ it. (27)
Here Eq. (25) and Eq. (26) are the nonlinear ordinary
differential equation for Id, Ike=, ,NF and Eq. (27) is
the algebraic equations for {¢iJi~i,, which is the constraints
for the above two equations.
3.3 Numerical Dispersion Relation.
If one solves the above variational formulation numeri
cally, one encounters the effect of numerical dispersion de
pending on the specific numerical schemes employed. How
ever, it is difficult to analyze the numerical dispersion erect
in the above nonlinear threedimensional problem. There
fore, we restrict ourselves to a linear problem and test the
numerical dispersion in the following. The linear version of
the Eq.'s (25~~~27) is given as
Tang = K° j~j (32)
K,°j Hi = 0, i ~ it (33)
It can be easily shown that the solution of the above Tkl Ail = Pal ~ (34)
discretized problem satisfies the conservation of mass and
total energy, i.e.
—(~ Pap) = 0 (28)
dt k,]
dd(~qkKijjj+~SkP.l`~) = 0 (29)
s,3 k,!
This property of the conservations is independent of the ten
sor Tk' if it satisfies
~ Ti' = ~ Pa (30)
k k
It should be also pointed out that the direct use of (26)
leads to some difficulty in the computations. This difficulty
arises from the first term in the righthand side which is the
derivative of the kinetic energy tensor with respect to the
wave elevation. We have avoided this difficulty by utilizing
the fact that (26) is equivalent to the condition of vanishing
of the righthand side in Eq. (18~. In Eq. (18) b; can be
regarded as test functions on SF. Then Eq. (26) can be
given as
Talkie =
—  Mk¢Z§t dS
SF
 ~  M~IV¢l dS
Peal .
If the integrals in Eq. (26) and (31) are evalutated ex
actly, they are equivalent. However, in the present compu
tation these integrals are calculated by integral quadrature
rules. Therefore, the conservation of energy may not be
satisfied exactly due to the error caused by the numerical
integration. This test result is given in the next section.
for k = 1, · · ·, NF where K,°j is the linearized kinetic energy
tensor evaluated in the undisturbed fluid domain.
We use the 8node brick element with linear interpolation
along each coordinate. Then by separation of variables the
velocity potential ¢(x, y, a, t) can be represented as a prod
uct of the following three sets of functions:
(x, t) = Hi (t) Xi~x), i = 1, · · ·, No,
(y,t) = ~~.(t)Yj~y), j = 1,···, Ny,
~z,t) = ~k~t)X~(z), k = 1,··., Nz
(35)
(36)
(37)
where Xi~x), Yelp), Z~(z) are the one dimensional basis func
tions and Nz, Ny' Nz are numbers of nodal points in x,y and
z coordinates, respectively. Taking Fourier transformation
from the time domain to the frequency domain, w, we have
the following set of eigenvalue problems.
.~(Xi,X'—kZXiXj~dx~jZ = 0 (38)
t(YiYJ  k2YiYj~dyy>1. = 0 (39)
to
<~2 bail  J Size + k2ZiZj~dz ~' = 0 ' ,
where ~' = LEO) and k2 = kz + k2. Here bit = 1 if i = 1
and zero otherwise. It is understood that the time depen
dent terms are now functions of frequencies without chang
(31) ing their notations.
Moreover, we can treat a twodimensional wave prop
agating in the x direction without loss of generality since
the linear wave is isotropic in OX plane as long as we use
the same mesh size in x and y directions. To obtain the
dispersion relation, we assume the type of solution for By as
~ (x, w) = Ae (41)
where A is a constant amplitude and
406
OCR for page 403
An = ,
n^z
n = 17 2, 3, · · · (42)
is the discrete wave numbers and /`x is the horizontal mesh
size. By solving Eq. (38) and (40) after substituting Eq.
(41) into Eq. (38), the dispersion relation can be obtained
as
where
~15.~
en
~ 10.
k2  6 1—cos(8n~x) <43~ 5
n '` =2 2 + COS(8n^X)
O) = W (en)
pNlGlF2 _ eNlG2F
§2 F1—,CNF2
(44)
2 Em/
Pi = 1 +~U /3  ,Bi(1  6 A, (46)
Gi = 1~2/6,Bi(l+ 3 ), i= li2 (47)
with the number of elements in zdirection N = Nz  1 and
= k"Az with vertical mesh size Liz = 1/N.
The wellknown exact linear dispersion relation corre
sponding to Eq. (43), (44) is
k2zact = (9n (48)
(~)c tact = k tanh k. (49)
The comparisons between these numerical and exact dis
persion relaions are shown in Fig. 1. It can be shown that
for Ax << 1 with en held fixed
kn = k2ZaC`(l+o(Ax2~), (50)
which means that W(8n) has an error of of, however
small values of fly we choose.
Go
..........................
0.5
~2
<1 k
0.0
1
O. ~
    ' // 2 ~
/ /N
'''''''' ''''' '''''''''' ' ''''' ' go'''''''''',? '' ~/~ ''
Exact
...~_...............................................................................
k
6. 8. 10.
Fig. 1 Comparison of numerical and exact dispersion
relations.
(a): Eq. (43) and Eq. (48)
(b): Eq. (44) and Eq. (49)
N: No. of elements in vertical direction.
Fig. 1 shows that the numerical dispersion relation pre
dicts higher values for co than the exact values. This inher
ent property of finite element method based on variational
formulation restricts the numerical stability in timedomain
analysis. A simple remedy to correct this inherent draw
back is the use of the socalled 'lumping' in the tensor To
which is originally equal to Pal as in Eq. (24~. The lumped
tensor Tk' is given as
TO = Ail ~ Pam (51)
m
where Ike is Kronecker's delta. Then the dispersion relation
= W ((7n) (52)
2 + COS(8n ~ 27) W `~9 ~
.
. , j ~ ;,
F.E.M/ /
Exact
.' /
. ~','
......................
0.0 0.2 0.4 0.6 0.8
Onyx
Or
(a)
The order of correction term in Eq. (52) is of) which is
same as the error in W(8n). For the shortest wave length
2~:z:, W(8n) of the lumped case predict one third of the
values of W(8n) of the unlumped case. This means that
with lumping we can take time step three times larger than
the standard finite element method without loss of the or
der of accuracy in the dispersion relation; The dispersion
correction given in Eq. (51) by the lumping also preserves
the conservation of the mass and energy, since Eq. (30) is
satisfied when Tk' is replaced by Tkl defined in Eq. (51~.
3.4 Treatment of the Convective Term.
We assume that a ship is moving with a Etroude number
Fh. Then the pure convection terms are added to free sur
face boundary conditions. Since the presence of convection
term due to the introduction of a moving coordinate does
407
OCR for page 403
not change the original dispersion and nonlinearity, we ana
lyze them separately and then add it to the original problem.
Therefore if suffices to consider a pure convection problem
separately here.
Let a convection problem for ~ be given as
~ =—Fang on SF. (53)
The numerical scheme based on a direct weak formulation
for the above equation is known to have a considerable phase
error for short waves. To reduce the contamination of the
phase errors in the computational domain caused by the
short wave components, an unwinding scheme by using asym
metric test function has been successfully employed as in
Hughes & Brooks (1982~. The discretized form of Eq. (53)
based on this scheme is obtained as
T,dn =—Fhii (Mi+ 2 a Ma) a Mi dS (54)
where ~ is the unwinding parameter and has a value between
O and 1. If c' = 0, it is equivalent to the standard Galerkin
method based on symmetric test functions. If c' = 1, the
righthand side of equation (54) is equivalent to the forward
difference formula in finite difference method.
~ .o
I <' k
3
0.5
3
............. ......................................
. ~
0.0 '/ ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''~''''''''''''''''''''''''''''''''''''~''''''''''''''~
0.0 0.2 0.4 0.6
And 2:
6.
4.
2.
Exact
(a) Real part
0.8 1.0
/
. ~
..... .......... .... ............
....... ....... . .
0.4 0.6 0.8 1.0
6'n~ Z:
Or
(b) Imaginary part
Fig. 2 Numerical dispersion relation due to the con
vective term for Fh = 1. (Eq. (55~)
In a similar procedure given in Eq. (41) through Eq.
(43), the following numerical dispersion relation is obtained
for Fh = 1:
+i
3 sink

Ax 2 + cos(8n~x)
3a 1—cost)
.
/`x 2 + cost)
(55)
The real and imaginary parts of the above equations are
plotted in Fig. 2a and 2b, respectively. In Fig. 2a the real
part of ~ is compared with the exact relation
~ = 6'n. (56)
In Eq. (55) it can be seen that the upwinding procedure
has no effect on the real part of w. The additonal term
in the imaginary part due to the unwinding behaves as a
diffusion term and plays a role of damping mainly to the
short wave component in time. The computations based
on Eq. (54) are made for several values of ~ to test the
effect of the numerical damping. The result is given in the
next section.
3.5 Treatment of the Radiation Condition.
In the present computations, the computation domain is
taken sufficiently large in the upstream direction so that the
generated solitons travelling upstream does not hit the up
stream boundary at the final time of computations. How
ever, along the downstream radiation boundary, we used a
simple boundary condition in a moving frame of reference
as
¢'n = 0 at x = XR (57)
where x = XR iS the radiation boundary at downstrea~rn.
This can be interpreted as the disturbance convected away
with the negative velocity of the moving coordinate. The
present computed flow cases are such that the Froude num
ber is around one. Accordingly, the velocity of the moving
frame of reference is the critical speed. If there are any
reflected waves from the downstream radiation boundary,
then the speed of the reflected waves will not be larger than
the critical speed. Therefore one may expect the reflected
wave cannot contaminate the solution in the computation
domain. Strictly speaking this is true only for a linear prom
fem. However, even in nonlinear case the speed of waves at
downstrem do not exceed the critical speed in most cases.
This simple numerical treatment of the downstream radia
tion boundary condition is tested successfully and the result
is given in the next section.
3.6 Time Integration on the Free Surface.
Up to the previous subsections, we presented the numer
ical treatments of the governing equations, the convective
terms and the radiation conditions as boundary value prob
lem. Once we discretize the computation domain by the
408
OCR for page 403
finite elements and obtain the discretized form of matrix
equation after integrating out all the space variables, then
the problem is reduced to a set of nonlinear ordinary dif
ferential equations to solve. In this reduced equations, the
time derivatives are present in the potential and the free sur
face elevation, both defined on the free surface. The reduced
set of ordinary differential equations becomes the following
form:
Tklp = —FhCt ~ (58)
+Ki~j~j + fin,
T`,` bib = —FhC~j~j—Pad (59)
—Hi Dsk'~— 2 Hi D2kj 45j
for k= 1,~.,NF,
where
Kiwi = f', i ~ it (60)
Hi, = ifs (Mi + 2 azMi)azM' dS
Ckj = ii (Me + 2 ~ Mi)3 Nj dS
Dsk' = ifs MiMia3Ni dS (61)
Dskj = /! M6VNi · VNj dS
fi = Fhii Nine dS.
so
It should be pointed out that Eq. (60) can be interpreted
as a constraint to Eq. (58) and (59~. Equation (60) is
obtained from the boundary value problem for Laplace equa
tion with an essential condition on free surface and a natural
condition on the body surface.
In the solution procedure for this problem, the constraint,
i.e., Eq. (60) is first solved by the Jacobi conjugate gradient
method. Then the matrix Tkl is inverted. Here the other
matrices, which are dependent on the free surface shape,
are treated as known from the previous time step. The fi
nal form in Eq. (58) and (59) is solved by the fourthorder
Runge Kutta method.
1
2b ~ ~ ~ ~ 2B
_ L
:\\\\\\W y = B A\\\\\\ ,
Fig. 3 Sketch of ship model and tank in horizon
tal plane view. Here the zaxis is against the
gravity.
4 Numerical Results and Discussions
The computed model as shown in Fig. 3, is a vertical
sided wedgeshaped ship extended from the free surface to
the bottom of the tank. The length of the wedge ship, the
length of the parallel middle body, and the beam are denoted
by L, Em, and 2b, respectively. The tank has the width of
2B and the mean waterdepth of h. The ship is assumed to
move with a constant velocity of—U along the centerline of
the tank. In presenting our computed results,all the phys
ical quantities are nondimensionalized by h, ph3 and
for the length, mass, and time, respectively as mentioned
previously.
Before we made the computations for the ship in the
tank, we tested the case of a twodimensional free oscillations
in a three dimensional rectangular tank as an intialvalue
problem. The results are given in Fig. 4. Here the motion
started from an initial hump on the free surface and let
the time increases. The normalized tank length is 40 and
the final time in nondimensional scale is 60 with the time
step of 0.2 . Fig. 4a and 4b show the results of 15 wave
elevations with time increment of 4. The computed wave
elevation and the speed of propagation agreed well with the
approximate theoretical results for the amplitude and speed
of the wave, which is given in Appendix A.
~ .
o
._
0.25
0.00
o
._

n an
t=0 t=36 /
V''.'''''..'.'''''''.''''''.'''''''''''''''''''''''''''''''''.'''''.''''''''''''''''''.'''''.''''''''.''''''''''''.'''.''''.'''''''''.~''.''''''.'..'''''''''''.'..'.'''.'''''''''
0 10.
(a)
20. 30. 40.
X
0.75
0.50
0.25
0.00,
0.25
t=60 · t t=40
1 ~
0. 10. 20.
X
(b)
30. 40.
Fig. 4 Numerical test of free oscillation in a tank
started from an initial hump at t=O. Tank
length = 40, Ax = 0.4 and At = 0.2
409
OCR for page 403
a, 0.25
ct
—0.2=
0.50
o
,,., 0.25
ct
0.00
.
(a)
t
0. 10. 20.
0.25
—20.  10. 0.
x
(b)
the formulation defined in the
. , . at, , . ~ .
inertial coordinates and the
results are shifted to—`,/ght to compare those given in Fig.
6b through Fig. Be. Fig. 6b through Fig. Be are obtained
from the results computed in the coordinates moving with
the velocity of1 for four values of the unwinding parame
ter, cz = 0.01, 0.03, 0.05, 0.1. These comparisons show the
dual roles of upwinding on the numerical damping and sta
bility. This means that too much unwinding would give too
much numerical damping effect in the numerical results. In
the moderate range of the values of a, the unwinding works
favorably to the stability.
Fig. 7 through Fig.12 are the results for the following
computation conditions:
Wedge model:
Computation Domain:
Finite Element Meshes:
Upwinding Parameter:
Mesh sizes:
L=8, b=0.4, Lm=4
x=30, 30),
y=~4, 0) (Symmetry is used)
150 x 10 x 1 elements
c'= 0.05
Ax=0.4, /`y=0.4, /`t=0.2
Fig. 7 throughout Fig. 10 shows the.results of wave
profiles as time increases with the time increment of 20, for
the cases of Fh =1.0, 0.9, 0.8, 0.7, respectively. The gen
eration of the upstream solitons are pronounced as the Fh
approaches to 1. For Fh = 0.7, the generation of upstream
solitons is barely noticeable. Fig. 11 shows the results of
supercritical Froude number, i.e., Fh=1.Q5. It is of interest
to note that the numerical instability shows up when t=60.
Further computations for t larger than 60 showed that the
water surface hit the bottom locally. At present, our pro
gram cannot incorporate with a local dry bottom. This can
be easily incorporated in the near future.
Fig. 12a and 12b show the computed wave resistance
for Fh=O.9 and 1.0, respectively. Fig. 13a and 13b show
the time record of wave elevations for Fh= 0.9 and 1.0, re
spectively. This shows that the timedependent hydrostatic
pressure is not negligible in the resistance computations.
Fig. 14a and 14b show the computed results of wave
resistance for two different tank conditions: Fig. 14a is for
the case of the tank and ship geometry being reduced to one
half in the direction of the ydirection. Fig. 14b is for the
case of those geometry reduced to onetenth in that direc
tion. Thus the wave resistances are shown by multiplied by
the factor of two and ten, in Fig.14a and 14b, respectively.
Table 1 shows the amplitude, speed and the generation
period of solitons for the blockage coefficient, S = 0.1. Here
the blockage coefficient, S is defined as in Ertekin et al.
(1986~. The amplitude A is taken from the first solitons.
The period is measured from the computed results at the
F.P. Also shown are the computed results by Ertekin et
al. who treated a pressure patch in Table 1.
Table 2 shows the amplitude, speed and the generation
period of solitons obtained in the tank reduced to one half of
the tank width for several values of the blockage coefficients
when Fh = 1. The amplitudes and the speed increases
whereas the generation period decreases as the blockage co
e~cient increases. This trend agrees well with the previous
experimental results of Ertekin.
Table 1. Amplitude, Speed and generation period of
solitons for the blockage coefficient, S=0.1.
The values in the parentheses are obtained
from Ertekin et al. (1986~.
U/~ A/h C/~b UT,/h
0.7 0.139 1.06 14.6
0.8 0.240 1.10 19.4
0.9 0.375 1.16 24.1
(0.5101) (1.224) (20.0)
1.0 0.553 1.24 30.0
(0.6248) (1.280) (29.6)
1.05* 0.677 1.27 34.7
1.1* 0.806 1.32 39.1
(0.7729) (1.390) 39.3
* The values are obtained from the narrow
tank reduced to the 1/10 in ydirection.
Table 2. Amplitude, Speed and generation period of
solitons obtained in the tank reduced to
one half of the tank width for U/~iZ = 1.
Block. coeff.
A/h
C/
UT§/h
0.06
0.08
0.1
0.12
0.404
0.475
0.545
0.604
1.18
1.22
1.24
1.26
44.2
36.2
29.8
27.0
As a concluding remark the accuracy of the present com
putations could be improved by employing finer meshes.
The computation time for a typical model of 1500 elements
was 260 seconds for each time step by IBM PC/XT with
T800 Monoputer. A typical number of the total time steps
was 500.
ACKNOWLEDGEMENTS.
This work has been supported by the Nonlinear Ship Hy
drodynamics Program supported by the Korean Science &
Engineering Foundation.
411
OCR for page 403
REFERENCES
Choi, H. S. & Mel, C. C. 1989 Wave resistance and squat
of a slender ship moving near the critical speed in re
stricted water. Proc. 5th Int. Conf. on Num. Ship
Nydrodynam., Hiroshima.
Dommermuth, D. G. & Yue D. 1988 The nonlinear three
dimensional waves generated by a moving surface dis
turbances. PTOC. 17th Symp. Naval. Hydrodynam.,
Hague.
Ertekin, R. C. 1984 Soliton generated by moving distur
bavces in shallow water: theory, computation and ex
periment. Ph.D. di~eration, University of California,
Berkeley.
Ertekin, R. C., Webster, W. C. & Wehausen, J. V. 1984
Shipgenerated solitons. Proc. 15th Symp. Naval Hy
drodynam., Hamburg, pp. 347364.
Ertekin, R. C., Webster, W. C. & Wehausen, J. V. 1986
Waves caused by a moving disturbance in a shallow
channel of finite width. J. Fluid Mech., 169 pp. 347
364.
Green, A. E. & Naghdi, P. M. 1976 A derivation of equations
for wave propagation in water of variable depth. J.
Fluid Mech. 78 pp. 237246.
Hess, J. 1977 Progress in the calculation of nonlinear free
surface problems by surfacesingularity techniques. Proc.
2nd Int. Conf. on Num. Ship Hydrodynam., Berkeley,
pp. 278284.
Hughes, T. J. R. & Brooks, A. 1982 A theoretical frame
work for PetrovGalerkin Methods with discontinuous
weighting Functions: Application to the Streamline
upwind procedure Finite Elements in Fluids, vol. 4,
John Wiley & Sons Ltd., pp. 4765.
Kim, Y. H. 1989 A numerical analysis of free surface wave
problems by source distribution method. M.S. Thesis,
Seoul Nat'1 Univ. (in Korean).
Korving, C. & Hermans, A. J. 1977 The wave resistance
for flow problems with a freesurface. Proc. 2nd Int.
Conf. on Num. Ship Hydrodynam., Berkeley, pp. 285
291.
Luke, J. C. 1967 A variational principle for a fluid with a
free surface. J. Fluid Mech. 27 pp. 395397.
412
Miles, J. W. 1977 On Hamilton's principle for surface waves.
J. Fluid Mech. 83 pp. 395397.
Mori, K. & Shin, M. 1988 Su~breaking wave: it's character
istics, appearing condition and numerical simulation.
Proc. 17th Symp. Naval. Hydrodynam., Hague.
Salvesen, N. & van Kerczek, C. 1978 Nonlinear Aspects of
Subcritical Shallow Water Flow Past TwoDimensional
Obstructions. J. Ship Res. vol. 22, No. 4, pp. 179
200.
Washizu, K., Nakayama, T. & Ikegawa, M. 1977 Applica
tion of the finite element method to some free surface
fluid problems. Finite Elements in Water Resources,
Pentech Press, London, pp. 4.2474.226
Wu, T. Y. 1981 Long waves in ocean and coastal waters. J.
Engng Mech. Div. ASCE 107 pp. 501522.
Wu, D. M. & Wu, T. Y. 1982 Threedimensional nonlinear
long waves due to moving surface pressure. Proc. 14th
Symp. Naval Hydrodynam., Ann Arbor, Mich., pp.
103129.
Yen, S. M., Lee, K.D. & Akai, T.J. 1977 Finiteelement
and finitedifference solutions of nonlinear freesurface
wave problems. Proc. 2nd Int. Conf. on Num. Ship
Hydrodynam., Berkeley, pp. 305318.
OCR for page 403
Fig. 7 Computed free surface for Fh =
413
~
OCR for page 403
Fig. 8 Computed free surface for Fh = 0 9
414
OCR for page 403
t—60
Fig. 9 Computed free surface for Fh = 0.8
~—60
Fig. 10 Computed free surface for Fh = 0 7
415
OCR for page 403
0.6
...
0.4
v
t=6G
Fig. 11 Computed free surface for Fh = 1.05
0.2 1
0.0
0.6
.......................
0.4
0.2
:
0.0 1
/ ~ ~
/
, ~
/ , .: ..... : ..... ........
..... ............................. ............ ..... ... .......
~ . ~
0.6
/
/''''''
..~ ~
. 1 _ . .~ ......
/ .. ..
:
~ ............................................
.................................. .............................. .. .... ....................... .........
...................... ' ~ ~
. ~ ......................................
.
~ ' . ~
, . . ~ 0.2
~ ' ~ ~
.... ....................... ..... .................. . .......... .... .. ... .. ............ .... .... ..... ........ .. .. ......... .......... ..... . .... ...
1 ' ~ . 0.0
0. 20. 40. 60. 80. 100. 0. 20. 40. 60. 80. 100.
t t
(a) (a)
................ ........ .. ..... .. ... .. .............. ..... ....... ..... ..... . .... . .. ....
.................. ..... ...... . ... ... ...... ... . . . ......
o
in, 0.4

o
., 0.4
C:

~ 0.2
Ct
0. 20. 40. 60. 80. 1 00. 0  
(b)
t
Fig. 12 Wave resistance for different Froude num
bers (B = 4, b = 0.4).
(a) Fh = 1
(b) Fh = 0.9
.. ... ....
'/ \~/' W`~ it_ _
0.6
2 f
0. 20. 40. 60.
(b)
t
Fig. 13 Wave elevation at F.P of a wedge.
(a) Fh = 1
(b) Fh = 0 9
416
80. 1 00.
OCR for page 403
~2
*
an
or 04
*
~ 0.2
0.4
0.2
/ . ~ '.
/ ..  ..................................
/ :
......... ~ ............................. ......... ..... .... .. .. ... ....... .
/
., ....... .... .......... . . . . . ....... ........ . ...... . ..
40. 60.
o.o 1
0. 20.
0.6
...........................................................
80. 100.
(a)
t
......................... ' ........................... ........................
. / ~
/
o.
GN, Wu :
L1
0.01 ~
so. loo. L2
20. 40. 60.
(b)
t
Fig. 14 Wave resistance for different tank and hull
geometry (Fh = 13.
(a) B = 2, b = 0.2
(b) B = 1, b = 0.1
A pp endix '~3 3
A Comparisons with Other Theories.
The comparisons are made on the linear dispersion rela
tion and the speed of the nonlinear wave of permanent form
among several different approximate theories. We present
the results obtained by the equations derived by Green &
Naghdi(1976)(hereafter the GNequation)andWu(1981~.
Since both equations are derived with the assumption of
shallowwater limit or equivalent assumption on the veloc
ity field, we compare them with the result of finite element
method with one element in depthwise direction. For the
finite element method we present two different approxima
tion in depthwise direction. In the first approximation a
linear interpolation is used. The results presented in the
section 4 were obtained by using this element. In the sec
ond, we use a quadratic interpolation which satisfies the
boundary condition on the bottom. The above two approx
imations will be refered to L1 and L2, respectively, hereafter.
The approximate velocity potentials of L1 and L2 schemes
in two dimensions are given as
L1: ¢' (x, z, t) = fit (z, t) + Z f2 (z, t) (A.1)
L2: Adz, z, t) = fit (x, t) + (z + 1~2 f2 (z, t). (A.2)
A.1 Linear Dispersion Relation
The linear dispersion relations of GN and Wu equations
are identical and can be found in Ertekin (1984~. The re
lations for L1 and L2 scheme can be derived from Eq. (40)
with appropriate basis functions. The results are
2 3k2
3 + k2
2 k2~12 + k2)
· ~ =
12 + 4k2
2 _ k2~15 + k2)
15 + 6k2
(A.3)
(A.4)
(A.5)
The above relations are plotted in Fig. 15 compared with
the exact relation. It can be found that the finite element
method gives the upper bound for the exact values of A,
whereas GN and Wu equation gives the lower bound. The
L2 scheme predicts more closely to the exact values of
whereas the L1 scheme deviates more than others.
6.
5.
4.
2.1...........
0.
O.
... .......... . ...... . . . . .. . . ..... . ......
... ... . .......
..~
............
. ~,_
. ~ ~ ~
..... .
1. 2. 3. 4.
k
Fig. 15 Linear dispersion relations in various schemes.
417
OCR for page 403
A.2 Speed of Permanent Wave Form
The speed of solitary waves for a given amplitude can
be found in Ertekin (1984) for GN equation. The speed
for Wu's equation is given only in limiting case of small
amplitude in Wu & Wu (1982~. But one can derive the
relation using the integral invariants of Wu's equation in
steady state. For L1 and L2 schemes the speed of permanent
wave form can be derived using the dynamic freesurface
condition and their two integral invariants, i.e., mass and
momentum flux. They are given below:
ON: C2 = 1 +A (A.6)
Wu : c2 = Am+ 2) '(log(1 + A)—l + A) (A.7)
L1 : C A + (1 + A)H _ A2 = 0 (A.8)
L2 : C A + (1 + A)H _ A2 = 0 (A.9)
where
CA
U =
1 + A
H = —~+A+~
(A.10)
and where C and A are the normalized speed and amplitude.
The above results are plotted in Fig. 16. It is of interest
to note that in the speedamplitude relation for L1 and L2
schemes two different speeds are present at an amplitude
(note the dotted lines), or vice versa. Presumably this is
due to the presence of waves other than the solitary waves.
Fig. 16 also shows that the L1 (or L2) scheme gives the
bounded amplitude and speed of permanent wave form with
A = 0.600 (or 0.714) for the maximum amplitude, and c2 =
1.473 (or 1.591) for the maximum speed.
0.8
~ 0.6
C~2
* 0.4
0.2
0.0
............................ ,'
................. , . ~ . . ... .
it, ~ at, ·,
, .
....:.. ~ ! ~
.... , // A. .......
:::::::::::1
0.0 0.2 0.4 0.6 0.8 1.0
A
Fig. 16 Speed of the permanent wave form in vari
ous schemes.
418
OCR for page 403
DISCUSSION
by R.C. Ertekin
I think this paper presents an excellent
example of how micro computers can be
effectively utilized to solve complicated
problems in ~ ~ ~ ~ __
solution to the problem is impressive. I
would like to comment on some points in the
paper.
time domain. Their method of
The GreenNaghdi equations that you refer
to are the simplest equations in a series of
equations (or theories) that can be obtained
by the Theory of Directed Fluid Sheets.
Shields (Ph. D. Thesis, U. C. Berkely, 1985)
have derived the theory II and III equations.
Their derivations (with Prof. W. C. Webster),
in some sense, resemble the Pohlhausen method
in laminar boundary layer theory. I think the
higher theory GreenNaghdi equations show
double valued phase speed (your Fig.16) like
Ll and L2. This doublevalue of c is not
very surprising since the highest wave is not
the fatest one as shown by Cokelet and
LonguetHiggins in a J. Fluid Mech. article in
70's. If you used a cubic interpolation
function, I believe the maximum amplitude
would become closer to A~0.78as shown by
Fenton in also a JFM article who used a very
high order perturbation expansion.
As you commented the GreenNaghdi
equations can be obtained by variational
methods. Solomon and Miles(1985) derived the
GN equations by using Hamilton's principle
(J. Fluid Mechanics). I will be happy to
furnish these references in a more complete
form if you wish.
I am not sure that "dry mode" is a good
idea. I am not even sure that it happens in
nature. We have experienced such dif
ficulties also. By changing x,Ay and At we
solved the problem. If you reduce your space
and time increments your problem of stability
may disappear.
In free oscillation tests (Fig.4), the
oscillations in the trailing part of the waves
are perhaps partially numerical. I remember
the thesis of H. Schember (Cal. Inst. of
Technology, Ph. D. thesis, 1982) who studied a
similar problem on these oscillations.
Author's Reply
Thank you for your nice comments. We have
admired your pioneering research on the Green
Naghdi method applied to the nonlinear free
surface flow problems.
We were the comment on the dry mode, we
partly agree with Professor Ertekin on the
possibility of dry bottom due to the numerical
instability. However, we have observed the
dry bottom in the experiments when we moved a
vertical flat plate in the direction normal to
the plate and along the centerline of the tank
in a shallow water in a wave tank (22cm X 20cm
X 110cm). The dry bottom was possible in
nature when the flat plate is moving so fast
that the water behind cannot fill the empty
space fast enough. Then there should be
locally dry bottom in nature.
On the last comment about Fig.4 for the
free oscillation test, we believe that the
trailing part of the smaller waves is due to
the dispersion of the initial freesurface
elevation (hump) into waves of different wave
numbers.
419
OCR for page 403