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

T?4SCHE UVER$ITfT

Laboratorium 'toor j

Sthchana

koMet ekelweg 2,2628 CD Deft

1SL O1-1-FsaOI-18135

To appear in

Nonlinear Water Wave Interaction, International Series on Advances in Fluid Mechanics, Computational Mechanics

Publications, Southampton, UK, 1998.

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-2145 USA

e-mail: rbeck@engin.umich.edu

Abstract

The use of a desingularized 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 ìn the Euler-Lagrange method is solved using singularities that are placed outside thefluid domain

-inside the body and above the free surface. The desingularization 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 modifying 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

forces due to incident waves.

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 on the order of io to iO9 so that viscous effects are confined to the boundary layer and 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 iterative 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 employ 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 consists 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}

[li], Comte, et aI. [12],

Sen [13], 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], Dommermuth}

& Yue [16], Kang & Gong [17], Zhou & Ou [18], Chan & Causal [19], Celebi, et al. [20], Xu & Yue _{[211, DiMascio, 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-three-dimensional

problems. Explicit Euler time-stepping schemes are unconditionally unstable, while implicit-like arid implicit Euler schemes and fourth-order Runge-Kutta schemes are conditionally stable. They defined a local stability index itg(tt)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 know 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 with the uniform stream. Von Karman [26] used the technique for the

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 desïngularized 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. However, _{in practice, 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 ¡Ir potential. The _{use of complex numerical quadratures to} integrate over panels is eliminated with an equivalent savings in computer time.

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.

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

4=U0(t)x+(x,y,z,t)

(1)

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

satisfy the Laplace equation in the fluid domain:

V24=O (2)

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 co the x, y, and z axis directions

respectively. _{There is also a kinematic condition applied on the} bottom:

=-UO(t)nl+VBn on SB

₍₄₎

where VB is the velocity of the bottom _{relative to the Oxyz system.}

For an infinitely deep _{ocean equation (4) reduces to}

as zco

₍₅₎

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}

ataz

on SF ₍₆₎

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 o zero. _{Using Bernoulli's equation the dynamic condition}

ax p 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-4cc (8)

In addition, the initial values of the potential and free surface elevation

must be specified such that

=O tO

influiddomain

(9)

i0 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 usedto} 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 ofa point moving} with a prescribed velocity _{y relative to the Oxyz coordinate system.} By adding y VT1 to both sides of (6) and y V4 _{to both sides of (7)}

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

on SF (lO)

&az

and on SF where 6 a

- +v.V

is the time derivative following _{the moving node. This is similar to the} usual material derivative of fluid mechanics except the velocity

_{is the}

given y rather than the fluid velocity.

If y is set equal to (U(t)1v(t).).

_{the node follows a prescribed}

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

with the free surface. Setting y =

(o.o.)

results in the x-y locations of the nodes remaining fixed in the Oxyz coordinate system and equations (10) and (li) reduce to:

on SF (13)

and

=

_gi_!

U0(t)

on SF

(14)

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)i+VØ} Dt

and

= -gT + - VO. V -

.-where XF = (XF(t),yF(t),ZF(t))

is the position vector of a fluid

particle on the free surface and

Dt&

(17)

is the usual material derivative.

The form of the free surface boundary conditions given above allows the value of the elevation 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 step; the spatial gradient of the potential can be determined analytically

