Fully nonlinear water wave computations using a desigularized Euler-Lagrange time-domain approach

ECHWSCfE tVESÇ1JT

'ahcí;rum 'ioor

&rehf Meke!wog 2, 2823 CD Oolft iaL 016 786873 Fa4. O1 7813

L

As published in Nonlinear Water Wave Interaction, Advances

in

Fluid Mechanics Series, WIT Press, 1999, pp. 1-58.

Fully nonlinear water wave computations using

a desingularized Euler-Lagrange time-domain

approach

Robert F. Beck

Department of Naval Architecture and Marine Engineering,

University of Michigan, Ann Arbor, MI 48109-2 145 USA

e-mail: rbeck@ engin. u,nich.edu

Abstract

The use of a desinguiarized Euler-Lagrange time stepping approach to solve fully

nonlinear water wave problems is reviewed. In the desingularized approach, the

mixed boundary value problem that arises at each time step in the Euler-Lagrange

method is solved using singularities that are placed outside the fluid domain

-inside the body and above the free surface. The desingulanzauon allows the normally continuous singularity distribution used in panel methods to be replaced

by isolated Rankine sources with the corresponding reduction in computational

complexity and little loss in accuracy.

Examples are given of the use of the method in both two and three

dimensions. Comparisons are made to analytic solutions, other nonlinear

computations, and experiments. A two-dimensional numerical wave tank is used to study highly nonlinear shallow water waves, the introduction of waves into a computational domain by modìfying the upstream boundary condition, and the

prevention of wave breaking using a local absorbing patch.

The diffraction of incident, shallow water waves by a circular cylinder are

investigated, and compared with other computations and experiments. The

extensively investigated mathematical Wigley hull form is used to demonstrate computations of wave resistance and calm water wave amplitudes along the hull. added mass and damping due to forced heave and pitch motions, and the exciting

i Introduction

As the state of the art in the design of_{ships and offshore structures has}

progressed, the demand for improved analysis _{tools has also increased.}

Marine hydrodynamic computations _{are very difficult because in any}

particular problem gravity waves, cavitation, surface tension, and

viscosity can all be important. Even _{the reduction to potential flow still}

leaves significant challenges because_{of the nonlinear nature of the free} surface boundary condition and the complex geometries of typical marine vehicles. This article will discuss recent _{progress in solving} fully nonlinear water wave problems.

We shall assume throughout that the fluid is inviscid, the density is

constant, and that the flow can be started from _{rest so that potential}

flow is applicable. The potential flow model is a major component of

the complete solution

including viscosity. Typical marine

hydrodynamics problems involve Reynolds _{numbers of the order of}

iû to lO so that viscous effects

_{are confined to the boundary layer}

arid viscous wake regions. While coupling between the flows in the

viscous region and the irrotational region are important, the majority of

wave effects are in the irrotational flow region. Since _{the potential flow}

problem is an order of magnitude easier to solve, it has been the basis

for much of the research.

_{At the present time, neither the high}

Reynolds number viscous flow around the body, nor the nonlinear waves can be adequately computed for all cases. Much more research needs to be done before a unified model of the entire flow field is available.

Fully nonlinear potential flow computations _{can be performed}

using a variety of methods. For steady _{forward speed problems, an}

iteratìve technique can be developed in _{which a series of linearized}

boundary value problems based

_{on the solution to the previous}

iteration and satisfied on the deformed free surface and exact wetted

surface of the body are solved. These _{iterative techniques have been}

successfully used by Jensen, et al. [1], Kim_{& Lucas [2] [3], Raven [4]}

and more recently Scullen & Tuck _{[5]. Bai, et al. [6] also used an}

iterative procedure but employed a localized finite-element method to

solve the boundary value problem at each _{iteration. For steady forward}

nonlinear solution using less computer time than a time-stepping method.

Longuet-Higgins & Cokelet [7] first introduced the mixed

Euler-Lagrange method to study two-dimensional fully nonlinear water waves near breaking. Their approach was a time-stepping procedure that

requires two major tasks at each time step. In the Lagrange phase,

individual nodes on the free surface are tracked by integrating the nonlinear free surface boundary conditions with respect to time; in the Euler phase, a mixed boundary value problem is solved to obtain the fluid velocities needed in the integration of the free surface boundary

conditions. Variations of this method have been applied to a wide

variety of two- and three-dimensional problems. The variations among the many published results consist primarily of different time-stepping methods and techniques to solve the resulting mixed boundary value

problem at each time step. Among the researchers who applied the

method to two-dimensional problems are Faltinsen [8], Vinje & Brevig [9], Baker, et al. [10], Grosenbaugh & Yeung [Il), Comte, el al. [12],

Sen f13], and Greaves, et

al. [14].

In three-dimensions the

computations become much more difficult because of the large

number of unknowns that are required. Results have been obtained by a number of researchers including Lin, et al. [15], Dornrnermuth & Yue [16], Kang & Gong [17], Zhou & Gu [18], Chan & Causal [19], Celebi, et al. [20], Xu & Yue [21], DLMascio, et al. [22], Wu & Taylor [23] and Wu, et al. [24], and the work to be reported on in this article done at the University of Michigan [25, 28, 29, 31, 32, 34, 36, 37, 38, 39, 40, 44, 45, 46].

The successful implementation of an Euler-Lagrange method requires a numerically stable time stepping scheme and a fast and accurate method to solve the mixed boundary value problem that results at each time step. In their original work, Longuet-Higgins & Cokelet [7] encountered a sawtooth instability of the free surface and

employed

a smoothing

technique to suppress its growth.

Dommermuth & Yue [16] also found

instabilities in solving axisymmetric problems. They postulated that the cause of the high wave number instability was the concentration of Lagrangian markers in the region of higher gradients that cause the local Courant stability

condition to be violated as the wave steepens. _{They eliminated the}

instability by regridding the free surface after each time step.

Park & Troesch [25] have investigated in detail the stability _of

time stepping for a variety of two- and three-dimensional _{problems and}

reached a number of conclusions.

_{They found that the stability}

depends upon the geometry of the specific problem, the closure _at

infinity, and the time integration method. They also note that

three-dimensional problems tend to be more stable than two-dimensional

problems. Explicit Euler ti me-stepping schemes are unconditionally

unstable, while implicit-like and implicit Euler schemes and

fourth-order Runge-Kutta schemes are conditionally stable. They defined a

local stability index Eg(t)2/x which if exceeded leads to an unstable

solution. The limit of the local stability index depends _{upon the}

geometry and the time-marching scheme used.

The solution of the mixed boundary value problem that results at

each time step in the Euler-Lagrange method can be solved by many

different approaches; in the work done at the University of Michigan,

we have used the desingularized boundary integral method. Similar_to

conventional boundary integral methods, it reformulates the boundary

value problem into a boundary integral equation. _{The integrals are}

discretized by the distribution of fundamental singularities _{over the}

