The Proceedings: Fifth International Conference on Numerical Ship Hydrodynamics (1990)

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.

Two-Dimensional Numerical Modelling of Large Motions of Floating Bodies in Waves D. Sen Memorial University of Newfoundland, Newfoundland, Canada J. S. Pawlowski NRC and Memorial University of Newfoundland, Newfoundland, Canada J. Lever and M. J. Hinchey Memorial University of Newfoundland, Newfoundland, Canada Abstract A numerical method is described, which s imulates in the time domain the propagation of steep two dimensional periodic waves and the large motions induced by the waves on free floating bodies. The method allows for mild transient phenomena. The algorithm is based on Green's formula for harmonic functions applied in a finite control fluid domain, together with the fundamental solutions of Laplace ' s equation. A lowest order boundary element discretization is used. In addition to several numerical results, computations of the sway forces and the roll and heave motions induced by steep periodic waves on a floating body restrained in the sway mode are presented and compared with the results of specially conducted model tests. 1. Introduction The increasing accessibility of computers of high capacity and the parallel development of computing techniques lead to the Leasability of algorithmic solutions of complicated hydrodynamic problems through a discretization of the corresponding governing equations in their fundamental form. In this paper, a solution of this type is discussed for the problem of motions induced by steep water waves on a floating body. The problem is formulated in two dimensions, in the time domain and potential fluid flow is assumed. A number of numerical models of the same problem have been developed s ince the appearance of the integral equation formulation combined with the time stepping scheme for the non- linear free surface conditions originally employed by Longuet-Higgins and Cokelet [1] in a study of the propagation of steep waves. Faltinsen [ 2 ] cons idered forced heave motions of a two-dimensional circular cylinder, as well as a related problem of sloshing [ 3 ] . Vinj e and co-workers extended their earlier work on breaking wave s imulation (Vinj e and Brevig [ 4 ] ) to include submerged and surface-piercing bodies in the fluid [ 5, 6 ], and next attempted 351 the problem of motions of floating bodies [ 7, 8 ] . Subsequently, the simulation of a capsizing of Salter' s duck, an ocean energy extracting device, was reported in [ 9 ] . For this latter study, experimental results supplemented the numerical simulation. Following largely the techniques of Vinj e and co-workers, with an important modification in the cons iteration of the body- free surface intersection point, Lin [ 10 ] simulated two - dimensional waves generated by a wave-maker in a f inite rectangular tank . This s tudy was followed by an extension to forced oscillations of axi-symmetric three-dimensional cylinders, [11] . Subsequently, Dommermuth and Yue [12,13] investigated the three-dimensional axi- symmetric problem and were able to s imulate forced heave oscillations of large amplitude of cylinders and inverted cones in an otherwise undisturbed free surface. Greenhow and co- workers also employed the approach of Vinje and applied the method in a study of the two- dimensional impact problem [14,15,16]. In the adaptation, specific improvements and developments of the algorithm were made to make it suitable for the particular application. Isaacson [ 17 ,18 ] reported on a similar method for studying two and three-dimensional fixed and floating body problems. In spite of those efforts, to the knowledge of the authors, a fully satisfactory solution of the problem has not been achieved. The study presented here was undertaken to investigate the possibility of developing a relatively simple but robust numerical model of the problem, which would be as suitable as possible for future generalizations to three dimens ions . This in particular lead to the requirement that the wave excitation in the control fluid domain and the radiation condition for the outgoing waves be implemented in forms which would not preclude an efficient modell ing of open water conditions . S imilarly, the modelling of wave breaking was excluded in order to ensure long periods of simulation for steep wave conditions. In this way the presented model acquired its characteristic features which distinguish it from the models

described above. The computer codes used in the study were written in FORTRAN and the computations were performed on DEC-VAX 8800 and -VAX 8530 cluster at Memorial University of Newfoundland. The study as presented here was completed in April 1988. 2. The Governing Equations and Their Boundary Element Formulation 2.1 The Governing Equations A finite two-dimensional control domain, D, containing fluid is considered. The domain is bounded by a piecewise smooth contour ~D, as shown in Fig. 1. The boundary ED is composed of the free surface IF, an impermeable bottom ADD, and the wetted contour DUB of a partially or fully submerged impermeable rigid body B. The remaining part of ED consists of two vertical geometrical control boundaries ~DC1 and ~DC2, separated by such a distance that D contains the submerged part of B. A Cartesian coordinate system Oxz is chosen, with the z axis directed vertically upwards and origin O coinciding with the intersection of ~DC1 with the undisturbed free surface. The fluid in D is considered to be inviscid, incompressible and homogeneous, and the flow is assumed to be irrotational. If B is completely submerged, thereby rendering D a multiply connected domain, the flow is assumed to be acyclic. Therefore the flow in D is described by a single valued velocity potential ¢(x,t), where x denotes the radius vector and t denotes time. The continuity of flow requires that the potential should satisfy Laplace's equation in D: a2¢ a2¢ + = 0 ~x2 ~z2 On the free surface IF the kinematic condition: 30 = 3¢ 3¢ 30 (2) Bt Liz ax ax and the dynamic condition: it ~ { P + go + 2 (V¢)2} (3) are imposed, where ~ denotes the free surface elevation, g denotes the acceleration of gravity, and V represents the gradient operator: V ~ B/8x; Pa is the applied pressure on the free surface, set equal to 0 in the sequel, and p signifies fluid density. On ODD the impermeability condition is applied Bi - O an where B/3n ~ n.a/8x, and n designates the unit normal vector on aD directed outwards of D. The same condition imposed on the wetted body surface aDB gives: In n (5) where Vn is the velocity component of aDB along its (outward to D) normal. For bodies fixed in space, Vn = 0. In addition, boundary conditions on control boundaries aDC1 and ~DC2 are assumed to be applied in such a form that either ~ or B¢/3n are known at all time instants t _ 0. The imposed conditions are described later in connection with specific applications. 2.2 The Boundary Element Formulation and Its Discretization 2.2.1 The Boundary Integral Formulation The application of Green's second identity to ~ and the fundamental solution In r(P,Q), with point P located in D or on ~D, results in the well known Green's formula for harmonic functions [19]: O(P)~(P) = lad [(A) an In r(P,Q) - -8 ¢(Q) In r(P,Q)] dS(Q) (6) ant) with point Q located on ~D. Here r(P,Q) ~ |X(P) - x(Q)| is the distance between points P and Q; the subscript of B/Bn indicates the point at which the differentiation is taken; n(P) = 2x for P located inside D but not on ~D. For P located on ~D, O(P) is the angle subtended inside D by the tangents to ED at P. which is equal to ~ when the normal to aD is continuous at P and the Cauchy principal value of the integral over ED is taken. Formula (6) expresses the potential ~ at any point P inside D in terms of its boundary values and those of its normal derivative. When P is placed on ED and values of ~ are prescribed on a part of ED and those of B¢/3n are imposed on the remaining part, (6) represents two coupled Fredholm's integral equations of the second kind with respect to ¢(P) as the unknown, and of the first kind with respect to B¢(P)/3n as the unknown. Formula (6) is valid at any instant of time, and therefore for solutions advancing in time, this relation can be applied at every consecutive time step. 2.2.2 Discretization (4) It is assumed that ED consists of Ns piece- wise smooth open contours: 352