after solving the boundary value problem and the wave elevation is known. The difficulty is the gradient of the free surface elevation (Vi 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 V4 _{reasonably well. In this case, the contribution of the} V1 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 _{V1 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 maybe be open

or Dirichlet.

_{In terms of a}

desingularized source distribution, the potential

_{at any point in the}

fluid domain is given by

= _dû

(x)

_Jf

i

I1xs,

Q

where û is the integration surface outside the fluiddomain.

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

fJa(x5) û Ixc_xsI a

_sI

Q E

where ; = a point on the integration surface, Q

= a point on the real boundary

= the given potential value on the free surface at x

= free surface on which F is

known

X = the given normal velocity on the boundary at

= surface on which _x is known

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 Jth direction of the body axis system is thus given by:

Fj =JJPnj (21)

SR

where n is the generalized unit normal into the hull defined as

(ni,n2,n3) = n (22)

(n4,n5,n6)= rxn

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

or

-U (t)_gz_fVV$+v.S7

p & (23b)

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.

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, (he 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. _{[291, 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

r = (Ly.)

Oy = body axis system

and

j = 1,2,3

_{corresponds to the directions of the 015!}

axis

respectively.

The pressure in the moving coordinate system is given by

Bernoulli's equation:

= - - _{- U0 (t) - gz -}

desinguIaiiz isolat soJrce point.s

/z

/

Ç- given on the free

surface

$

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 =tdO)m) (24)

where 1d 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 1d = 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, another body, or the bottom. Using an array of isolated Rankine

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:

(x) =

j jx

-

+

- 4 I

NF 0F 'B

_aB

(25)

where NF is the number of the _{isolated sources above the free surface}

x _{is the location of the th source and}

is the strength of the th _{source at}

Similarly, NB, z

and a

_{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

NF

aF

N8

_a8

'ç, J

i':iIi

_4:I ijx

z

(i=1...NF)

or

(26)

where x and x _{are the node points on the free surface and body} surface respectively and _{N = NF + NB being the total number of} unknowns. Equations (26)_{can also be written in matrix form as}

AI=B

₍₂₇₎

1PBa1

i

ii'.'J

_xIJ

an1

U4

where A2)YLW

AAE2)

F si i=l...NF, j=L..NF F B i=L..NF, j=(NF+1)...N Ixc _xsjI

Ai21[B

']

i=(NF+l)...N. j=I...NF

A2)=_[B

']

i=(NF+l)...N, j=(NF+l)...N I)_{= aF} i=L..NF

z2

B _i=(NF+I)...N

B»

= ¡=L..NF

B2=x(xfl

i=(NF+I)...N

Equations (27) are a system of N = (NF + NB) equations for the unknown source strengths X1 . Once

Z

are known, the velocity

potential and fluid velocity can be analytically calculated anywhere in the fluid domain.

Domain Decomposition

Equations (27) can be. solved by direct, _{iterative, or multipole imthods}

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 CTh{RES

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 thistechnique _{(1) and Z2 are}

determined separately through

_{an iterative procedure.}

_Rewrite

equation (28) as

A(11) . (l) _{= B1 - AU2).} (2) ₍₂₉₎

A(22) (2)= B(2) - A21 . (1)

(30)

(1) _{is determined by solving (29) with}

(2) known (or guessed), then (1) is substituted into (30) to determine (2)_{. The new (2) is then}

used in (29) to update

_X'

. _{The procedure repeats until both}

equations are satisfied to the given error tolerance. _{The individual} matrices AU') and 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}

computationaily 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, _{Vr1 , are not required.}

The use of generalized nodes moving with prescribed velocity y

requires the spatial derivatives of TI 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 hydrodynamic 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, 84/& or must be determined at

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 ö4)/6t or d$iat must be known at the present time step. The conventional technique is to compute &416t or

a/at 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 aiat 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 arid 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].

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,

A1. The condition, and

hence the convergence rate of an iterative solver, can be greatly

improved if a good approximation for the inverse of A- is _readily

available and not too expensive

_{to compute. Let A be the}

approximation to

Aj'

_{, then the iterative solver can be preconditioned}

by solving

AjjAjksk =b1 ₍₃₁₎

for Sk, where _5k is related to the _{source strength by Z,= 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 technique, the inverse to}

A1 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 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)}

Multipole Acceleration

Multipole acceleration is a method to accelerate 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 N logN

_{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:

a\

The solution of this equation _{by separation of variables results in a}

series of spherical harmonic terms,

--ir

i

a(2atì

_i

a(.

O ₍₃₂₎

(34)

1L'r

(33)

n=Om=-n r

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 a1 _), i = 1,..., k and strengths . 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 at any point, P(r,G.9) , outside of

the sphere, in the form

0 Mm

V

0+1 0

n=Om=-n r

where the coefficients of the expansion are given by,

M =

(35)

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 lay 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) lay. The potential due to these sources at

any point, P(r,O,ç) , inside of the sphere is given by,

(b(P) =

where the coefficients are given by,

Ltm =_D i=1 Pi

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 sorne 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. lt}

_{is the careful}

arrangement, or clustering, of sources into groups that leads to the N logN 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 to 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-structured 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.}