"integration surface" and satisfying the boundary conditions _{on the}

"control surface." _{In traditional panel methods, the integration and}

control surfaces coincide so that the resulting kernels are singular. In

the desingularized method, the integration surface is _{separated from the}

control surface by moving it outside the fluid domain, resulting in a

"desingularized" kernel that no longer contains _{a singularity.}

Techniques to represent bodies in potential flow using singularities outside the fluid domain have been in use since the earliest days of

hydrodynamic computations. The Rankine ovoid is probably the _best

known example. A source and sink of equal strength are combined

with a uniform stream to yield a closed stream surface_{surrounding the}

two singularities. The combination thus represents the flow about _an

elongated body in a free stream.

_{The Rankine ovoid can be}

generalized to the flow about an arbitrary axisymmetric _{body by}

distributing a greater number of singularities along the _{axis aligned}

calculation of pressure distributions on airship hulls. The strength of

the source distribution is determined by the kinematic boundary

condition on the body surface. More recently, Webster [27] used

triangular patches of sources with linearly distributed strengths

"submerged within the body surface to study the steady flow past

arbitrary three-dimensional bodies in an infinite fluid. Schultz &

Hong [28] showed the effectiveness and accuracy of the desingularized

method for two-dimensional problems. Convergence rates and error

limits for simple three-dimensional flows, including a source-sink pair

traveling below a free surface, were given by Cao, et al. [29].

All panel methods, whether desingularized or not, eventually require the solution of a system of linear algebraic equations in order to determine the unknown singularity strength, given the boundary

conditions. In first order panel methods, the matrix to be inverted

expresses the influence of each panel or singularity on each collocation point where the boundary conditions are to be satisfied. Thus, for a uniquely determined system with N panels or singularities, there must be N collocation points and the resulting influence matrix has N2 elements. In traditional panel methods, care must be taken because the kernel is singular, but this does lead to large diagonal elements in the influence matrix. In desingularized methods, such care is not needed

because the kernel

is

no longer singular.

Simple numerical quadratures can be used in computing the influence matrix with a

corresponding reduction in computational effort, particularly by

avoiding transcendental functions. The desingularization does lead to a loss of diagonal dominance, which can cause numerical difficulties in the inversion process. In practice, however, if the desingularization is "not too large", the loss of diagonal dominance is not a problem. We have found that because of the desingularization, source distributions over panels are not necessary, and the panels can be replaced by simple

isolated sources Higher-order singularities such as dipoles can easily

be incorporated.

The isolated Rankine sources allow direct

computation of the influence matrix since in three dimensions it is simply a hr potential. The use of complex numerical quadratures to

integrate over panels is eliminated with an equivalent saving in

In the next section, the fully nonlinear boundary value problem

that must be solved to predict the hydrodynamic loads on ships and

offshore structures in the time domain will be set up. The following

section will discuss the desingularized method developed at the

University of Michigan using isolated Rankine

_{sources, and its}

numerical solution. Finally, numerical results for both two- and

three-dimensional problems will be presented.

2 Fully nonlinear problem formulation

The fluid is assumed to be inviscid, homogeneous and incompressible.

Surface tension on the free surface is neglected. The problem is started

from rest so that the flow remains irrotational; this implies the existence

of a velocity potential such that the fluid velocity is given by its

gradient.

A right-hand coordinate system, Oxyz, is translating in the

negative x _{direction relative to a space fixed frame. The}

time-dependent velocity of translation is given by U0(t) . The Oxyz axis

system is oriented such that the z = O plane corresponds to the calm water level with z positive upwards. The total velocity potential of the

flow can then be expressed as

U0 (t)x + Ø(x. y, z, t) ₍₁₎

where (x,y,z,t) is the perturbation potential. Both and must

satisfy the Laplace equation in the fluid domain:

V2D-O ₍₂₎

Boundary conditions must be applied on all surfaces surrounding

the fluid domain: the free surface (SF) , the body surface (SH) , the

bottom (SB) and the surrounding surface at infinity (S)

. The

kinematic body boundary condition is applied on the instantaneous

position of the body wetted surface:

where n = (n1, n2, n3) is the unit normal vector into the surface (out

of the fluid domain) and VH is the velocity relative to the moving

coordinate system of a point on the body surface including rotational

effects. The subscripts 1,2,3 refer to the x, y, and z axis directions

respectively.

There is also a kinematic condition applied on the

bottom:

=-U0(t)n1+Vgn on SB

(4)

where _B is the velocity of the bottom relative to the Oxyz system.

For an infinitely deep ocean, equation (4) becomes

as

Z-4-

(5)

Finite depth will increase the computational time because of the

additional unknowns necessary

to meet the bottom boundary

condition, but there is no increase in computational difficulty. Flat

bottoms are very easy to compute since an image system can be added

to the influence matrix for relatively small computational cost. This

contrasts with the typical Green function approach where a finite depth Green function is significantly harder to compute than an infinite

depth Green function and a nonflat bottom can not in general be

accommodated.

On the instantaneous free surface both the kinematic and dynamic

conditions must be satisfied. The kinematic condition is

& on SF (6)

where z = (x, y, t) is the free surface elevation. The dynamic

condition requires that the pressure everywhere on the free surface

equals the ambient pressure, _a The ambient pressure is assumed

known and may be a function of space and time; normally it would be

set equal to zero. Using Bernoullis equation, the dynamic condition

ax on SF (7)

where p is the fluid density and g the gravitational acceleration. Because we are solving an initial value problem with no incident

waves, the fluid disturbance must vanish at infinity

as

R=jx2+y2

400 ₍₈₎

In addition, the initial values of the _{potentiai and free surface elevation}

must be specified such that

=O tO

influiddomain

(9)

fl0 tO

In the Euler-Lagrange approach, the

_{mixed boundary value}

problem that must be solved at each time step is obtained by time marching (or integrating) the boundary conditions from the previous

time step. _{In our computations, a fourth-order Runge-Kutta method is}

typically used for the time marching. _{At each step, the value of the}

potential is given on the free surface and the value of the normal derivative of the potential is known on the body surface, the bottom

surface, and the surrounding surfaces _{at infinity. The free surface}

potential and elevation are determined by integrating the free surface

boundary conditions (6) and (7). _{The normal velocity on the body}

surface is prescribed for the forced motion _{problem. In the case of a}

freely floating body, the equations of motion must be integrated with

respect to time in a manner similar to the free surface_{conditions. In}

computations to date, the bottom body condition _{has been satisfied}

using an appropriate image system for finite depth or an open

boundary in the case of infinite depth. The _{surrounding surfaces at}

infinity typically have prescribed values of the normal velocity, either

zero for walls or a precomputed value in order _{to introduce incident}

On the free surface, the kinematic condition (6) is used to time step the free surface elevation and the dynamic condition (7) is used to update the potential. Many different approaches are possible to time march the free surface boundary conditions. The most common is a material node approach in which the nodes or collocation points follow