Ns ED ~ U ask k=1 (7) The open contours ask are further subdivided into a finite number of segments, each approximated by a straight line between its end points: Mk Ok ~ U ~ S , it (8) A collocation point Qk is chosen on each of the segments bSj. The functions ~ and B¢/3n on AD are approximated by constant values on each segment and the values are determined at the corresponding collocation points Qk . Following the usual practice, the collocation points are placed at the centre of each segment. However, in principle Qk need not be located centrally in Ask. The discretization scheme maintains a consistent order of approximation (see e. g. [20] ) and represents an application of the lowest order boundary element method or b.e.m., [ 21 ] . A system of linear algebraic equations is thus obtained from (6), (7) and (8), with respect to the unknown ok or B¢k/3nj values at collocation points Qk. The influence coefficients in this system depend only on the geometry of the boundary contour and can be determined explicitly (see e . g. [ 3 ] ) for the present type of straight line segments. The bottom condition (4) permits exclusion of ADD from the contour of integration in the integrand in (6) if ADD is horizontal, by augmenting the fundamental singularity with its symmetric image with respect to ADD. Thus, when the bottom is horizontal at a depth d, In r(P,Q) in (6) is replaced by In r(P,Q) + In r(P,Q'), and ADD is discarded from aD. Here Q' is the symmetric image of Q with respect to aDD. This results in a reduction of the system of linear equations by the number of segments used in discretizing ADD. The influence coefficients then contain additional integrations over the image segments. It is possible to apply a higher order formulation of the b . e . m., for example, an isoparametric representation of Ok and B¢k/3nj over the segments, {Sk (see e . g. [21] ) . Such refinements are achievable at the expense of additional algorithmic complexities. In particular, complications arise in the treatment of the point at which the normal to AD is discontinuous or at the intersection of two parts of the boundary on each of which different conditions (i. e. ~ or B¢/3n) are imposed, which in the present application are the features of the body- free surface intersection point. In addition to a discontinuity of the normal at this point, conditions posed on the body and the free surface are different (Neumann condition on the body boundary and Dirichlet condition on the free surface) . In this respect, the central collocation discretization scheme adopted here 353 avoids an explicit occurrence of the problems and as a result is the most straightforward for numerical implementation. In addition, higher order applications of b. e .m., which lead to a better resolution with lesser number of segments (see e.g. [22]), are known to be more susceptible to numerical instability [ 24 ] which is a maj or concern in the considered application. The present choice of discretization scheme is founded in the belief that a 'workable ' model can be developed based on this simplest scheme, for the final task of the simulation of large motions of floating bodies in steep waves, since it is supposed that several of the anticipated problems may not necessarily be remedied by applying more refined discretizations but initially may even be obscured by the inherent difficulties of such applications. More refined discretizations can, in principle, be adopted latter. 2.3 The Basic Algorithm of Solution The general discretization scheme described above is applied to the boundary value problem ( 1 ) to ( 5 ) by identifying the open contours Dk with the free surface aDF, the bottom ADD, the body contour DUB and the control boundaries aDC1 and ~DC2. The simulation sought for is accomplished through a time marching procedure. Initial conditions provide the boundary data on ED at t e 0 ~ on IF, B¢/8n on ADD and DUB ~ and either ~ or B¢/8n on ~DC1 and ~Dc2. From the solution of the appropriately rearranged system of linear equations, B¢/3n on IF and on DUB for that time instant are determined. The boundary conditions are then invoked to obtain the required boundary data for the next time level. In particular, appropriate evolution equations for the free surface, deduced from ( 2 ) and ( 3 ), and for the control boundaries, if applicable, are integrated in time to determine the updated free surface contour (i.e. the configuration of IF for the advanced time level) and the values of ~ or In the updated boundaries. On the body, B¢/3n on DUB is related to the body velocity by virtue of (5), which in turn is related to ~ through the body equations of motion and Bernoulli's equation. The boundary contour AD as well as the boundary data for the advanced time level are now established and the solution process can be repeated. The potential ~ at any desired location in D can be calculated from the discretized form of relation (6). Other information, e.g. fluid velocity and pressure are calculable from ~ by utilizing Bernoulli's equation and employing numerical difference techniques in space and time. The evolution of the free surface and the motion of the body, which constitute necessary information for advancing the solution in time, are extracted as the simulation proceeds.

The system of linear equations to be solved, in general corresponds to a full coefficient matrix and thus benefits admissible in solutions of matrices with special features are not available. In the present algorithm, a standard IMSL (International Mathematical and Statistical Library) routine is utilized which employs a Gaussian elimination technique for matrix inversion [25]. The evolution equations for the free surface can be put in the following general form: dot ~ f(y,t) In the present algorithm, a 4th order implicit Adams-Bashforth-Moulton (A-B-M) scheme is adopted and is found to be convergent for all required integrations. To achieve convergence to the third significant figure, usually not more than one corrector step is necessary. This scheme requires information at the preceding four steps. In the initiation of the solution, the first three steps are therefore treated by means of successively lower order schemes with a greater number of iterations. A variety of other schemes exists for the integration of equations of the form (9), e.g. Runge-Kutta schemes, Hamming's method, etc. Fourth order Runge-Kutta starters are popular for analogous initial-value problems (e.g. [1,2 and 13]. However, the starter scheme employed here is found adequate for the applications considered. Limited numerical experiments with other schemes were also performed and the algorithm appeared insensitive to the choice of any particular scheme. 3. Applications to Linear Free Surface Flow Problems 3.1 General Considerations A simple means of partly testing the effectiveness and reliability of the algorithm is to apply it to problems which involve small amplitude free surface elevations. The simulation of the propagation of small amplitude waves allows one to linearize the free surface conditions and therefore reduces the possible sources of numerical errors. It also allows us to examine and solve the remaining numerical problems as discussed below. In addition, solutions of linearized flows are usually available in closed forms thereby providing an exact basis for comparison. Upon applying the usual approximations, the free surface conditions (2) and (3) take the following linearized forms: an an _ = _ at az (10) 354 an Bt ~ ~ (11) and are applied on the undisturbed free surface z ~ 0. In addition, for the applications considered in this section, the fluid domain is represented by the rectangular area depicted in Fig. 1, with the body contour removed. The bottom surface is taken to be at a constant depth d. The free surface part of the boundary on which the integrand in (6) is to be evaluated remains undisturbed at all time instants, and B/3n ~ 3/8z on aDF. The entire boundary aD is therefore independent of time and consequently the influence coefficients in the resulting system of linear equations remains unchanged in time. 3.2 The Simulation of Airy Waves As a test case, the method was applied to simulate the propagation of steady Airy waves in the control domain. In the simulation, the initial values of the potential on the undisturbed free surface z ~ 0 were specified according to the Airy wave potential: /(x,t) ~ H2T co hi[h (Zd/)/ ] sin 2> (x - at) (12) with t ~ 0. This corresponds to a wave of height H. length ~ and period T. progressing in the positive x direction with celerity c. For the linearized simulations, either ~ or a¢/3n, computed from (12), were applied on the control boundaries. The following notation is used in the discussion of computed results: Ax denotes the length of the segments (or the spatial grid size), suffixed appropriately to indicate the part of aD on which they are chosen, viz. AXF, xc1, Axc2 are the segment sizes on aDF, aDC, and aDC2 respectively. The time step size is denoted by At. Nt respresents the time step level of computation: Nt = t/^t. The distance between the control boundaries aDC1 and aDC2, i.e. the horizontal extent of the free surface, is denoted by L. In the linearized simulations, the spatial grid sizes are kept constant on each part of the boundary, and At is constant over the entire time of simulation. As an example, in Figs. 2(a) and (b), the free surface elevations are shown for a simulation where Neumann conditions were imposed on both the control boundaries. The control domain was relatively long, L = 7A, with water depth d = 0. 4) . Other parameters were: AXF, ~XC,, AXC2 = A/2 5 and At = T/40 . The comparison with theoretical free surface contours clearly demonstrates that the algorithm is capable of following the wave motion with acceptable degree of precision over long periods of simulation time. Computations were also performed for a wide variety of combinations of the spatial and

temporal grid sizes, for different values of L/) and d/A, and for different initial distributions of the free surface potential (i.e. initial values of ~ on IF given by (12) with values of t different than 0). In all computations, the quality of agreement between the numerical and theoretical results was similar to the presented examples. The numerical solutions did not exhibit any discernible evidence of degeneration even after, for example, 400 time steps or up to 10 wave periods. 3.3 The Unsteady Wave Propagation The simulation of the propagation of unsteady waves is achieved by specifying an exciting wave potential on one of the control boundaries. The fluid in D is initially at rest with z ~ O as the initial contour of IF. The potential given by (12), which corresponds to an Airy wave propagating in the positive x direction, is applied in a modulated form on aDC1 at all time instants thy. Therefore, a disturbance is introduced at one end of the control domain to excite a fluid motion in the initially unperturbed fluid in D. For this simulation, the boundaries aDC1 and aDC2 can be referred to as upstream and downstream boundaries respectively. For the initially unperturbed state of fluid ¢(x,t)|t=o = 0 in the entire of D, including ED (the value of ~ could strictly be any constant, but it is convenient to make this constant O by redefining ¢, see e.g. [26]). What is not so apparent is the requirement that B¢(x,t)/3tlt=o = 0 be maintained simultaneously. Examining eqn. (3), it can be noticed that I| t=0 ~ O and ¢|~=o ~ O imply a¢/3t|~=o ~ O on aDF. It follows that B¢/3t on aDC1 must have a zero value at t = 0 for the compatibility of the initial boundary data, at the intersection of aDC1 and aDF. The potential given by (12) maintains To ~ 0 on aDC1 but ai/at | t=0 has a finite value. Simulations which did not impose a¢/at | t=0 on aDC1 were not successful due to a numerical instability which started at the origin (at ~Dc1n8DF) and slowly spread over the entire domain. Although this instability is of a weak type, in the sense that the solution can still be continued, the free surface contour shows undesired 'zigzag' patterns and eventually diverges from the required Airy wave profile. In the present formulation, a modified excitation potential ~ is used, defined by introducing a modulation function M(t): ~ (t) = M(t)~(t) (13) with: M(t) = ~ 10.5(1 - cOS(xt/~)) t < ~ (14 where ~ > 0. This function has the property that M(t) It=o ~ O and BM(t)/3t|~=o ~ 0. Therefore, regardless of the form of ~ on aDC1, the initial values of B~ /at are equal to zero by virtue of (13) and: 3~*(t) ~ M(t) halt) + aMatt) i(t) (15) The time span over which the excitation potential is modulated can be controlled by selecting an appropriate value of the parameter p. When (13) is applied with a sufficiently large value of ~ in M(t), the instability disappears. With p/T = 1.0, the existence of some undesired behaviour was still detected. With further increase of p/T to 2, waves evolved smoothly although no numerical smoothing was applied. By inserting (13) in Bernoulli's equation it is found that the modulation excludes an impulsive application of pressure on ~Dc, at too. These considerations appear to parallel a recent study of the impulsive wave-maker problem, [27]. 3.4 Compupted Examples Examples of computed results in terms of the free surface elevations are shown in Fig. 3. For these computations, the downstream boundary was placed at the distance of 2) from the upstream boundary. The discretization parameters were: AXE ~ A/24 and ~t/T ~ 1/36, where ~ and T refer to the length and period of the excitation wave. The downstream boundary was considered to be a rigid wall, thus the condition posed on ~DC2 was B¢(t)/8n = 0 at all time instants. The water had a depth of d = 0.5A, and the excitation potential was modulated over 2T, i.e. p/T = 2. Fig. 3. shows plots of the evolutions of free surface elevations at four stations situated at x = 0.26N, 0.50A, 0.74A and 0.98A, together with the theoretical Airy wave evolutions computed from (12) at the corresponding stations. For comparative purposes, the Airy wave evolutions are also modulated by the same modulation function. It is clear that for t/T < 5, the reflected waves do not reach the locations x/A = 0.98. At locations x/A = 0.24 and 0.50, more than two periods of steady state results are achieved. The algorithm was also applied to simulate wave generation by piston and flapper types of wave-makers. The corresponding Neumann boundary conditions were imposed on a ~DC1 fixed at the mean position of the wave-maker board. In both cases, the normal velocity values were modified by the modulating function (14). Stable propagating waves were simulated. The gain functions (wave amplitude to half stroke ratios) were found to be within 0.2% and 1% relative error when compared to the linear theory values [28]. The presented results exemplify the robustness of the numerical time domain 355

