Twenty-Fourth Symposium on Naval Hydrodynamics (2003)

Below is the uncorrected machine-read text of this chapter, intended to provide our own search engines and external engines with highly rich, chapter-representative searchable text of each book. Because it is UNCORRECTED material, please consider the following text as a useful but insufficient proxy for the authoritative book pages.

A Spectral-Shell Solution for Viscous Wave-Body Interactions R. W. Yeung and ]. A. Hamilton (University of California at Berkeley) ABSTRACT To solve viscous wave-body interaction problems, a domain decomposition technique is used in which a viscous flow in the near field is matched to an outer inviscid flow. A pseudo-spectral solution of the inte- gral equation describing the outer flow is formulated in such a way that the outer solution can be applied to a variety of interior-specific inner problems. This outer flow effectively absorbs outgoing waves and can independently supply incident waves. Application of this outer flow to viscous and inviscid interior flows is demonstrated and the specifics and limitations of the matching scheme are described. The interior prob- lems are solved by a pseudo-spectral finite-difference technique (Yeung and Yu, 1994) which requires the geometry to be cylindrical but does not restrict in any way the interior and exterior flow to be axisymmetric. INTRODUCTION Unsteady wave-body interaction problems are typ- ically solved using a velocity-potential formulation that neglects viscosity of the fluid. The inclusion of the effects of viscosity in these problems has always been a challenge because of the extensive computa- tional requirements associated with solving a viscous- flow problem over the large domain required to cap- ture wave effects, e.g., see Alessandrini and Delhom- meau (1996) and Chen and Huang (1998~. This paper describes a formal methodology of matching a near-field solution of a viscous wave-body interaction problem to the solution of an inviscid wave-like flow in the outer-field. The advantages of this domain decomposition method are particularly evident in the case of zero forward speed wave-body interactions. In these problems, vorticity is generated primarily by the body without a substantial amount of convective transport. The wave energy, however, travels relatively fast and if not accounted for, reflects off any domain-truncation boundary. This was in fact the restriction found by the very accurate viscous- flow work of Yeung and Yu (1994~. The domain of the inviscid solution is assumed to lie outside of a cylindrical "shell" which extends from the bottom to the free-surface (see figure 1~. The axisymmetric geometry does not restrict the general three-dimensionality of the problem and it is helpful to imagine this cylindrical surface as a flexible wave- maker whose motion couples the inner and outer so- lutions. This outer solution alone is formulated in terms of a time-dependent integral equation utilizing the unsteady free-surface Green function. To solve the necessary boundary integral equation, a pseudo- spectral technique is employed in which orthogonal basis functions are used to discretize the solution. Such use of the "shell-function method" (Hamilton and Yeung, 1997), which was carried out for a sub- merged spherical shell, allows the solution to be char- acterized by a set of universal coefficients that are integrals in time and space of the Green function. These integrals are typically computationally ex- pensive to evaluate but the methodology allows their pre-computation and subsequent usage for any prob- lem inside the shell. Additionally, with the axisym- metric shell geometry, the integration and evaluation of the free-surface Green functions are carried out us- ing semi-analytical techniques that take advantage of the orthogonal properties of the basis functions. The resulting generalized solution can be applied as an outer boundary condition for a variety of problems and solution methods in the inner region. Boundary- integral methods have previously been used in the in- ner region to solve nonlinear body motion problems (Hamilton and Yeung, 2000~. Even if one were to solve a purely potential-flow problem in the time domain, the outer boundary con- dition at an arbitrary truncation surface cannot be

stated as a simple condition since generated waves points on the bounding surface. This global property are dispersive and there are memory effects relating cannot be ignored because the outer field, despite the the current solution to the past solution. Difficul- presence of waves, is elliptic in property. A predictor- ties encountered in the early treatments of the open- corrector scheme is thus needed to find the solution boundary condition that are "non-reflective" can be on the matching surface that satisfies the boundary found in Yeung (1982a). A condition relating the conditions on the inner and outer flow-fields simul- velocity potential ~ and the normal derivative of po- taneously. This new treatment is described and its tential ~~ in a global (rather than point-wise) sense was first used by Nestegard and Sclavounos (1984) for water-wave problems in the frequency domain. Global interaction means all points on a matching surface are interacting, resulting in a matrix rela- tion. The ~ to (Pn relation was used simultaneously in three-dimensional time-dependent problems by Lee (1985) and Lin et al. (1985), the former used what we now call a "shell", the latter used a vertical line- dipole, to represent the outer flow. In a review pa- per, Yeung (1985) introduced the concept of "shell functions" to be used in conjunction with a "general shell." The work defined how a set of universal coe~- cients could be computed only once and re-used sub- sequently for any three-dimensional time-dependent interior problems. The idea was pursued further by Yeung and Cermelli (1993) and Hamilton and Ye- ung (1997, 2000) for a submerged shell with success in two and three dimensions, respectively. In a re- lated context, Keller and Givoli (1989) conducted stability studies of the Dirichlet-to-Neumann (DIN) formulation for acoustic problems. However, only recently has the shell-method's effectiveness been fully realized, particularly with the use of a pseudo- spectral outer representation and the introduction of a surface-piercing shell. Lin et al. (1999) reported computations for ship-motion problems using similar ideas, but provided no details on the treatment of the matching procedure, particularly some of the dif- ficulties associated with "surface-piercing" matching surfaces. In this paper, the exterior solution provided by the shell-function method is used as an outer bound- ary condition for two types of internal problems, in- viscid and viscous fluid. Inviscid internal problems have been applied to this type of domain decom- position before (Lee, 1985; Dommermuth and Yue, 1987), primarily for the case where the interior prob- lem is also solved by a boundary-integral equation method. The matching to a finite-difference solu- tion of the interior problem presents a new challenge. Finite-difference methods require point-wise bound- ary conditions while the integral equation solution of the outer flow provides only global condition. As explained earlier, a global condition is a relation in which the flow characteristics at one point are de- scribed in terms the flow characteristics at all other 2 accuracy verified here. Somewhat surprisingly, when tested in the frequency domain, the technique de- scribed here does not seem to demonstrate the usual irregular frequency effects associated with the use of a free-surface Green function. The predictor-corrector matching scheme can be extended to provide an outer boundary condition for the solution of viscous flow in the interior region. In this work, the viscous problem in the interior re- gion is solved using a spectral Navier-Stokes equa- tion solver for cylindrical geometries, developed by Yeung and Yu (1994) and previously presented for the case of a closed domain in the 20th ONR Sym- posium. The inviscid-viscous matching technique is applied to two types of scenarios: (a) an initial-value problem with waves in the viscous region passing out through the matching boundary without reflection, and (b) a transient-wave problem with waves, exist- ing initially in the outer inviscid region, propagating into the viscous interior region and diffracting about a body. INVISCID OUTER SOLUTION In this section, a general solution of the linear, invis- cid, outer flow is found by solving a boundary-integral equation for a velocity potential. A representation of this potential in terms of Chebyshev-Fourier basis functions allows the general solution to be expressed in terms of integrals of the free-surface Green func- tion. Numerical evaluation of these coefficients is per- formed differently than those in classical boundary- element method applications. The semi-analytical integration methods used are described below. The goal of this outer solution is to provide a general solu- tion to the outer flow which can be used as a bound- ary condition for a variety of interior flow problems. This work is a "'surface piercing" extension of con- cepts demonstrated in Hamilton and Yeung (1997~. Mathematical Problem for Outer Flow As illustrated in figure 1, the outer fluid domain VO is bounded by the matching surface SS, the free sur- face SFO' the bottom SHOP and a far-field surface A, located at infinity. All quantities in the problem are non-dimensionalized by a characteristic length L, a characteristic velocity U. and fluid density p, where

Outer h Flow z, , a: ,~6 '-) x So ·P = (r'E'z') Q = (a, 8.Z) ~ss Figure 1: Schematic of shell geometry. the tilde indicates dimensional variables. For each ~ problem studied, U and L can be chosen as appropri- ate and will be specified as needed. The velocity of the fluid is given by the gradient of a velocity potential that is initially zero and satisfies the Laplace equation at all time in the fluid domain, V2~(P,t) = 0 P ~ Vo (1) Linearized conditions are applied on all boundaries of the fluid domain, ptt (P. t) + (LIZ (P. t) = 0 P ~ SFO (2) /)Z(P, t) = 0 P ~ SHO (3) AMP, t) + B~l/(P, t) = F(P, t) P ~ SS (4) where the subscript I' indicates differentiation with respect to the surface normal, directed in the positive r direction on SS, and U is taken as A. The con- stants A, B. and the function F(P, t) are considered known, their exact form depends on the particular interior problem. Additionally, initial values of sb(P) and ¢~(P) must be provided and in the far field (as P- ~ x), the fluid velocities must vanish. This unsteady formulation is an initial-value prob- lem in which the boundary conditions advance the so- lution in time and a boundary-value problem is solved at each time-step. This requires initial conditions on the potential and on the free-surface elevation. Integral Equation Solution The hydrodynamic problem defined above is best solved by using a Green function that satisfies some of the boundary conditions to convert equation (1) to an integral equation, (Wehausen and Laitone, 1960~. An unsteady Green function G(P; Q,t—r) is used which represents the velocity potential at time t and position P due to the introduction of a source at point Q. in the presence of a free surface, at time a. The Green function must satisfy the following field equa- tion, boundary conditions, and initial condition. V G(P, Q. t - 7) = 5(P—Q) P ~ Vo Gtt(P,Q,t - T) + Gz(P'Q't - T) = 0 P ~ SFO Gz(P'Q,t-~) = 0 P ~ SHO Additionally, homogeneous initial conditions for G and Go must be satisfied on the free surface and the velocity resulting from the source must vanish in the far field. Such a function is given by Wehausen and Laitone (1960) as the sum of a singular part and a regular time-dependent part G(P, Q. t—'a) = Or (P. Q) + H(P, Q. t—r) t ~ ~ (6) By applying Green's theorem to the time derivative of the potential and the time-dependent Green func- tion, and integrating with respect to time, a time- dependent integral equation can be formulated, see Yeung (1982b). Because the Green function satisfies the linearized free-surface boundary conditions, the bottom boundary condition, and the far-field condi- tions, only integrals over the matching surface remain in the integral equation. If the geometry of Ss is taken as a cylinder of radius a and depth h (as shown in figure 1) and the variable of surface integration Q is expressed in cylindrical coordinates, Qtr = a ~ z' the integral equation can be shown to be, - 2~(P, t)+ `~27r `~0 at J ~~e z t)Gz,fP;6,z,t=0)dzdb— O —h All ~ '.27r `~0 - ~ at ~ ~6 z T)H~(P;6,Z, t—r~dzdd d7= O O —h '.27r `~0 at ~ ~ ¢§ z t)G(P;6,z,t=0)dzdd— O —h fit - rem r° a: ~ (PIJ(§ z 7)HT(P;6,Z, t—7)dzd~ dr O O —h (7) Where the field point P lies on the cylinder and H(P; Q. t—irk is the time-dependent part of the free- surface Green function, G(P; Q. t—r). Time Integrations The integrations with respect to time in equation (7) represent a convolution in time of the solution and the free surface Green function. The initial step in per- forming these integrations is an integration by parts. 3

As an example, consider the last time integral in (7~. where, I = Jt - t2~ to aJ ~ (P (A z T)HT(P;6,Z, t—r)dzdb do O O —h - 27r 0 t = a | | it z T)H(P; A, Z. t—T) O —h O t2~ r° pt at ~ ~ (PI7(6 z T)H(P;6,Z, t-'r)d~dzdb O —h O (8) To perform these integrations in time, the solution up to the current instant of time to is discretized by the time sequence Elk = i(tk), (tk = kAt; k = 0, . . ., K), and the potential is assumed to vary in a prescribed way between time-steps. A linear variation is used here, which removes ~~ from the time integrals, leav- ing only time integrals of the Green function which can be evaluated analytically. This formulation shows clearly how a higher order variation (quadratic, cu- bic, etc.) of potential between time-steps can be im- plemented. Previous work of Young (1982b) did not utilize this integration by parts but rather assumed the potential to be constant and equal to the average between time-steps. That constant-potential formu- lation is the lowest order approximation possible and can be recovered from (8) by recognizing the contri- bution from the infinite rate of change in potential due to the jumps at the beginning and end of each time-step. Each higher order of approximation intro- duces a new term into the integral equation which is another time integral of the Green function. Because this outer solution is developed for a fixed cylindrical geometry and a "compute once, use many times" ap- proach is being used, this is not a penalty as it would be in the case where the boundary integral equation method is applied directly to a moving body. With the approximation of linearly varying poten- tial, ~,~ is removed from the time integrals and the remaining time integrals of H can be defined as: ~K,k~p;§ z) = ~ | H(P;6,z, tK—T)dr (9) The first integral on the right hand side of (8) is zero because the potential at time zero is vanishes and H(P;0,z,O) = 0. The integration in time may be written as a sum, K—1 27r 0 I ~ Am, a / / A, Z,tk) INK kelp; 6, z~dzd§- ~ (P; §' Z) = {>K k+1 INK k k K By defining ~K-k~p; A, Z) as an analog to INK k ~ p; a, Z) with Ho (P; 6, z, tK—T) in place of H. the integral equation (7) for the potential outside of the cylinder at time to can be written in terms of these "shell functions": —2~r~ (P. tK ~ + r2~ r0 a; ~ ale, z, tic) [GT/(P; §, z, O)+~°(P; 6, z)]dzd0- O —h r2~ r0 a) ~ ~,,(6, z, tu) [G(P; §, z, O) +~°(P; 6, z)] dzdb = O —h K—1 27r 0 ~ al | (p(0 t ~ K—kelp ~ id do k=1 0 —h K—1 27r 0 ~ a; | BUNS z tk)AK k(P;§ z~dzd~ (12) k=1 0 —h The time-discretization procedure has split the in- tegral equation into terms that involve the velocity potential at the current time to and terms that in- volve the velocity potential at past time-steps. These "memory" terms have been collected on the right hand side of equation (12) as known quantities. Pseudo-Spectral Solution of Integral Equation To solve the integral equation (12) numerically, it is necessary to represent the solution discretely, enforce the integral equation at a set of collocation points, and solve the resulting linear system. In this work, the geometry of the integration surface is exploited by representing the potential over the cylindrical shell in terms of global orthogonal functions, Chebyshev polynomials [Tj~x) = cos(j cos-1 x)] are used in the vertical direction and Fourier components are used in the circumferential direction. This choice of global basis functions is a deviation from the usual piece- wise approximation of earlier works (Hamilton and Young, 1997) and superior accuracy of the solution can be demonstrated. The representation of ~ on Ss by orthogonal poly- nomials is summarized by the following equations J—1 N/2 - 1 (~(§,z,t) = ~ ~ ~nj~t)T2j~h + Lena (13) k—1 JO J—h j=o n=—N/2 27r 0 J—1 N/2 - 1 a; | cute z tu'>°(P;e zy~z(lc (joy (p~(§,z,l) = it, ~ ~ovni(~)T2j(—h + 1)e (14) O —h j=o n=—N/2 4

The prime on the first summation indicates that a factor of one half should be included when j = 0. In- serting these decompositions into the integral equa- tion removes ~ from inside the integrals and the ex- pansion coefficients cPnj become the unknowns. The integrals that remain in (12) are defined as Fourier- Chebyshev coefficients of the Green functions, h ~ ~ ~ l;2~/~0G¢P;~9 z O>~T2j(~h + 1~einedzd~ (~15~) r2~r r ° j( ) 2~Jo J_GL,(P;~,z'O)T2j(h + l)ein~dzdd (16) 1 r27r rO to J ~ (P; 0, Z)T2j(h + 1) (17) r27r r° nj (P)= 2~) / BY (P; §' Z)T2j(h + l) (18) A A Note here that G and G`, are coefficients of the im- pulsive (or Rankine) part of the Green function (i.e. the Green function evaluated at t = 0) and fly and ~ are coefficients related to the time-dependent part of the Green function as defined by equations (9) and (11). Denoting sonKj = inj(tK) and inserting these co- efficients into the integral equation (12), one obtains the following integral equation for ~ at time to, 2a a it= ' ~ Q(a,O,z) p(r',~',z') : ~ \ / R=2asin 1 2 1 (ifr' = a) Figure 2: Polar coordinate system for evaluation of Green function coefficients on shell surface. COMPUTATION OF SHELL COEFFI- CIENTS Because the formulation presented here is intended to be a universal solution of the outer flow that can be applied to a variety of internal problems, it is impor- tant to determine the most efficient way to compute and store the various coefficients for re-use. Regardless of the Green function in question, when P and Q are expressed in cylindrical coordinates, the coefficients (15) - (17) are defined by integrals of the form, 1 r27r Ib0 2~1o ;-h T2j~h + 1)dzein~d~ ~—1 N/2 - 1 (p(P,lKi)+aE, A, (pnj (Gv~~j(~P)+7nj(~P)] =e ins Gnj~z Are) (21) j=0 n=—N/2 J—1 N/2 - 1 a ~ ~ (PVK j tGnj~p) + ~nj(P)] = j=0 n=—N/2 K—1 "T_1 N/2 - 1 a ~ ~ ~ Nicks—kelps_ Thus, all required coefficients can be found easily from evaluations of Gnj (z', r', 0). The unsteady free-surface Green function that sat- isfies all boundary and initial conditions in (5) is given by Wehausen and Laitone (1960), k=l j=0 n=—N/2 K-l J-l N/2-l G(<P;Q,t—r) = Gr(`P;Q) + H(<P;Q,t—r) (22) a ~ ~ ~ (~,lj ~nj (P) (19) where Gr is the "Rankine" part of the Green function, k=l j=0 n=—N/2 Numerical solution of this equation can be achieved Gr(P; Q) = 1 + 1 - by a collocation technique, resulting in a linear sys- r tern that relates (pKj and <,OKnj in terms of the shell coefficients which are integrals in space and time of the unsteady, free-surface Green function, tAI(pK + ABACK = EN (20) where EN contains the convolution of the solution previous to to with the Green function. 5 {°° kh cosh k(z' + h) cosh k(z + h) Jo(kR)dk ~ J0 ~ cosh kh (23) Here r is the distance between the field and source point, R is the radial component of r, and r" is the distance between the field point and an image of the source about the bottom. The time-dependent part

H is given by, H(P; Q. t—7) = J2 °° cosh k(z' + h) cosh k(z + h) 0 sinh kh cosh kh |1—cosELw(k, h)(t—a)] Jo(kR)dk (24) the Fourier series representing the singular part is needed. This series is found analytically to be, lntR] = In [2a sin 2~] = lnta] + ~ k costly (27) Although non-singular, the time-dependent part of G is expensive to compute for a different reason. The last integrand in (25) is highly oscillatory with in- creasing k. The spatial integrations in (18) and (17) can be performed analytically in the ~ direction with the help of the 'addition formula' for the Bessel func- tion, where w(k, h) = Ok tanh kh. Although efficient methods have been developed to evaluate the unsteady free-surface Green function (Beck and Liapis, 1987; Newman, 1985), the total computation is still intensive in panel methods be- cause the integrations in space over the elements are typically done numerically. This can require many Jo(kR)= Jo(kr'jJO(ka)+ evaluations of the Green function, especially at large times when the Green function varies rapidly. Be- cause the shell surface is a fixed cylindrical surface, and global basis functions are used on this surface, analytic integrations can be employed to efficiently and accurately compute the shell coefficients. Integrals of (24) with respect to ~ are needed in the formulation. This integration is performed ana- lytically, Jtk Auk—k = H(P; Q. tK—T)d'T = tk—1 cosh Liz' + h) cosh k'z + h 2(tk—tk_1 ) L sinh kh cosh kh Jo (kR)dk Too cosh k(z' + h) cosh k(z + h) +2 ~ Jo w (k, h) sinh kh cosh kh tk sin[w(k, h) (tK—T)] Jo (kR)dk (25, tk—1 Integrals in space of these Green functions need to be evaluated as defined by equations (15) - (17) and figure 2. The singular terms in the Rankine part of the Green function (Gr) are integratable but direct numerical integration fails due to the singularity. To overcome this, the singular part must be subtracted and accounted for analytically. For the Rankine part, the technique developed here is to perform the inte- gration in the vertical direction first, resulting in a function G. of §. The remaining integration with respect to ~ is equivalent to finding the Fourier coef- ficients. However, Gas is singular at ~ = 0 (R = 0~. The strength of the singularity can be shown to be, —2T2j(h +l~lntR] (26, Subtracting this singular part from Gil results in a regular function, the Fourier coefficients of which can be found by FFT. To complete the evaluation of G', oo 2 ~ Jm (kr/) Jm (ka) cos(mb) (28) m=1 Multiplication of both sides by cos rid and integration from 0 to 2~ provide the needed integrals, 1 r2~ IJo(kR) cos ruddy = Jn(kr') Jn(ka) (29) 2~ o The introduction of a second Bessel function in (29) increases the rate of decay of the integrand by a factor of 1/~ for large k, making numerical integration more tractable. The remaining integration in the ver- tical direction is well behaved and can be done numer- ically with a change of variable and Filon quadrature. APPLICATION TO INVISCID FLOWS Matching of a linear outer-flow field has been done before for the specialized case when the interior prob- lem is also solved by a boundary-integral equation method (Lee, 1985; Dommermuth and Yue, 1987; Hamilton and Young, 1997~. In these studies the boundary of the interior domain is discretized in the same manner as the outer region and the relation between ~ and A', represented by (20) is included implicitly in the solution of the interior problem. In this technique, the entire flow field (inner and outer) is solved for simultaneously at each time-step. Recently, volume-discretization methods have been found to be competitive with the boundary-integral methods because of the sparseness of the resulting linear systems (Ma et al., 2001~. However, finite- difference and finite-element methods for solving the interior potential-flow problem require a point-wise boundary condition in which the relation between po- tential and normal velocity is specified at each point on the surface. Unfortunately, the shell method for the outer flow admits only a global relation between 6

z ~ ~ ~r-_~ C ~ ~— ' it= Figure 3: Schematic of an inviscid interior problem. the shell-function method, the general outer flow so- lution from the previous section is used to provide a boundary condition that exactly mimics the wave- like outer flow. In the inner region, numerical inte- gration of the free-surface boundary condition is used to advance the potential and wave elevation from one time-step to the next. The hydrodynamic forces and moments on the body are defined as the integral of pressure over the body: F(t) = / P(<P, t) n dS (37) SB (t) r , M(t) = / P(P, t) (OP xn) dS (38) SB (t) where the pressure p(P t) is given by either the lin- and ~L,, he. ~ at any point on Ss depends on ~~ . ~ ~ ~ . earlzed or exact Euler's integral depending on the everywhere on as. lo overcome Ants a prea~ctor- ' assumption of the body boundary condition. Consis- corrector Iterative process Is used to find the veloc- . ~ tent with the non-dimensionalization process these lty potential on the matching surface which satlsnes - the boundary conditions for both the inner and outer flows. The technique is described below. Formulation of Inner Inviscid Flow The mathematical problem for the inner flow is well known and has been studied extensively for both lin- ear and non-linear boundary conditions. A velocity potential o (separate from the outer potential ~) is defined and boundary conditions are specified, again all quantities are non-dimensionalized by some char- ~ ~ acteristic L, U. and p. The velocity potential must satisfy the Laplace equation in the interior region, V2<b(P, t) = 0 ~ ~ Vitt) (30) where Vi (t) is bounded by SB (<t) U SFi (t) U SHi U Ss With reference to the schematic of the inner flow in figure 3, the field equation and linearized boundary conditions are, Vo(P, t) n = Vn(P't) Cop, t) = -BLIP, t) It (P. t) = oz (`P, t) Oz(P, t) = 0 >r(<P, t) = Fn(P, t) - or - P ~ SBO (31) P ~ Ski (32) P ~ SFi (33) P ~ SHi (34) P ~ Ss (35) ¢(P, t) = Ed (<P. t) P ~ Ss where Vntt) is the body velocity in the direction of the body surface normal and ~ is the wave elevation. On the outer boundary Ss, if the normal velocity denoted by Fn(P' t) is chosen to be zero, the bound- ary condition becomes a no-leak wall condition. In forces and moments are normalized by pU2L2 and pU2L3, respectively. Inner Inviscid-Flow Solution Yu and Yeung (1995j develop a highly accurate pseudo-spectral finite-difference (PSFD) technique for solving the Poisson equation in an annular region. V2<b(P) = S(P) (~39) where S(P) is the "source term." This method is used here with a homogeneous right-hand side to solve for the inner velocity potential in conjunction with the outer shell solution which provides a bound- ary condition on the matching surface Ss. Numer- ical advancement of the free-surface conditions, to- gether with specified normal velocities on the body surface, provide a boundary-value problem for the interior potential, which is solved by the PSFD tech- nique at each time-step. Details of this procedure are explained in the reference. Matching Because the finite-difference techniques require point- wise boundary conditions, the shell function solution of the outer flow cannot be applied directly. In- t36; stead, the procedure is to predict the normal velocity on Ss at the next time-step, solve the outer-region flow using the shell function boundary-integral equa- tion with the predicted Neumann boundary condi- tion, and use this result as a Dirichlet boundary con- dition for the interior region. A follow-up corrector step improves the solution of the normal velocity at 7

no. on n1 n ? 0 20 40 60 time Figure 4: Horizontal force on cylinder as function of time, co = ~r/4, A = 0.05. the new time-step and stabilizes the method. Specif- ically, for the predictor step, we write lit at I (40) tK—1/2 tK—3/2 or equivalently, in differencing form: ~u = 2¢K—1 _ oK—2 (41) which provides a new value of ~K. The outer solu- tion represented by equation (20) now provides OK on SS, which is used with the matching conditions as a boundary condition for the inner flow. Note that the elliptic property of the outer flow is preserved in this manner. A corrector step uses the predictor-step solution to obtain a better estimate of OK, ant = Veil`' ~ (42) tK—1/2 tK—1 OK = OK This OK value is again used with equation (20) to provide a boundary condition on the inner flow to complete the advancement. This procedure can be seen as using the outer solution to provide the pres- sure on Ss (time-derivative of velocity potential), ap- propriate to the outer flow, subject to the predicted normal velocity on Ss. An alternate matching procedure is to predict the velocity potential ~ on Ss (instead of ~~) at the next time-step and use the outer solution to find the re- sulting normal velocity, thus providing a Neumann boundary condition for the inner problem. This al- ternate procedure is found to be unstable because of l Added Mass Damping 2 3 5 6 7 8 Figure 5: Non-dimensionalized sway added mass and damping as a function of wave number. an amplification of errors near the free-surface/shell- surface intersection. It appears that a slight error in the advancement of ~ leads to an error in the ver- tical velocity, which is coupled with Mitt through the free-surface condition. Another consideration is that equation (19) is a Fredholm integral equation of the second kind if ~ is considered unknown instead of in, thus offering more stability. The resulting linear sys- tem (20) has a diagonally dominant matrix A, which makes solving ~ from (pn more stable than solving i)n from A. The instability associated with a Dirichlet condition was found to be weaker if extremely small time-steps are taken. This is the first documented three-dimensional work involving a surface-piercing matching that provides a consistent and perfectly transmissive outer condition. Campana and Iafrati (2001) reported some success of the Dirichlet match- ing but for a two-dimensional convective flow and a submerged configuration. Inviscid-fluid Results This explicit matching procedure is applied to the un- steady swaying of a vertical cylinder in finite depth water. Characteristic length is chosen as the radius of the cylinder ri and U is defined as >/~. The cylin- der is initially at rest and for t > 0 the prescribed horizontal velocity is given by U(t) = Aw sin wt. The potential flow in the interior region is solved by the PSFD method. In this case, the interior region is truncated at a radius of five times the cylinder radius and the shell solution accounts for the wave behav- ior outside of this region. Figure 4 shows the result- ing horizontal force, non-dimensionalized by pgri3, for A = 0.05 and ~ = ~/4. After a few periods of os- 8

~ = 15.0 ~ ' i , ~r -2~ ,',: :~2' ~4,~-0-1 x t=30.0 _ , ~s ~ i__ _ _ _ __ ~ ~_ _ ~ I, . _ ~ 4~,,0~: t=20.0 1 ~ ' i - , l W_ l . W_ 1 y non \ ~ 0.1 O ~ at, . , ~ ~-4' 'I -0.1 'to, x Figure 6: Snapshots of incident waves diffracting about cylinder. cillation, the force on the cylinder achieves a steady state, and importantly, this steady state persists for many periods after the waves generated by the body oscillation have passed out of the domain of the in- terior solution, having been effectively accounted for by the outer shell solution and the explicit matching technique. The amplitude and phase lag (thus the added mass and damping coefficients) of the steady- state solution agree well with the analytical solution (Young, 1981) in the frequency domain. In order to verify the performance of the match- ing technique across a range of frequencies, this sim- ulation is repeated for many frequencies of oscilla- tion corresponding to a non-dimensional wave num- ber range of k = (O. 8~. At each frequency, the added- mass and damping in sway motion is computed from the steady state part of the force response. Figure 5 shows the computed results along with the analytical solution of the frequency-domain problem achieved by separation of variables. Excellent agreement is ev- ident but more significantly, there appears to be no irregular-frequency effects that are usually associated with the use of a free-surface Green function. The normal breakdown of the solution for a vertical cylin- der with the same radius of the shell occurs at the roots of the Bessel function of order one: Jo (krO) = 0. A very dense grid of computations is carried out near the first root k = 0.7663412... and no irregular be- havior is observed. Neither is there irregular behav- ior occurring at the first root of J~(kri) = 0 (or k = 3.831706...), corresponding to an irregular fre- quency based on the dimension of the physical cylin- der (rather than the matching surface). This amazing property of the time-dependent shell warrants further study. Incident-wave problems can be studied with the shell-function method by superposing the incident wave in the outer region only and solving for the to- tal potential in the inner region. This departure from the usual technique of modifying the body boundary condition to reflect the presence of incident waves is done with the aim of including nonlinear effects in the inner region, providing an outer boundary condition that not only absorbs outgoing waves but that can in- dependently supply incoming waves. The flow about a vertical cylinder in a transient incoming wave field using this technique is computed and the free sur- face elevation is shown in figure 6. To create a strong transient effect, the incident wave is generated by the collapse of a two-dimensional hump of fluid outside of the shell surface. The elevation at t = 0 is given by: 710(X)= xy~e-(X+7.5)2/2~ (44~ with ~ = 0.5. Note that transient diffracted waves are generated and pass out of the interior domain as the free surface returns to a quiescent state. The dominant diffracted wave is a ring wave that is clearly seen to be passing out of the interior region after the incident wave has passed. The time history of the 9

Inviscid Flow 7 Figure 7: Schematic of viscous-flow problem. The incompressible Navier-Stokes equations with their associated boundary conditions are solved using x primitive variables: velocity V and pressure p. The velocity vector in cylindrical coordinates has the com- ponents a, v, and w in the radial r, circumferential 0, and vertical z directions, respectively. Kinematic conditions are required on all boundaries and on the free-surface, stress-continuity relations appropriate to the wave behavior are also included. The unsteady Navier-Stokes equations in cylindri- cal coordinates are: horizontal hydrodynamic force associated with this transient event is shown in figure 10 and is plotted kit + us, + - ,~ + wig with the wave slope at ~ = 0. 1 [v2u _ u _ 2 dv] _ dP APPLICATION TO VISCOUS FLOWS This section applies the shell solution for the outer inviscid flow as an outer boundary condition to a viscous-flow problem in the interior region. The matching technique is similar to the case of an in- viscid interior flow described above. The interior viscous-flow problem remains challenging in three dimensions. In the present validation, a spectral fractional-step method for solving the Navier-Stokes equations in cylindrical coordinates is employed. This technique was developed by Young and Yu (1994, 2001~. A limitation of this technique is that it requires axisymmetric geometries, although there are no limitations on the flow itself, as can be seen from the more recent convective-flow computations of Ye- ung and Yu (2001~. The high accuracy and efficiency of the method make it attractive for developing a test of the viscous-inviscid matching techniques presented here. Formulation of Viscous Inner Flow MU MU V 8~ 811 112 TV TV V0V TV US +U + - +W — = At Or r 00 Liz r Re [V v - r2 - r2 06 I dw dw v dw bw at +u`, + - ,~ + w .~ = _ 1 HIP (46) Re [V w] - ,~~ (47) Here, the Laplacian operator in cylindrical coordi- nates is, V2 = 02/0r2 + 0/rbr + 02/0r232 + 02/~z2. The continuity equation in cylindrical coordinates completes the Navier-Stokes equations. 1 0(w) + ~ i~33v + jaw = o. (48) The quantity P used in this section, which is not to be confused with the field point variable in the ear- Figure 7 illustrates the physical geometry of the inner lier sections, represents the non-dimensional dynamic flow being considered. A vertical cylindrical strut Of pressure. It is related to the total pressure p by radius ri extends from a flat bottom to the free sur- face. Outside of the cylinder is an annular region filled with a viscous fluid, the outer cylinder shown can either be a rigid wall, or the more interesting scenario, a matching surface, on which the bound- ary condition provided by the shell method which behaves the same as a wave-like outer flow. As before, all quantities in figure 7 are non- dimensionalized by a characteristic length L, taken as ri, a characteristic velocity U. and fluid density p. The physical constants of viscosity and grav- ity will appear as the non-dimensional parameters of Ffoude number Fr = U/~, and Reynolds number — — Re= UL/z/. P=p+ F2' r (49) On the cylindrical body, prescribed velocities (U(t), V(t), W(t)) are enforced: u = U(t), v = V(t), w = W(t), at r = ri (50) There are no limitations on V(t) and W(t), but U(t) should be such that the resulting motion of the inner cylinder is small enough to be well modeled by the linearized conditions applied on r = ri. On the free surface, the boundary conditions can be linearized from the exact kinematic and stress- continuity relations (see e.g., Wehausen and Laitone, 10

1960). The linearized dynamic boundary conditions sor are flu + bw = 0 Liz Or p+ ~ + 2 bw Fr2 Re Liz Rev 1 bw _ + __ = 0 0, atz=0 (51) They provide the appropriate boundary conditions for velocities and pressure on SFO- The kinematic condition is ,~77 = w, at z = 0, (52) which determines the free-surface elevation A. Note the boundary conditions in equations (51) and (52) are satisfied on the mean free-surface z = 0, in order to be consistent with the linearization procedure of the outer flow. In the original work of Yeung and Yu (1994), no- slip boundary conditions were applied to the outer cylinder at r = rO. This wall condition effectively limited the length of simulations and also made the inclusion of an incident-wave field impossible. In this work, the rigid wall condition is replaced by condi- tions which couple the pressure and velocities to the outer flow through the shell functions. Two types of boundary conditions can be applied at the bottom of the viscous flow. A no-slip wall condition is the most realistic, as was carried out in Yeung and Yu (2001~. u= 0, v = 0, w= 0, at z =—h. Alternatively, a free-slip boundary condition is useful for alleviating the incompatibilities of the boundary conditions at the intersection of the bottom and the outer cylinder at r = rO Fx = / t~rTcosy ore sin04 dS SB JO (27r =—ri dz) d6Pcos0+ —h O Ri | demo do [2,' cost—,9 sine] = ~ ~rrzdS SB = R | do; d§,09W (56) My = J [err cos ~ ore sin d] zdS SB loo P27r = - ri J zdz; d§P cost+ —h O Ri I Adz Jo do [2~ cost—~ sine] (55) (57) These quantities are non-dimensionalized by pU2ri2 and pU2r3 for forces and moments respectively. Note that Fx and My consist of two terms, one due to pressure and the other due to viscous stresses. Inner Viscous-Flow Solution A time-stepping method for solution of equations (45) - (48) is used to solve for the velocity and pressure field in the viscous domain at a sequence of time- steps t = knot, (k = 1, 2, . . .~. The basic method fol- lows Chorin (1968) and is a fractional-step method in which an intermediate velocity field is found that satisfies the Navier-Stokes equation with the pressure term removed. To complete each step, a pressure field is then found that corrects the intermediate velocity field to form the velocity and pressure field at the new time-step. Considering a time-difference scheme (53) between the old (K—1) and new (K) time-steps, of equations (45~- (48~: ~ (uK _ uK—1) = ~ t_(U . V)U +—V2U~ —VP (58) V uK = 0 (59) du dv where v2 is the differential operator inside the brack- `, = O. ,~ = O. w = O. at z =—h. (54) ets of equation (45) and Q is a suitable difference operator. To solve these equations numerically, an intermediate velocity field u is introduced which sat- After the hydrodnamic problem is solved for a, v, isfies the momentum equation (58) without the pres- w, and P. the hydrodynamic forces and moments on sure terms. the cylinder are found as the sum of the pressure forces and viscous stresses contained in the stress ten- 11 MENU- UK-1) = 5! [ - (U V)U+ R Mu] (60)

Subtracting this from (58J, an equation for the pres- Dirichlet boundary conditions may be used in con- sure field results. junction with the Poisson equation for pressure at the new time-step pa. In the case of free-surfaces, the velocities at the new time-step are not known, Yeung and Yu devel- oped an algorithm in which the physical velocities and pressures from the previous (K—1) time-step are used in the right hand side of (63~. The boundary condition on pressure is found from the free-surface boundary condition and a predictor-corrector scheme in which the pressure equation is solved twice at each time-step. The matching technique developed below follows a similar procedure, using the velocity of the K—1 and K—2 steps on the matching surface to provide boundary conditions on the intermediate ve- locity field, then using the pressure at the new time step supplied by the outer flow as a boundary condi- tion on the pressure equation. ,,`~`UK _ Uy = _vpK (61) Taking the divergence of this equation, and using the continuity equation (59), one obtains a Poisson equa- tion for the pressure field. v2pK = i\t YOU, (62) Use of (62) ensures mass conservation in the numer- ical scheme without iteration. The algorithm to compute the velocity and pres- sure field at the new time-step K is to solve equation (60) first for the intermediate auxiliary velocity field, u. This provides a right-hand side for the pressure equation (62~. Note that this equation for pa iS a Poisson equation in the exact form of (39) and can be solved by the same algorithm as that used earlier. Finally, equation (61) provides the velocity field at the new time-step. The solution of (60) for u must be done accurately and is accomplished by a spec- tral collocation method which is also developed and presented in Yeung and Yu (1994~. The high effi- ciency and accuracy of the techniques demonstrated in these works stems from the cylindrical geometry and the decomposition of the solution into spectral modes. In the above solution algorithm for solving the vis- cous flow field, it is necessary to provide boundary conditions for the intermediate velocity field u = (u, v, w). The boundary conditions for the interme- diate velocity can be found in terms of the boundary conditions on the physical velocity and pressure. Ye- ung and Yu (1994) carefully develop the boundary conditions on the intermediate velocity in terms of physical velocities, these results may be summarized as: U = fl(/\t~uK,vK WE pa) V = f2~\t,uK,vK WE pa) - = f3~/\t, uK, vK wK pK) with the resulting ODE's advanced by an ADI method. When Dirichlet conditions are applied to the ve- locity field, these boundary conditions can be used directly. When free-slip conditions are applied as in the case of the bottom surface (54), derivatives of (63) provide the needed expressions. Boundary conditions on the pressure field must also be specified, either Neumann boundary conditions or 12 Matching of Inner and Outer Flows The matching of the viscous-flow solution to the inviscid-flow shell solution can be carried out in a manner similar to the inviscid-interior case described above. Again, the radial velocity on the shell surface from the interior solution is used to predict the new radial velocity at the new time-step, which is used as a Neumann boundary condition for the outer-flow solution. Since a velocity potential no longer exists in the interior domain, the new outer potential implies a Dirichlet condition on the pressure of the interior flow. Again a predictor-corrector sequence is required for stability. Hence, by analogy to equation (41), the predictor-step normal derivative of the outer poten- tial is OK = 2nK—~ _ ok—2 (64) This is used as a Neumann condition for equation (20), which provides a new outer potential ~K. Pres- sure on the shell surface at the new time-step is ob- tained from a backward-difference form of the lin- <63y earized Euler integral. _ OK _ OK—~ (65) The quadratic terms in velocities are not needed for reason of consistency. Application of this pressure as a boundary condition on the interior problem pro- vides a better estimate for the radial velocity at the new time-step, K K_i + UK—us 2 (66)

y o -1 -2 -4 _5 4 3 2 Y n -1 -2 3 4 -5— a -2 n 4 -0.020 -0.018 -0.016 -0.014 -0.01 1 -0.009 -0.007 -0.005 -0.003 -0.001 0.001 0.004 0.006 0.008 0.010 3 2 Y O -1 -2 -3 4 -2 0 x 2 4 -4 -2 X 2 4 Figure 8: Wave-elevation contours ~(x, y) and velocity vectors for a non-axisymmetric Cauchy-Poisson prob- lem in a viscous fluid near a cylinder, Re = 5, 000. 13

0.2 0.1 ~ o -0.1 -0.2 o It=7.6~ 3 1 2 .......... ............... - - -- r0 = 10 Calculation r0 = 5 Calculation \ , . ;~ ~ , 0.500 0.409 0.318 0.227 0.136 0.045 -0.045 -0.136 ·0.227 -0.318 -0.409 -0.500 6 7r Figure 9: Comparison of velocity vectors and ~ component of vorticity in ~ = 0 plane for a non-axisymmetric Cauchy-Poisson problem with the shell located at r0 = 10 (top) and r0 = 5 (bottom). which gives the final pressure boundary condition for the interior region after (20) is used again: K INK _ INK—~ P — /\t In viscous inner flow, care must be exercised in the tangential velocities. While derivatives of the outer potential could also be used as tangential velocity boundary conditions for the inner flow, consistent re- sults are achieved by implementing a free-slip veloc- ity boundary condition on the shell surface. This ap- proach is based on the physical idea that in a viscous- invisid matching, the inviscid outer flow is expected to be incapable of supporting shear stresses on the matching surface. The predicted radial velocity could be used as the final, required boundary condition for the inner flow but good results with less computa- tion are achieved by simply using the most current time-steps radial velocity in the expressions for the boundary conditions on the auxiliary velocity A. In this way, the flow is driven by an applied pressure on the shell matching surface, analogous to the free- surface treatment used by Yeung and Yu (1994) in the original development of the PSFD solution of the viscous flow problem. This matching procedure differs from that used by ~ ~ v - v Campana and Iafrati (2001), as noted earlier, their Figure 9 shows the flow velocity and vorticity in the approach is the opposite, using pressure from the in- ~ = 0 plane of the interior viscous flow. Evidently, terior flow to advance the outer velocity potential through the Euler integral. Although they demon- strate success and good results, that approach has `67' been found unstable when used in conjunction with a surface-piercing shell. Results for a Viscous Inner Flow To demonstrate the effectiveness of the viscous- inviscid matching procedure, a Cauchy-Poisson type problem is solved in which a Gaussian-shaped hump of fluid near the body is allowed to collapse, generat- ing waves which impinge upon the cylindrical piling and pass out of the viscous domain. The hump is given by the following form at t = 0: POOPS = 0.1e-2 6r (68) where r- is a polar coordinate system centered at (2.3,0~. Figure 8 shows the wave elevation at several time- steps after the release of this hump in the viscous region. These plots show clearly that the shell- boundary condition at r0 = 5 is effective in absorbing the wave energy without reflection. For a consistency check, the problem is solved for two radial locations of the shell matching boundary. rig = 5 and rig = 10. 14

0.15 0.1 0.05 o -0.05 -0.1 -0.15 Horizontal Pressure Force, Re = 1000 - ' . Horizontal Viscous Force, Re = 1000 . Horizontal Force, Inviscid Flow Slope of Incident Wave at Origin ,1: \ 1 0 5 10 15 20 25 30 35 time Figure 10: Incident-wave slope at origin and resulting force on a vertical cylinder. the viscous-inviscid shell has little effect on the flow strafed and explored in this work. This is the first details, especially close to the body. For each of these three-dimensional free-surface piercing shell which two test cases, the characteristic velocity U is cho- sen as/ and viscosity is chosen to give a Reynolds number of 5, 000. has proven to be most efficient when used in con- junction with a pseudo-spectral representation of the outer flow in finite water-depth. The shell is fully The transient incident-wave computations from the transmissive in wave properties in both incoming inviscid section are repeated here with the viscous in- terior problem fully coupled to the inviscid outer flow, which supplies transient incident waves. The total wave elevations that result have a similar appearance to the inviscid case, with a ring wave being the dom- inant diffracted wave which is transmitted properly outwards by the shell boundary condition. Figure 10 compares the resulting horizontal force on the cylin- der with the inviscid-fluid computations. The pres- sure forces resulting from the viscous flow are similar to the inviscid results in the early part of the simu- lation but ultimately have a smaller magnitude and decay faster as the transient incident-wave passes. A viscous force due to shear stress is also present but remains only a few percent of the pressure force since flow separation has not developed in this short time interval. In this solution, the flow stays attached to the cylinder. For separation to occur, the frequencies must be much lower than those in this example, re- quiring a non-transient incident wave field that starts from rest and achieves a steady state wave behavior; but this is simply a matter of having an appropriate incident wave for the problem. CONCLUDING REMARKS A highly effective outer boundary condition based on the use of shell functions for an inviscid outer flow, as initially introduced by Yeung (1985), is demon- 15 and outgoing directions. It is characterized by a set of pre-computed coefficients based on the time- dependent free-surface Green function of Finkelstein (1957~. The pre-computed coefficients can be stored and re-used for a variety of internal wave-body inter- action problems. A new "explicit" point-wise matching procedure is presented for coupling the interior solution to the fully transmissive shell. The procedure enables the advancement of both interior and shell solutions in a stable manner. The basic reasoning relies on the fact that the velocity potential in the outer field must re- ceive a normal-velocity condition based on a matched advancement of the pressure field on the shell. The el- liptic property of the flow is preserved through the use of time-dependent Green's theorem which relates the potential and its normal derivative. This simple pro- cedure is validated for several Cauchy-Poisson prob- lems with the presence of a vertical cylinder as a body in the interior region. Solutions in the interior region are obtained by highly accurate and efficient pseudo- spectral finite-difference (PSFD) methods for invis- cid flow (Yu & Yeung 1995) and full Navier-Stokes flow (Yeung & Yu, 1994~. When used in combination with the shell theory and the new matching condi- tion, these PSFD methods can uncover interesting fluid physics associated with the viscous interaction of surface waves and a body heretofore not possible.

References Alessandrini, B. and Delhommeau, G., "A multigrid velocity-pressure-free surface elevation fully cou- pled solver for calculation of turbulent incompress- ible flow around a hull," Proceedings of the 21st ONR Symposium on Naval Hydrodynamics, liond- heim, Norway 1996, pp. 40-55. Beck, R. F. and Liapis, S. J., "transient motions of floating bodies at zero forward speed," Journal of Ship Research, 313:164-176, 1987. Campana, E. F. and Iafrati, A., "Unsteady free sur- face waves by domain decomposition approach," Proceedings of the 16th International Workshop on Water Waves and Floating Bodies, Hiroshima, Japan, 2001. Chen, H. C. and Huang, E., "Validation of a chimera rans method for transient flows induced by a full- scale berthing ship," Proceedings of the 22nd ONR Symposium on Naval Hydrodynamics, Washington, D.C., 1998, pp. 948-963. Chorin, A. J., "Numerical solution of incompressible flow problems," Studies on Numerical Analysis, 2: 64-71, 1968. Dommermuth, D. G. and D.K.P. Yue, D. K. P., "Numerical simulations of nonlinear axisymmetric flows with a free surface," Journal of Fluid Me- chanics, 178:195-219, 1987. Finkelstein, A. B., "The initial value problem for transient water waves," Communications of Pure and Applied Mathematics, 10:511-22, 1957. Hamilton, J. A. and Yeung, R. W., "Shell-function solution for 3-d nonlinear body-motion problems," Schiffstechnik, 44:62-70, 1997. Hamilton, J. A. and Yeung, R. W., "Motion of a submerged body in incident waves," Proceedings of the 15th International Workshop on Water Waves and Floating Bodies, Ceasarea, Israel, 2000. Keller, J. B. and Givoli, D., "Exact non-reflecting boundary conditions," Journal of Computational Physics, 82:172-192, 1989. Lee, M. Y., Unsteady Fluid-Structure Interaction in Water of Finite Depth, Ph.D. thesis, University of California, Berkeley, 1985. Lin, W. M., Newman, J. N., and Yue, D. K., "Nonlin- ear forced motions of floating bodies," 15th Sym- posium on Naval Hydrodynamics, Hamburg, Ger- many, pp. 33-49, 1985. Lin, W. M., S. Zhang, S., Weems, K., and Yue, D. K. P., "A mixed source formulation for nonlin- ear ship-motion and wave-load simulations," Pro- ceedings 7th International Conference on Numeri- cal Ship Hydrodynamics, Paris, Ffance, 1999. Ma, Q. W., Wu, G. X., and Eatock Taylor, R., "Fi- nite element simulation of fully non-linear inter- action between vertical cylinders and steep waves. partl: Methodology and numerical procedure," In- ternational Journal for Numerical Methods in Flu- ids, 3:265-285, 2001. Nestegard, A. and P.D. Sclavounos, P., "A numer- ical solution of deep-water wave-body problems," Journal of Ship Research, 1984. Newman, J. N., "Algorithms for the free-surface green function," Journal of Engineering Mathe- matics, 19:57-67, 1985. Wehausen, J. V. and Laitone, E. V. "Surface waves," In S. Flugge, editor, Handbuch der Physik, vol- ume 9, pp. 446-778, Springer-Verlag, 1960. Yeung, R. W., "Added mass and damping of a verti- cal cylinder in finite-depth waters," Applied Ocean Research, 33:119-133, 1981. Yeung, R. W., "Numerical methods in free-surface flows," Annual Review of Fluid Mechanics, 14:395- 442, 1982a. Yeung, R. W., "liansient heaving motions of floating cylinders," Journal of Engineering Mathematics, 16:97-119, 1982b. Yeung, R. W., "A comparative evaluation of numer- ical methods in free-surface flows," IUTAM Sym- posium on Hydrodynamics of Utilization of Ocean- Wave Energy, Lisbon, Portugal, 1985, pp. 325-356. Yeung, R. W. and Cermelli, C., "Shell functions: A global method for computing time-dependent free- surface flows,", Proceedings of the 8th International Workshop on Water Waves and Floating Bodies, St. Johns, Newfoundland, 1993. Yeung, R. W. and Yu, X., "Three-dimensional flow around a surface-piercing body," Proceedings of the 20th ONR Symposium on Naval Hydrodynamics, Santa Barbara, CA, 1994. Yeung, R. W. and Yu, X., "Three-dimensional free- surface flow with viscosity: A spectral solution," Hydrodynamics in Ship and Ocean Engineering, RIAM, Kyushu University, pp. 87-114, April 2001. Yu, X. and Yeung, R. W., "Interaction of tran- sient waves with a circular surface-piercing body," ASME, Journal of Fluids Engineering, 117:382- 388, September 1995. 16