DAVID W. TAYLOR NAVAL SHIP; RESEARCH AND DEVELOPMENT CENTER Bethesda, Maryland 20084-5000 DTNSRDC-85/086 CALCULATION OF TWO-DIMENSIONAL NONLINEAR FLUID FLOW RESULTING FROM LARGE-AMPLITUDE FORCED HEAVE MOTION OF A U-SHAPED CYLINDER INA FREE SURFACE by John G. Telste APPROVED FOR PUBLIC RELEASE; DISTRIBUTION IS UNLIMITED. COMPUTATION, MATHEMATICS AND LOGISTICS DEPARTMENT RESEARCH AND DEVELOPMENT REPORT ATION OF TWO-DIMENSIONAL NONLINEAR FLUID FLOW RESULTING FROM LARGE- UDE FORCED HEAVE MOTION OF A U-SHAPED CYLINDER IN A FREE SURFACE December 1985 MAJOR DTNSRDC ORGANIZATIONAL COMPONENTS —_— DTINSROC COMMANDER TECHNICAL DI Ree OFFICER-IN-CHARGE | ANNAPOLIS OFFICER-IN-CHARGE CARDEROCK SHIP SYSTEMS INTEGRATION DEPARTMENT 12 AVIATION AND SURFACE EFFECTS DEPARTMENT SHIP PERFORMANCE DEPARTMENT 15 COMPUTATION, MATHEMATICS AND LOGISTICS Pepe en STRUCTURES DEPARTMENT PROPULSION AND AUXILIARY SYSTEMS DEPARTMENT SHIP ACOUSTICS DEPARTMENT SHIP MATERIALS CENTRAL INSTRUMENTATION DEPARTMENT ENGINEERING DEPARTMENT GPO 866 993 NDW-DTNSRDC 3960/43 (Rev. UNCLASSIFIED ee e — ——E —————e SECURITY CLASSIFICATION OF THIS PAGE REPORT DOCUMENTATION PAGE Ja. REPORT SECURITY CLASSIFICATION 1b. RESTRICTIVE MARKINGS UNCLASSIFIED Pee ragd seat/s he oe a 2a. SECURITY CLASSIFICATION AUTHORITY 3 DISTRIBUTION/ AVAILABILITY OF REPORT a APPROVED FOR PUBLIC RELEASE; \2b. DECLASSIFICATION / DOWNGRADING SCHEDULE DISTRIBUTION IS UNLIMITED. 4 PERFORMING ORGANIZATION REPORT NUMBER(S) 5 MONITORING ORGANIZATION REPORT NUMBER DINSRDC-85/08 = Ce 6a NAME OF PERFORMING ORGANIZATION David Taylor Naval Ship Research and Development Center 6b OFFICE SYMBOL (If applicable) Code 1843 7a. NAME OF MONITORING ORGANIZATION 6c ADDRESS (City, State, and ZIP Code) 7b. ADDRESS (City. State, and ZIP Code) Bethesda, MD 20084-5000 Ba. NAME OF FUNDING/ SPONSORING ORGANIZATION 8b. OFFICE SYMBOL (If applicable) } 9. PROCUREMENT INSTRUMENT IDENTIFICATION NUMBER Be ADDRESS (City, State, and ZIP Code) 10 SOURCE OF FUNDING NUMBERS PROGRAM PROJECT TASK ELEMENT NO NO NO 61152N WORK UNIT Ast 38083 DN 8 TITLE (Include Secur CALCULATION Of CRA BitiNs TONAL NONLINEAR FLUID FLOW RESULTING FROM LARGE-AMPLITUDE FORCED HEAVE MOTION OF A U-SHAPED CYLINDER IN A FREE SURFACE 13b. TIME COVERED 14 DATE OF REPORT (Year, Month, Day) |15 PAGE COUNT FROM TO 1985, December Sl ZROLLO201 11 12 PERSONAL AUTHOR(S) 18 SUBJECT TERMS (Continue on reverse if necessary and identify by block number) Free Surface Effects, Heaving Cylinder, Perturbation Theory, Ship Motions, Ship Hydrodynamics, Water Waves 19 ABSTRACT (Continue on reverse if necessary and identify by block number) Progress in developing a tool to compute large amplitude ship motions is reported. In particular, a method to calculate transient two-dimensional potential flow about a body moving in a free surface is described. The flow problem is formulated as an initial boundary value problem in which the velocity potential along the free surface and the positions of the moving boundaries are sought as solutions of a coupled system of differential equations. An implicit finite-difference method is used to advance the solution of the coupled system of equations in time. The auxiliary problem of computing the velocity potential inside the fluid region is solved by a method which is based on boundary-fitted coordinates and is directly extensible to three-dimensional flows. Results from calculating the potential flow about a body in forced heave motion are presented. The hydrodynamic force on the body has been obtained and compared with the hydrodynamic force predicted from second-order perturbation theory. 21 ABSTRACT SECURITY CLASSIFICATION UNCLASSIFIED D2b TELEPHONE (include Area Code) ees (202) 227-1929 Code 1843 83 APR edition may be used until exhausted SECURITY CLASSIFICATION OF THIS PAGE All other editions are obsolete ON CHEANSIS ET ICT ID)” : a 20 DISTRIBUTION/ AVAILABILITY OF ABSTRACT (J UNCLASSIFIED/UNLIMITED (1) SAME AS RPT 22a NAME OF RESPONSIBLE INDIVIDUAL John G. Telste DD FORM 1473, 84 MAR C1 oric USERS UNCLASSIFIED SECURITY CLASSIFICATION OF THIS PAGE UNCLASSIFIED OO SECURITY CLASSIFICATION OF THIS PAGE TABLE OF CONTENTS EST OF FIGURES) (ey joy cp coh pow ellen teb © ved lojilon toyiel ¥ cease) ays RABIGE ie Gs) © 0: se) 6) (0) fob teltiet fo} cmiclhiciiet (6 Vell lol Jou io oh folie NOTATION) \) js: “s) (0) 10) 0) teltionte: fefttel lot elt toy) Wes felh ol fel! (oNNell le) te IABSIERAGT. (7s) (ou se op ote) fon olbioly ow letlioMiielajejaseuiceseNelhie! kop torptoll fe ADMINISTRATIVE INFORMATION . «© « © « «© « © © © © © e@ LNERODUGEEON circ) le) ie) fotbton tolate) feltfe) follis) (ollie) .omhoMlen teh (01) (0 MATHEMATICAL FORMULATION . 2. « « «© «© «© © © © © © @ NUMERICAL SCHEME . . « «© « «© «© © © © © © © © «© © © « BOUNDARY FUNCTIONS . . « « « © «© « © © © © © © © «© LAPLACE SOLVER . « «© @ « © « eee « «© = & © « RESULTS ie) ie) 0, ie) ce) fotetot NeW otMolntoliden (oom oh loNbiolh oh Kon, oMlts GCONCIUSTON, ‘eo je) lo: io Tomiotstettiol ©. fotiteision lon tefl ok lofnte) fobiton te REFERENCES 0, ye! 6)" tee ce He! (0) el cer sel) Tepe: 8) ken fe) er 18, eo e ee LIST OF FIGURES 1 - U-Shaped Body in the Free Surface at Time t = 0.0 2 - Geometric Description of the Moving Body Contour 3 - Fluid Region for Flow Symmetric about the y—Axis 4 — Subregions into Which the Fluid Region is Divided 5 - Mesh Size in Each Subregion of the Computational Region 6 — Initial Mesh Near the Body .... +--+ « « « « 7 - Free-Surface Elevations Near Body for t 54, 5-6, 8 - Free-Surface Elevations Near Body for t = 9.0, 9.2, @eey @ ee, 9 — Free-Surface Positions in a Frame of Reference Fixed to t = Oler25 64. eee, 9.0 cee to ib ett ese eh 0 e. 0 Ve strlal Wil o& 6 e e the Body for Page 25 26 27 28 29 30 Sill 32 33 Page 10 - Free-Surface Positions in a Frame of Reference Fixed to the Body for t = oes Oa eeey No? es e e e e e e e e e e e e es e e e e e e es e e 34 11 - Comparison of Computed Nonlinear Force with Linear Force and with Second-Order Force as Functions of Time t . . « »« » »« «© © «© ee »© © © © 39 12 - First- and Second-Order Forces versus Frequency Number for U-Shaped Cylinder e e es e e e e e e e e e e e e e e e e e e e e e e e s e e e e 36 13 - Phase Angle of Forces of First- and Second-Order versus Frequency Number for U=-Shaped Gyliinder <6 < 5 © «6 © « \s) 0) s) 1) «lo 0) os eon TABLE 1 - First- and Second-Order Force Coefficients . . »« « « » « » « » « e 22 iv AB BC By, (a,t),By,(a,t) b CD ED j,k,2 NOTATION Constant used in describing the body contour Contour of moving body Symmetry boundary of fluid region Angle parameterizing the body contour A function used to denote the initial relation between the parameters a and e Value of the parameter a at the intersection of the free surface and the body contour Free-surface boundary of fluid region Functions describing the moving body contour Half-beam of the body Far vertical boundary of fluid region Bottom boundary of fluid region Parameter describing free surface and wetted body contour; perturbation parameter = hg/b Dimensionless force = F'/L302p Dimensional force First-order vertical force on heaving body, normalized by 2 2pgb ho/b Second-order sinkage on heaving body, normalized by 2pgb“(ho/b) Second-order vertical force on heaving body, normalized by 2pgb*(hg/b)? Dummy variable used to describe a convergence criterion Gravitational acceleration Depth of the rectangular fluid region Amplitude of forced heave motion Location of a grid point xp(e,t),yp(e,t) xp(e,t),yp(e,t) X0»Y0 u,V Vy(a,t) Length scale Number of body grid points Time levels Number of free-surface grid points contour Outward directed unit normal on (By(a,t),By(a,t)) body Origin of Cartesian coordinate system Dimensionless pressure = p'/po2L2 Dimensional pressure Arclength along the body contour Dimensionless time = t'o Dimensional time Dimensionless Cartesian coordinates = (x',y')/ L Dimensional Cartesian coordinates Functions describing position of body contour Functions describing free-surface position Coordinates of moving center of body Velocity of a fluid particle; u = ox, V = by Normal component of body-contour velocity outward) at (xp(a,t), yp(a,t)) Half-width of the rectangular fluid region Functions of (&,n) describing the mapping of the region onto the computational region Laplace operator Time increment Phase angle of first-order force on heaving body Phase angle of second-order force on heaving body vi at fluid particles along (positive fluid Parameters used to test for convergence of procedure Coordinates in computational region Dimensional wavelength Fluid density Frequency of heave motion Dimensionless velocity potential = $'/oL2 Dimensional velocity potential Dimensionless velocity potential along the free vii iterative surface svtjaip2! " ie, Vt “ saa tine io se | 1 I | : it 7 - 7 i il) teoach echle - humbe pr of body qetd potore " gun agi wwnoy 3 sa07° 02 nee #189 Aan 7 Ai fine Lavelle ‘SHms5074 alia Peete tengdseamgnen pent x20 1390 bag tn sg ig taindb bhai Oecreh WE ERT COA TAR 588M Set Pat pr Moma? $00 484 HS FGF Eo Pad * ie nok awamia on verd frente Ay (A t), My Casbd) hime hatane | APC! GRETA sheets see Ito 1 (EWA m0 ocak Pineneiontibea 1 lar = 0'@ wens iene) {yee % ~ = p : hii pean 3 i? B2Sn ct r LTate t? im ' area af , Ki iphone i te Loine i y of a hi : ' ’ His nog tha ler ' ‘yi : t i } l x? y } iy il Of ; } L] i vy t ENG | j ( i wo : | f ‘1 sin ; } t ye ae a ‘ t ‘ ae { \ edi eet ar Bore wi hha PEISNG2 6 Cisulev San LAr ty % aT, and ’ uf ee ia jig body ABSTRACT Progress in developing a tool to compute large amplitude ship motions is reported. In particular, a method to calculate transient two-dimensional potential flow about a body moving in a free surface is described. The flow problem is formulated as an initial-boundary value problem in which the velocity potential along the free surface and the positions of the moving boundaries are sought as solutions of a coupled system of differential equations. An implicit finite-difference method is used to advance the solution of the coupled system of equations in time. The auxiliary problem of computing the velocity potential inside the fluid region is solved by a method which is based on boundary-—fitted coordinates and is directly extensible to three- dimensional flows. Results from calculating the potential flow about a body in forced heave motion are presented. The hydrodynamic force on the body has been obtained and compared with the hydrodynamic force predicted from second-order perturbation theory. ADMINISTRATIVE INFORMATION This work was supported by the Numerical Ship Hydrodynamics Program at the David W. Taylor Naval Ship Research and Development Center. This Program is jointly sponsored by the Office of Naval Research under contract NOOO1484AF00001, NR-334-001, Work Unit 1542-019, and by DTNSRDC under its Independent hesearch Program, Task Area ZRO140201, Program Element 61152N, and Work Unit 1843-045. INTRODUCTION In recent years attention of naval architects and ocean engineers’ has focused on how vessels and offshore structures react to large-amplitude ocean waves. Th attention has been motivated by the capsizing of vessels in large breaking waves and structural failure due to the slamming forces associated with such waves. It is therefore of interest to have a method available to determir the forces on a floating body and how the body will react in these extreme conditions. Since little is known of how a ship reacts even to non- breaking waves, it would be extremely valuable to have a method available for predicting ship motions in the presence of large-amplitude nonbreaking waves. The method would be of practical value in a systematic study of how ship design changes would add stability in rough seas. Researchers have spent much effort in devising methods for computing large-amplitude ship motions. Their work is usually based on the assumptions that the fluid is incompressible and the fluid motion is irrotational. The assumptions lead to the existence of a velocity potential, which simplifies the problem formulation. But, since the equations describing the free surface are nonlinear and cannot be linearized for large-amplitude waves, great difficulties arise in the computation of solutions to free-surface potential flow problems. When a body is present in the free surface, additional difficulties related to the intersection of the free surface and the body Occur. For instance, the potential flow in the region is known to be singular. Lin et al. [1]* have recently described some aspects of the singularity. Dagan and Tulin [2] and Fernandez [3] discuss nonlinearities in fluid flow about blunt bodies. Because of the formidable difficulties, most of the work on nonlinear free surface flows has been for two-dimensional (2-D) flows. Many authors have formulated the 2-D problem as an initial-boundary value problem whose solution is obtained from an integral equation for functions defined along the boundaries of the fluid. Longuet-Higgins and Cokelet [4] used such a method combined with a time-stepping procedure to calculate free- surface heights with no body present. Faltinsen [5] and Vinje and Brevig [6-8] *A complete listing of references is given on pages 39 and 40. have used the integral-equation approach in their studies of fluid motion in the presence of bodies. They make the restrictive assumption that the fluid domain is periodic and use complex-variable techniques that cannot be extended to three-dimensional problems. Greenhow et al. [9] have applied the method of Vinje and Brevig to the capsizing of a body in the free surface. Baker et al. [10] have developed a generalized vortex integral-equation technique that has been used for a body under the free surface. It is not clear whether the generalized vortex method is suitable for numerical computations when a _ body intersects the free surface or whether it will be computationally efficient when it is extended to three dimensions. Thus, even for computing nonlinear two-dimensional free-surface potential flows, full generality has not been attained. This report describes progress made in developing a tool to compute large-amplitude ship motions. The method discussed is based on an initial- boundary value formulation, and it is a method that is directly extensible to three dimensions. The velocity potential along the free surface and the posi- tions of the moving boundaries are sought as solutions of a coupled system of differential equations. An implicit finite-difference method is used to march the solution of the coupled system of equations forward in time. The auxiliary problem of computing the velocity potential inside the fluid region is solved by a finite-difference method based on boundary-fitted coordinates. Haussling [11] has presented a review of such techniques used for fluid flow problems. Results from calculating the potential flow about a cylinder in forced heave motion are presented. The hydrodynamic force on the body has been ob- tained and compared with the hydrodynamic force predicted from ‘second-order perturbation theory. MATHEMATICAL FORMULATION The physical flow problem is to compute transient two-dimensional flow about a body moving in a free surface. It is formulated mathematically as a potential flow problem in which the velocity potential along the free surface and the position of the moving free surface are sought as the solution to a nonlinear initial-boundary value problem. In particular, the physical problem is to determine the fluid motion caused by the prescribed movement of a body partially submerged in a fluid and the resulting hydrodynamic force on the body. The prescribed motions are forced harmonic heave motions never so large that the body becomes completely submerged in or rises out of the fluid. The body considered in this report is a closed U-shaped cylinder. Gravity is the only body force acting on the fluid which is inviscid, incompressible, and initially at wrest. The fluid motion is irrotational and thus a velocity potential »' is assumed to exist. Surface tension is neglected. All variables are nondimensionalized. Lengths are scaled by a length L characterizing the size of the body. Time is scaled by 1/o where o is the frequency of the body motion in radians per second. Velocities are scaled by OL; the velocity potential ', by oL2; pressure, by po2L2; and force, by po2L3 . Here p is the fluid density. Thus, for example: xe ES Ih Sep We Ib ye ENS e/e OO Agi pl eS edtiga,, PY Ss polis, where (x,y) are variables representing the coordinate system, t is time, p is pressure, and F is force. The primed variables represent dimensional quantities; the nonprimed variables, nondimensional quantities. The fluid region is described in terms of a fixed (x,y)-coordinate system chosen so that the y-axis points vertically upward and the undisturbed free surface is at y = 0 (Fig. 1). A fluid region of infinite depth and infinite lateral extent is modeled by a rectangular tank so deep that the effect of the bottom boundary is insignificant and so wide that no waves reflect from the side boundaries during the time for which the fluid motion is modeled. The rectangular tank is bounded by the lines y=-h, x =w, and x = -w. The contour of the body moving in the free surface is given as a function of time t and a parameter, a, by the equations 56 = LE((igis) = IN (eos (a) = Ocil cos (3 a))) (la) y = By(a,t) = SAN (sin Ca) RTO isin @)a)))— hip) cos! "Ce) (1b) where A is a measure of the size of the body, a is the angle measured counterclockwise from the direction of the positive x-axis (Fig. 2), and hg is the amplitude of the heave motion. The position of the free surface is given in terms of a parameter e and the time t by x = xp(e,t) and y = yp(e,t). The functions xp(e,t) and yp(e,t) are to be calculated. Since the flows considered in this paper are symmetric about the y-axis, only the half of the fluid region where x > 0 is considered (Fig. 3). The region is bounded by five curves. Across the boundaries AE (x = 0), CD (x = w), and ED (y = -h) there is no flow. The curved line AB, given by Equations (la,b), is the contour of the moving body. BC is the free surface, whose location must be computed. The velocity potential satisfies the Laplace equation Ad = 0 (2) in the fluid region and is subject to certain boundary conditions. (Fluid velocity in the x-direction is given by u= 9,3; fluid velocity in the y- direction, by v = bys) At a solid boundary, the normal velocity of the fluid must equal the normal velocity of the solid boundary since fluid cannot penetrate the boundary and no cavities are assumed to form in the fluid. In particular, at stationary boundaries the normal velocity must vanish. At the right vertical boundary of the flow domain, about 16 half-beams away from the body, this condition is given by dx = O at x =w (3) Similarly, at the bottom boundary, about six half-beams below the free surface, vanishing normal velocity is specified by the equation Oy ae Oaays = sal (4) At the boundary AE directly beneath the body, where a symmetry condition is specified as a wall condition, the velocity potential must satisfy the equation dy = 0 at x = 0 (5) Along the body contour, the normal velocity is known from the _ prescribed motion of the body. In fact, the normal velocity of the body at (B,(a,t), By(a,t)) is vp(a,t) = nydBy/dt + nydBy/at where n(a,t) = (ny ny) is the unit normal directed into the fluid given by (n, ny) = (OBy/2a, 8B, /3a)/[(aB,/3a)? + (aB,/aa)2) 1/2 The required boundary condition for » at (B,(a,t), By(a,t)) on the body contour is thus ~ = oxNy + dyny = nydBy/dt + nydBy/dt (6) The free-surface coordinates xp(e,t) and yr(e,t) have been parameterized in terms of e. The parameterization e is chosen such that for fixed e the functions xp(e,t) and yp(e,t) describe the path of a fluid particle. In other words, e is a Lagrangian variable. The velocity potential on the free surface is also parameterized in terms of the Lagrangian variable e by > = op(e,t). An equation to be satisfied by $p(e,t) can be obtained from Bernoulli's equation, which can be expressed as p + d¢/at + ($2 + $2)/2 + (g/o2L)y = 0 (7) x y where p is pressure, g is the acceleration of gravity, o is the frequency of the forced harmonic heaving, and L is the characteristic length. Bernoulli's equation is valid on the free surface and throughout the fluid region. At the free surface, the pressure is assumed to be zero and a particular case of Bernoulli's equation, the dynamic free-surface boundary condition, results: d D ey (ex) = = (xp(e,t),yp(e,t),t) = (4% + 62)/2 - (g/o7L) yp(e,t) (8) dt Dt x y Here D/Dt = $,0/dx + byd/ay + 0/dt is the derivative following the motion of a fluid particle. The kinematic free-surface boundary condition, which states that no fluid particle on the free surface can leave the free surface, is expressed by the two equations dxy, Dx = (e,t) = Se (xp(e,t), yp(e,t), t) = $,, (xp(e,t), yp(e,t), t) (9) dt Dt SVE (ai, 4) te Gea(S, ON FMC), 6) pitty) = Galen (ey Oe yA(GUOL amEGIED ae ’ a F ’ r) YFP ’ ) y \*F ’ 9 YF i 9 At t = 0, the velocity potential on the free surface is given by op(e,t=0) = 0 (11) and the free surface is such that yp(e,t=0) = 0 Gi2)) The parameter e and the function xp(e,t) can be arranged so that xp(e,t=0) =e (13) It is convenient to parameterize the fluid at the body contour by a Lagrangian variable e so that x = xp(e,t) and y =yp(e,t) along this boundary. This is possible since fluid particles along the solid boundary can never leave that boundary, except possibly to become free-surface particles at the intersection of the body and the free surface. Thus d D B (e,t) = — (xple,t), yg(est), ¢) dt Dt Ut $,, (xp(e,t), yp(e,t), t) (14) d D —E (e,t) = — (aplest)s yg(ert)» £) = dy(xp(est), yp(ert), ©) (15) dt subject to the initial conditions xp(e,0) = By (a(e),0) (16) yp(e,0) = By(a(e),0) (7) where a(e) is a prescribed function. In summary, we seek the solution of an initial-boundary value problem for tal I xp(e,t), y = yr(e,t), = op(e,t) on the free surface and x = xp(e,t), y = yp(e,t) on the body contour, in which e is a Lagrangian variable that parameterizes the free surface and the body contour. These five functions obey the evolution Equations (8), (9), (10), (14), and (15) subject to the initial conditions (11), (12), (13), (16), and (17). The vellocity potential in these equations must satisfy Equation (2) subject to the Neumann boundary conditions (3), (4), (5), and (6) and a Dirichlet boundary condition along the free surface governed by Equation (8). The force on the body is calculated by integrating the pressure over the wetted surface of the body. Because the flow problem is symmetric about x = 0, the x-component of the force on the body vanishes: Fy = 0 (18) The y-component of the force, positive upward, is given by Fy = 7) p(a,t) ny(a,t) (ds/da) da (19) where (n, »Ny ) is the unit normal at the body contour directed into the fluid and s, increasing in the counterclockwise direction along the body, represents arclength. Because of symmetry, the pressure is integrated over the half the wetted length of the body given by x = B,(a,t) and y = By(a,t) where - 1/2 < a ¢ ap(t). (The function ap(t) depends on the position of the intersection of the free surface and the body contour.) NUMERICAL SCHEME The functions xp(e,t), yple,t), op(e,t) along the free surface and xp(e,t), yple,t) along the body contour obey five coupled first order differential equations in time with specified initial conditions. In addition, the velocity potential must satisfy specified boundary conditions at all times. To solve these equations a finite-difference method is used. BOUNDARY FUNCTIONS Each of the five boundary functions is discretized with respect to time and space using a fixed time step At. The discretized forms of the functions are denoted by “i = xp(ej, ndt) (20a) me = yp(e,;, ndt) (20b) iy = $p(ej, nat) (20c) fortes Te. Nee add aa = xu(e,, nAe) (20d) 10 ne = yp(e,, nAt) (20e) for k=1, ... ,M. (The subscript j will be used for a free surface variable and the subscript k for a variable along the body contour.) Thus, the boundary functions are discretized into boundary grid points that are now assumed to move as if they were associated with fluid particles. The evolution Equations (8), (9), and (10) for op(e,t), xp(e,t), yp(e,t) and the evolution Equations (14) and (15) for xp(e,t), yple,t) are applicable to these particles and are replaced by finite difference equations based on the Euler- modified method. The finite-difference equations are given by x (n+1) = x(n) + At[o(™) + p(ntl) ]/2 (21a) Fj Fj xj xj yntl) = y(m) + arg) + g(mt1)]/2 (21b) Fj Fj yj yj g(nt1) = ¢(n) + At p(n) D, + p(n) 2 + p (ntl) 2 + g (ntl) 2 2 Fj Fj xj yj xj yj - (g/10?) "s ‘ oD f2 (als) Fj Fj x(nt+1) = x() + At ee + ue /2 (21d) Bk Bk xk xk y(ntl) = y(n) + at Pe 2 al /2 (2le) Bk Bk yk yk where o(m) = > (x(m), y(m), mt) xh x GL GL o(m) = > (x(m), y(m), mat) yh y GL GL 11 for 2 = j or k, m =n or ntl, and G = B or F. The initial conditions become x(0) = oe (22a) Fj j y69) = 0 (22b) Fj (0) = 0 (22c) Fj x(0) = x (a(e ), 0) (22d) Bk B k yee vacate): 0) (22e) In Equations (2la-e) and Equations (22a-e), the subscript j runs over all possible free surface grid points, and the subscript k runs over all possible body grid points. The intersection of the free surface and the body is treated as a body point obeying Equations (2ld,e) subject to an initial condition given by Equations (22d,e). The numerical scheme is implicit. An initial estimate for the five functions at the (n+l)-st time step is obtained by linearly extrapalating from two previous time steps. (For the first time step, the initial estimate is the initial condition.) The functions are corrected iteratively by using Equations (2la-e). The iterative procedure is stopped when the x- and y- values have satisfied an absolute error criterion of the form le ee) a g(n+],£+1) | < Ey @23)) and the ¢-values along the free surface have satisfied the relative error criterion given by lee pont, 2) 7g (ntl, 241) | < €5 whenever [gS oly | > €3 12 The superscripts (n+1,2) and (n+1,2+1) in Equations (23) and (24) refer to the 2-th and (2+1)-st corrected solutions of the system of differential equations at time step n+l. Generally, the stopping criteria for the iterative procedure in the Euler-modified method at each time step are that the maximum absolute change in the x- and y-values on the free surface and the body is €; = 0.001 and the maximum relative change in the 4-values on the free surface is €7 = 0.001 for d-values greater in absolute value than €3 = 0.000001. Each of the five discrete evolution equations has 6, or dy on the right side. In particular, to compute the right sides of Equations (2la-e), 6, and dy along the moving boundaries of the fluid region must be computed from the solution of the Laplace equation for 6. The solution method for this equation is described in the next section. To prevent numerical instabilities on the computed free surface from arising, a linear filtering scheme due to Shapiro [12] is used. The filtering scheme has been used successfully by Ohring and Telste [13], Haussling and Coleman [14], and several other researchers. It has been applied at fixed intervals of time, and has been especially helpful near the intersection of the free surface with the body contour. Unless some method is used to maintain a reasonable grid spacing along the free surface and the body contour, grid points will congregate in some areas and become sparse in other areas. A method is used to keep a uniform distribution of points along the body contour and a prescribed distribution of grid points along the free surface. The prescribed free-surface distribution is such that the free-surface length between the body and the j-th free- surface grid point is a constant fraction of the total free-surface length 13 between the body and the outer boundary. The scheme allows one to follow the movement of boundary grid points within a time step as if they were fluid particles and to shift the grid points to other fluid particles at the end of each time step. To shift the grid points, cubic spline interpolation is used to fit the arclength as a function of grid point number along the free surface and the body contour. The positions of the grid points along the free surface and the hull are shifted and the values of all pertinent functions interpolated, using the cubic splines, to their new values. The redistribution of grid points, of course, affects the initial guess for the Euler-modified method at time step n+l. However, if the shifting is done every time step and the time step is sufficiently small, the redistribution scheme has been found to proceed smoothly. Because of numerical errors, grid points cannot be expected to remain exactly on the body as the solution of the initial-boundary value problem is advanced in time. Numerical errors arise from the redistribution scheme and from replacing the differential equations by finite difference equations. To correct for such errors, grid points that move off the body are shifted back to the _ body. This is accomplished by computing the counterclockwise angle about the center of the body from the direction of the positive x-axis (iio Do Every body grid point at a location slightly off the body surface at a certain value of that angle is relocated to the point on the body having the same value of that angle. At the intersection of the free surface and the body contour difficulties arise. At this point both the free-surface boundary conditions and the boundary conditions associated with the solid boundary apply at a _ particular instant of time. However, the fluid particle at the intersection may move to 14 a position on the free surface or may move to a position on the body below the surface. In other words, the intersection point may not move with the fluid. Since the free-surface and hull boundary conditions all involve time derivatives following the fluid motion, they cannot, in general, be directly applied over a time interval to predict quantities at the intersection point. Special methods for handling this point must be developed. However, for heaving motions of an almost wall-sided body, the intersection will be essentially a fluid particle. Therefore for the current study, the body Equations (14) and (15) are applied directly at this point. Inaccuracies in this approach will become apparent in the form of a deviation from zero of the pressure at the intersection point. In fact, such pressure deviations might be used in a method to more accurately follow the intersection point as would be necessary for more complicated body shapes. One such scheme has_ been tested but has proved to be numerically unstable. The pressure at all the grid points along the body including the intersection is calculated from a finite-difference version of Bernoulli's (n+1/2) = (n)]} 2 (n)} 2 (n+1)} 2 (n+1)} 2\ 74 an je i Psa] ‘ ts 4 E , - (g/o2L) oe + “| /2 (25) equation: k The pressure is the pressure at the k-th fluid particle on the body at time t = (mtl/2)) At. The force on the body at this time is calculated by numerically integrating the pressure along the body contour. | Trapezoidal quadrature is used. 15 LAPLACE SOLVER To solve the Laplace equation for 9, a finite difference method based on boundary fitted coordinates, due to Thompson et al. [15], is chosen. The finite-difference method involves mapping the time-dependent fluid region onto a fixed computational region. The coordinates € and n in the computational region are such that they obey the Laplace equation with x and y as dependent variables: I oO Exx at Eyy (26a) ise oe iy = (26b) The boundary conditions for & and n along a given boundary are Dirichlet if a particular mesh distribution along the boundary is prescribed. The boundary condition of one coordinate is Dirichlet and that of the other coordinate is Neumann if mesh orthogonality near a particular baundary is desired. Since all computations are to be done in the fixed (&,n)-computational region, it is convenient to interchange the independent and dependent variables. When this is done, the Equations (26a,b) for & and n become ax - 2Bx + yx = 0 (27a) Bis En nn CHV Nee SZ) EO) (27b) ata En nn where a=x2 + y2 (27c) n nN Sy St ass PSL SY (27d) ew En y = x2 + y2 (27e) E & 16 To obtain a mesh that wraps around the body and conforms to the other boundaries, the physical fluid region is divided into several subregions (Fig. 4). Each subregion is mapped onto a rectangle of computational space. In each rectangle of computational space, the inverted Laplace equation is solved subject to Dirichlet boundary conditions, transformed Neumann boundary conditions, or matching boundary conditions where rectangles overlap. The Laplace equation for the velocity potential @ transforms exactly as Equations (26a,b). The transformed Laplace equation is given by oe) > A) ap yf) = (27£) BE En nn where a, 8, Y, are defined by Equations (27c-e). Wherever possible, central differencing is used to discretize Equations (27a-£). At the boundaries of the fluid region where a Neumann boundary condition is specified, second-order one-sided finite differences replace some of the derivatives in Equations (27a-f). (See Coleman and Haussling [16] for more details.) The resulting system of quasilinear equations for x, y, and 9 at the grid points is solved by successive overrelaxation. Lin et al. [1] cite the works of various researchers who show that a logarithmic singularity exists in the velocity potential of free-surface potential flow near the intersection of the free surface with a vertical wavemaker in horizontal motion. But, since the body contour for the problem considered is nearly vertical at the intersection with the free surface and since the horizontal velocity component of the body is zero, the logarithmic singularity in the velocity potential may be relatively weak. Thus special numerical treatment of the singularity may not be critical. In fact, nothing special has been included in the numerical method to accommodate such a singularity if it should arise. 17 RESULTS Several forced harmonic heave motions of the U-shaped body in the free surface have been considered. The shape of the body has been fixed by setting the parameter A in Equations (la,b) to 0.7407 in all cases. With A set to this value, the half-beam of the body and the draft of the body are both 0.6667. Amplitudes considered were ho/b = 0.05, 0.3, and 0.4, in which hg is the amplitude of the motion and b is the half-beam of the body. Most researchers consider values of the frequency parameter bo2/g that lie between 0.0 and 2.0. In this study the frequency parameter is restricted to lie in the interval from 1.5 to 2.0 since this interval contains the frequencies for which nonlinear effects are greater. Linear theory for the problem of an oscillating body in a fluid of infinite depth and lateral extent predicts that the wavelength far from the body will approach h/L = 2ng/Lo2 = 2n(b/L)g/bo2 asymptotically in time [17]. The dimensions of the rectangular tank are taken to be about one such wavelength deep and four wavelengths in half-length. Thus, if the frequency parameter bo2/g is no smaller than 1, the depth h should be about 4 and the length should be about 16. This region is long enough that no waves will reach the side boundary during the first few periods of forced motion. Since the region is more than half a wavelength deep, the effects of finite depth can be ignored. The fluid region has been divided into five time-dependent subregions. Each subregion has been mapped onto a fixed rectangle of computational space (Figs. 4, 5). The size of the mesh in each of these subregions has _ been depicted in Figure 5. A total of 2916 grid points, counting twice those 18 points where two subregions overlap, has been used for the grids that cover the computational and fluid regions. There are 34 equally spaced grid points on the half-body contour. At the far right vertical boundary, at the bottom horizontal boundary, and at the vertical boundary below the body, grid points position themselves in a way that makes the mesh near these boundaries orthogonal. At the far right boundary there are 25 grid points; at the bottom boundary, 102 grid points; on the vertical boundary directly beneath the body, where a symmetry condition is specified as a wall condition, 18 points. The initial distribution of grid points on the free surface is arranged so that the mesh near the intersection of the free surface and the body is approximately uniform in both the x- and y-directions (Fig. 6), and the mesh near the intersection of the free surface and the far right boundary is also approximately uniform in both directions. There are 90 grid points along the free surface. The iteration for the solution to the Laplace equations for the changing mesh and for the velocity potential is done simultaneously. It is found that the mesh does not change greatly from time step to time step and hence a _ test for the convergence of 6 suffices. This convergence test is that the square root of the sum of the squares of the residuals at the 34 body points’ should be less than 0.001. A spot check of the solution iterates indicated that the relative change in near the body was about 0.1] percent when the criterion was satisfied. For the solution of the Laplace equations, an overrelaxation factor of 1.4 was chosen for grid points inside the computational region. For grid points on the body, where a Neumann boundary condition is specified, such a large relaxation factor caused instabilities in the solution process for }. 9 The instabilities first appeared in areas where the body curvature was large. No completely satisfactory explanation of the problem was found, but numerical experimentation led to the use of underrelaxation factors between 0.75 and 1.0 on the _ body. Such underrelaxation factors on the body eliminated the instabilities. Figures 7-10 depict some results of computing the free-surface position during the first two periods of forced motion with the amplitude such that ho/b = 0.4 and the frequency such that bo 2/g = 2. JA timeysisiz:ep ys OL At = ORO? resolved one period of motion into about 314 time steps. Figures 7 and 8 show the free surface and the body with a fixed coordinate system. Figure 7, corresponding to times between 5.4 and 8.8, shows a rising body and a_ sinking free surface near the body; Figure 8, for times between 9.0 and 11.4, shows a descending body and a rising free surface. Figures 9 and 10 show similar results, but the coordinate system in these figures is fixed to the heaving body. The time-dependent details of the free surface near the body are clearer in this frame of reference. Figure 9 shows a rising body between the times 6.2 and 9.0; Figure 10, a descending body between the times 9.2 and 12.2. It is interesting to note that the slope of the free surface near the body is largest at about the time the body has attained its maximum height. Results for computing the vertical force on the heaving U-shaped body are depicted in Figure 11. From the figure it is seen that the force has _ become periodic in less than one period of forced motion. Also shown on the figure is a curve of the force versus time predicted from the second-order theory of Papanikolaou and Nowacki [18]. The other curve is a linear magnification of the results from calculating the force when the amplitude of motion is eight times smaller (ho/b = 0.05) but the frequency of the motion is the same. 20 Qualitative agreement among the curves appears to be _ good. The deviation between the computed nonlinear force and either the linear or the second-order force represents aspects of force which could not be predicted from linear or second-order theory. Another comparison of the computed force with the force predicted by second-order theory is obtained from a Fourier analysis of the curve of force versus time for the last period of forced motion. According to second-order theory, the total dimensional force on the body when its center (xg,yo) oscillates vertically according to the formula y = hg sin (ot) is given by F = 2 pg A (0.97 1/4) + 2 og b2 e F, sin (ot + 61) + 2 og b2 e2 Foy + 2 pg b2 e2 F,, sin (2ot + 65) (28) where e = ho/b is the perturbation parameter. The first term represents the hydrostatic force on the body at its neutral water position. The second term represents the first-order force on the body, and the last two terms represent second-order modifications to the first-order force on the body. A Fourier analysis of the computed nonlinear force versus time curve will produce coefficients of the Fourier expansion in the orthogonal functions {1, sin (nt), cos (nt) }. The coefficients are obtained from the usual integrals by numerical integration. In the case of the computed nonlinear force, the center of the body (xg,yg) oscillates vertically according to the formula yg = hg cos (ot) instead of y = hg sin (ot). After considering this difference, the nonlinear results can be used to arrive at the five parameters describing the force in Equation (28). The procedure, of course, assumes that the contribution to the force from third-order effects is insignificant. Table 1 shows how the computed results compare with those of Lee [17], and Za Papanikolaou and Nowacki [18]. The computed results seem to be good except for the second-order amplitude F9;. The discrepancy in the second-order force amplitude is probably due to insufficient accuracy in the computed nonlinear force. Since first-order contributions to the force are dominant, the TABLE 1 — FIRST- AND SECOND-ORDER FORCE COEFFICIENTS Papanikolaou Computed and Nonlinear Nowacki Results 0.67 relative accuracy of the remaining force when the first-order predictions are subtracted is small. The error is compounded since values calculated from the remainder are divided by the square of the perturbation parameter, a small number, to get FQ and F9,. The computed second-order phase angle is, however, excellent. 22 A rough estimate of the accuracy of the force computations can be obtained from two sets of computations in which only the time step is different. One such numerical test was conducted for this problem with the frequency parameter o2b/g = 1.7 and the amplitude of the motion such that ho/b = 0.3. The two time steps used were At = 0.02 and At = 0.01. After about t = 2, the average difference in the computed force was less than one percent. Thus the error seems to be quite small. Various frequencies of forced motion have been considered for the amplitude corresponding to ho/b = 0.3. Predictions of the coefficients for second-order theory were made and compared with those of Papanikolaou and Nowacki [18], Potash [19], and Lee [17]. The results are presented in Figures 12 and 13. The computed first-order coefficients agree well with the previously computed results. The computed second-order phase angle agrees well with the results of Papanikolaou and Nowacki. The magnitude of the sinkage (the third term in Equation (28)) agrees well with the previously computed results, but the magnitude of the oscillatory part of the second- order force does not agree so well. This is a reflection of the relative error in the nonlinear force when the first-order force, accounting for most of the force, has been subtracted. It is doubtful that there is more than about one digit of accuracy in the computed second-order results. To obtain these results, the fluid motion resulting from about two periods of forced harmonic oscillation was computed, and the force for the last period of the motion was decomposed into its Fourier components. 23 CONCLUSION A computer program to compute the fluid motion resulting from forced harmonic heaving of a U-shaped body in a free surface has been produced. The program computes forces in reasonable agreement with the results of previous researchers who used second-order perturbation theory. The forces have been computed for cases of significant nonlinear fluid motion. The computer program will be modified to handle a variety of body shapes and motions. Unlike the complex-variable techniques of Faltinsen [5] and Vinje and Brevig [6-8], the method chosen to solve the problem is extensible to the problem of 3-D nonlinear ship motions. The fluid domain has not been assumed periodic, as was assumed by Vinje and Brevig [6-8]. The location of the intersection of the free surface and the body has not been found from extrapolation, as was used by Vinje and Brevig [7]. Additional research needs to be performed before this method can be applied to the general problem of ship motions in two or three dimensions. Satisfactory results have been obtained for a heaving U-shaped cylinder, but heaving wedge-shaped bodies or blunt bodies in horizontal motion in the free surface are known to have more singular fluid flow behavior. Before such behavior can be treated properly, a more sophisticated treatment of the flow near the hull/water surface intersection will be required. Such a _ treatment should include both the dynamic free-surface equation and the kinematic condition associated with the solid boundary of the cylinder. 24 Figure 1 - U-Shaped Body in the Free Surface at Time t = 0.0 25) (xg.VQ) IN (xg ap 0.9A, Yo) XO =0 Yo = -hg cos (t) x = Xg + A(cos (a) - 0.1 cos (3a)) Y = Yo + Alsin (a) + 0.1 sin (3a)) Figure 2 - Geometric Description of the Moving Body Contour 26 Figure 3 - Fluid Region for Flow Symmetric about the y-Axis 27 Figure 4 —- Subregions into Which the Fluid Region is Divided 28 18x9 86 x9 Figure 5 - Mesh Size in Each Subregion of the Computational Region 31 obs Doo oe Figure 7 — Free-Surface Elevations Near Body for t 32 0°6 Coos Oo) C7206) = al Bios Apog oy} 02 pexXTy soUeTSFOY Jo owWeAg © UL SUOTITSOg 90eFANS-9erg - 6 OINBTY x GLL GL Gc l L SLO sO Sc0 0) Sc 0 0) (1) soo Oy + WE'n + A SLO Sc lL 33 CAG Toe GG on = et Os) APOg etd 0] pexXTy sVoUeTeTey Fo oWleAIg e UT SUOTITSOd 90BFANS-Sei1W — OL ean3ty GLL GL Gel L SL0 S'0 Sc 0 0) Sc'0- sc 0 (3) soo Oy + wero + A SL0 34 q OWT], FO SUOTJOUNYA se 9d10q ATepig—puodes UJTM pue 2d10q APOUTT YITM 9DI0g APOUTTUON peyndwog Jo uostaedmog -— [TT eansty I ug uy ug UG up ug uz u fe) LO a x, | a 7 . t i ¥SZL°O th +4 i tH. \ | Sco 4 SZe'0 02 = 6/20q © = AONANODAYS £999'0 = q = INV38-41VH (YVANITNON) SJWIL SASH3SA 3D0HO4 0 —e—4+— { (Y3GYO GNODAS) AWIL SNSHSA 3DHOSA —+—+— °u = 3qGNLINdNV £99¢°0 (YVANIT) AWIL SNSY3A ADHOS’ = —+—+—— vo 35 Zeputts9 pedeyg-p Tox Aequmy AoUSnbeirg snsteA soo10g TepzQ-—puooss pue -4SsATq - ZL eunSty 6/20q bid? ¢€ (Gl Sse orl St vy Sb eh bt tL 60 80 40 90 GO v0 €0 20 10 O. a ee | ee pe eel aie] a eee 80 GaLNdWOD soso FER) So oS SS HSVLOd --------. L IISVMON GNV NOVIONINVavd ————— 36 AJeputTA9 pedeys-n 4oJ Joequny Aduenberqg snsieA Tepig—puodes pue -}SATg™ JO sed10q Jo sTZuy sseug —- ET 2An3Ty 6/20q abe 4b Cb ih A ei ee tL 60 80 £40 90 GSO v0 €0 20 10 O “] | I | I Pls os- OSL OO0€ d3aLndWwos 431 ose HSV1Od DIOVMON GNV NOVIONINVdVd oor (Bap) &9 ‘bg 37 “Yo2 tsa -yonsupety auatery gitod sisi