simulation algorithm for fluid flow problems that include a free surface. Computations were performed for a number of other combinations of parameters, and showed a similar quality of agreement with theoretical solutions. No rigorous rule could be established for the minimum size of /\x. As a rough guide, a size of '\x ~ A/12 was found to describe adequately the evolution of the free surface for most of the simulations. Further relaxation resulted in a lack of resolution, although the fluid motion could still be followed (which means the solution did not break down). For temporal grid s ize, the Courant - Friedrichs - Lewy ( C - F- L) type of condition, [ 29 ]: NCFL ~ IC I\XI < 1 (16) with c representing the wave celerity, was used for guidance, although no formal stability analysis of the algorithim was carried out. In the present computations, for most of the cases, a value of NCFL between 0.5 and 0.7 was used. However, the solutions were found to exhibit a tendency towards numerical instability upon successive refinements of the spatial mesh sizes. When a collocation point was located very close to a corner where the boundary undergoes a sharp change in curvature, such as the intersections of ~DF with ~DC1 and ~DC2, relatively large errors occurred in the computed velocities at these locations, in comparison with points far from such corners. These non-uniform differences appeared to introduce numerical instability when the grid size was very fine, typically when Ax/) < 1/100. A similar behaviour of solutions near corners in applications of boundary element methods is documented in literature [24]. This instability is thought not to be a serious limitation in the applications of the algorithm but serves to indicate a lower bound on the grid sizes. In addition to the 4th order A- B -M scheme, other schemes for the integration of eqns. (10) and (11) were also tried. First order implicit schemes were found to be inadequate in that the solution showed poor convergence characteristics and contained large errors. In contrast, 2nd order schemes lead to substantial improvements. Further improvements were achieved by us ing 3rd and 4th order schemes, although the relative improvements between those two latter orders were practically ins ignificant ( the numerical values differed only in the sixth significant figure). In no case was more than one corrector level required for a convergence to the third significant figure. It is observed here that, from a computational point of view, higher order schemes do not necessarily require additional computational effort. To start the 4-th order A-B-M schemes, lower order A-B-M schemes were found to be adequate. 4. The Unsteady Propagation of Steep Waves 4.1 The Evolution Equations for The Free Surface For the s imulation of the propagation of steep waves, it is necessary to consider the non-linear free surface conditions (2) and (3) without the linearizing approximations. These equations are to be satisfied on the exact location of the free surface and therefore the evolution of the free surface within the control domain must be followed. As a result the fluid domain must be redefined at every consecutive time instant. In addition the evolution of the boundary data on the free surface must also be determined. In the present work, an Eulerian formulation of the free surface conditions is used. Assuming the free surface elevation ~ to be a single valued function of coordinate x, the evolutions of ,7 and ~ on the instantaneous free surface are followed at image points of the undisturbed free surface, obtained by the proj ection along the vertical axis . The variation of the potential at the image points on the free surface, which undergo vertical displacements, is determined by (see e.g. [3]): di ~ ~~¢ at + ~¢ d,7 ( 17 ) since ~ ~ ¢(z, t) at these points. Here d,7 is the differential of the vertical displacement of the image point: d'7 = ~ ,7 dt (18) From (2), (3), (17) and (18), the evolution equation for ~ is obtained as: di = -gt7 ~ 21[ (~) ~ (~¢ ) ] 8¢ 3¢ 8?1 (19) which defines the rate of change of the potential at the image points. In the present method, therefore, eqns. (2) and (19) respectively are the evolution equations to be integrated in time in order to determine the instantaneous free surface contour and potential . The above method of following the evolution of free surface is different from the Lagrangian method utilized in most of the previous investigations of non- linear water waves, which were also based upon boundary integral formulations (e . g. [ 1, 5, 11, 13, 30 and 31 ] ) . The attractiveness of the Lagrangian method follows from its ability to describe multivalued free surface contours. In contrast, the applicability of the present method is restricted to single-valued free surface contours. The possibility of simulating extreme wave conditions, pertinent 356

to wave breaking, is therefore excluded. However, the present method provides several computational advantages. The image points on the free surface, which through the discretizaton are identified with collocation points, are not allowed to cluster. So, the scheme avoids the adverse numerical effects often inherent in the Lagrangian methods in which the particles tend to concentrate in some regions. From the previous experience of other workers, it is known that some form of control of the Lagrangian points is necessary to prevent them from clustering, e.g. the introduction of a 'tangential' velocity component as discussed by Baker in [30] or a regrinding of the free surface points at every time step, as employed in [13]. The present method of following the free surface is free from such complications. In addition, the collocation points cannot leave the computational domain at any time, therefore the additional task of tracing such points is avoided. Computational experience gained from performed simulations suggests that in the presented method numerical difficulties arise when a collocation point is situated very close to one of the vertical control boundaries. By preventing horizontal displacements of the collocation points, such problems are also . . . mlnlmlze' .. Yet another point with regard to the applicability of the present method needs to be mentioned. The ultimate objective is to be able to simulate the motion of a floating body in steep waves, for a sufficiently long time, preferably over a number of periods of oscillation. It must be noted that in the Lagrangian method the simulation cannot be extended much beyond the time when the wave breaks. Similar restrictions in applications are typical of most finite-difference algorithms (see e.g. [32]). 4.2 The Simulation Procedure and Wave Excitation The simulation of the unsteady propagation of steep waves is accomplished by a procedure similar to the one described in §3.3. A wave potential, representing an oncoming wave travailing in the positive x direction, is imposed on ~DC1 at t>0. This applied potential is hereinafter called the excitation potential, since it provides the necessary excitation of fluid motion in D. It was found through numerical results that the form of the excitation potential has negligible influence on the generated numerical wave at points in the interior of D, provided they are sufficiently far from the upstream boundary. An alternative way of simulating waves is to provide a moving wave-maker at one end of the control domain. Such a procedure was applied in earlier works, e.g. [10,11,14]. In the present method of following the evolution of the free surface, such an approach would necessitate either a redistribution of the free surface grid or a successive introduction and deletion of collocation points, since the wave- maker enters and withdraws from the free surface grid. It could cause collocation points to come very close to the wave-maker and thus generate numerical difficulties, as it was explained above. 4.3 The Non-reflective Downstream Boundary The application of the impermeability condition at the downstream boundary ~DC2' as it was done in the applications presented above, is not satisfactory for simulations of long duration. An appropriate open boundary, or radiation condition must be specified, which makes the boundary sufficiently transmissive to allow wave patterns generated in D to pass through the boundary without causing appreciable reflection effects. In the present algorithm, a simple open boundary condition is adopted on ~DC2, which assumes that the potential at this boundary can be written as a wave form of the same celerity as that of the applied excitation potential on aDC,: ¢(x,t) ~ ¢(x - at) (20) where c represents the celerity of the excitation wave (cf. eqn. (12)). This results in the following relation: ai ~ - c a~ Bt an (21) in which B/3n - B/3x on ~DC2 was used. Eqn. (21) has a form similar to Orlanski's radiation condition but its application here is not strictly equivalent. In [33] and in many finite-difference algorithms (see e.g. [34]), the value of c is taken as the celerity of the local waves approaching the downstream boundary, and c is determined by a numerical differentiation at the neighbouring grid points. In [35], a similar simple form is adopted with c determined from: c ~ j(gd) where d denotes the local water depth at the downstream boundary. Eqn. (22) represents the shallow water approximation of the phase velocity and is therefore different from the condition applied in the present method (both methods become equivalent in the limiting situation of d/) << 1). (22) The evolution of ~ is now easily determined by the time-integration of eqn. (21), using the same scheme as the one used for integrating (2) and (19). Simple as it appears, in the considered applications this procedure results in minimal reflection effects on ~DC2, as the results presented below illustrate. 4.4 Specific Considerations 4.4.1 Spatial Derivatives on The Free Surface 357

