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 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 finite-element 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
wedge-shape 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, two-dimensional 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 Korteweg-de Vries, Boussinesq
and Kadomtsev-Petviashivili equations and Green-Naghdi
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
three-dimensional problems.
In the present paper, we initially planned to investigate
two different computational methods, i.e., the finite-element
method and the fundamental singularity-distribution 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 wedge-shape 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 404
OCR for page 405
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 z-axis, 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 z-direction 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 406
The stationary condition on J = f Ldt is equivalent to
the following Euler-Lagrange 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 three-dimensional 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 right-hand 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 right-hand 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 8-node 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 eigen-value 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 two-dimensional 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 407
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)
pN-lGlF2 _ eN-lG2F
§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 z-direction N = Nz - 1 and
= k"Az with vertical mesh size Liz = 1/N.
The well-known 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 time-domain
analysis. A simple remedy to correct this inherent draw-
back is the use of the so-called '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 408
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
right-hand 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 409
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 fourth-order
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 z-axis is against the
gravity.
4 Numerical Results and Discussions
The computed model as shown in Fig. 3, is a vertical-
sided wedge-shaped 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 water-depth 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 two-dimensional free oscillations
in a three dimensional rectangular tank as an intial-value
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 410
OCR for page 412
OCR for page 413
OCR for page 414
OCR for page 415
OCR for page 416
OCR for page 417
OCR for page 418
OCR for page 419
OCR for page 420
Representative terms from entire chapter:
dispersion relation
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 of-1 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 time-dependent 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 y-direction. Fig. 14b is for the
case of those geometry reduced to one-tenth 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 y-direction.
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
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 non-linear 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
Ship-generated solitons. Proc. 15th Symp. Naval Hy-
drodynam., Hamburg, pp. 347-364.
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. 237-246.
Hess, J. 1977 Progress in the calculation of nonlinear free-
surface problems by surface-singularity techniques. Proc.
2nd Int. Conf. on Num. Ship Hydrodynam., Berkeley,
pp. 278-284.
Hughes, T. J. R. & Brooks, A. 1982 A theoretical frame-
work for Petrov-Galerkin Methods with discontinuous
weighting Functions: Application to the Streamline-
upwind procedure Finite Elements in Fluids, vol. 4,
John Wiley & Sons Ltd., pp. 47-65.
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 free-surface. 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. 395-397.
412
Miles, J. W. 1977 On Hamilton's principle for surface waves.
J. Fluid Mech. 83 pp. 395-397.
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 Two-Dimensional
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.247-4.226
Wu, T. Y. 1981 Long waves in ocean and coastal waters. J.
Engng Mech. Div. ASCE 107 pp. 501-522.
Wu, D. M. & Wu, T. Y. 1982 Three-dimensional nonlinear
long waves due to moving surface pressure. Proc. 14th
Symp. Naval Hydrodynam., Ann Arbor, Mich., pp.
103-129.
Yen, S. M., Lee, K.D. & Akai, T.J. 1977 Finite-element
and finite-difference solutions of nonlinear free-surface
wave problems. Proc. 2nd Int. Conf. on Num. Ship
Hydrodynam., Berkeley, pp. 305-318.
Fig. 7 Computed free surface for Fh =
413
~-
Fig. 8 Computed free surface for Fh = 0 9
414
t—60
Fig. 9 Computed free surface for Fh = 0.8
~—60
Fig. 10 Computed free surface for Fh = 0 7
415
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.
~2
*
an
or 0-4
*
~ 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
shallow-water limit or equivalent assumption on the veloc-
ity field, we compare them with the result of finite element
method with one element in depth-wise direction. For the
finite element method we present two different approxima-
tion in depth-wise 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
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 free-surface
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 speed-amplitude 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
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 Green-Naghdi 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 Green-Naghdi equations show
double valued phase speed (your Fig.16) like
Ll and L2. This double-value of c is not
very surprising since the highest wave is not
the fatest one as shown by Cokelet and
Longuet-Higgins 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 Green-Naghdi
equations can be obtained by variational
methods. Solomon and Miles(1985) derived the
G-N 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 free-surface
elevation (hump) into waves of different wave
numbers.
419
