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 351
TwoDimensional Numerical Modelling of Large Motions
of Floating Bodies in Waves
D. Sen
Memorial University of Newfoundland, Newfoundland, Canada
J. S. Pawlowski
NRC and Memorial University of Newfoundland, Newfoundland, Canada
J. Lever and M. J. Hinchey
Memorial University of Newfoundland, Newfoundland, Canada
Abstract
A numerical method is described, which
s imulates in the time domain the propagation
of steep two dimensional periodic waves and the
large motions induced by the waves on free
floating bodies. The method allows for mild
transient phenomena. The algorithm is based
on Green's formula for harmonic functions
applied in a finite control fluid domain,
together with the fundamental solutions of
Laplace ' s equation. A lowest order boundary
element discretization is used. In addition
to several numerical results, computations of
the sway forces and the roll and heave motions
induced by steep periodic waves on a floating
body restrained in the sway mode are presented
and compared with the results of specially
conducted model tests.
1. Introduction
The increasing accessibility of computers
of high capacity and the parallel development
of computing techniques lead to the Leasability
of algorithmic solutions of complicated
hydrodynamic problems through a discretization
of the corresponding governing equations in
their fundamental form. In this paper, a
solution of this type is discussed for the
problem of motions induced by steep water waves
on a floating body. The problem is formulated
in two dimensions, in the time domain and
potential fluid flow is assumed.
A number of numerical models of the same
problem have been developed s ince the
appearance of the integral equation formulation
combined with the time stepping scheme for the
non linear free surface conditions originally
employed by LonguetHiggins and Cokelet [1] in
a study of the propagation of steep waves.
Faltinsen [ 2 ] cons idered forced heave motions
of a twodimensional circular cylinder, as well
as a related problem of sloshing [ 3 ] . Vinj e
and coworkers extended their earlier work on
breaking wave s imulation (Vinj e and Brevig [ 4 ] )
to include submerged and surfacepiercing
bodies in the fluid [ 5, 6 ], and next attempted
351
the problem of motions of floating bodies
[ 7, 8 ] . Subsequently, the simulation of a
capsizing of Salter' s duck, an ocean energy
extracting device, was reported in [ 9 ] . For
this latter study, experimental results
supplemented the numerical simulation.
Following largely the techniques of Vinj e
and coworkers, with an important modification
in the cons iteration of the body free surface
intersection point, Lin [ 10 ] simulated two 
dimensional waves generated by a wavemaker in
a f inite rectangular tank . This s tudy was
followed by an extension to forced oscillations
of axisymmetric threedimensional cylinders,
[11] . Subsequently, Dommermuth and Yue [12,13]
investigated the threedimensional axi
symmetric problem and were able to s imulate
forced heave oscillations of large amplitude
of cylinders and inverted cones in an otherwise
undisturbed free surface. Greenhow and co
workers also employed the approach of Vinje and
applied the method in a study of the two
dimensional impact problem [14,15,16]. In the
adaptation, specific improvements and
developments of the algorithm were made to make
it suitable for the particular application.
Isaacson [ 17 ,18 ] reported on a similar method
for studying two and threedimensional fixed
and floating body problems. In spite of those
efforts, to the knowledge of the authors, a
fully satisfactory solution of the problem has
not been achieved.
The study presented here was undertaken to
investigate the possibility of developing a
relatively simple but robust numerical model
of the problem, which would be as suitable as
possible for future generalizations to three
dimens ions . This in particular lead to the
requirement that the wave excitation in the
control fluid domain and the radiation
condition for the outgoing waves be implemented
in forms which would not preclude an efficient
modell ing of open water conditions . S imilarly,
the modelling of wave breaking was excluded in
order to ensure long periods of simulation for
steep wave conditions. In this way the
presented model acquired its characteristic
features which distinguish it from the models
OCR for page 351
described above.
The computer codes used in the study were
written in FORTRAN and the computations were
performed on DECVAX 8800 and VAX 8530 cluster
at Memorial University of Newfoundland. The
study as presented here was completed in April
1988.
2. The Governing Equations and
Their Boundary Element Formulation
2.1 The Governing Equations
A finite twodimensional control domain, D,
containing fluid is considered. The domain is
bounded by a piecewise smooth contour ~D, as
shown in Fig. 1. The boundary ED is composed
of the free surface IF, an impermeable bottom
ADD, and the wetted contour DUB of a partially
or fully submerged impermeable rigid body B.
The remaining part of ED consists of two
vertical geometrical control boundaries ~DC1
and ~DC2, separated by such a distance that D
contains the submerged part of B. A Cartesian
coordinate system Oxz is chosen, with the z
axis directed vertically upwards and origin O
coinciding with the intersection of ~DC1 with
the undisturbed free surface.
The fluid in D is considered to be inviscid,
incompressible and homogeneous, and the flow
is assumed to be irrotational. If B is
completely submerged, thereby rendering D a
multiply connected domain, the flow is assumed
to be acyclic. Therefore the flow in D is
described by a single valued velocity potential
¢(x,t), where x denotes the radius vector and
t denotes time. The continuity of flow
requires that the potential should satisfy
Laplace's equation in D:
a2¢ a2¢
+ = 0
~x2 ~z2
On the free surface IF the kinematic
condition:
30 = 3¢ 3¢ 30 (2)
Bt Liz ax ax
and the dynamic condition:
it ~ { P + go + 2 (V¢)2}
(3)
are imposed, where ~ denotes the free surface
elevation, g denotes the acceleration of
gravity, and V represents the gradient
operator: V ~ B/8x; Pa is the applied pressure
on the free surface, set equal to 0 in the
sequel, and p signifies fluid density. On ODD
the impermeability condition is applied
Bi  O
an
where B/3n ~ n.a/8x, and n designates the unit
normal vector on aD directed outwards of D.
The same condition imposed on the wetted body
surface aDB gives:
In n
(5)
where Vn is the velocity component of aDB along
its (outward to D) normal. For bodies fixed
in space, Vn = 0.
In addition, boundary conditions on control
boundaries aDC1 and ~DC2 are assumed to be
applied in such a form that either ~ or B¢/3n
are known at all time instants t _ 0. The
imposed conditions are described later in
connection with specific applications.
2.2 The Boundary Element Formulation and Its
Discretization
2.2.1 The Boundary Integral Formulation
The application of Green's second identity
to ~ and the fundamental solution In r(P,Q),
with point P located in D or on ~D, results in
the well known Green's formula for harmonic
functions [19]:
O(P)~(P) = lad [(A) an In r(P,Q) 
8 ¢(Q) In r(P,Q)] dS(Q) (6)
ant)
with point Q located on ~D. Here r(P,Q) ~
X(P)  x(Q) is the distance between points
P and Q; the subscript of B/Bn indicates the
point at which the differentiation is taken;
n(P) = 2x for P located inside D but not on
~D. For P located on ~D, O(P) is the angle
subtended inside D by the tangents to ED at P.
which is equal to ~ when the normal to aD is
continuous at P and the Cauchy principal value
of the integral over ED is taken.
Formula (6) expresses the potential ~ at any
point P inside D in terms of its boundary
values and those of its normal derivative.
When P is placed on ED and values of ~ are
prescribed on a part of ED and those of B¢/3n
are imposed on the remaining part, (6)
represents two coupled Fredholm's integral
equations of the second kind with respect to
¢(P) as the unknown, and of the first kind with
respect to B¢(P)/3n as the unknown. Formula
(6) is valid at any instant of time, and
therefore for solutions advancing in time, this
relation can be applied at every consecutive
time step.
2.2.2 Discretization
(4)
It is assumed that ED consists of Ns piece
wise smooth open contours:
352
OCR for page 351
Ns
ED ~ U ask
k=1
(7)
The open contours ask are further subdivided
into a finite number of segments, each
approximated by a straight line between its end
points:
Mk
Ok ~ U ~ S ,
it
(8)
A collocation point Qk is chosen on each of
the segments bSj. The functions ~ and B¢/3n on
AD are approximated by constant values on each
segment and the values are determined at the
corresponding collocation points Qk . Following
the usual practice, the collocation points are
placed at the centre of each segment. However,
in principle Qk need not be located centrally
in Ask. The discretization scheme maintains
a consistent order of approximation (see e. g.
[20] ) and represents an application of the
lowest order boundary element method or b.e.m.,
[ 21 ] . A system of linear algebraic equations
is thus obtained from (6), (7) and (8), with
respect to the unknown ok or B¢k/3nj values at
collocation points Qk. The influence
coefficients in this system depend only on the
geometry of the boundary contour and can be
determined explicitly (see e . g. [ 3 ] ) for the
present type of straight line segments.
The bottom condition (4) permits exclusion
of ADD from the contour of integration in the
integrand in (6) if ADD is horizontal, by
augmenting the fundamental singularity with its
symmetric image with respect to ADD. Thus,
when the bottom is horizontal at a depth d, In
r(P,Q) in (6) is replaced by In r(P,Q) + In
r(P,Q'), and ADD is discarded from aD. Here
Q' is the symmetric image of Q with respect to
aDD. This results in a reduction of the system
of linear equations by the number of segments
used in discretizing ADD. The influence
coefficients then contain additional
integrations over the image segments.
It is possible to apply a higher order
formulation of the b . e . m., for example, an
isoparametric representation of Ok and B¢k/3nj
over the segments, {Sk (see e . g. [21] ) . Such
refinements are achievable at the expense of
additional algorithmic complexities. In
particular, complications arise in the
treatment of the point at which the normal to
AD is discontinuous or at the intersection of
two parts of the boundary on each of which
different conditions (i. e. ~ or B¢/3n) are
imposed, which in the present application are
the features of the body free surface
intersection point. In addition to a
discontinuity of the normal at this point,
conditions posed on the body and the free
surface are different (Neumann condition on the
body boundary and Dirichlet condition on the
free surface) . In this respect, the central
collocation discretization scheme adopted here
353
avoids an explicit occurrence of the problems
and as a result is the most straightforward for
numerical implementation. In addition, higher
order applications of b. e .m., which lead to a
better resolution with lesser number of
segments (see e.g. [22]), are known to be more
susceptible to numerical instability [ 24 ] which
is a maj or concern in the considered
application.
The present choice of discretization scheme
is founded in the belief that a 'workable '
model can be developed based on this simplest
scheme, for the final task of the simulation
of large motions of floating bodies in steep
waves, since it is supposed that several of the
anticipated problems may not necessarily be
remedied by applying more refined
discretizations but initially may even be
obscured by the inherent difficulties of such
applications. More refined discretizations
can, in principle, be adopted latter.
2.3 The Basic Algorithm of Solution
The general discretization scheme described
above is applied to the boundary value problem
( 1 ) to ( 5 ) by identifying the open contours
Dk with the free surface aDF, the bottom ADD,
the body contour DUB and the control boundaries
aDC1 and ~DC2. The simulation sought for is
accomplished through a time marching procedure.
Initial conditions provide the boundary data
on ED at t e 0 ~ on IF, B¢/8n on ADD and DUB ~
and either ~ or B¢/8n on ~DC1 and ~Dc2. From
the solution of the appropriately rearranged
system of linear equations, B¢/3n on IF and
on DUB for that time instant are determined.
The boundary conditions are then invoked to
obtain the required boundary data for the next
time level. In particular, appropriate
evolution equations for the free surface,
deduced from ( 2 ) and ( 3 ), and for the control
boundaries, if applicable, are integrated in
time to determine the updated free surface
contour (i.e. the configuration of IF for the
advanced time level) and the values of ~ or
In the updated boundaries. On the body,
B¢/3n on DUB is related to the body velocity by
virtue of (5), which in turn is related to ~
through the body equations of motion and
Bernoulli's equation. The boundary contour AD
as well as the boundary data for the advanced
time level are now established and the solution
process can be repeated.
The potential ~ at any desired location in
D can be calculated from the discretized form
of relation (6). Other information, e.g. fluid
velocity and pressure are calculable from ~ by
utilizing Bernoulli's equation and employing
numerical difference techniques in space and
time. The evolution of the free surface and
the motion of the body, which constitute
necessary information for advancing the
solution in time, are extracted as the
simulation proceeds.
OCR for page 351
The system of linear equations to be solved,
in general corresponds to a full coefficient
matrix and thus benefits admissible in
solutions of matrices with special features
are not available. In the present algorithm,
a standard IMSL (International Mathematical and
Statistical Library) routine is utilized which
employs a Gaussian elimination technique for
matrix inversion [25].
The evolution equations for the free surface
can be put in the following general form:
dot ~ f(y,t)
In the present algorithm, a 4th order implicit
AdamsBashforthMoulton (ABM) scheme is
adopted and is found to be convergent for all
required integrations. To achieve convergence
to the third significant figure, usually not
more than one corrector step is necessary.
This scheme requires information at the
preceding four steps. In the initiation of the
solution, the first three steps are therefore
treated by means of successively lower order
schemes with a greater number of iterations.
A variety of other schemes exists for the
integration of equations of the form (9), e.g.
RungeKutta schemes, Hamming's method, etc.
Fourth order RungeKutta starters are popular
for analogous initialvalue problems (e.g. [1,2
and 13]. However, the starter scheme employed
here is found adequate for the applications
considered. Limited numerical experiments with
other schemes were also performed and the
algorithm appeared insensitive to the choice
of any particular scheme.
3. Applications to Linear Free
Surface Flow Problems
3.1 General Considerations
A simple means of partly testing the
effectiveness and reliability of the algorithm
is to apply it to problems which involve small
amplitude free surface elevations. The
simulation of the propagation of small
amplitude waves allows one to linearize the
free surface conditions and therefore reduces
the possible sources of numerical errors. It
also allows us to examine and solve the
remaining numerical problems as discussed
below. In addition, solutions of linearized
flows are usually available in closed forms
thereby providing an exact basis for
comparison.
Upon applying the usual approximations, the
free surface conditions (2) and (3) take the
following linearized forms:
an an
_ = _
at az
(10)
354
an
Bt ~ ~
(11)
and are applied on the undisturbed free surface
z ~ 0. In addition, for the applications
considered in this section, the fluid domain
is represented by the rectangular area depicted
in Fig. 1, with the body contour removed. The
bottom surface is taken to be at a constant
depth d. The free surface part of the boundary
on which the integrand in (6) is to be
evaluated remains undisturbed at all time
instants, and B/3n ~ 3/8z on aDF. The entire
boundary aD is therefore independent of time
and consequently the influence coefficients in
the resulting system of linear equations
remains unchanged in time.
3.2 The Simulation of Airy Waves
As a test case, the method was applied to
simulate the propagation of steady Airy waves
in the control domain. In the simulation, the
initial values of the potential on the
undisturbed free surface z ~ 0 were specified
according to the Airy wave potential:
/(x,t) ~ H2T co hi[h (Zd/)/ ] sin 2> (x  at)
(12)
with t ~ 0. This corresponds to a wave of
height H. length ~ and period T. progressing
in the positive x direction with celerity c.
For the linearized simulations, either ~ or
a¢/3n, computed from (12), were applied on the
control boundaries.
The following notation is used in the
discussion of computed results: Ax denotes the
length of the segments (or the spatial grid
size), suffixed appropriately to indicate the
part of aD on which they are chosen, viz. AXF,
xc1, Axc2 are the segment sizes on aDF, aDC,
and aDC2 respectively. The time step size is
denoted by At. Nt respresents the time step
level of computation: Nt = t/^t. The distance
between the control boundaries aDC1 and aDC2,
i.e. the horizontal extent of the free surface,
is denoted by L. In the linearized
simulations, the spatial grid sizes are kept
constant on each part of the boundary, and At
is constant over the entire time of simulation.
As an example, in Figs. 2(a) and (b), the
free surface elevations are shown for a
simulation where Neumann conditions were
imposed on both the control boundaries. The
control domain was relatively long, L = 7A,
with water depth d = 0. 4) . Other parameters
were: AXF, ~XC,, AXC2 = A/2 5 and At = T/40 .
The comparison with theoretical free surface
contours clearly demonstrates that the
algorithm is capable of following the wave
motion with acceptable degree of precision over
long periods of simulation time.
Computations were also performed for a wide
variety of combinations of the spatial and
OCR for page 351
temporal grid sizes, for different values of
L/) and d/A, and for different initial
distributions of the free surface potential
(i.e. initial values of ~ on IF given by (12)
with values of t different than 0). In all
computations, the quality of agreement between
the numerical and theoretical results was
similar to the presented examples. The
numerical solutions did not exhibit any
discernible evidence of degeneration even
after, for example, 400 time steps or up to 10
wave periods.
3.3 The Unsteady Wave Propagation
The simulation of the propagation of
unsteady waves is achieved by specifying an
exciting wave potential on one of the control
boundaries. The fluid in D is initially at
rest with z ~ O as the initial contour of IF.
The potential given by (12), which corresponds
to an Airy wave propagating in the positive x
direction, is applied in a modulated form on
aDC1 at all time instants thy. Therefore, a
disturbance is introduced at one end of the
control domain to excite a fluid motion in the
initially unperturbed fluid in D. For this
simulation, the boundaries aDC1 and aDC2 can be
referred to as upstream and downstream
boundaries respectively.
For the initially unperturbed state of fluid
¢(x,t)t=o = 0 in the entire of D, including ED
(the value of ~ could strictly be any constant,
but it is convenient to make this constant O
by redefining ¢, see e.g. [26]). What is not
so apparent is the requirement that
B¢(x,t)/3tlt=o = 0 be maintained
simultaneously. Examining eqn. (3), it can be
noticed that I t=0 ~ O and ¢~=o ~ O imply
a¢/3t~=o ~ O on aDF. It follows that B¢/3t on
aDC1 must have a zero value at t = 0 for the
compatibility of the initial boundary data, at
the intersection of aDC1 and aDF. The
potential given by (12) maintains To ~ 0 on
aDC1 but ai/at  t=0 has a finite value.
Simulations which did not impose a¢/at  t=0 on
aDC1 were not successful due to a numerical
instability which started at the origin (at
~Dc1n8DF) and slowly spread over the entire
domain. Although this instability is of a weak
type, in the sense that the solution can still
be continued, the free surface contour shows
undesired 'zigzag' patterns and eventually
diverges from the required Airy wave profile.
In the present formulation, a modified
excitation potential ~ is used, defined by
introducing a modulation function M(t):
~ (t) = M(t)~(t) (13)
with:
M(t) = ~ 10.5(1  cOS(xt/~)) t < ~ (14
where ~ > 0. This function has the property
that M(t) It=o ~ O and BM(t)/3t~=o ~ 0.
Therefore, regardless of the form of ~ on aDC1,
the initial values of B~ /at are equal to zero
by virtue of (13) and:
3~*(t) ~ M(t) halt) + aMatt) i(t) (15)
The time span over which the excitation
potential is modulated can be controlled by
selecting an appropriate value of the parameter
p. When (13) is applied with a sufficiently
large value of ~ in M(t), the instability
disappears. With p/T = 1.0, the existence of
some undesired behaviour was still detected.
With further increase of p/T to 2, waves
evolved smoothly although no numerical
smoothing was applied. By inserting (13) in
Bernoulli's equation it is found that the
modulation excludes an impulsive application
of pressure on ~Dc, at too. These
considerations appear to parallel a recent
study of the impulsive wavemaker problem,
[27].
3.4 Compupted Examples
Examples of computed results in terms of the
free surface elevations are shown in Fig. 3.
For these computations, the downstream boundary
was placed at the distance of 2) from the
upstream boundary. The discretization
parameters were: AXE ~ A/24 and ~t/T ~ 1/36,
where ~ and T refer to the length and period
of the excitation wave. The downstream
boundary was considered to be a rigid wall,
thus the condition posed on ~DC2 was B¢(t)/8n
= 0 at all time instants. The water had a
depth of d = 0.5A, and the excitation potential
was modulated over 2T, i.e. p/T = 2. Fig. 3.
shows plots of the evolutions of free surface
elevations at four stations situated at x =
0.26N, 0.50A, 0.74A and 0.98A, together with
the theoretical Airy wave evolutions computed
from (12) at the corresponding stations. For
comparative purposes, the Airy wave evolutions
are also modulated by the same modulation
function. It is clear that for t/T < 5, the
reflected waves do not reach the locations x/A
= 0.98. At locations x/A = 0.24 and 0.50, more
than two periods of steady state results are
achieved.
The algorithm was also applied to simulate
wave generation by piston and flapper types of
wavemakers. The corresponding Neumann
boundary conditions were imposed on a ~DC1
fixed at the mean position of the wavemaker
board. In both cases, the normal velocity
values were modified by the modulating function
(14). Stable propagating waves were simulated.
The gain functions (wave amplitude to half
stroke ratios) were found to be within 0.2% and
1% relative error when compared to the linear
theory values [28].
The presented results exemplify the
robustness of the numerical time domain
355
OCR for page 351
simulation algorithm for fluid flow problems
that include a free surface. Computations were
performed for a number of other combinations
of parameters, and showed a similar quality of
agreement with theoretical solutions. No
rigorous rule could be established for the
minimum size of /\x. As a rough guide, a size
of '\x ~ A/12 was found to describe adequately
the evolution of the free surface for most of
the simulations. Further relaxation resulted
in a lack of resolution, although the fluid
motion could still be followed (which means the
solution did not break down). For temporal
grid s ize, the Courant  Friedrichs  Lewy ( C  F L)
type of condition, [ 29 ]:
NCFL ~ IC I\XI < 1 (16)
with c representing the wave celerity, was used
for guidance, although no formal stability
analysis of the algorithim was carried out.
In the present computations, for most of the
cases, a value of NCFL between 0.5 and 0.7 was
used.
However, the solutions were found to exhibit
a tendency towards numerical instability upon
successive refinements of the spatial mesh
sizes. When a collocation point was located
very close to a corner where the boundary
undergoes a sharp change in curvature, such as
the intersections of ~DF with ~DC1 and ~DC2,
relatively large errors occurred in the
computed velocities at these locations, in
comparison with points far from such corners.
These nonuniform differences appeared to
introduce numerical instability when the grid
size was very fine, typically when Ax/) <
1/100. A similar behaviour of solutions near
corners in applications of boundary element
methods is documented in literature [24]. This
instability is thought not to be a serious
limitation in the applications of the algorithm
but serves to indicate a lower bound on the
grid sizes.
In addition to the 4th order A B M scheme,
other schemes for the integration of eqns. (10)
and (11) were also tried. First order implicit
schemes were found to be inadequate in that the
solution showed poor convergence
characteristics and contained large errors.
In contrast, 2nd order schemes lead to
substantial improvements. Further improvements
were achieved by us ing 3rd and 4th order
schemes, although the relative improvements
between those two latter orders were
practically ins ignificant ( the numerical values
differed only in the sixth significant figure).
In no case was more than one corrector level
required for a convergence to the third
significant figure. It is observed here that,
from a computational point of view, higher
order schemes do not necessarily require
additional computational effort. To start the
4th order ABM schemes, lower order ABM
schemes were found to be adequate.
4. The Unsteady Propagation of Steep Waves
4.1 The Evolution Equations for The Free
Surface
For the s imulation of the propagation of
steep waves, it is necessary to consider the
nonlinear free surface conditions (2) and (3)
without the linearizing approximations. These
equations are to be satisfied on the exact
location of the free surface and therefore the
evolution of the free surface within the
control domain must be followed. As a result
the fluid domain must be redefined at every
consecutive time instant. In addition the
evolution of the boundary data on the free
surface must also be determined.
In the present work, an Eulerian formulation
of the free surface conditions is used.
Assuming the free surface elevation ~ to be a
single valued function of coordinate x, the
evolutions of ,7 and ~ on the instantaneous free
surface are followed at image points of the
undisturbed free surface, obtained by the
proj ection along the vertical axis . The
variation of the potential at the image points
on the free surface, which undergo vertical
displacements, is determined by (see e.g. [3]):
di ~ ~~¢ at + ~¢ d,7 ( 17 )
since ~ ~ ¢(z, t) at these points. Here d,7 is
the differential of the vertical displacement
of the image point:
d'7 = ~ ,7 dt
(18)
From (2), (3), (17) and (18), the evolution
equation for ~ is obtained as:
di = gt7 ~ 21[ (~) ~ (~¢ ) ] 8¢ 3¢ 8?1
(19)
which defines the rate of change of the
potential at the image points. In the present
method, therefore, eqns. (2) and (19)
respectively are the evolution equations to be
integrated in time in order to determine the
instantaneous free surface contour and
potential .
The above method of following the evolution
of free surface is different from the
Lagrangian method utilized in most of the
previous investigations of non linear water
waves, which were also based upon boundary
integral formulations (e . g. [ 1, 5, 11, 13, 30 and
31 ] ) . The attractiveness of the Lagrangian
method follows from its ability to describe
multivalued free surface contours. In
contrast, the applicability of the present
method is restricted to singlevalued free
surface contours. The possibility of
simulating extreme wave conditions, pertinent
356
OCR for page 351
to wave breaking, is therefore excluded.
However, the present method provides several
computational advantages. The image points on
the free surface, which through the
discretizaton are identified with collocation
points, are not allowed to cluster. So, the
scheme avoids the adverse numerical effects
often inherent in the Lagrangian methods in
which the particles tend to concentrate in some
regions. From the previous experience of other
workers, it is known that some form of control
of the Lagrangian points is necessary to
prevent them from clustering, e.g. the
introduction of a 'tangential' velocity
component as discussed by Baker in [30] or a
regrinding of the free surface points at every
time step, as employed in [13]. The present
method of following the free surface is free
from such complications. In addition, the
collocation points cannot leave the
computational domain at any time, therefore the
additional task of tracing such points is
avoided. Computational experience gained from
performed simulations suggests that in the
presented method numerical difficulties arise
when a collocation point is situated very close
to one of the vertical control boundaries. By
preventing horizontal displacements of the
collocation points, such problems are also
. . .
mlnlmlze' ..
Yet another point with regard to the
applicability of the present method needs to
be mentioned. The ultimate objective is to be
able to simulate the motion of a floating body
in steep waves, for a sufficiently long time,
preferably over a number of periods of
oscillation. It must be noted that in the
Lagrangian method the simulation cannot be
extended much beyond the time when the wave
breaks. Similar restrictions in applications
are typical of most finitedifference
algorithms (see e.g. [32]).
4.2 The Simulation Procedure and Wave
Excitation
The simulation of the unsteady propagation
of steep waves is accomplished by a procedure
similar to the one described in §3.3. A wave
potential, representing an oncoming wave
travailing in the positive x direction, is
imposed on ~DC1 at t>0. This applied potential
is hereinafter called the excitation potential,
since it provides the necessary excitation of
fluid motion in D. It was found through
numerical results that the form of the
excitation potential has negligible influence
on the generated numerical wave at points in
the interior of D, provided they are
sufficiently far from the upstream boundary.
An alternative way of simulating waves is
to provide a moving wavemaker at one end of
the control domain. Such a procedure was
applied in earlier works, e.g. [10,11,14]. In
the present method of following the evolution
of the free surface, such an approach would
necessitate either a redistribution of the free
surface grid or a successive introduction and
deletion of collocation points, since the wave
maker enters and withdraws from the free
surface grid. It could cause collocation
points to come very close to the wavemaker and
thus generate numerical difficulties, as it was
explained above.
4.3 The Nonreflective Downstream Boundary
The application of the impermeability
condition at the downstream boundary ~DC2' as
it was done in the applications presented
above, is not satisfactory for simulations of
long duration. An appropriate open boundary,
or radiation condition must be specified, which
makes the boundary sufficiently transmissive
to allow wave patterns generated in D to pass
through the boundary without causing
appreciable reflection effects. In the present
algorithm, a simple open boundary condition is
adopted on ~DC2, which assumes that the
potential at this boundary can be written as
a wave form of the same celerity as that of the
applied excitation potential on aDC,:
¢(x,t) ~ ¢(x  at) (20)
where c represents the celerity of the
excitation wave (cf. eqn. (12)). This results
in the following relation:
ai ~  c a~
Bt an
(21)
in which B/3n  B/3x on ~DC2 was used. Eqn.
(21) has a form similar to Orlanski's radiation
condition but its application here is not
strictly equivalent. In [33] and in many
finitedifference algorithms (see e.g. [34]),
the value of c is taken as the celerity of the
local waves approaching the downstream
boundary, and c is determined by a numerical
differentiation at the neighbouring grid
points. In [35], a similar simple form is
adopted with c determined from:
c ~ j(gd)
where d denotes the local water depth at the
downstream boundary. Eqn. (22) represents the
shallow water approximation of the phase
velocity and is therefore different from the
condition applied in the present method (both
methods become equivalent in the limiting
situation of d/) << 1).
(22)
The evolution of ~ is now easily determined
by the timeintegration of eqn. (21), using the
same scheme as the one used for integrating (2)
and (19). Simple as it appears, in the
considered applications this procedure results
in minimal reflection effects on ~DC2, as the
results presented below illustrate.
4.4 Specific Considerations
4.4.1 Spatial Derivatives on The Free Surface
357
OCR for page 351
The evolution equations ( 2 ) and ( 19 ) require
the evaluation of spatial derivatives of ,7 and
at the collocation points. To determine
B,7/3x, ,7 as a function of x is approximated by
a cubic spline. From this approximation, the
components of the outward normal vector can be
calculated. For the spatial derivatives of ¢:
8¢ an
~ n
as Z ax
(23)
since 1/(ds/dx) ~ nz, where B/3s denotes the
tangential derivative . To determine B¢/3x, in
turn ~ as a function of x is approximated by
a cubic spline. From B¢/3s and B¢/8n and the
outward normal vector, other components of the
spatial derivatives can be determined. In the
software, an IMSL routine for cubic splines
with natural end conditions is used in which
no conditions are prescribed at the end points
and continuity of second derivatives is
enforced at the penultimate points [ 36 ] .
4.4. 2 The Instability at The Intersection of
~DC1 and aDF
When a timemodulated excitation Airy
potential is applied on ~DC1 in the described
algorithm, an instability is found to originate
at the intersection of ~DC1 and IF. The form
of the instability is qualitatively similar to
that observed in the analogous linear
application when the modulation function was
not used. The solution breaks down, typically
within 10 time steps irrespective of the step
size .
The occurrence of these oscillations can be
attributed to the incompatibility of the
excitation wave potential and elevation applied
externally on ~DC1 with the wave potential and
free surface elevation induced in D in the
vicinity of IDA. In other words, the free
surface conditions implicitly satisfied on the
left of ~DC1 (i. e. by the exciting wave
potential) are inconsistent with the conditions
on ~DF immediately to the right of IDA. On
the basis of conducted numerical experiments,
this discontinuity is believed to cause large
velocity gradients across the vertical boundary
and these initiate the instability. The
application of other, 'non linear' excitation
potentials, e.g. Stokes second order potential,
were tried without success. Difficulties
originating from analogous discontinuities were
known earlier, e. g. similar problems were
discussed in [ 37 ] and in [ 11 ], the
discontinuity was identified with the
difficulty of the matching of non linear
interior with linear exterior solutions,
reported in [ 8 ] .
In order to achieve a smooth variation of
the free surface elevation and potential across
~DC1, the matching procedure, described below,
was developed. Another vertical boundary ED c1
in the interior of the control domain D is
introduced at a short distance 1 from IDA.
The free surface elevation and potential
evaluated without the matching procedure are
represented by f2(x). A transfer function g(x)
is introduced to redefine f2(x) as f,?(x) in the
'matching zone' between aDC1 and ED cl:
f2(x) ~ g(x) f1 (x) (24)
where f2(x) is defined in x1 < x < x1 + 1, f1 (x)
indicates the theoretical wave elevation or
potential corresponding to the excitation
potential applied on aDC1, and x1 represents
the x coordinate of aDC1. A quadratic
polynomial is chosen for g(x) whose
coefficients are determined from the
conditions:
f2 (X1 ) f 1 (X1 )
f2(x1 + 1) f 2(X1 + 1) (25)
ax 2(X1 + 1) = aX f2(X1 + 1)
The above procedure requires the evaluation of
a [ f2 (X1 + 1 ) ] /3x and this is determined from a
second order central difference scheme. A
linear form of g(x) was found not as
satisfactory. The quadratic function for g(x)
is however very effective in suppressing the
instability and enables the fluid motion to be
followed without further problems originating
at the intersection point.
On the downstream boundary ~DC2 it is
necessary to determine the intersection of aDF
with ~DC2. This is determined via a second
order Lagrangian scheme us ing the data at the
three preceding collocation points on aDF.
4.4.3 The Instability on The Free Surface
In addition to the instability originating
at the upstream end of the free surface,
another instability was found to develop on the
entire free surface as the solution progressed.
Sawtooth instabilities of this type have been
reported by earlier investigators. Numerical
experiments with various combinations of the
spatial and temporal grid sizes were performed
with the hope of identifying a stability
criterion related to these discretization
parameters. However no such criterion could
be established. In the present formulation,
in which collocation points on the free surface
are restricted to move vertically, the arc
lengths between the adjacent collocation points
never reduce below the horizontal grid spacing.
Consequently, the relation between the time
step size and the horizontal projection of the
free surface segments is easily controlled.
In all computations, the CFL type condition,
see eqn. (16), was satisfied in the entire
fluid domain and throughout the whole
s imulation period. A form of stability
criterion based upon a linear van Neumann
stability analysis for a 4th order RungeKutta
scheme was provided in [ 13 ]:
358
OCR for page 351
8 ~XF
(26)
This condition was also maintained in the
discussed computations.
The present computational experience
indicates that the instability is closely
associated with the shape of the free surface
elevation. It becomes more pronounced as the
wave steepens. It should be observed that in
the analogous linear application, no such
problem was encountered. Computations with
successively higher levels of iteration in the
timeintegration of the free surface conditions
and a close examination of the computed free
surface profiles and boundary data suggest an
insensitivity of the instability to the time
integration schemes. Therefore, the violation
of a stability condition of the above type
might not be the root mechanism in the
initiation of the instability.
To suppress this instability, the smoothing
scheme originally employed by LonguetHiggins
and Cokelet [1] was adopted. In view of their
computationalexperience, the fivepoint scheme
was used instead of a sevenpoint scheme. The
formula provided in [1] is inapplicable for the
special cases of the first and last two
collocation points on the free surface and a
modified scheme had to be applied there. The
smoothing scheme is applied at regular time
intervals. Usually, the application at every
4th time step was found effective in
eliminating the unwanted oscillations. It is
possible to employ a variety of other available
smoothing schemes. A third degree five point
least squares smoothing scheme was also tested
and was found to be equally effective.
Although the application of artificial
smoothing is known to cause a loss of the local
accuracy of the solution, the global solution
fortunately remained within acceptable limits
of accuracy, which was demonstrated by the
computed results. This feature of the
artificial smoothing was also noted by earlier
workers (e.g. [11, 14 and 38]).
4.5 Computed Results
boundary are generally kept fixed in space,
except for the uppermost segment. Depending
on the length of aDC2, a segment is deleted or
an additional segment is introduced so that the
length of the segment in comparison with the
length of the adjacent segment on aDF maintains
a ratio between 0.5 and 2.0.
In the present context, the notation
described above applies with the exception that
xF denotes the spacing of the free surface
collocation points along the x axis, instead
of the actual lengths of the segments. Unless
otherwise specified, the applied excitation
potential on aDC1 is the Airy potential. The
normalizing parameters for horizontal and
vertical length scales and time scale are
respectively the length A, height H and period
T of the Airy wave corresponding to the
excitation potential (cf. eqn. (12)). In all
presented computations, Axe, AXc1, ~Xc2 and At
are constants. However, due to wave run up,
aDC2 continuously changes in length. Since the
left hand side of eqn. (21) is an Eulerian time
derivative, the collocation points on this
Numerical experiments carried out to
investigate the effectiveness of the matching
procedure showed that the number n of
collocation points in the matching zone is the
dominant factor in comparison with the length
1 of the zone. The choice of n ~ 4 was very
effective in removing the oscillations,
regardless of the grid sizes and wave heights.
The subsequent results are all computed with
this value of n = 4.
The convergence characteristics of the
solution in the entire fluid domain was
studied. The chosen fluid domain extended
horizontally over L ~ 2A and vertically over
d ~ 0.5A, and the applied excitation
corresponded to a large nominal wave steepness
of H/> ~ 0.10. For these computations, the
segment sizes were: AXF, AXC1 , AXC2 = A/16,
A/20, A/24, A/28, A/32 and A/40. The time step
size was ~t/T `1/40 for the first five values
of N. while for N ~ 120, it was halved to At/T
= 1/80. This was necessary because of the
reduced segment sizes. Otherwise for At/T =
1/40, which corresponds to NCFL = 1 (cf. eqn.
(16) ), the solution exhibited an instability.
Free surface smoothing was applied at every
fourth step, p/T = 1 was taken in the
modulation function, and n = 4 was used for the
matching region. Except for the value of N
48, the results demonstrated good convergence
characteristics. For this value of N. the
results were affected by comparatively poorer
resolution. These computations (and many
others) indicated that the value of AxF ~ A/24
and comparable values for ~Xc1, AXc2 were
adequate for describing the fluid motion
without appreciable effects of the lack of
resolution.
The effectiveness of the open boundary
condition (21) was examined by selecting a
range of values of the celerity of the outgoing
waves. Taking c in eqn. (21) as c', with c'
= ac, where c represents the celerity of the
exciting wave at aDC1 (as in eqn. (12)),
computations covered a variation of ~ from 0
to 1, with specific values of ~ = 0, 0.25,
0.50, 0.75, 0.90 and 1.00. The other
parameters were chosen as: L = 2~; d/A ~ 0.5;
H/) = 0.10; ~xF, ~xc, and ~xc2 = A/24, and ~t/T
= 1/40. The free surface elevations at the
simulation time of t/T ~ 8.75 are shown in Fig.
4. It is apparent that the reflection effects
at the downstream boundary increase with the
difference between ~ and 1. Results for
1 indicate that the wave passes through ~DC2
with imperceptible reflection. From these and
other computations the effectiveness of
choosing ~ ~ 1.0 for making the downstream
boundary transmissive was evident, although
3S9
OCR for page 351
values slightly less than 1 also appeared to
work well. Computations were also attempted
for values of c' greater than 1, but even for
a value moderately greater than one, e.g.
1.05, the solution broke down after about t
ST, which was approximately the time for the
wave to grow fully at the downstream boundary.
This breakdown resulted from an instability
originating at the downstream intersection of
IF and ~DC2. In view of the success of
1.0 in making the boundary sufficiently non
reflective, this aspect was not pursued any
further.
A comparison of evolutions of wave
elevation, for H/A ~ 0.10 at x/A  0.98, with
theoretical profiles for the Airy wave given
by (12), Stokes second order wave and Miche's
second order theoretical profile [39] are shown
in Fig . 5 . The numerically s Emulated wave
compares well with the second order profiles,
but displays stronger non linear
characteristics . The comparatively more peaky
crest and shallower trough of the computed wave
are visible.
To investigate the influence of the
excitation potential on the interior solution,
computations were performed with a Stokes
second order potential as the excitation
instead of the Airy wave potential. The
applied excitations had a value of H/) ~ 0.10,
for which the second order correction in wave
amplitude is almost 10% of the first order
amplitude, but both excitations have the same
average energy dens ity . The fluid domain
chosen and the discretization parameters were:
L ~ 2>, d ~ 0. 5A, Axe ~ A/24, At ~ T/40 . The
evolutions of wave elevation in time at x/)
0.48 and 1.48 are shown in Figs. 6 (a) and (b).
The differences between the two simulations are
undetectable .
The results presented above show that the
simulation of unsteady steep wave propagation
can be achieved by imposing an excitation
potential on one of the vertical control
boundaries encompass ing a rectangular fluid
domain. The interior solution apppears to be
not sensitive to the exact form of the
potential. The simulated wave profile displays
typical non linear characteristics of
relatively more peaky crest and shallower
trough in comparison with linear waves. As
expected, the non linearities are more
pronounced for s teeper waves . Very s teep waves
were simulated for time durations of over 20
wave periods, for which a steady state
behaviour occurred in the entire domain.
The simple outgoing wave condition (21)
produced good results over the entire period
of each computation as well as for all
combinations of H. T. L and d for which
computations were performed. The interior wave
was not observed to be contaminated by
numerical reflection effects even after long
times of simulation and at locations close to
the downs beam boundary . This demons bated the
effectiveness of (21) in the modelling of the
propagation of nonlinear periodic waves.
5. The Floating Body Problem
5.1 The Governing Equations
In order to simulate motions of a floating
body in steep waves, a floating body B is
introduced in the fluid, such that its wetted
contour DUB is completely contained in D ( see
Fig. 1) . The obj ective is to expose B to an
incident steep wave train and subsequently to
follow the induced motion of B. A propagating
steep wave is generated in D in the manner
described above. For such simulations, it is
necessary to know the exact location of DUB at
every time instant . In addition, a relation
between ~ and B¢/3n on DUB must be established
such that the evolution of boundary data on DUB
can be followed. The required information is
obtained by invoking the equations of rigid
body motion and the impermeability condition.
For the following cons iterations an
additional coordinate system fixed with the
body is introduced. A rectangular Cartes fan
coordinate system Gx' z ' is used with its origin
G located at the body centre of gravity CG, and
Gz' axis directed vertically upwards in the
undisturbed position of the body. The axis
through G perpendicular to the x' and z ' axes
is assumed to coincide with a principal axis
of inertia of the body. The body geometry is
invariant in the coordinate sys tem and
therefore the instantaneous wetted contour DUB
is completely defined by the location and
orientation of Gx'z' system with respect to the
Oxz system and wave elevation. The components
in the Gx'z' system denoted by {x} ', of the
radius vector of a point fixed with the body,
are related to the components in the fixed in
space system, of the radius vectors of the same
point and of CG, denoted by {x) and {XG)
respectively, by the following relation:
_ ~ _
{X  XG} ~ [R]' {x)' (27)
where matrix [R] represents the tensor of
rotation and the superscript T indicates a
transpose. [R] is given by:
[ ] [sin ~ cos ~ ] (28)
where ~ denotes the angular displacement of
Gx'z' system with respect to Oxz system,
measured as positive counterclockwise. The
equations of motion for the body are written
in the Newtonian form:
{Fx, Fz, My) ~ {M} XG' ME ZG, IS 8) (29)
where Fx and Fz are the components of the force
F acting upon the body, whereas XG and ZG are
360
OCR for page 351
those of XG, in x and z directions
correspondingly; Me represents the moment of
force F about the axis passing through G and
orthogonal to Gx' and Gz'; MB denotes the body
mass and IF denotes the body moment of inertia
about the axis with respect to which M} is
defined, and dots indicate differentiation with
respect to time.
The external forces and moment exerted on
B can be obtained by a direct integration of
fluid pressure p on ADS:
FX ~ T3D B PnXdS (30)
Fz = T3D pnzdS  g MB ( 3 1 )
M} ~ T3D B P[ (Z ZG)nX+(X XG)nZ]dS (32)
where nx and nz denote the components of the
normal vector on DUB in the fixed in space
system of reference. The pressure p is
computed from unsteady Bernoulli's equation:
= P{gZ + B~ + 21 (Vi)23 (33)
On the body surface, the fluid normal
velocity B¢/3n is equal to the normal component
of the body velocity by virtue of
impermeability condition (5). Therefore, from
(5) and (27), the following form of the
impermeability condition is obtained:
`~¢ = nX XG + nz ZG + ~ [  (ZZG)nX+(XXG)nZ]
(34)
to be satisfied on ODE. The above expression
provides B¢/3n at any point on DUB in terms of
the body configuration, velocity and geometry.
5.2 The Algorithm for The Computation of
Hydrodynamic Forces
The solution algorithm described in §2. 3 can
now be adapted for the simulation of motions
of B. At any instant, presuming B¢/8n to be
known on ADD, the other boundary data are
determined from the solution of the integral
relation (6). From this, the fluid pressure
p exerted on the body is determined by
utilizing Bernoulli's equation (33).
Subsequently, the fluid excitation loads on B
are determined by a direct integration of
pressure over the wetted body surface, (30) to
(32). To establish the boundary data for the
next time level, the equations of motion are
invoked. By integrating (29), XG ~ VG ~ ~ and ~
at the advanced time are determined. The
necessary boundary data (in) on the body
surface are then established from (34) and the
computation for the next time step can begin.
The use of Bernoulli's equation (33)
requires the evaluation of B¢/3t on the body.
To this end the following relation:
di Bi + ~ 3¢ (35)
at at ax
is applied, where d/dt signifies the rate of
change of any quantity following a material
point of the rigid body and v denotes the
velocity of the point. The quadratic term in
Bernoulli's equation is readily obtainable
from:
(V¢)2 = (~n)2 + (~52 (36)
At a collocation point i the tangential
derivative B/3s of ¢, is determined in the
form:
(pi = (~¢i)i/~)i = ASP (~¢i)i (37)
since for the straight line segments, (Bs/Bi)
= ASP which denotes the length of the segment.
To determine (/i), appropriate secondorder
difference formulae are employed. This is
permissible despite large variations in DUB,
since ~ on the surface is in general a smoothly
varying continuous function, which was
confirmed by plotting he against i for several
conditions. In the present algorithm,
Simpson's rule was applied for the evaluation
of force and moment expressions (30) to (32).
To carry out the computations on the wetted
body surface aDB, it is convenient to describe
the body geometry with respect to the Gx'z'
system, in which it is invariant. Denoted by
D'B, it can be subdivided into segments once
and for all. To determine the wetted contour
DUB (BOB ~ ~D'B), the intersection points ~DF n
DUB need to be found. This is achieved using
an extrapoltion scheme in which a second order
polynomial is fitted to the three points on ~DF
adjacent to DAB However, the application of
a fixed discretization of DAB through the
determination of DAB n ~DF and the subsequent
consideration of only those segments of DAB
which belong to BOB, produced an instability in
the force computation and a resulting
divergence of the solution. This was caused
by the introduction or deletion of segments on
the body near the intersection point with the
free surface, which in turn produced a large
variation of pressure in the computation of the
dynamic part of the pressure distribution.
This problem is overcome by redistributing the
collocation points on the body at every time
step such that the segment sizes vary smoothly
in time. Because of the redistribution of the
collocation points, a spatial interpolation of
becomes necessary in the computation of
d¢/dt. In general, this spatial interpolation
introduces very small approximation errors,
since the changes of collocation points between
two consecutive time steps are very small.
361
OCR for page 351
input conditions. The comparison of these two
records provided the comparison of the oncoming
wave conditions. For presentation, these two
records were synchronized. The synchronization
was standardized by matching the peak of the
simulated wave in the time interval 3
The numerical results obtained by the
described above procedure showed a significant
improvement of the predicted roll motion. The
changes of computed sway force and heave motion
were marginal. Figs. 7 and 8 show the
comparisons of the computed and measured
records, obtained without and with the
inclusion of viscous roll damping respectively.
The results are for the largest wave steepness
H/> = 0.028 at which the worst agreement
between the experiment and computation was
observed. For all three wave steepnesses the
differences between computed and measured peak
topeak roll values reamined below 6% for the
simulations in which eqn. (39) was applied.
6.4 Summary
Taking into account the experimental errors
and inaccuracies of the comparison resulting
from the unsteady characteristics of the
experimental and computed data, the predictions
by the numerical model show good agreement with
experimental records. The observed differences
between computed and measured sway forces and
heave motions result, at least partly, from
the absence of viscous effects in the numerical
model. A similar remark applies to the roll
motions were the inclusion of the viscous roll
damping is semiempirical.
It should be observed that the experimental
data used in the comparison included three
kinds of nonnegligible nonlinear phenomena.
For the shorter waves with higher steepness,
the waves near the body approached the breaking
limit. In a number of tests, a foam formation
was observed. Large heave and relative motions
of the body with respect to waves were observed
in several tests. For these tests, the video
records showed that the runup profiles
resembled closely those generated by the
simulation. In addition, the tests included
the occurrence of roll amplitudes up to 20 deg.
combined with large motions relative to waves.
The comparison confirmed the validity of the
numerical method in the presence of the
described phenomena. It also demonstrated that
the method can provide realistic estimates of
roll motions within the specified range if a
semiempirical model of viscous roll damping
is included in the algorithm.
Conclusion
The numerical results presented above show
that the propagation of steep longcrested
periodic waves, which may include mild
transients, can be modelled numerically in the
time domain by means of a relatively simple,
low order boundary element method algorithm.
The efficiency of a radiation condition of
Orlanski's type is also demonstrated. The use
of the condition makes possible the avoidance
of other means of eliminating the reflection
of waves at the downstream boundary, such as
e.g. the introduction of an artificial wave
damping in the free surface conditions, which
are not compatible with open water conditions.
Similar comments apply to the wave excitation
in the control fluid domain, achieved in the
algorithm by a matching of an imposed exciting
wave potential with the flow in the control
domain, in the vicinity of the upstream control
boundary. The procedure is different from the
more usual simulation of a physical wavemaker,
and corresponds better to open water
conditions. However, in the procedure the
exciting wave potential must be modulated in
time to enforce the compatibility of the
initial conditions at the boundary and in the
control domain. The problem parallels similar
difficulties encountered in the modelling of
the wavemaker.
Another characteristic feature of the
algorithm constitutes the use of an Eulerian
form of the nonlinear free surface conditions
based on the assumption of a singlevalued wave
elevation. At present this form of the
conditions appears to be more suitable for the
extended in time modelling of motions of
floating bodies in steep waves, than its
Lagrangian counterpart, although the latter is
capable of modelling wave breaking. The
application of the Eulerian free surface
conditions makes possible a strict observance
of local stability condition of Courant type,
without the use of procedures (such as e.g. a
regridding of collocation points) which may
introduce a numerical smoothing. The results
presented above suggest that the violation of
a Courant type stability condition may not be
the reason of the occurrence of the typical
instability in the free surface data.
The extension of the basic steep wave
propagation algorithm, which includes the
presence of a free floating body in the control
fluid domain, provided time domain simulations
of body motions. In the simulations steady
state motions of the body excited by steep
periodic waves, preceded by short transients,
were achieved. The computed records compare
well with experimental data, thus illustrating
the applicability of the extended algorithm.
The experimental data were obtained from
specially performed model tests, with a body
model restrained in the sway mode. In several
of those tests significant nonlinear phenomena
related to the wave propagation and interaction
with the body were observed. In the numerical
simulations no significant wave reflection at
the downstream boundary was detected. However,
for the form of the algorithm presented here,
simulation times are limited by the reflection
of waves scattered by the body from the
upstream control boundary. To ensure the
stability of the computation of hydrodynamic
forces and of the integration of the quations
of body motion, special procedures were
developed which include a smoothing of the
hydrodynamic forces in time.
365
OCR for page 351
References
1. LonguetHiggens , M. S . and Cokelet , E . D .,
"The deformation of steep surface waves on
water: I. a numerical method of
computation." Proc. R. Soc. Land., Series
A, 350, pp. 126 (1976) .
2. Faltinsen, O.M., "Numerical solution of
transient nonlinear free  surface motion
outside or inside moving bodies. " Proc. 2nd
Int. Conf. Num. Ship Hydro., Berkeley, CA,
pp . 347  357 (1977) ~
3. Faltinsen, O . M ., "A numerical nonlinear
method of sloshing in tanks with two
dimensional flow. " J. Ship Res., 22, 3, pp.
193  202 (1978) .
4. Vinj e , T . and Brevig, P ., "Numerical
simulation of breaking waves. " Adv. Water
Resources, 4, pp. 7782 (1981) .
5. Vinj e, T. and Brevig, P., "Numerical
calculation of forces from breaking waves. "
Int. Sym Hydro. Ocean Eng., Norwegian Inst.
Tech., Trondheim, pp. 547565 (1981). 19
6. Brevig, P., Greenhow, M. and Vinje , T.,
" Extreme wave forces on submerged energy
devices." Appld. Ocean Res., 4, 4, pp. 219
225 (1982).
7. Vinj e, T. and Brevig, P., "Nonlinear ship
motions. " Proc. 3rd Int. Conf. Num. Ship
Hydro., Paris, pp. 257268 (1981).
8. Vinj e , T ., Maogang, X. and Brevig, P ., "A
numerical approach to nonlinear ship
motion. " Proc. 14th ONR Symp. Naval Hydro.,
National Academy Press, pp. 245278 (1982).
9. Greenhow, M ., Vinj e, T ., Brevig, P . and
Taylor, J., "A theoretical and experimental
study of the capsize of Salter's duck in
extreme waves. " J. Fluid Mech., 118, pp.
221 239 (1982) .
10. Lin, W. M., "Nonlinear motion of the free
surface near a moving body. " Ph.D. Thesis, 23
Dept . Ocean Eng ., MIT , 127 p . (1984) .
11. Lin, W. M., Newman, J. N. and Yue, D.K.P.,
"Nonlinear forced motions of floating
bodies. " Proc. 15th ONR Symp. Naval
Hydro ., National Academy Press, Washington,
D. C., pp . 3349 (1989) .
12. Dommermuth, D.G. and Yue, D.K.P., "Study of
nonlinear axisymmetric bodywave
interactions. " Proc. 16th ONR Symp. Naval
Hydro., National Academic Press,
Washington, D . C ., pp . 116  136 (1986) . 25.
13. Dommermuth, D . G . and Yue, D . K . P .,
"Numerical simulations of non linear
axisymmetric flows with a free surface. "
J. Fluid Mech., 178, pp. 195219 (1987) .
14. Greenhow, M. and Lin, W. M ., "Numerical
simulation of free surface flows generated
by wedgeentry and wavemaker motions. "
Proc. 4th Int. Conf. Num. Ship Hydro.,
National Academy of Sciences, Washington,
D.C. pp. 94106 (1985).
15. Greenhow, M., "Water entry and exit of a
horizontal cylinder." Proc. Second
International Workshop on Water Waves and
Floating Bodies, Rep. No. AM876, Univ.
Bristol, pp. 2528 (1987) .
16. Greenhow, M., "Wedge entry into initially
calm water. " Appld. Ocean Res., At press
(1988) .
17. Isaacson, M. St . Q ., "Nonlinear wave 
effects on fixed and floating bodies. " J.
Fluid Mech., 120, pp . 267  281, also
Corrigendum, 133, pp. 469 (1982).
18. Isaacson, M. St. Q., "Steep wave forces on
large offshore structures." Soc. Petrol
Eng. J., 23, 1, pp. 184190 (1983).
Krzyzanski, M., "Partial Differential
Equations of Second Order. " Vol. 1, Polish
Scientific Publishers (1971).
20. Hess, J. L., "Higher order numerical
solution of the integral equation for the
twodimensional Neumann problem. " J. Comp.
Meth. Appld. Mech. Eng., 2, pp . 1  15
(1973) .
__ . Wardle, L. J ., "An Introduction
Boundary Element Method. " Numerical
solution of partial differential equations,
J . Noye ea., NorthHolland, pp . 289  312
(1982) .
22. Hess, J. L., "The use of higher order
surface singularity distributions to obtain
improved potential flow solutions for the
twodimensional lifting airfoils. " J.
Comp. Meth. Appld. Mech. Eng., 5, pp. 11
35 (1975).
Breit, S. R., Newman, J. N. and Sclavounos,
P. D ., "A new generation of panel programs
for radiation diffraction problems . " Proc .
3rd Int. Conf. Beh. Off. Str. (BOSS),
Elsevier Sc. Pub. B.V., pp. 531544 (1985).
24. Schultz, W.W., "A complexvalued integral
method for free surfaces with intersecting
bodies. " Proc. Second International
Workshop on Water Waves and Floating
Bodies, Rep. No. AM8706, Univ. Bristol,
pp . 101  103 (1987) .
Forsythe, G. and Moler, C. B., "Computer
solution of linear algebraic equations. "
Chapter 9, PrenticeHall, Inc., Englewood,
Cliff ., N . J . (1967) .
06. Lamb, Sir H., "Hydrodynamics. " 6th
Printing, Dover Pub ., (1945) .
366
OCR for page 351
27. Cointe, R., Jami, A. and Molin, B.,
"Nonlinear impulsive problems." Proc.
Second Int. Workshop on Water Waves and
Floating Bodies, Rep. No. AM8706, Univ.
Bristol, pp. 1316 (1987).
28. Biesel, F., "Les Appareils Generateures de
houle en Laboratoire." La Houille Blanche,
6, pp. 147165 (1951).
29. Roache, P.,"Computational Fluid Dynamics."
Hermosa Pub. (1972).
30. Baker, J. R., Meiron, D. I. and Orszag,
S.A., "Applications of a generalized vortex
method to nonlinear freesurface flows."
Proc. 3rd Int. Conf. Num. Ship Hydro.,
Paris, pp. 179191 (1981).
31. New, A. L., McIver, P. and Peregrine, D.
H., "Computation of overturning waves." J.
Fluid Mech., 150, pp. 233251 (1985).
32. Telste, J. C., "Calculation of fluid motion
resulting from large amplitude forced heave
motion of a twodimensional cylinder on a
free surface." Proc. 4th Int. Conf. Num.
Ship Hydro., National Academy of Sciences,
pp. 8193 (1985).
33. Orlanski, J., "A simple boundary condition
for unbounded hyperbolic flows." J. Comp.
Phy., 21, pp. 251269 (1976).
`r. Chan, R.KC. and Chan, F.WK., "Numerical
solution of transient and steady free
surface flows about a ship of generalized
hull shape." Proc. 13th ONR Symp. Naval
Hydro., Tokyo, pp. 257280 (1980).
35. Wu, DM. and Wu, T.Y., "Three dimensional
nonlinear long waves due to a moving
surface pressure." Proc. 14th ONR Symp.
Naval Hydro., National Academy Press, pp.
103125 (1982).
40. Kyozuka, Y., "Experimental study on second
order forces acting on cylindrical bodies
in waves." Proc. 14th Symp. Naval Hydro.,
Tokyo, pp. 319382 (1982).
41. Muggeridge, D.B. and Murray, J.J.,
"Calibration of a 58 m. wave flume." Can.
J. Civil Eng., 8, 4, pp. 449455 (1981).
42. Himeno, Y., "Prediction of ship roll
damping a state of the art." University of
Michigan Dept. of Naval Arch. Rep. No. 239,
65 pp (1981).
Z
UNDISTURBED
FREE SURFACE
aDF r— —%\
~ aDF
D
~DD
aDB
Fig. 1 Definition diagram.
~ F  ~  T HCTHO0. t'T" .. O
· PRC3CHS MCTH=, t/T. ~ O. O
~ T  ~rIC". t/T  . O. I~ O
0.00 1.00
aDc2
~5. Ahlberg, J., Nilson, E. and Walsh, H., "The .
theory of splines and their applications." (a) Free surface elevat~ons
Academic Press, New York. ( 1967 ) .
. Han, F. S. and Stansby, P.K., "On the
application of the boundary element method
to twodimensional free surface
interactions with bodies." Proc. Second
International Workshop on Water Waves and
Floating Bodies, Report. No. AM8706,
Univ. of Bristol, pp. 39 42 ( 1987 ) .
38. Dold, J.W. and Peregrine, D.H., "An
efficient boundaryintegral method for
steep unsteady water waves." In Numerical
Methods for Fluid Dynamics II.Ed.K.W.
Morton and M.J. Baines, Oxford Univ.
Press., pp. 671 679 ( 1986 ) .
39. Miche, R., "Mouvements ondulatoires de la
mer en profondeur croissante ou
decroissante." Annales des Ponts et
Chausses (1944).
~ PR£SENT METHOD,
o ~ THEORET ~ CRL,
_ _ _ ~ ~ ~ ~ ~
0.00 2.00 4.00 6.00 8.00 10.00
t/T
(b) Evolution of the free surface in time at the center of the domain ~ = 3.5)
Propagation of a linear steady
progessive wave: (L/A = 7,
d/) = 0.4, ~XF/) = 1/25 and
At/T = 1/40).
367
OCR for page 351
x present method, x/A = 0.26
theoretical, =/A = 0.26
present method, x/)t—0.50
theoretical, :~/A = 0.50
_~
, . . .
t.00 2.0~ 3.~0
0.00
x present method, x/) = 0.74
o .~. theoretical, z/) = 0.74
_ ~ present method, ~/A = 0.98
~ theoretical, :z /) = 0.98
o
0
"L'
_ ·~
o'
t/T
I ~_~
o
o
, , ~ ~ ~ , ~
0.00 1.00 2.00 3.00 ~ OQ 5.~0
t/T
1. 00 2. 00 3. 00
Fig. 3 The evolution in time of a free surface elevation started from rest.
o
_ 0
_ ~
_ o~
o
~ =
_ I ~ J I I
0.0
~ 0.2s
x O.SO
0.7s
O. 90
1. 00
\~
//~'
f::
0.00 0.2s 0.50 0.75 1.00 1.25 1.50 1.75 . 2.00
z/)
Fig. 4 Free surface elevations for different values of a:
(L2A, d/A0.5, H/A'O.1, AxF/A ~ 1/24 and ~t/T:1/40).
368
OCR for page 351
Fig. 5
Of (a)ati= 12(~/) = 0.48)
_~
o~
lo
o—
X PRESENT METHOD
MICHE
~ STOKES 2 Ad OROEfl
tI) RIRY
an
_
to
a,
q 1 1
4.se 4.7s s.oo
t/T
boo
5.25 5.50
Comparison of theoretical and numerical time histories of wave
elevation.
Airy wave pot.on ~Dc1
Stokes 2nd order potion dice

¢=o~it~l~M ~AM/
o
o
i I ~ ~ I ~ ~ ~ ~
0.00 1.00 2.00 3.00 4.OQ 5.00 6.00 7.00 8.00 9.00
t/T
To
° (b)ati= 36(~/) = 1.48)
A ADAM
. . .
0.00 1.00 2.~O 3.00 4.00 5.00 6.00 7.00 8.00 9.00
t/r
Fig. 6 The evolution in time of free surface elevations induced by different
excitation potentials.
369
OCR for page 351
~ incident WAVE
_ ~
h: $'
Fig. 7 Comparison of experiment and theory for w~b/2g 2 O. 54 and H/,'`
0.028: (T1~14.5 sec and T2~21.0 see).
370
OCR for page 351
~ incident .~
 16 16 1~ 1~ Ik. ^ ^ ^ ~
t [ "
hew motlen
S . ~ 16 IL~ 1~ ~ 1~. ^ ~ ~ ~
t ~ a ~
Fig. 8 Comparison of theory and experiment for 07b/2g ~ 0.54 and H/) = 0.028
when roll viscous damping is included in the simulation:
(T1~14.5 sac and T2~21.0 see).
371
OCR for page 351
DISCUSSION
by R. Cointe
I would like to congratulate the authors
for a very interesting and extensive work. We
are now working on a similar problem and we
find it very difficult to accurately compute
the free motions of the body using finite
differencesin time to evaluate the 3~/at
term in Bernoullie's equation. An alternative
method consists in evaluating this term by
solving the corresponding boundary integral
equation. Would the authors discuss their ex
perience on this problem?
Author's Reply
The work presented by Dr. Cointe at this
conference[A1], addresses several problems
which are either identical or analogous to the
ones we have addressed in our paper. The main
differences between the two works result from
the assumed directions of approach, with the
"numerical flume" approach taken by Dr.
Cointe, and "open water" approach chosen by
us. As a result certain comparisons between
the achieved results may be instructive.
In particular I believe, it is worthwhile
to notice that the wave damping condition
employed in [Al] at the down stream boundary,
is tuned to the dominant wave frequency, and
our application of the radiation condition
depends on the use of an appropriate wave
celerity.
It could also be interesting to see how
an introduction of the compatibility of the
initial condition on the exciting boundary and
free surface, which we impose through a
modulation function, would affect the
performance of Dr. Cointe's algorithms.
Referring directly to Dr. Cointe's
question, we could not avoid the necessity to
solve the boundary value problem twice in
order to advance the simulation of motion of a
floating body in time. The technique used is
described in section 5.3, and it relies on
using a central difference scheme to obtain
time derivatives of the potential on the body
surface.
We very __ _~
question and wish him success in finalizing
the development of his model.
milch ~nnrmm; ~~  Do Hi of  ~=
[Al] R. Cointe, "Nonlinear Simulation of
Transient Free Surface Flows", 5th
International Conference on Numerical
ShipHydrodynamics,Hiroshima, 1989.
DISCUSSION
by R.C. Ertekin
The authors should be commended for a very
careful and detailed analysis carried out in
their paper. It was a pleasure to read it.
By using the plane speed obtained from the
potential for Airy Waves, you appear to be
assuming, at least implicitly, that downstream
waves are linear. And by setting c'= ac and
varying a, aren't you basically determining
the value of the phase speed which you would
have if you numerically calculated it?
We have recently solved the problem of
diffraction of nonlinear waves by 2
dimensional submerged objects by using the
BEM. The results will be presented in the
next OMAE conference in Houston (1990, 9th
meeting). We have used CrankNicholson scheme
in time stepping and found no sawtooth waves
that you have seen in your results by using
the ABM method. We have not filtered our
results. This makes we believe that your time
stepping algorithm is causing the problem. It
is not very unheard of that a multistep method
is associated with such difficulties (Burden &
Faires, "Numerical Analysis", 1985, Prindle &
Weber). One way of handling this difficulty
is by using adaptive time stepping in the
algorithm.
I would like to also point out that
Shapiro, R. ("Linear filtering", Mats.
Comput., Vol. 29, 1975) presented filtering
formulas for any number of points earlier.
Finally, could you explain why sway forces
could not have been repeated in the
experiments?
Author's Reply
The comments by Professor Ertekin are
appreciated. considering the formulation
presented in the paper, the upstream boundary
can be used as a permeable boundary for
oncoming waves which satisfy fully the non
linear free surface conditions. However,
applications on that boundary, of wave
potentials which do not conform to the
conditions, are equivalent to wave excitations
by a wave maker. As much as the kinematics of
a wave maker board does not imply the
linearity of the generated waves, since they
satisfy the nonlinear free surface
conditions, the application of an Airy wave
potential on the upstream boundary does not
imply the linearity of the wave generated in
the control fluid domain. A wave excitation
372
OCR for page 351
of this type, by a wave potential which in
cludes a wave celerity c in its definition,
makes it convenient to use the same celerity
in the radiation condition. This procedure
appeared to be effective in the performed
simulations. Otherwise the celerity would
have to be determined numerically from the
solution in the control fluid domain.
The source of the sawtooth instability in
the wave propagation simulation is not
entirely clear. In our simulations it was
related to the use of the nonlinear free
373
surface conditions and to wave steepness, and
therefore the accuracy of the computation of
the spatial derivatives in the free surface
conditions seems to be the probable source of
the instability.
Finally, the measurements of sway forces
were repeated in all the tests carried out
twice. However the peaktopeak values from
those measurements showed a variation of at
least +4% between the original and repeated
test records. The cause of these variations
was not fully determined.
OCR for page 351