The evolution equations ( 2 ) and ( 19 ) require the evaluation of spatial derivatives of ,7 and at the collocation points. To determine B,7/3x, ,7 as a function of x is approximated by a cubic spline. From this approximation, the components of the outward normal vector can be calculated. For the spatial derivatives of ¢: 8¢ an ~ n as Z ax (23) since 1/(ds/dx) ~ nz, where B/3s denotes the tangential derivative . To determine B¢/3x, in turn ~ as a function of x is approximated by a cubic spline. From B¢/3s and B¢/8n and the outward normal vector, other components of the spatial derivatives can be determined. In the software, an IMSL routine for cubic splines with natural end conditions is used in which no conditions are prescribed at the end points and continuity of second derivatives is enforced at the penultimate points [ 36 ] . 4.4. 2 The Instability at The Intersection of ~DC1 and aDF When a time-modulated excitation Airy potential is applied on ~DC1 in the described algorithm, an instability is found to originate at the intersection of ~DC1 and IF. The form of the instability is qualitatively similar to that observed in the analogous linear application when the modulation function was not used. The solution breaks down, typically within 10 time steps irrespective of the step size . The occurrence of these oscillations can be attributed to the incompatibility of the excitation wave potential and elevation applied externally on ~DC1 with the wave potential and free surface elevation induced in D in the vicinity of IDA. In other words, the free surface conditions implicitly satisfied on the left of ~DC1 (i. e. by the exciting wave potential) are inconsistent with the conditions on ~DF immediately to the right of IDA. On the basis of conducted numerical experiments, this discontinuity is believed to cause large velocity gradients across the vertical boundary and these initiate the instability. The application of other, 'non- linear' excitation potentials, e.g. Stokes second order potential, were tried without success. Difficulties originating from analogous discontinuities were known earlier, e. g. similar problems were discussed in [ 37 ] and in [ 11 ], the discontinuity was identified with the difficulty of the matching of non- linear interior with linear exterior solutions, reported in [ 8 ] . In order to achieve a smooth variation of the free surface elevation and potential across ~DC1, the matching procedure, described below, was developed. Another vertical boundary ED c1 in the interior of the control domain D is introduced at a short distance 1 from IDA. The free surface elevation and potential evaluated without the matching procedure are represented by f2(x). A transfer function g(x) is introduced to redefine f2(x) as f,?(x) in the 'matching zone' between aDC1 and ED cl: f2(x) ~ g(x) f1 (x) (24) where f2(x) is defined in x1 < x < x1 + 1, f1 (x) indicates the theoretical wave elevation or potential corresponding to the excitation potential applied on aDC1, and x1 represents the x coordinate of aDC1. A quadratic polynomial is chosen for g(x) whose coefficients are determined from the conditions: f2 (X1 ) f 1 (X1 ) f2(x1 + 1) f 2(X1 + 1) (25) ax 2(X1 + 1) = aX f2(X1 + 1) The above procedure requires the evaluation of a [ f2 (X1 + 1 ) ] /3x and this is determined from a second order central difference scheme. A linear form of g(x) was found not as satisfactory. The quadratic function for g(x) is however very effective in suppressing the instability and enables the fluid motion to be followed without further problems originating at the intersection point. On the downstream boundary ~DC2 it is necessary to determine the intersection of aDF with ~DC2. This is determined via a second order Lagrangian scheme us ing the data at the three preceding collocation points on aDF. 4.4.3 The Instability on The Free Surface In addition to the instability originating at the upstream end of the free surface, another instability was found to develop on the entire free surface as the solution progressed. Saw-tooth instabilities of this type have been reported by earlier investigators. Numerical experiments with various combinations of the spatial and temporal grid sizes were performed with the hope of identifying a stability criterion related to these discretization parameters. However no such criterion could be established. In the present formulation, in which collocation points on the free surface are restricted to move vertically, the arc lengths between the adjacent collocation points never reduce below the horizontal grid spacing. Consequently, the relation between the time step size and the horizontal projection of the free surface segments is easily controlled. In all computations, the C-F-L type condition, see eqn. (16), was satisfied in the entire fluid domain and throughout the whole s imulation period. A form of stability criterion based upon a linear van Neumann stability analysis for a 4th order Runge-Kutta scheme was provided in [ 13 ]: 358