the individual fluid particles. Another technique is to prescribe the

horizontal movement of the node but allow the node to follow the vertical displacement of the free surface. The prescribed movement may be zero such that the node locations remain fixed in the x-y plane. Depending on the problem, one of the techniques may be easier to apply than the others.

It is convenient to rewrite the free surface boundary conditions, equations (6) and (7), in terms of the time derivative of a point moving with a prescribed velocity y relative to the Oxyz coordinate system.

By adding y Vr1 to both sides of (6) and v V to both sides of (7)

and after some algebraic manipulation, the kinematic and dynamic conditions can be put in the form

on _SF (10) & az ax and on where a

E+vV

is the time derivative following the moving node. This is similar to the

usual material derivative of fluid mechanics except that the velocity is

the given y rather than the fluid velocity.

If y is set equal to (U(t)V(t)iJ

the node follows a prescribed

path with velocity (U(t), V(t)) in the x-y plane and moves vertically

with the free surface. Setting y=

_(oo.1J

_{results in the x-y}

locations of the nodes remaining fixed in_{the Oxyz coordinate system}

and equations (10) and (11) reduce to:

on SF (13) at az and and aî1a on _{SF (14)}

atar

In the material node approach, the velocity is set equal to the fluid

velocity such that y =U0(t)i +

_{V, resulting in}

DXF(t)

=U0(t)iV

Dt

=grl + - V

-where XF

= _{(XF(t),yF(t),ZF(t))}

_{is the position vector of a fluid}

particle on the free surface and

= --+v V

Diat

is the usual material derivative.

The form of the free surface boundary _{conditions given above}

allows the value of the elevatìon and potential to be stepped forward in

time. The left hand sides of equations (10) _{- (16) are the derivatives}

with respect to time of the potential and wave elevation moving with a node. The quantities on the right hand side are all known at each time

after solving the boundary value problem and the wave elevation is known. The difficulty is the gradient of the free surface elevation (Vii in equations (10) and (13)) that must be evaluated numerically. This

leads

to increased computer time and numerical inaccuracies.

However, this term is only needed in the prescribed horizontal node

movement approach; in the

material node approach no extra

derivatives need to be evaluated. This probably explains why the

material node approach is most often used. On the other hand, the

unconstrained nature of material nodes makes regridding in three dimensions difficult and one must always be concerned that the nodes

do not penetrate the body surface between time steps. In zero forward

speed problems, material nodes or fixed nodes seem to be the most

appropriate. In problems with forward speed, the material node

approach has difficulties near the body because nodes tend to pile up near the bow and stern stagnation regions. The prescribed horizontal node movement approach does not have this difficulty since the node

movement is constrained. An appropriate choice of y is one that

approximates Vb reasonably well. In this case, the contribution of the Vii term to the right hand side of (10) will be small, and numerical inaccuracies will be minimized. Consequently, fast, simple numerical derivatives can be used to evaluate the Vil term.

At each time step a mixed boundary value problem must be

solved; the potential is given on the free surface and the normal

derivative of the potential

is known on the body surface and the

bottom. The boundary conditions on the surrounding surfaces at

infinity are usually of the Neumann type, but depending on the

problem they may be be open or Dirichiet.

In

terms of a

desingularized source distribution, the potential at any point in the fluid domain is given by

=

(x) fJa(x)

Q

where ) is the integration surface outside the fluid domain.

Applying the boundary conditions, the integral equations that must be solved to determine the unknown source strengths a(x) are (18)

5Jc(x)

d14F(xC)

Q _1xcxsI

'

dç2=:X(Xc)

Q

_ncI1cXsI

E

where x _{= a point on the integration surface, f}

x _{= a point on the real boundaiy}

= the given potential value on the free surface_at

= free surface on which _{F IS knOwn}

x = the given normal velocity on the boundary at

surface on which _{x is known.}

2.1 Hydrodynamic forces

The hydrodynamic forces acting on the body are found by integrating the pressure over the instantaneous wetted surface. The generalized

force acting on the body in the th _{direction of the body axis system is}

thus given by:

=f5 Pn ds (21)

SR

where ni _{is the generalized unit normal into the hull defined as}

(n1,n2,n3)= n ₍₂₂₎

(n4,n5,n6)= rxn

n = unit normal to body surface (out of fluid) E rd

or

r = (i,y,i)

O57 = body axis system

and j = 1,2,3 corresponds to the directions of the

O7

axis

respectively.

The pressure in the moving coordinate system is given by

Bernoulli's equation:

= - - U0 (t) gz

--U

p & °

where is the time derivation of the potential following a moving

node on the body and y

is the velocity of the node relative to the Oxyz system.

3 Numerical techniques

In the usual manner, the integral equations at each time step (19) and (20) may be solved by a collocation method in which the integrals are discretized to form a system of linear equations with the boundary conditions satisfied at collocation or node points on the boundaries of

the fluid domain. In the desingularized method, the sources are

distributed on an integration surface outside the fluid domain. Since

the integration surface never coincides with the collocation points, the resulting kernels are nonsingular.

As discussed in Cao, et al. [29], the desingularization allows the replacement of the surface distribution of sources by a distribution of isolated Rankine sources with little loss in accuracy as long as the desingularization is sufficient. This greatly reduces the complexity of

a desingu1ariz isolated Source points

/z

a _t a

t

a a a * a s t

/

Ç - given on the free surface

- given on the body

Figure 1: Schematic of isolated source and node locations the form of the influence coefficients that make up the elements of the

kernel matrix.

As shown in figure

_{1, the collocation points are}

distributed on the free surface and body surface; the isolated _sources

are distributed a small distance above each of the free surface nodes

and inside the body. _{The nondimensional desingularized distance is}

given by

Ld=td(Dmf' (24)

where td and u are constants to be chosen by the user and Dm is a measure of the local mesh size (typically the local panel area in

three-dimensional problems). Cao, et al. [29] found values of

td = 1.0 and

u = .5 to be about optimum. The accuracy and convergence of the solutions are relatively insensitive to the choices of td and u.

As an example of the numerics involved in desingularized

computations, consider the case of just a body and free surface as

shown in figure 1. _{Further source distributions and collocation points} can easily be added for other domains such as surfaces at infinity,

or

NF X i=1 f

NFFa

i

NBBa

F j=i

a,

_\1C1

-XJ

j=1 an, (i=l...NB)

where x

and x

are the node points on the free surface and body

surface respectively.

N = NF + NB being the total number of

unknowns. Equations (26) can also be written in matrix form as

AE=B

(27)

w

sources above the free surface and inside the body surface to construct the velocity potential, equation (18) reduces to a simple summation of the influence of each source:

NF

cF

NB

J=IiX_xJ

J=11XJ

where NF is the number of the isolated sources above the free surface

x is the location of theJth free surface source and

is the strength of the th free surface

source at ¿

Similarly, NB, x

and c

are the number, location and strength of

the sources inside the body surface.

The integral equations (19) and (20) are satisfied at the nodes on the free surface and body surface such that

cF

NB

aB

¿i_X! i=IXl_IJ!

(i=l...NF) = (25) (26)

(A(' 1)

_{A'2 '(z0) )}

_{(B0) '} I1A21 _{A2)JI,Z(2)f_1B(2))} where

A''

Jx xJ

i=l...NF, j=L..NF

F B

i=l...NF, i(NF+1)...N

xc xs

f a i' LB F

tci

XS. "I J J

i=(NF+I)...N, j=L..N

(28)

A22)=[

xj

i=(NF+1)...N, j=(NF+1)...N i=L..NF

(2)B

_i=(NF1)...N

B=o(x)

i=I...NF

B2=x(xfl

i=(NF+1)...N

Equations (27) are a system of N _{= (N} + _{NB) equations for the}

unknown source strengths

_Z

. Once

I

_{are known, the velocity}

potential and fluid velocity can _{be analyticaliy calculated anywhere in}

3.1 Domain decomposition

Equations (27) can be solved by direct, iterative, or multipole methods depending on the size of the matrix. For most of the two-dimensional problems considered in this paper, we use a LU decomposition solver.

For the three-dimensional problems, an iterative solver, such as GMRES

(Saad & Schultz [30]), or a multipole method (Scorpio, et al. [31]) are

used.

For problems where different

grid

spacings and/or

desingularization distances are used in different parts of the domain, solving equation (27) as a whole system may not be efficient because

the A1

matrix is poorly conditioned.

A domain-decomposition

technique may often work better. In this technique

(1) and X2 are

determined separately through an iterative procedure.

Rewrite

equation (28) as

1(1) is determined by solving (29) with (2) known (or guessed) then

(1) is substituted into (30) to determine (2) _{The new is then}

used in (29) to update

E'

. The procedure repeats until both

equations are satisfied to the given error tolerance. The individual

matrices A(11 arid A(22) will have better conditioning because they

have relatively constant grid densities. For most water wave problems, the majority of the nodes are used to discretize the free surface. The

number of nodes on the body is often small enough to use direct

factorization.

The domain-decomposition technique

will be

computationally faster than solving the entire matrix at once as long as the number of iterations for (29) and (30) to converge are not too large.

After the are determined, the fluid velocity on the free surface can be computed analytically so that the right hand sides of the free surface conditions (10) and (11) can be computed. The free surface

elevation and the potential are then updated by a time stepping

procedure. In the free surface updating, if material nodes are used, the

spatial derivatives of the free surface elevation, V'ri , are not required.

A(11) .

(l)= B(1)- A'2) .

(2) (29)

The use of generalized nodes moving with prescribed velocity _y

requires the spatial derivatives of _{and central-differencing or cubic}

splines are used. _{The time-stepping of the free surface is done by}

using a fourth-order Runge-Kutta _scheme.

As given in (21), the hydrodyriamic forces acting on the body are

found by integrating the pressure over the body surface at each time

step.

_{The pressure may be evaluated by either (23a) or (23b)}

depending on whichever form is _{more convenient. In either case, the}

time derivative of the potential, &/St or 4Iat, must be determinedat

each time step. For forced body motion _{problems, the forces acting on}

the body do not enter the _{time-stepping procedure. Therefore, the}

calculation of the hydrodynamic forces can be post-processed and any

numerical time derivative, such as central differences, can be used.

However, for problems involving free _{body motions, the equations}

of motion must be integrated _{at each time step in order to determine}

the wetted surface and normal velocities for the next time step. A

critical component of the integration process is the accurate evaluation

of the hydrodynamic forces, hence E4/6t or _{It must be known at the}

present time step. The conventional technique is to compute 641& or

aiat using backward differencing,

_{but at times this might not be}

accurate enough. Beck, et al. [32] presented _{a technique that solves an}

auxiliary problem for aqiat _{directly. Only the right hand side of the}

original problem is altered in the _{auxiliary problem so that minimal}

additional computational effort is needed as long as the inverse of the

influence matrix is known. _{Unfortunately, for iterative solvers such} as

GMRES the inverse matrix is not computed and the problem must

essentially be resolved for each right _{hand side. Similar techniques}

have been used by other researchers such as Vinje & Brevig [9], Yeung

[33], Kang & Gong [17], and Comte, _{et al. [12].}

3.2 Preconditioning

An alternative to the domain decomposition cited above for problems

with different grid spacings and/or desingularizations distances is to use

preconditioning of the influence matrix, _A. . The condition, and

improved if a good approximation for the inverse of A1 is readily

available and not too expensive to compute.

Let A

be the

approximation to A' , then the iterative solver can be preconditioned

by solving

A1JAJk Sk = b1 (31)

for sk, where sk is related to the source strength by = A s For

the computations being done at the University of Michigan, Scorpio [34] developed an overlapping block preconditioner based on the work

of Nabors, et al. [35].

In this technìque, the inverse to

is

computed by solving overlapping subproblems using the lowest level boxes in the tree structure developed for the multipole algorithm described below. The sub-influence matrices are inverted and only the inverted coefficients which represent the influence of the box on itself

are saved and loaded into the appropriate locations in the A1 matrix.

The resulting matrix

is sparsely populated and contains

information about the local inverse problems. The work and storage

requirements for the preconditioner are of 0(N) 3.3 Multipole acceleration

Multipole acceleration is a method of accelerating the evaluation of the potential (or velocity) at N points due to N sources. Normally, this would be an order N2 process; multipole acceleration can reduce this to

order N or NiogN depending on the specific algorithms. As

discussed in Scorpio [34] or Scorpio, et al. [31], multipole acceleration

starts with multipole expansions that are developed from the solution to the Laplace equation in spherical coordinates:

i

(2 a

i i

ac

=0 . (32)

-Tr

_--)

_{r2 sine ao1}

1n8.)+

r2sin2S 2

The solution of this equation by separation of variables results in a series of spherical harmonic terms,

'=

n

lLr°-'-

( ₍₃₃₎

n=Om=n

where M and L

are the moments of the expansion and Y(8,)

are the spherical harmonics. Using the series solution _{to the Laplace}

equation, the far field potential due to a collection of near field sources

can be expressed in a multipole expansion. _{Consider a sphere of}

radius a containing the k sources with coordinates

(p1,a1 _), i = l,...,k and strengths y1 Let the center of a multipole

expansion be defined as the center of this sphere. We can then express

the potential due to the k sources atany point, P(r,8,p) , outside of

the sphere, in the form

O

M

Z

.!_ym(P)

_o+I n

n=Om=n I

where the coefficients of the expansionare given by

M =

ap0y_m(a)

1=1

Likewise, the near field potential due to a collection of far field sources

can be expressed in a local expansion. Suppose the k _{sources now lie}

outside of the sphere of radius a centered at the origin. The center of

a local expansion is defined by the center of the sphere inside which_all

of the evaluation points (P) lie. The potential due to these sources at

any point, P(r.8,p) ,inside of the sphere is given by

(P) ₍₃₆₎

where the coefficients are given by

L ₍₃₇₎

Several theorems are used in the multipole algorithm which involve

shifting the origins of multipole and local expansions and the

conversion of multipole expansions into local expansions. Greengard

[36] gives details of these theorems and their error bounds.

The underlying idea of multipole acceleration is that the potential due to a group of sources may be evaluated at some distant point x by first accumulating the influences into a multipole expansion and then evaluating the single expansion rather than having to compute

each of the individual influence functions.

It is the careful arrangement, or clustering, of sources into groups that leads to the NiogN efficiency. The fast multipole algorithm used by Scorpio, et al. [31] reduces the cost yet further, to order N , by a complementary arrangement of the evaluation points into groups so that accumulated multipole expansions may be transformed into local expansions

centered around each group which are then evaluated instead. More

specifically, the use of multipole and local expansions is orchestrated by a tree-stnictured hierarchy of source groups and evaluation point

groups. As shown in figure 2, multipole expansions for groups of

sources are accumulated from the leaves of the tree to the root, and

Figure 2: Multipole and local expansion manipulation

mbe

c2usz

p

du

local expansions

_{are distributed from the root to the leaves for}

evaluation at the collocation points. This is accomplished in order N

operations while maintaining

_{a uniform precision as shown by}

Greengard [36].

The hierarchical partitioning of _{the domain into a tree structure}

ensures that the potential is accurately approximated everywhere in the

problem domain. _{The partitioning is accomplished by successive}

divisions of the fluid domain into finer and finer cubes. Since the free surface lies more or less in the plane of the calm water level, Scorpio

[34] replaces the three-dimensional

_{tree structure with}

_its

two-dimensional equivalent in order to simplify some of the logic required in partitioning the free surface. As presented in figure 2, the largest box is referred to as the root, or level O , box. The smallest boxes (the

leaves) are designated level L . The number of levels is usually

selected so that the finest level boxes contain no more than some small fixed number of sources.

Setting up the

tree structure requires a

_{fair amount of}

computational overhead; thus, multipole acceleration is only profitable

for large numbers of unknowns. _{Since the majority of unknowns in}

desingularized water wave problems are the free surface sources,

Scorpio [34] only applies multipole _{acceleration to the influence of}

the free surface sources. When used inconjunction with a fast iterative

solver such as GMRES, the solution

_{can be arrived at in order N}

operations. Since the influence matrix is _{no longer needed, the N2}

storage requirement is also eliminated.

_{This makes a significant}

difference in the computer time and _{memory requirements when N}

becomes large. Scorpio [34] or Scorpio, et al. [31] give graphs of the

savings in computer time and

_{memory for simple mathematical}

boundary value problems

_{versus the number of unknown source}

strengths. On a linear

_{processor work station, the trade-off between}

using GMRES or the multipole acceleration occurs around 1000

unknowns. As expected, the CPU time and memory requirements grow

as N2 for the iterative solver and as N when using

_{the multipole}

expansions. _{For the problem of a source-sink pair below a =0 free}

surface with 5000 unknowns, Scorpio [34] finds the multipole method

is approximately 4 times faster and_{requires 10 times less storage. With}

speeded up by a factor of approximately 5. The memory requirements

with preconditioning are increased for both methods since the 0(N2) preconditioner matrix must be stored.

4 Numerical results

Cao, et al. [29] examined various aspects of the convergence and errors in using the desingularized method with a direct application of the Green theorem, a panel distribution of sources, and a distribution of

isolated sources. The findings for the case of a dipole below a = O

infinite flat plane are shown in figures 3 and 4. Figure 3 presents a

schematic of the location of the desingularized control points for the direct application of the Green theorem and the isolated source points in the indirect method. The desingularization distance is based on the

square root of the panel area times the constant . The integrals in

the direct method are evaluated using a 2x2 Gaussian quadrature.

X'

4.

4.

Figure 3: Schematic diagram of a dipole below a = O

infinite fiat plane (from Cao et al. [29])

Otec1 Meod Method

-Co1r Pi - Corvos P1

Nod. Poul - So,a Pi

z

4. 4. 4. 4. 4. 4. 4.

4.

In figure 4, the numerical results for _{/an on the z = O plane are} compared with the exact solution, which can be written down for this

simple problem. The results are presented _{as the RMS error for the two}

calculation methods versus the desingularization distance for three

different numbers of nodes. As expected, the _{error reduces for both}

methods as the number of nodes is increased. The direct method

initially has a smaller error, but for larger

_{desingularizations the}

indirect method using isolated _{sources performs better. In either case,}

there is a wide range over which the errors are relatively insensitive to

the desingularization distance. _{At zero desingularization distance, the}

u E

0 1 2 3 4 5 6

Figure 4: Effects of desingularization on the RMS

_{error -, direct}

application of the Green theorem;

_{-, indirect method}

using isolated sources; O, N = 231; N = 496; +,N = 861

isolated source solution blows up because the sources are distributed

directly above the nodes. Cao, et al. [29] did examine staggeredgrids

and recommended against them because staggering strategies are

difficult to choose and only marginal improvement was obtained.

They also examined the results using a desingularized panel

distribution of sources (similar to Webster [27]) rather than isolated

sources. The panel distribution did allow a larger range of 'd for

minimum error, but the optimum error was the same for both methods. However, the condition number of the system of equations to be solved using isolated sources is improved by a factor of 100 over a panel distribution. In addition, the panel distribution requires significantly

more computer time to formulate the influence matrix. The final

conclusion was that nonstaggered, desingularized, isolated Rankine sources were the preferred method of computation.

The results of the desingulanzed method have been compared to a

number of other published fully nonlinear computations. Figure 5,

from Beck [37], is one example of such a comparison. It shows the time history of the vertical force acting on a heaving three-dimensional cylinder in finite water depth equal to the radius. Also shown in the figure are the linear and nonlinear results computed by Dommermuth

& Yue [16] using entirely different numerical methods. The

agreement for the nonlinear results is very good; the only difference is a slight phase shift on the upward part of the cycle. Because the actual numerical results were not available for comparison, the phase shift could be the result of the xerox copy used to acquire the Dornrnermuth

& Yue results.

The force time histories clearly show the nonlinear effects due to

the extremely large amplitude of motion of half the draft. The linear

results have zero mean and are basically sinusoidal.

The fully nonlinear results have force amplitudes approximately the same as the

linear results during the down cycle (or peak positive force), but

significantly larger suction forces when the bottom is near the free surface. Consequently, there is a mean shift in the nonlinear results. This is one of the advantages of nonlinear computations. The mean

shifts and higher-order forces are determined naturally from the

4.0

H'H

Figure 5: Time history of the vertical force for a heaving

three-dimensional cylinder in finite water depth equal to the radius (R = l-O, T= O-5), heave

amplitude aIT = 0-5. _Motion:

- - nonlinear calculation (present method);

-.. non-linear calculation (Dornmerrnuth

& Yue [16]);

- - linear calculation

(Dommermuth & Yue [16]) (after Lee _[38])

-4.0

_i

_i

0.0

1.0

2.0

_3.0

4.1 Two-dimensional numerical wave tank

A two-dimensional numerical wave tank has been used for many

different calculations. The tank has a wavemaker at one end and a

beach or end wall at the other end. We have used various types of wavemakers including pneumatic, a wedge plunger, a paddle, or a piston in shallow water. The best type depends on the problem and the physical system that is being modeled. Cao, et al. [391 found that the

most effective beach appears to be an energy absorbing pressure

distribution that is added to the dynamic free surface boundary

condition. The strength of the distribution builds up gradually from the start of the beach and is proportional to the value of the potential or normal velocity on the free surface.

As an example of the application of the method to highly

nonlinear shallow water waves, Cao, et al. [40] studied solitary waves in

shallow water generated by moving disturbances. All computations

were in the time domain started from rest. Figure 6 shows the free

surface elevation at different times for a cylinder of radius r=.15

traveling at depth d = .6 below the free surface. The water has a depth

of h = 1.0 and the Froude number based on depth,

Fn = uo/ji

, is

1.0. The generation of runaway solitons ahead of the disturbance are

clearly visible. A moving computational window was used for the

calculations; the cylinder is located underneath the arrow for each

realization of the free surface.

For the upstream boundary, the elevation and potential on the free surface are precomputed using some extra upstream points outside the window that are convected by the

computed velocity from the already-solved source strength. On the

downstream side no special treatment is used since the method does not require spatial derivatives of the free surface elevation. Each time the

computational window is shifted, one node point is dropped. In figure

6, a nonphysical reflection from the open downstream boundary

appears to build up with time.

Using an open boundary on the downstream end of the computational window is expedient because it

saves computer time; it might or might not work well, however,

depending on the problem. If the computations in figure 6 were to

continue much longer, it would be necessary to have either a better downstream boundary condition or a longer computational window.

One of the difficulties with nonlinear computations to determine the seakeeping performance of ships is the incident waves. Normally,

seakeeping computations _{use an analytic form for the incident}

waves,

either linear or Stokes higher

_{order waves. For fully nonlinear}

computations, however, the incident _{waves must also be nonlinear and}

therefore need to be computed. In order to get a sufficiently _long

12.0

10.0

8.0

6.0

4.0

2.0

0.0

V-_

F=1.O

AW

Jw-I I I 0 20 40 60 80 100 x + Frt

Figure 6: Waves generated _{by a submerged circular cylinder traveling}

at a depth Froude number of 1.0. _{rIb = .15, d/r = 4,}

= .2, free surface nodes = 251, body nodes_{= 30 (from}

Cao, et al. [39])

2c0

160

120

time history, the wavemaker must be far upstream; this leads to a great deal of "wasted" computational effort in three-dimensional problems because of all the free surface nodes between the wavemaker and the ship that do not effect the flow around the ship. Scorpio et al. [31]

present a method to overcome this difficulty by using relatively

inexpensive computations in a two-dimensional wave tank to set up the input for a much smaller three-dimensional computational window

around the ship. As shown in figure 7, waves are made in a

two-dimensional tank in fixed coordinates X0-Y0. The computational

window containing the ship is then moving through the incident wave

field by prescribing 4,I?in on the vertical upstream boundary and 4,

and i

on the upstream free surface. The same boundary conditions

can be used on the downstream boundary of the computational

window, but for cases with a body present the downstream waves are not known a priori. A radiation condition or some type of numerical beach could be applied, but because of the forward speed an open boundary seems to work.

To verify the technique, the waves in the moving computational window without a body present were compared with the original fixed

coordinate waves. Since no body is present, the wave elevations, r, should be identical. Figure 8 shows the results of the computations for

Direction of window mouon Direction of

wave propagation

Figure 7: A two-dimensional wave tank with a moving computational window (from Scorpio, et al. [31])

-4 -o o 4 2 o -2 -4 -4 4 2 o -2 -4 -4

dizvcjoo or mc idcnt wv a 4 dut,c. of udo*_nù to..

\

= 00 Fr -Fr=0.27 -Fr =0 40 - - -2 _{s I} I o 2 3 'o o 2 3 4 S S a XoA

Figure 8: Incident waves in a moving computational _{window. Solid}

line - wave elevations in the computational _window. Dashed line - wave elevations in the stationary wave tank

(from Scorpio, et al. [31])

2 s 7 a a Io

s

2

o -2

fout different translation speeds. The incident waves have wavelength X and the speed of the computational window is defined by a Froude

number based on the wavelength, Fn = Uo/J Because the

downstream boundary condition is left open, the agreement while very good is not perfect. Specifying the downstream boundary in the same manner as the upstream boundary results in agreement to within graphical accuracy.

Another difficulty with fully nonlinear water wave computations is

breaking waves and spray. In linear or higher-order expansion

methods, the wave amplitudes can get unrealistically large without

causing numerical difficulties. In fully nonlinear computations,

however, wave breaking or spray at the body/free surface intersection will cause the computations to stop even though it may be confmed to

a small area such as the breaking bow wave.

To overcome this difficulty we have been working on the local suppression of wave

breaking and spray by absorbing energy using the free surface

boundary conditions.

As reported by Subramani, et

al.

[41], a

technique is proposed to absorb energy locally from waves that are about to break, thereby suppressing wave breaking. Two steps needed in the 'local absorbing patch" model are: 1) the detection of the likely occurrence of the wav&s breaking, and 2) determining the appropriate amount of local damping so as to render reasonably realistic waves in the post-breaking regime.

There are many proposed criteria for wave breaking, see for

example Griffin, et al. [42] or Wang, et al. [43]. Unfortunately, most

of these criteria are hard to invoke andlor are computationally intensive for use in Euler-Lagrange computations. For this reason, Subramani, et al. [41] developed a criterion based on the fact that the wave profile

has a sharp crest with infinite curvature (ic) near breaking and there is

a limiting steepness (assumed to be HA 1/12). Combining the above

with extensive testing of fully nonlinear computations of regular,

two-dimensional waves resulted in a proposed threshold of Kai = .35. For

values greater than this number, it is assumed that the waves are about

to break and a local wave damper is applied to the wave crest.

The actual damper is developed by applying an artificial pressure on the free surface whose phasing is such that it always takes energy

out of the waves. Thus, the dynamic free surface boundary condition (11) can be rewritten as:

}'d;nP

Uo(t)

on SF (38)

where _damp is the damping term. The damping is confined to a small region around the crest and in zero-speed computations has the form:

damp= y

u(x) IV2 sgu(Ian)

(39)

where y = constant coefficient that defines the magnitude of the

damping and the sgn function ensures the correct sign on the damping term. The shape function, u(x) , is chosen to ensure that the damper is a smoothly varying patch:

u(x) = .5{l + co4r(x - xo)/LoJ}

where x0 = location where Kai = .35 is exceeded

L0 = half-length, centered at x0 over which the damper acts. L0 is presently assumed to be a0/.35, a number chosen so that energy is extracted from approximately the crest portion of the wave between

zero-crossings.

Figure 9 from Subramani, et aI. [41] is an example of the variation

in the curvature for a typical wave train. Note the extremely large

negative curvature of the leading crest which subsequently broke 4.5

seconds later. Figure 10 demonstrates the effectiveness of the local

absorbing patch model for shallow water waves. The waves were made by a piston type wavemaker sweeping over the frequency range so that the waves coalesce at a point down the tank. As can be seen, without

the damper the waves break at about t = 11.4.

With the damper activated, the wave breaking is suppressed and the wave continues to propagate.

10 20 30 40

Distance along tank (m)

10 20 30 40

Distance along tank (m) 50

Figure 9: A representative _{snapshot at t=1 &8 of: (a) the surface}

displacement, and (b) the curvature of the surface, for waves generated by a wedge wave-maker of motion

amplitude 0.12m and frequency _{0.559Hz (nominal}

A.=5m). The leading _{wave front broke at about t=23.3} in this simulation (from Subramani, _{et al. [41])}

4.2 Three-dimensional computations

The Euler-Lagrange approach _{using desingularized isolated}

sources to

41 = IIn

code denominated UM-DELTA for

University of Michigan

Desingularized Euler-Lagrange Time-domain Approach. The code

distributes nodes on the body surface, free surface, and bounding surfaces at infinity. Finite depth calculations for a flat bottom are done using image sources below the bottom. The free surface nodes are

Distance along tank

Distance along tank

Figure 10: Time history of shallow water wave breaking caused by the coalescing of waves of different frequencies generated by

a piston wavemaker. Wave breaks in (a) at t=11.4 sec. In (b) a "local absorbing patch" has been implemented to prevent breaking (after Subramani [41])

(b) Direciioiio(wan

.lt=ii.s1

:

4= 11.41->

ti -t=iiol 41 _{r I i}

constrained to move along tracks with a velocity that is the derivative of

the potential in the direction of the track. _{Equations (10 ) and (11) are}

thus the appropriate free surface boundary _{conditions. The tracks are}

altered to keep the appropriate spacing as the body shape changes due

to non-wailsidedness. The prescribed tracks allow _{easy cubic spline}

interpolation for the regridding that occurs at the end of each major

time step. No numerical smoothing has _{been found to be necessary.}

The free surface elevations are either stable and smooth or the results quickly blow up if the wrong gridding or time step size has been used.

The primary difficulty with the code has been the geometric modeling

necessary to fmd the Unit normals at the arbitrary node locations on the

body. Unlike panel methods that_{use flat quadrilateral panels to define}

the surface and the unit normal, UM-DELTA _{distributes individual}

nodes over the body surface. Determining the surface location and

unit normal at arbitrary locations

_{on a typical ship hull has been}

difficult. Celebi & Beck [44] use a NURBS surface to define the hull but

_{that technique appears to be impractical}

_{because of the}

computational effort required at each time step. At present we are

using a cubic spline technique that is much faster but at times has

difficulties near the ends. What is needed is a very fast, robust routine

that given two coordinates (x and z) will return the third coordinate (y)

and the unit normals for an arbitrary point on the hull surface.

UM-DELTA has been used to solve a variety of three-dimensional problems; several examples are given below:

4.3 Incident waves diffracted by a vertical cylinder

Regular, shallow water waves incident upon a vertical cylinder have

been investigated by Scorpio [34] and Scorpio _{& Beck [45]. Waves}

were introduced at the upstream boundary by a piston wavemaker,and

wall boundary conditions _{were used on the other three truncation}

boundaries.

All of the wave energy that reaches the

_truncation

boundaries is reflected back into the basin; thus, _{the simulation is}

limited to the time it takes the reflected _{waves to return to the cylinder.}

The bottom and centerplane are accounted for Using images.

Three incident wave cases have been investigated. In _{all cases, the}

to water depth ratio (AIH) was 0.1. The wavelength was varied so that the wave number times the radius (kR) equaled 1.324, 1.094 and 0.825

for cases 1, 2, 3 respectively. For each case the non-dimensional

horizontal force, defined as F*

F/pgR2A

, is computed by

integrating the pressure over the cylinder surface.

An isometric view of the waves for case (1) at maximum run-up

on the

front of the cylinder

is shown ru figure 11. Contours of the

free surface elevations are shown at the base of the cylinder. The waves

are generated on the right side of the domain and propagate

downstream in the positive x-direction past the cylinder. The wave run-up for this case (cf. Beck & Scorpio [46]) was found to be 2.23 times the incident wave amplitude.

An example of the horizontal force time history in case (1) for both a coarse and a fine grid is shown in figure 12. The coarse grid has 900 nodes on the free surface, 55 on the cylinder, 75 on the piston wavernaker and 275 on the truncation walls. Similarly, the fine grid has 2136 nodes on the free surface, 55 on the cylinder. 140 on the piston wavemaker and 520 on the truncation walls. As can be seen, the

Figure 11: Isometric view and contour plot of the free surface elevation for incident waves diffracted by a cylinder

(HJR = 1.16, A/H = .1, kR = 1.324) (from Scorpio

two force histories agree very well. The force time histories are

predominately linear (the actual harmonic _{amplitudes are presented in}

Table 2 of Scorpio & Beck _{[45]) because the wave amplitudesare}

relatively small. _{The small amplitudes}

_{were chosen in order to}

compare the UM-DELTA predictions with the experimental results of

Chakrabarti [47] and calculations by Yang & Ertekin [48]. Yang &

Ertekin used a boundary element method with source panels

distributed on the free surface and body surface; a second order Stokes

wave was used to create the incident

_{waves and an extended}

Sommerfeld radiation condition developed by the authors was applied

on the outer boundaries. Table I (from _{[45]) compares the maximum}

non-dimensional horizontal force amplitude

_{found by the three}

methods. The agreement between the two very different computational

methods and the experiments is typically _{very good.}

s

Finegrid,

N =2851

Coarse grid, N = 1305

4 o .4

Figure 12: Convergence of the horizontal force time histories for

incident waves diffracted by a cylinder (HIR = 1.16,

A/H = .1, kR = 1.324) Coarse grid, N_{= 1305; Fine}

grid, N = 2851

o 2

B(

(2x

y(x, z) = .- '

-

LiT

Table 1: Comparison of the nondimensional maximum horizontal

force on a vertical cylinder

4.4 Wigley hull computations

Because of its simple mathematical form, many of the computations using UM-DELTA have been made for the Wigley hull form and are

reported in Beck, et al. [32] and Scorpio [34]. In this article, examples

will be given of calm water resistance, added mass and damping in heave and pitch, and the exciting forces in regular head seas.

The Wigley hull is a mathematical form with the hull surface defined by the equation:

(2x 2'\

(Ji+a2_)

J

where L = model length

B = modelfulibearn

T = modeidraft

a2 = coefficient for bow fullness = 00, standard hull

= .2 for modified Wigley hull ifi.

For all the computations presented in this article, L/B = IO and

B/T = 1.6 for both the standard and modified hull m.

Experiments Chakrabarti

Case

Calculations Yang

& Ertekin Present

from Yang & Ertekin % difference from experiment 1 3.07 3.14 2.98 5.1 2.9 2 3.45 3.41 3.42 0.3 0.9 3 4.14 4.09 4.00 2.2 3.4

The standard hull form has been extensively _{tested in model}

basins around the world. The experimental _{results presented in figures}

15 and 16 are from Noblesse & McCarthy _{[49] using the data taken at}

the University of Tokyo on a 2.5 meter model fixed in sinkage and

trim. A series of four modified Wigley _{models were tested at Delfi}

University of Technology by Gemtsma and _{Journée. Journée [50]}

presents the complete heave and pitch results for the _{added mass and}

damping, exciting forces, and motion amplitudes and phases for a variety of frequencies and Froude numbers.

As explained in Beck, et al. [32], the sharp bow and stern regions

of a ship hull can present problems for _{the desingularized method}

since the sources can cross the centerplane

_{and the perturbation}

velocities are large due to the stagnation _{points. From numerical}

experiments with the Wigley hull and the _{Karrnan-Trefftz airfoil, for}

which the analytic solution is known, techniques to handle these difficulties have been found. As seen in figure 13, a source point that

crosses the centerplane is moved back to the

_{centerplane without}

changing the x-coordinate. _{For a source corresponding to a leading}

edge node, the desingularization distance _{is greatly reduced in orderto}

prevent an unrealistic large tangential velocity on the body surface.

The leading edge _{source must be strong enough to cancel the free}

stream velocity at the leading edge node which _{is also a stagnation}

point. _{The farther from the leading edge node this source is placed,}

the larger its strength has to be. This large _{source strength will in turn}

induce a spike in the tangential velocity. As the leading edge source is

placed closer to the node, _{its required strength is reduced and the}

induced tangential velocity is likewise reduced. _{The net result is that a}

spike in the tangential velocity does _{not develop. The narrower the}

entrance angle, the greater the effect. Beck, et al. [32] show that with this source placement the agreement with the Karman-Trefftz airfoil analytic solution is within graphical accuracy along the entire surfacé

of the airfoil including the stagnation points.

The computed wave patterns for the standard Wigley hull at a Froude number of .3 are shown in figure 14 (results for other Froude

numbers are presented in Scorpio [34]). In figure 14, contours and

collocation points

Figure 13: Desingularization near the leading edge of a Karman-Trefftz airfoil (from Beck, et al. [32])

contour lines indicate z>0 while dashed lines are for z<0. Only half the grid (y>O) was used in the actual computations: the lower half is simply a mirror image included for illustrative purposes. The grid

spacing on the free surface is shown

in the figure.

In the

computations, 3267 nodes were distributed on the free surface and 663 nodes on the body wetted surface. The nodes were concentrated near the bow and stern regions in order to resolve the high flow gradients. The computational grid extends 1.2 ship lengths downstream, .5 ship lengths upstream, and .85 ship lengths in the y-direction. Upstream the potential and wave elevation are required to be 0.0, while the sides and downstream boundaries are open. The Kelvin wave pattern shows no signs of wave reflection at the open boundaries.

The experimental (from Noblesse & McCarthy [491) and

numerical predictions of the wave elevations along the side of the hull

for four Froude numbers are shown in figure

15. The numerical results were linearly interpolated from the even Froude numbers used in the calculations to the slightly different Froude numbers used in the experiments and shown in figure 15. Wave cut data at various distances from the hull can be found in Scorpio [34]. The agreement between

0-S

o

.05

Figure 14: Free surface contours and isometric view of the standard

Wigley hull advancing at Fr=O.30 throughcalm water.

In the isometric view, the z-axis is scaled_{by a factor of}

five. 3267 nodes are distributed on the free surface and

663 nodes on the body wetted surface_{(from Scorpio}

[34]).

.15

45

45

Figure 15: Comparison between computations and experimental

measurements for the wave elevation along the Wigley

hull. Experiments were conducted at the University of Tokyo on a 2.5 meter model fixed to sink and trim (cf.

Noblesse & McCarthy [49]). (from Scorpio [34])

the experiments and the numerical predictions is very good over the entire range of Froude numbers, especially the prediction of the bow

wave height that is typically underpredicted by linear theories.

A comparison of the experimental and numerical wave resistance

coefficients is given in figure 16. The wave resistance coefficients are

of the form C,,,, = F/(l/2pUO2S) , _{where F is the surge force and S is}

55

O F.. S5 U?T

F IS l-5 ISIS

h.Ll.flPT

T F l. .DflT

-ç-i

L

w

Lits -- _{O---- F,.L315}

?4T

k__

L.

T&I

the calm water wetted surface. The numerical result, CFI , is computed

by a pressure integration _{over the wetted surface after the model}

has

come to a steady state. The total resistance of the model is dominated

by the skin friction of the hull.

The wave resistance

_{can not be}

experimentally measured directly, but must be inferred from other measurements. Different measurements lead to different values for the

experimentally predicted _{wave resistance. Cpr. C , and C are the}

experimental coefficients determined by pressure integration over the hull, the residual resistance found in a calm water resistance test, and

linear wave-cut analysis respectively. _{As can be seen in the figure, the}

00020 00019 0.0018 0.0W 7 0.0016 0.0015 0.0014 00013 -. 0.0012 00011 ) 0.0010 o0009 0.000$ o 000i o cxxo o 0.250 0 260 0.270 0.280 0 290 0.300 Froude Number

Figure ¡6: _{Comparison of the wave resistance coefficients versus}

Froude number. Experiments _{were conducted at the}

University of Tokyo on _{a 2.5 meter model fixed in} sinkage and trim (cf. Noblesse & McCarthy [49]. C, C,,,,,, and C,, are the experimentalresults

obtained by pressure integration over the hull wetted surface (Cpr) residual resistance from the resistance

test (C), and wave pattern analysis (C,,,,,). CFI is the

computed force using the UM-DELTA _method.

o--c

- DELTA n

o--- c,

O--- C,.

MA

-'

J 0.310 0.320 0.330