muTapo4e

M

chiva

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 in conjunction 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 [341 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 =O 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

wìth preconditioning _{are increased for both methods since the 0(N2)} preconditioner matrix must be stored.

Numerical Results

Cao, et al. [29J 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 _{desmgularization 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.}

t

t

Figure 3: Schematic diapam of _{a dipole below a $ = O}

inimite flat plane (from Cao et al. [29]). Oiec* Method Method Coiroi PoUI .CO(l/VP*1 Hod. Pcirl _{j, . Swc* PI}

r

r

-

r

e. a q. q. q. q. q. X, _t

In figure 4, the numerical results for d4'/dn 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 desingulariz.ation distance, the

a

E

Figure 4: Effects of desingularization on the RMS error ---, direct application of the Green theorem; -, indirect method

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

isolated source solution blows up because the sources are distributed directly above the nodes. Cao, et al. [29) did examine staggered grids 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 1271) 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 desingularized _{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 Dornmermuth}

& 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 determine naturally from the}

4.0

2.0-Q

0.0

-2.0--4.0

0.0

1.0

2.0

Time/period

4"

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

three-dimensional cylinder in finite water depth equal to the radius (R = lO, T = O5), heave

amplitude a/T = 05. Motion;

- nonlinear calculation (present method);

- - - non-linear calculation (Dommermuth

& Yue [16]); - - - linear calculation

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

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. [39] 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/J}