8 ~XF (26) This condition was also maintained in the discussed computations. The present computational experience indicates that the instability is closely associated with the shape of the free surface elevation. It becomes more pronounced as the wave steepens. It should be observed that in the analogous linear application, no such problem was encountered. Computations with successively higher levels of iteration in the time-integration of the free surface conditions and a close examination of the computed free surface profiles and boundary data suggest an insensitivity of the instability to the time- integration schemes. Therefore, the violation of a stability condition of the above type might not be the root mechanism in the initiation of the instability. To suppress this instability, the smoothing scheme originally employed by Longuet-Higgins and Cokelet [1] was adopted. In view of their computationalexperience, the five-point scheme was used instead of a seven-point scheme. The formula provided in [1] is inapplicable for the special cases of the first and last two collocation points on the free surface and a modified scheme had to be applied there. The smoothing scheme is applied at regular time intervals. Usually, the application at every 4th time step was found effective in eliminating the unwanted oscillations. It is possible to employ a variety of other available smoothing schemes. A third degree five point least squares smoothing scheme was also tested and was found to be equally effective. Although the application of artificial smoothing is known to cause a loss of the local accuracy of the solution, the global solution fortunately remained within acceptable limits of accuracy, which was demonstrated by the computed results. This feature of the artificial smoothing was also noted by earlier workers (e.g. [11, 14 and 38]). 4.5 Computed Results boundary are generally kept fixed in space, except for the uppermost segment. Depending on the length of aDC2, a segment is deleted or an additional segment is introduced so that the length of the segment in comparison with the length of the adjacent segment on aDF maintains a ratio between 0.5 and 2.0. In the present context, the notation described above applies with the exception that xF denotes the spacing of the free surface collocation points along the x axis, instead of the actual lengths of the segments. Unless otherwise specified, the applied excitation potential on aDC1 is the Airy potential. The normalizing parameters for horizontal and vertical length scales and time scale are respectively the length A, height H and period T of the Airy wave corresponding to the excitation potential (cf. eqn. (12)). In all presented computations, Axe, AXc1, ~Xc2 and At are constants. However, due to wave run up, aDC2 continuously changes in length. Since the left hand side of eqn. (21) is an Eulerian time derivative, the collocation points on this Numerical experiments carried out to investigate the effectiveness of the matching procedure showed that the number n of collocation points in the matching zone is the dominant factor in comparison with the length 1 of the zone. The choice of n ~ 4 was very effective in removing the oscillations, regardless of the grid sizes and wave heights. The subsequent results are all computed with this value of n = 4. The convergence characteristics of the solution in the entire fluid domain was studied. The chosen fluid domain extended horizontally over L ~ 2A and vertically over d ~ 0.5A, and the applied excitation corresponded to a large nominal wave steepness of H/> ~ 0.10. For these computations, the segment sizes were: AXF, AXC1 , AXC2 = A/16, A/20, A/24, A/28, A/32 and A/40. The time step size was ~t/T `1/40 for the first five values of N. while for N ~ 120, it was halved to At/T = 1/80. This was necessary because of the reduced segment sizes. Otherwise for At/T = 1/40, which corresponds to NCFL = 1 (cf. eqn. (16) ), the solution exhibited an instability. Free surface smoothing was applied at every fourth step, p/T = 1 was taken in the modulation function, and n = 4 was used for the matching region. Except for the value of N 48, the results demonstrated good convergence characteristics. For this value of N. the results were affected by comparatively poorer resolution. These computations (and many others) indicated that the value of AxF ~ A/24 and comparable values for ~Xc1, AXc2 were adequate for describing the fluid motion without appreciable effects of the lack of resolution. The effectiveness of the open boundary condition (21) was examined by selecting a range of values of the celerity of the outgoing waves. Taking c in eqn. (21) as c', with c' = ac, where c represents the celerity of the exciting wave at aDC1 (as in eqn. (12)), computations covered a variation of ~ from 0 to 1, with specific values of ~ = 0, 0.25, 0.50, 0.75, 0.90 and 1.00. The other parameters were chosen as: L = 2~; d/A ~ 0.5; H/) = 0.10; ~xF, ~xc, and ~xc2 = A/24, and ~t/T = 1/40. The free surface elevations at the simulation time of t/T ~ 8.75 are shown in Fig. 4. It is apparent that the reflection effects at the downstream boundary increase with the difference between ~ and 1. Results for 1 indicate that the wave passes through ~DC2 with imperceptible reflection. From these and other computations the effectiveness of choosing ~ ~ 1.0 for making the downstream boundary transmissive was evident, although 3S9

values slightly less than 1 also appeared to work well. Computations were also attempted for values of c' greater than 1, but even for a value moderately greater than one, e.g. 1.05, the solution broke down after about t ST, which was approximately the time for the wave to grow fully at the downstream boundary. This breakdown resulted from an instability originating at the downstream intersection of IF and ~DC2. In view of the success of 1.0 in making the boundary sufficiently non- reflective, this aspect was not pursued any further. A comparison of evolutions of wave elevation, for H/A ~ 0.10 at x/A - 0.98, with theoretical profiles for the Airy wave given by (12), Stokes second order wave and Miche's second order theoretical profile [39] are shown in Fig . 5 . The numerically s Emulated wave compares well with the second order profiles, but displays stronger non- linear characteristics . The comparatively more peaky crest and shallower trough of the computed wave are visible. To investigate the influence of the excitation potential on the interior solution, computations were performed with a Stokes second order potential as the excitation instead of the Airy wave potential. The applied excitations had a value of H/) ~ 0.10, for which the second order correction in wave amplitude is almost 10% of the first order amplitude, but both excitations have the same average energy dens ity . The fluid domain chosen and the discretization parameters were: L ~ 2>, d ~ 0. 5A, Axe ~ A/24, At ~ T/40 . The evolutions of wave elevation in time at x/) 0.48 and 1.48 are shown in Figs. 6 (a) and (b). The differences between the two simulations are undetectable . The results presented above show that the simulation of unsteady steep wave propagation can be achieved by imposing an excitation potential on one of the vertical control boundaries encompass ing a rectangular fluid domain. The interior solution apppears to be not sensitive to the exact form of the potential. The simulated wave profile displays typical non- linear characteristics of relatively more peaky crest and shallower trough in comparison with linear waves. As expected, the non- linearities are more pronounced for s teeper waves . Very s teep waves were simulated for time durations of over 20 wave periods, for which a steady state behaviour occurred in the entire domain. The simple outgoing wave condition (21) produced good results over the entire period of each computation as well as for all combinations of H. T. L and d for which computations were performed. The interior wave was not observed to be contaminated by numerical reflection effects even after long times of simulation and at locations close to the downs beam boundary . This demons bated the effectiveness of (21) in the modelling of the propagation of nonlinear periodic waves. 5. The Floating Body Problem 5.1 The Governing Equations In order to simulate motions of a floating body in steep waves, a floating body B is introduced in the fluid, such that its wetted contour DUB is completely contained in D ( see Fig. 1) . The obj ective is to expose B to an incident steep wave train and subsequently to follow the induced motion of B. A propagating steep wave is generated in D in the manner described above. For such simulations, it is necessary to know the exact location of DUB at every time instant . In addition, a relation between ~ and B¢/3n on DUB must be established such that the evolution of boundary data on DUB can be followed. The required information is obtained by invoking the equations of rigid body motion and the impermeability condition. For the following cons iterations an additional coordinate system fixed with the body is introduced. A rectangular Cartes fan coordinate system Gx' z ' is used with its origin G located at the body centre of gravity CG, and Gz' axis directed vertically upwards in the undisturbed position of the body. The axis through G perpendicular to the x' and z ' axes is assumed to coincide with a principal axis of inertia of the body. The body geometry is invariant in the coordinate sys tem and therefore the instantaneous wetted contour DUB is completely defined by the location and orientation of Gx'z' system with respect to the Oxz system and wave elevation. The components in the Gx'z' system denoted by {x} ', of the radius vector of a point fixed with the body, are related to the components in the fixed in space system, of the radius vectors of the same point and of CG, denoted by {x) and {XG) respectively, by the following relation: _ ~ _ {X - XG} ~ [R]' {x)' (27) where matrix [R] represents the tensor of rotation and the superscript T indicates a transpose. [R] is given by: [ ] [-sin ~ cos ~ ] (28) where ~ denotes the angular displacement of Gx'z' system with respect to Oxz system, measured as positive counterclockwise. The equations of motion for the body are written in the Newtonian form: {Fx, Fz, My) ~ {M} XG' ME ZG, IS 8) (29) where Fx and Fz are the components of the force F acting upon the body, whereas XG and ZG are 360

those of XG, in x and z directions correspondingly; Me represents the moment of force F about the axis passing through G and orthogonal to Gx' and Gz'; MB denotes the body mass and IF denotes the body moment of inertia about the axis with respect to which M} is defined, and dots indicate differentiation with respect to time. The external forces and moment exerted on B can be obtained by a direct integration of fluid pressure p on ADS: FX ~ T3D B PnXdS (30) Fz = T3D pnzdS - g MB ( 3 1 ) M} ~ T3D B P[ (Z ZG)nX+(X XG)nZ]dS (32) where nx and nz denote the components of the normal vector on DUB in the fixed in space system of reference. The pressure p is computed from unsteady Bernoulli's equation: = -P{gZ + B~ + 21 (Vi)23 (33) On the body surface, the fluid normal velocity B¢/3n is equal to the normal component of the body velocity by virtue of impermeability condition (5). Therefore, from (5) and (27), the following form of the impermeability condition is obtained: `~¢ = nX XG + nz ZG + ~ [ - (Z-ZG)nX+(X-XG)nZ] (34) to be satisfied on ODE. The above expression provides B¢/3n at any point on DUB in terms of the body configuration, velocity and geometry. 5.2 The Algorithm for The Computation of Hydrodynamic Forces The solution algorithm described in §2. 3 can now be adapted for the simulation of motions of B. At any instant, presuming B¢/8n to be known on ADD, the other boundary data are determined from the solution of the integral relation (6). From this, the fluid pressure p exerted on the body is determined by utilizing Bernoulli's equation (33). Subsequently, the fluid excitation loads on B are determined by a direct integration of pressure over the wetted body surface, (30) to (32). To establish the boundary data for the next time level, the equations of motion are invoked. By integrating (29), XG ~ VG ~ ~ and ~ at the advanced time are determined. The necessary boundary data (in) on the body surface are then established from (34) and the computation for the next time step can begin. The use of Bernoulli's equation (33) requires the evaluation of B¢/3t on the body. To this end the following relation: di Bi + ~ 3¢ (35) at at ax is applied, where d/dt signifies the rate of change of any quantity following a material point of the rigid body and v denotes the velocity of the point. The quadratic term in Bernoulli's equation is readily obtainable from: (V¢)2 = (~n)2 + (~52 (36) At a collocation point i the tangential derivative B/3s of ¢, is determined in the form: (pi = (~-¢i)i/~)i = ASP (~¢i)i (37) since for the straight line segments, (Bs/Bi) = ASP which denotes the length of the segment. To determine (/i), appropriate second-order difference formulae are employed. This is permissible despite large variations in DUB, since ~ on the surface is in general a smoothly varying continuous function, which was confirmed by plotting he against i for several conditions. In the present algorithm, Simpson's rule was applied for the evaluation of force and moment expressions (30) to (32). To carry out the computations on the wetted body surface aDB, it is convenient to describe the body geometry with respect to the Gx'z' system, in which it is invariant. Denoted by D'B, it can be subdivided into segments once and for all. To determine the wetted contour DUB (BOB ~ ~D'B), the intersection points ~DF n DUB need to be found. This is achieved using an extrapoltion scheme in which a second order polynomial is fitted to the three points on ~DF adjacent to DAB- However, the application of a fixed discretization of DAB through the determination of DAB n ~DF and the subsequent consideration of only those segments of DAB which belong to BOB, produced an instability in the force computation and a resulting divergence of the solution. This was caused by the introduction or deletion of segments on the body near the intersection point with the free surface, which in turn produced a large variation of pressure in the computation of the dynamic part of the pressure distribution. This problem is overcome by redistributing the collocation points on the body at every time step such that the segment sizes vary smoothly in time. Because of the redistribution of the collocation points, a spatial interpolation of becomes necessary in the computation of d¢/dt. In general, this spatial interpolation introduces very small approximation errors, since the changes of collocation points between two consecutive time steps are very small. 361

Computations using linear and second order interpolation rules produced practically indistinguishable results. In the present algorithm, the second order rule is applied. integrates the equations of motion using the explicit scheme. In this way, it is possible to use a stable difference formula for the dynamic pressure term in the force calculation, in exchange for an explicit scheme of the integration of the equations of motion. On the free surface, a similar redistribution of the collocation points becomes necessary. Except for a wall-sided body in heave motion, any other combination of the body geometry and modes of motion causes a change in the size of the segments adjacent to the body, eventually leading to a deletion or introduction of a collocation point. For the same reason as on the body, this destabilizes the force computations. Therefore the original locations of the collocation points cannot be retained, and collocation points must be redistributed at every time step. This necessitates an interpolation of and ~ in space for the integration of the free surface evolution equations. Instead of storing and interpolating between ~ and ~ values, the values of the right hand sides of (2) and (19) are stored for the required number of past time steps (four steps for the integration scheme employed) and interpolated. This results in a slight reduction of the computation. A cubic spline interpolation is employed. With regard to the associated- numerical inaccuracies, the above comments on the discretization of DUB apply. 6.1 The Experimental Program 5.3 The Integration of The Equations of Motion In order to integrate them in time, the equations of motion (29) are transformed into six ordinary differential equations of the first order: · · · · · — {UG, XG' VG, ZG, A, ~ } it, UG, it, VG, ~i, A) (38) A number of standard techniques are available for integrating the above system of equations. For convenience, as well as to be consistent with the integration of the free surface conditions, a 4th order A-B-M scheme was originally used. However, the solution was divergent, and this could not be remedied by increasing the number of iterations per time step. The problem was found to be caused by the computation of the d¢/dt term. When an implicit scheme is used for the equations of motion, this term can only be computed from a backward difference scheme in time in the corrector part of the algorithm and this leads to the instability. In the present algorithm, the instability is avoided by computing the predicted values by using a backward difference scheme for d¢/dt in the force evaluation, and an explicit scheme for integrating the equations of motion. For the second and higher iterations, the scheme returns to the preceding time step and corrects the forces, this time utilizing a central difference scheme for d¢/dt, and once more Although the adopted procedure was more stable and not divergent at the initial time steps, after sufficient time, in particular when surge and/or roll motions were involved, another 'saw-tooth' type instability was found in the time variation of the forces, velocities and displacements. These appeared first in the force computations and gradually contaminated the velocities and displacements. The mechanism which initiates this instability was not identified, but the singularity at the body and free surface intersection point, associated with larger horizontal velocity components in these modes, is probably a factor. This problem is circumvented by smoothing the forces in time. The smoothing formulae applied are the same as those used on the free surface (notice that the smoothing here is in time). For the chosen integration rules, forces over only the past four steps need to be smoothed. 6. The Evaluation Against Experimental Data An experimental program was undertaken to validate the numerical model. This was considered necessary due to the inadequacy of published analytical, numericalor experimental data on two-dimensional motions of floating bodies in steep waves. To the authors knowledge, no systematic two-dimensional experimental data are readily available in open literature, in which a floating body is subjected to an incident wave train such that the motions and/or waves contain significant non-linear characteristics. The only exception appears to be the experiment by Kyozuka [40] who conducted an experiment on a body of Lewis form by subjecting it to oncoming waves and presented the results in the frequency domain. The experiments were performed in the Memorial University wave tank which has interior dimensions of 54.74 m x 4.8 m x 3.04 m and is equipped with a piston type wave-maker and a parabolic beach (for more details on the tank, see [41]). The arrangement where the body is completely free to float was considered not to be favourable for experimental purposes, since it would be difficult to prevent the body from undergoing motions in the transverse plane (yaw and pitch). This lead to the choice of experiments in which the body was restricted from swaying. The body chosen for testing was of rectangular 40 cm x 40 cm cross section, and its length was 120 cm. To avoid sharp corners and minimize resulting flow separation, a bilge radius of 2.5 cm was provided. The body was ballasted to a draft of 20 cm. To achieve two- dimensionality of the flow, a channel within the wave tank was constructed by erecting 362

vertical walls. The channel length was 6.1 m. For the period for which test results were collected, reflected waves from either end of the tank did not reach the test site. Experiments were performed in a water depth of 0.9 m. The wave field was monitored using standard twin wire wave probes of resistor type. The collected data were therefore the wave heights, the horizontal force exerted on the body and roll and heave displacements of the body. Prior to the testing of body motions, a series of preliminary tests were performed in which waves were generated and wave heights were measured at four locations along the centre line of the channel. One of them coincided with the designed location of the body. The purpose of these tests was to determine the range of frequencies and heights for which waves of acceptable quality could be generated in the channel. Also, the data generated were needed for the comparison with the simulation. The tests with the model were conducted in two series. In one the model was free to heave but restrained from rolling. The roll moment was not measured in these tests. In the other, the model was free to heave and roll but restrained from swaying. The steepest waves generated in the preliminary wave test series could not be applied in the main series of the tests because they caused significant flooding. Altogether 39 tests were performed covering the range of wave periods 0.054 ~ w j(b/2g) 1.08, with ~ ~ 2~/T and b denoting the breadth of the body cross-section, and of wave steepness 0.068 < H/> < 0.013. The tests with the body free to roll were carried out for two radii of gyration 0.29b and 0.36b, and for GM/b = 0.054. For all body motion tests, data were recorded which included the transient information. Except for the horizontal force measurement, all other measurements (wave heights, heave and roll motions) had insignificant amount of noise content. The noise was removed by a five point averaging technique. For the force measurement, however, the time records contained a relatively large proportion of high frequency noise. The noise was subsequently removed by applying a digital filtering technique. A majority of the tests was conducted twice to verify repeatability of the tests. The results showed very good repeatability, with the exception of the sway forces. 6.2 The Comparison of Experimental and Numerical Results 6.2.1 General Comments In order to establish a proper basis for the comparison of the experimental data with the corresponding results of computation, it is necessary to replicate the experimental wave conditions in the numerical simulation as close as possbile. Previous computed results (cf. §4) showed that the height of the generated waves inside the control domain is closely comparable to the height (H) of the Airy wave describing the excitation potential. Also, the fundamental period of the generated wave was shown to be very close to the excitation period. Ideally, the excitation potential should be selected such that the waves simulated at the location of the body match the test wave conditions but this would lead to a trial-and-error search for the correct excitation potential. Considering the simulation time and the number of experimental conditions that were to be simulated, this process would have been prohibitively expensive and time consuming. Instead, the application of the Airy excitation potential with H and T taken from the test conditions was found to the give close enough simulations of experimental waves. All of the numerical results, used in the comparisons with the experimental data, unless specifically mentioned, were computed for a standardized control domain extending over L = 4.0~. The CG of the body at rest was located at L1 ~ 2.5) from the boundary ~DC1. The two values of time variable T1 and T2 indicate approximately the time at which the fully developed wave pattern reached the location of the body (T1), and the time at which the reflected from the body wave pattern reached the excitation boundary (T2). The values provide a guidance for the time interval within which the comparisons are meaningful. For the size of the domain chosen, about 3 to 5 wave periods could be obtained within the interval between T1 to T2 for most of the tests. Considering the wide range of wavelengths for which numerical results were to be generated, the grid sizes were not standardized. The spatial grid sizes were chosen such that a reasonably good description of the entire boundary can be obtained, At for each computation was chosen such that it satisfied the condition given by eqn. (16). The other relevant parameters were: §/T ~ l and n = 4 for all computations (this corresponded to the matching length between 0.1) to 0.13~). The shape of the body in the numerical model was rectangular. The comparison between the numerical and experimental time records were presented in the following manner. The numerical wave elevation was monitored at a station inside the control domain, located at a distance of 0.5) from the excitation boundary. The station was about 0.4) outside of the matching region on the average, therefore the time history for the free surface elevation remained uninfluenced by reflections from the body for a relatively long time. For the experimental input conditions, the wave elevation records from the preliminary test series without the presence of the body were taken as the experimental 363