, 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; however, depending on the problem it might or might not work well. If the computations in figure_{6 were to continue} much longer, it would be _{necessary to either have 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 forni 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

o = 1.0 t I 20 40 60

x+Frt

80 100 200 160 120 80 40

Figure 6: Waves generated by a submerged circular cylinder traveling

at a depdi Froude number of 1.0. r/b = .15, d/r =4,

at = 2, free surface nodes = 251, body nodes = 30 (from Cao, et al. [39]).

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 a4/an on the vertical_{upstream boundary and}

and T on the upstream free _{surface. The same boundary}

conditions can be used on the

_{downsteam 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, i,} should be identical. Figure 8 shows the results _{of the computations for}

ComputaxionaJ

window

Diiection o(

window ncxi.

Xo

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

o

Direction o(

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

dutctjoo of ,nC8dCM Ø - dirt1ios of t,do

o s I

Figure 8: Incident waves in a moving computational window. Solid line - wave elevations in the computational window. Dashed line - wave elevations in the stationazy wave tank (from Scorpio, et al. (31]).

-a '

:

\

Fr 13

---Fror7 8

---FrOlO -o 3 I I 80 2 o -2 -4 -4 o 3 s s s to 3 s s I 10 8 IQ

four different translation speeds. The incident waves have wavelength

and the speed of the computational window is defined by a Froude number based on the wavelength, Fn = uo/J . Because the

downsteam boundary condition is left open, the agreement while very good is _{not perfect. Specifying the downsteam 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 confined 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 wave 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 criterion for wave breaking, see for example Griffin, et al. [42] or Wang, _{et al. [431. Unfortunately, most}

of these criterion are hard to invoke

and/or 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 (K) near breaking and there is a limiting steepness (assumed _{to be H/X 1/12). Combining}

the above with extensive testing of fully nonlinear _{computations of} regular, two-dimensional waves resulted in _{a proposed threshold of}

IKal = .35. For values greater than this number, it is assumed 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:

U0(t)

on SF (38)

where is the damping term. The damping is confined to a small

region around the crest and in zero-speed computations has the form:

damp yu(x) IVI2 sgn(/&n) (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 the damper is a smoothly varying patch:

= .5{l + co4x(x

-

xo)/Loj}

where x0 = location where ka] = .35 is exceeded

L0 half-length, centered at x0 over which the damper acts.

L0 is presently assumed to be a3/.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 al. [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.

04

-' o

.1r V > -0 2 J .1 O IO 20 30 40 50

Distance along tank (m)

-4' I

0 _10. 20 30 40

Distance along tank (m)

Figure 9: A representative snapshot _{at t=l8.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 X=5m). The leading wave front broke _{at about t=23.3} in this simulation (from Subramani, et al. [41]).

Three-Dimensional Computations

The Euler-Lagrange approach using _{desingutarized isolated sources to} solve three-dimensional problems has been_{developed into a computer}

j

code denominated UM-DELTA for University of M ich igan

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]),

Di000(wsve ._-4 =1111 t) _{.__4tii.J.}

-'-s--> -°'

.t.

C o t) t) t) t=II4 ° ___1

ii

.o---\_-4'.

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

the potential in the direction of the track. Equations (IO ) and (il) 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-walisidedness. 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 find 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:

Incident Waves Diffracted by a Vertical Cylinder

Regular, shallow water waves incident upon a vertical cylinder have

been investigated by Scorpio [34] and Scorpio,_{et al. [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 (AJH) 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 in 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 (I) for both a coarse and find grid are shown in figure 12. The coarse grid has 900 nodes on the free surface, 55 on the cylinder, 75 on the piston wavemaker and 275 on the truncation walls. Similarly, the fine grid has 2136 nodes on the free surface, 55 on the cylinder, 14.0 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. (14fR = 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 [45j) because the_{wave amplitudes are} 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.

Fine giid,

N = 2851

Coarse grid, N

= 1305

4 o -4 i 4

Figure 12: Convergence of the horizontal_{force time histories for}

incident waves diffracted by a cylinder (H/R = 1.16, AIH = .1, kR = 1.324). Coarse grid. N = 1305; Fine

Table 1: Comparison of the nondimensional maximum horizontal

force on a vertical cylinder.

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'\

BI

Y(x.z)=_l()Jl()Il+a2r)

_J where L = model length

B = modelfull beam

T = modeidraft

a2 = coefficient for bow fullness

= 0.0, standard hull

= .2 for modified Wigley hull Ill

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

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

Case

Experiments Chakrabarti

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 Wigtey models _{were tested at Deift} University of Technology by Gemtsma and Journ6e. _{Journee [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 Karman-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 order to 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 further 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. [321 show that with this source placement the agreement with the Karman-Trefftz _airfoil analytic solution is within graphical accuracy along the entire surface 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_{ocher Froude}

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

/ s'

'--

-s. s' 5 coI1ocaion points 01 L.d

Figure 13: Desingularization near the leading edge of a

Karman-Trefftz airfoil (from Beck, e ai. [32]).

contour lines indicate z>0 while dashed lines are for z<O. Only half the grid (y>0) 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 arid 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 [49]) 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

0.5

o

.05

Figure 14: Free surface contours and _{isometric view of the standard} Wigley hull advancing at Fr=O.30 through_{calm 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)).

ai

.

I.t.D1

-LI

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 m figure 16. The wave resistance coefficients are

of the form C.1, = F/(lt2pUO2S), where F is the surge force and S is

.1 s_ils C P. Lw. WU.sVT

-t--

L.

_____

h Lw. ll_.._s_a -'a's sali

C

-4-s u LI 4-s si I-5

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 _{measure directly, but must be}

inferred from other measurements. Different _{measurements lead to different values}

for the

experimentally predicted _{wave resistance. Cpr}

.

C,

, and 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}

o---A1--- C

.--*---- c_p

O-- C,. - D ELIA

J 0310 0.3w 0.330

Figure 16: _{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 and McCarthy [49].

C, C.,,, and C,

_{are the experimental results}

obtained by pressure integration over the hull _wetted surface (Cpr) _{residual resistance from the}

resistance

test (Ce), and wave pattern

_{analysis (C. C}

_{is the}

computed force using the_{UM-DELTA method.}

o

0000 0.0019 o00is 00017 00016 00015 0.0014 00013 0001 0.0011 00010 o cx*» 00007 o 0006 0260 0.270 _{0210 090 0.300} Froude Number 0.250