input conditions. The comparison of these two records provided the comparison of the oncoming wave conditions. For presentation, these two records were synchronized. The synchronization was standardized by matching the peak of the simulated wave in the time interval 3<t/T<4 with the peak of the experimental wave record. The primary outputs of the experimental and numerical model were the two displacements, heave and roll, and the horizontal (sway) force. These results were plotted after their records were synchronized with respect to the undisturbed wave pattern at the horizontal location of the body CG at rest. The synchronization was achieved by first presenting the experimental records relative to the undisturbed wave pattern at the CG location. The undisturbed wave pattern was obtained from a reference probe record and the phase difference between the reference probe and the probe at the location of the body CG, measured in the absence of the model. Subsequently the records of the responses determined by the numerical model were adjusted for the phase difference between the experimental and simulated wave patterns at the location of the reference probe. The error inherent in this procedure was estimated within + 15 deg. 6.2.2 The Comparison of Sway Forces and Heave Motions A detailed review of the comparison of the experimental and computed records is beyond the scope of the present discussion, only a brief summary and some examples can be presented here. With the exception of three instances over the whole experimental range, the sway force differences of peak-to-peak values were contained within +10% with respect to the experimental data, for the steady responses. However, relatively large, on the average +20°, phase differences were observed, larger for the model restrained from rolling, on the average +26°. It should be noticed that the experimental peak-to-peak values were assessed to be repeatable within not less than +4% error. The differences of steady peak-to-peak heave responses, between the experiment and computation were well below +5% for the model free to roll and occassionally above that value for the model restrained from rolling. The phase differences were within +15 deg. with the exception of three tests. In the range of heave resonance the numerical values consistently overpredicted the experimental data. The repeatability of the experimental data was estimated within +2%. 6.3 The Comparison of Roll Motions With the exception of the lowest experimental frequency in the vicinity of roll natural frequency, the body displayed very small roll displacements in all other tests. The roll amplitudes were mostly less than 4 deg. The numerical method predicted a similar behaviour. Large roll amplitudes were obtained experimentally for the frequency w7(b/2g) = 0.54 with the radius of gyration 0.033b. At this frequency, experimental data was gathered for three different wave steepnesses, H/) = 0.013, 0.023 and 0.028. A larger steepness could not be achieved due to the limitation of the dynamometer (maximum allowable roll of + 30 deg.) as well as due to water spilling inside the body. For these tests, although the agreement of measured and computed phases was very good, with phase differences within +15 deg. the roll amplitudes were considerably over-predicted, with peak-to-peak values differing by between 17% to 32%. The peak-to- peak sway forces were under-predicted on the average by 19%. The heave motions correlated well, with peak-to-peak values over-predicted on the average by 5% and approximately zero phase difference. The over-predictions of roll amplitudes is believed to result from effects of fluid viscosity which were not accounted for in the numerical model. The significant influence of fluid viscosity on the damping of roll is well documented in literature (see e.g. [42]). Therefore, the incorporation in the numerical model of the viscous effects upon roll damping is expected to improve the predictions. In order to include the viscous damping of roll in the numerical model, the equation of roll motion was modified to the form: I`8 = Me - bed) (39) where the second term on the right-hand side of (39) represents the contribution of viscous damping to the roll moment and is considered as a function of angular speed 8. It should be observed that moment Me includes the effects of roll damping resulting from wave scattering. The expression for the viscous damping moment was identified from free rolling numerical and phys ical experiments in the form: . . . . bib) = BIB + B21818 ( ) where B1 and B2 denote constant coefficients derived from the identification. In the identification the assumption was made that the damping contribution of wave scattering to the roll moment is linear with respect to 8. 364

The numerical results obtained by the described above procedure showed a significant improvement of the predicted roll motion. The changes of computed sway force and heave motion were marginal. Figs. 7 and 8 show the comparisons of the computed and measured records, obtained without and with the inclusion of viscous roll damping respectively. The results are for the largest wave steepness H/> = 0.028 at which the worst agreement between the experiment and computation was observed. For all three wave steepnesses the differences between computed and measured peak- to-peak roll values reamined below 6% for the simulations in which eqn. (39) was applied. 6.4 Summary Taking into account the experimental errors and inaccuracies of the comparison resulting from the unsteady characteristics of the experimental and computed data, the predictions by the numerical model show good agreement with experimental records. The observed differences between computed and measured sway forces and heave motions result, at least partly, from the absence of viscous effects in the numerical model. A similar remark applies to the roll motions were the inclusion of the viscous roll damping is semi-empirical. It should be observed that the experimental data used in the comparison included three kinds of non-negligible non-linear phenomena. For the shorter waves with higher steepness, the waves near the body approached the breaking limit. In a number of tests, a foam formation was observed. Large heave and relative motions of the body with respect to waves were observed in several tests. For these tests, the video records showed that the run-up profiles resembled closely those generated by the simulation. In addition, the tests included the occurrence of roll amplitudes up to 20 deg. combined with large motions relative to waves. The comparison confirmed the validity of the numerical method in the presence of the described phenomena. It also demonstrated that the method can provide realistic estimates of roll motions within the specified range if a semi-empirical model of viscous roll damping is included in the algorithm. Conclusion The numerical results presented above show that the propagation of steep long-crested periodic waves, which may include mild transients, can be modelled numerically in the time domain by means of a relatively simple, low order boundary element method algorithm. The efficiency of a radiation condition of Orlanski's type is also demonstrated. The use of the condition makes possible the avoidance of other means of eliminating the reflection of waves at the downstream boundary, such as e.g. the introduction of an artificial wave damping in the free surface conditions, which are not compatible with open water conditions. Similar comments apply to the wave excitation in the control fluid domain, achieved in the algorithm by a matching of an imposed exciting wave potential with the flow in the control domain, in the vicinity of the upstream control boundary. The procedure is different from the more usual simulation of a physical wave-maker, and corresponds better to open water conditions. However, in the procedure the exciting wave potential must be modulated in time to enforce the compatibility of the initial conditions at the boundary and in the control domain. The problem parallels similar difficulties encountered in the modelling of the wave-maker. Another characteristic feature of the algorithm constitutes the use of an Eulerian form of the nonlinear free surface conditions based on the assumption of a single-valued wave elevation. At present this form of the conditions appears to be more suitable for the extended in time modelling of motions of floating bodies in steep waves, than its Lagrangian counterpart, although the latter is capable of modelling wave breaking. The application of the Eulerian free surface conditions makes possible a strict observance of local stability condition of Courant type, without the use of procedures (such as e.g. a regridding of collocation points) which may introduce a numerical smoothing. The results presented above suggest that the violation of a Courant type stability condition may not be the reason of the occurrence of the typical instability in the free surface data. The extension of the basic steep wave propagation algorithm, which includes the presence of a free floating body in the control fluid domain, provided time domain simulations of body motions. In the simulations steady state motions of the body excited by steep periodic waves, preceded by short transients, were achieved. The computed records compare well with experimental data, thus illustrating the applicability of the extended algorithm. The experimental data were obtained from specially performed model tests, with a body model restrained in the sway mode. In several of those tests significant nonlinear phenomena related to the wave propagation and interaction with the body were observed. In the numerical simulations no significant wave reflection at the downstream boundary was detected. However, for the form of the algorithm presented here, simulation times are limited by the reflection of waves scattered by the body from the upstream control boundary. To ensure the stability of the computation of hydrodynamic forces and of the integration of the quations of body motion, special procedures were developed which include a smoothing of the hydrodynamic forces in time. 365

References 1. Longuet-Higgens , M. S . and Cokelet , E . D ., "The deformation of steep surface waves on water: I. a numerical method of computation." Proc. R. Soc. Land., Series A, 350, pp. 1-26 (1976) . 2. Faltinsen, O.M., "Numerical solution of transient nonlinear free - surface motion outside or inside moving bodies. " Proc. 2nd Int. Conf. Num. Ship Hydro., Berkeley, CA, pp . 347 - 357 (1977) ~ 3. Faltinsen, O . M ., "A numerical nonlinear method of sloshing in tanks with two dimensional flow. " J. Ship Res., 22, 3, pp. 193 - 202 (1978) . 4. Vinj e , T . and Brevig, P ., "Numerical simulation of breaking waves. " Adv. Water Resources, 4, pp. 77-82 (1981) . 5. Vinj e, T. and Brevig, P., "Numerical calculation of forces from breaking waves. " Int. Sym Hydro. Ocean Eng., Norwegian Inst. Tech., Trondheim, pp. 547-565 (1981). 19 6. Brevig, P., Greenhow, M. and Vinje , T., " Extreme wave forces on submerged energy devices." Appld. Ocean Res., 4, 4, pp. 219- 225 (1982). 7. Vinj e, T. and Brevig, P., "Nonlinear ship motions. " Proc. 3rd Int. Conf. Num. Ship Hydro., Paris, pp. 257-268 (1981). 8. Vinj e , T ., Maogang, X. and Brevig, P ., "A numerical approach to nonlinear ship motion. " Proc. 14th ONR Symp. Naval Hydro., National Academy Press, pp. 245-278 (1982). 9. Greenhow, M ., Vinj e, T ., Brevig, P . and Taylor, J., "A theoretical and experimental study of the capsize of Salter's duck in extreme waves. " J. Fluid Mech., 118, pp. 221- 239 (1982) . 10. Lin, W. M., "Nonlinear motion of the free surface near a moving body. " Ph.D. Thesis, 23 Dept . Ocean Eng ., MIT , 127 p . (1984) . 11. Lin, W. M., Newman, J. N. and Yue, D.K.P., "Nonlinear forced motions of floating bodies. " Proc. 15th ONR Symp. Naval Hydro ., National Academy Press, Washington, D. C., pp . 33-49 (1989) . 12. Dommermuth, D.G. and Yue, D.K.P., "Study of nonlinear axisymmetric body-wave interactions. " Proc. 16th ONR Symp. Naval Hydro., National Academic Press, Washington, D . C ., pp . 116 - 136 (1986) . 25. 13. Dommermuth, D . G . and Yue, D . K . P ., "Numerical simulations of non- linear axisymmetric flows with a free surface. " J. Fluid Mech., 178, pp. 195-219 (1987) . 14. Greenhow, M. and Lin, W. M ., "Numerical simulation of free surface flows generated by wedge-entry and wave-maker motions. " Proc. 4th Int. Conf. Num. Ship Hydro., National Academy of Sciences, Washington, D.C. pp. 94-106 (1985). 15. Greenhow, M., "Water entry and exit of a horizontal cylinder." Proc. Second International Workshop on Water Waves and Floating Bodies, Rep. No. AM-87-6, Univ. Bristol, pp. 25-28 (1987) . 16. Greenhow, M., "Wedge entry into initially calm water. " Appld. Ocean Res., At press (1988) . 17. Isaacson, M. St . Q ., "Nonlinear wave - effects on fixed and floating bodies. " J. Fluid Mech., 120, pp . 267 - 281, also Corrigendum, 133, pp. 469 (1982). 18. Isaacson, M. St. Q., "Steep wave forces on large offshore structures." Soc. Petrol Eng. J., 23, 1, pp. 184-190 (1983). Krzyzanski, M., "Partial Differential Equations of Second Order. " Vol. 1, Polish Scientific Publishers (1971). 20. Hess, J. L., "Higher order numerical solution of the integral equation for the two-dimensional Neumann problem. " J. Comp. Meth. Appld. Mech. Eng., 2, pp . 1 - 15 (1973) . __ . Wardle, L. J ., "An Introduction Boundary Element Method. " Numerical solution of partial differential equations, J . Noye ea., North-Holland, pp . 289 - 312 (1982) . 22. Hess, J. L., "The use of higher order surface singularity distributions to obtain improved potential flow solutions for the two-dimensional lifting airfoils. " J. Comp. Meth. Appld. Mech. Eng., 5, pp. 11- 35 (1975). Breit, S. R., Newman, J. N. and Sclavounos, P. D ., "A new generation of panel programs for radiation- diffraction problems . " Proc . 3rd Int. Conf. Beh. Off. Str. (BOSS), Elsevier Sc. Pub. B.V., pp. 531-544 (1985). 24. Schultz, W.W., "A complex-valued integral method for free surfaces with intersecting bodies. " Proc. Second International Workshop on Water Waves and Floating Bodies, Rep. No. AM-87-06, Univ. Bristol, pp . 101 - 103 (1987) . Forsythe, G. and Moler, C. B., "Computer solution of linear algebraic equations. " Chapter 9, Prentice-Hall, Inc., Englewood, Cliff ., N . J . (1967) . 06. Lamb, Sir H., "Hydrodynamics. " 6th Printing, Dover Pub ., (1945) . 366

27. Cointe, R., Jami, A. and Molin, B., "Nonlinear impulsive problems." Proc. Second Int. Workshop on Water Waves and Floating Bodies, Rep. No. AM-87-06, Univ. Bristol, pp. 13-16 (1987). 28. Biesel, F., "Les Appareils Generateures de houle en Laboratoire." La Houille Blanche, 6, pp. 147-165 (1951). 29. Roache, P.,"Computational Fluid Dynamics." Hermosa Pub. (1972). 30. Baker, J. R., Meiron, D. I. and Orszag, S.A., "Applications of a generalized vortex method to nonlinear free-surface flows." Proc. 3rd Int. Conf. Num. Ship Hydro., Paris, pp. 179-191 (1981). 31. New, A. L., McIver, P. and Peregrine, D. H., "Computation of overturning waves." J. Fluid Mech., 150, pp. 233-251 (1985). 32. Telste, J. C., "Calculation of fluid motion resulting from large amplitude forced heave motion of a two-dimensional cylinder on a free surface." Proc. 4th Int. Conf. Num. Ship Hydro., National Academy of Sciences, pp. 81-93 (1985). 33. Orlanski, J., "A simple boundary condition for unbounded hyperbolic flows." J. Comp. Phy., 21, pp. 251-269 (1976). -`r. Chan, R.K-C. and Chan, F.W-K., "Numerical solution of transient and steady free surface flows about a ship of generalized hull shape." Proc. 13th ONR Symp. Naval Hydro., Tokyo, pp. 257-280 (1980). 35. Wu, D-M. and Wu, T.Y., "Three dimensional nonlinear long waves due to a moving surface pressure." Proc. 14th ONR Symp. Naval Hydro., National Academy Press, pp. 103-125 (1982). 40. Kyozuka, Y., "Experimental study on second- order forces acting on cylindrical bodies in waves." Proc. 14th Symp. Naval Hydro., Tokyo, pp. 319-382 (1982). 41. Muggeridge, D.B. and Murray, J.J., "Calibration of a 58 m. wave flume." Can. J. Civil Eng., 8, 4, pp. 449-455 (1981). 42. Himeno, Y., "Prediction of ship roll damping -a state of the art." University of Michigan Dept. of Naval Arch. Rep. No. 239, 65 pp (1981). Z UNDISTURBED FREE SURFACE aDF r— —%\ ~ aDF D ~DD aDB Fig. 1 Definition diagram. ~ F - ~ - T HCTHO0. t'T" .. O · PRC3CHS MCTH=, t/T. ~ O. O ~ T - ~rIC". t/T - . O. I~ O 0.00 1.00 aDc2 ~5. Ahlberg, J., Nilson, E. and Walsh, H., "The . theory of splines and their applications." (a) Free surface elevat~ons Academic Press, New York. ( 1967 ) . . Han, F. S. and Stansby, P.K., "On the application of the boundary element method to two-dimensional free surface interactions with bodies." Proc. Second International Workshop on Water Waves and Floating Bodies, Report. No. AM-87-06, Univ. of Bristol, pp. 39 -42 ( 1987 ) . 38. Dold, J.W. and Peregrine, D.H., "An efficient boundary-integral method for steep unsteady water waves." In Numerical Methods for Fluid Dynamics II.Ed.K.W. Morton and M.J. Baines, Oxford Univ. Press., pp. 671- 679 ( 1986 ) . 39. Miche, R., "Mouvements ondulatoires de la mer en profondeur croissante ou decroissante." Annales des Ponts et Chausses (1944). ~ PR£SENT METHOD, o ~ THEORET ~ CRL, _ _ _ ~ ~ ~ ~ ~ 0.00 2.00 4.00 6.00 8.00 10.00 t/T (b) Evolution of the free surface in time at the center of the domain ~ = 3.5) Propagation of a linear steady progessive wave: (L/A = 7, d/) = 0.4, ~XF/) = 1/25 and At/T = 1/40). 367

x present method, x/A = 0.26 theoretical, =/A = 0.26 present method, x/)t—0.50 theoretical, :~/A = 0.50 _~ , . . . t.00 2.0~ 3.~0 0.00 x present method, x/) = 0.74 o .~. theoretical, z/) = 0.74 _ ~ present method, ~/A = 0.98 ~ theoretical, :z /) = 0.98 o 0 "L' _ ·~ o' t/T I ~_~ o o , , ~ ~ ~ , ~ 0.00 1.00 2.00 3.00 ~ OQ 5.~0 t/T 1. 00 2. 00 3. 00 Fig. 3 The evolution in time of a free surface elevation started from rest. o _ 0 _ ~ _ o~ o ~ = _ I ~ J I I 0.0 ~ 0.2s x O.SO 0.7s O. 90 1. 00 \~ //~' f:: 0.00 0.2s 0.50 0.75 1.00 1.25 1.50 1.75 . 2.00 z/) Fig. 4 Free surface elevations for different values of a: (L-2A, d/A-0.5, H/A'O.1, AxF/A ~ 1/24 and ~t/T:1/40). 368

Fig. 5 Of (a)ati= 12(~/) = 0.48) _~ o~ lo o— X PRESENT METHOD MICHE ~ STOKES 2 Ad OROEfl tI) RIRY an _ to a, q 1 1 4.se 4.7s s.oo t/T boo 5.25 5.50 Comparison of theoretical and numerical time histories of wave elevation. Airy wave pot.on ~Dc1 Stokes 2nd order potion dice - ¢=o~--it~l~M ~AM/ o o i I ~ ~ I ~ ~ ~ ~ 0.00 1.00 2.00 3.00 4.OQ 5.00 6.00 7.00 8.00 9.00 t/T To ° (b)ati= 36(~/) = 1.48) -A ADAM . . . 0.00 1.00 2.~O 3.00 4.00 5.00 6.00 7.00 8.00 9.00 t/r Fig. 6 The evolution in time of free surface elevations induced by different excitation potentials. 369

~ incident WAVE _ ~ h: $' Fig. 7 Comparison of experiment and theory for w~b/2g 2 O. 54 and H/,'` 0.028: (T1~14.5 sec and T2~21.0 see). 370

~ incident .~ - 16- 16- 1~- 1~- Ik. ^- ^- ^- ~- t [ " hew motlen S .- ~ 16- IL~ 1~- ~ 1~. ^- ~- ~- ~- t ~ a ~ Fig. 8 Comparison of theory and experiment for 07b/2g ~ 0.54 and H/) = 0.028 when roll viscous damping is included in the simulation: (T1~14.5 sac and T2~21.0 see). 371

DISCUSSION by R. Cointe I would like to congratulate the authors for a very interesting and extensive work. We are now working on a similar problem and we find it very difficult to accurately compute the free motions of the body using finite differencesin time to evaluate the 3~/at term in Bernoullie's equation. An alternative method consists in evaluating this term by solving the corresponding boundary integral equation. Would the authors discuss their ex- perience on this problem? Author's Reply The work presented by Dr. Cointe at this conference[A1], addresses several problems which are either identical or analogous to the ones we have addressed in our paper. The main differences between the two works result from the assumed directions of approach, with the "numerical flume" approach taken by Dr. Cointe, and "open water" approach chosen by us. As a result certain comparisons between the achieved results may be instructive. In particular I believe, it is worthwhile to notice that the wave damping condition employed in [Al] at the down stream boundary, is tuned to the dominant wave frequency, and our application of the radiation condition depends on the use of an appropriate wave celerity. It could also be interesting to see how an introduction of the compatibility of the initial condition on the exciting boundary and free surface, which we impose through a modulation function, would affect the performance of Dr. Cointe's algorithms. Referring directly to Dr. Cointe's question, we could not avoid the necessity to solve the boundary value problem twice in order to advance the simulation of motion of a floating body in time. The technique used is described in section 5.3, and it relies on using a central difference scheme to obtain time derivatives of the potential on the body surface. We very __ _~ question and wish him success in finalizing the development of his model. milch ~nnrmm; ~~ - Do Hi of - ~= [Al] R. Cointe, "Nonlinear Simulation of Transient Free Surface Flows", 5th International Conference on Numerical ShipHydrodynamics,Hiroshima, 1989. DISCUSSION by R.C. Ertekin The authors should be commended for a very careful and detailed analysis carried out in their paper. It was a pleasure to read it. By using the plane speed obtained from the potential for Airy Waves, you appear to be assuming, at least implicitly, that downstream waves are linear. And by setting c'= ac and varying a, aren't you basically determining the value of the phase speed which you would have if you numerically calculated it? We have recently solved the problem of diffraction of nonlinear waves by 2- dimensional submerged objects by using the BEM. The results will be presented in the next OMAE conference in Houston (1990, 9th meeting). We have used Crank-Nicholson scheme in time stepping and found no saw-tooth waves that you have seen in your results by using the A-B-M method. We have not filtered our results. This makes we believe that your time stepping algorithm is causing the problem. It is not very unheard of that a multistep method is associated with such difficulties (Burden & Faires, "Numerical Analysis", 1985, Prindle & Weber). One way of handling this difficulty is by using adaptive time stepping in the algorithm. I would like to also point out that Shapiro, R. ("Linear filtering", Mats. Comput., Vol. 29, 1975) presented filtering formulas for any number of points earlier. Finally, could you explain why sway forces could not have been repeated in the experiments? Author's Reply The comments by Professor Ertekin are appreciated. considering the formulation presented in the paper, the upstream boundary can be used as a permeable boundary for oncoming waves which satisfy fully the non- linear free surface conditions. However, applications on that boundary, of wave potentials which do not conform to the conditions, are equivalent to wave excitations by a wave maker. As much as the kinematics of a wave maker board does not imply the linearity of the generated waves, since they satisfy the non-linear free surface conditions, the application of an Airy wave potential on the upstream boundary does not imply the linearity of the wave generated in the control fluid domain. A wave excitation 372

of this type, by a wave potential which in- cludes a wave celerity c in its definition, makes it convenient to use the same celerity in the radiation condition. This procedure appeared to be effective in the performed simulations. Otherwise the celerity would have to be determined numerically from the solution in the control fluid domain. The source of the saw-tooth instability in the wave propagation simulation is not entirely clear. In our simulations it was related to the use of the non-linear free 373 surface conditions and to wave steepness, and therefore the accuracy of the computation of the spatial derivatives in the free surface conditions seems to be the probable source of the instability. Finally, the measurements of sway forces were repeated in all the tests carried out twice. However the peak-to-peak values from those measurements showed a variation of at least +4% between the original and repeated test records. The cause of these variations was not fully determined.