Slender Ship Procedures That Include the Effects ofYaw, Vortex Shedding and Density StratificationByflaw LWongB.Sc. (Hons), Naval Architecture/Shipbuilding, University of Newcastle Upon Tyne, 1985M.A.Sc., Mechanical Engineering, University of British Columbia, 1990A THESIS SUBMITFED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDOCTOR OF PHILOSOPHYinTHE FACULTY OF GRADUATE STUDIESDEPARTMENT OF MECHANICAL ENGINEERINGWe accept this thesis as conformingto the required standardOctober 1994© Haw L Wong, 1994In presenting this thesis in partial fulfilment of the requirements for an advanceddegree at the University of British Columbia, I agree that the Library shall make itfreely available for reference and study. I further agree that permission for extensivecopying of this thesis for scholarly purposes may be granted by the head of mydepartment or by his or her representatives. It is understood that copying orpublication of this thesis for financial gain shall not be allowed without my writtenpermission.(Signature)________________________Department of QAThe University of British ColumbiaVancouver, CanadaDate 0c%r .99‘DE-6 (2/88)AbstractThe accurate determination of hydrodynamic loads on moving ships is important for hullform design and optimization and structural design purposes. This is especially true at thepreliminary design stage during which time quick predictions of the forces and momentsacting on a ship advancing steadily with, and without, yaw would be extremely useful. hiview of this, simple numerical cross-flow algorithms has been developed. The numericalprocedures are based on slender body theory, which is used to convert the threedimensional problem into a series of two dimensional wavemaker problems in the plane oftransverse sections, marching in small steps from the bow section towards the stem.Fluid density stratification, vortex shedding, finite water depth and nonlinear freesurface effects can be allowed for in the algorithms. A procedure for handling densitystratified flow has been developed and successfully used for the calculation of surface andinterfacial waves created by a prolate spheroid. Vortex shedding is modelled using thediscrete vortex method. A hybridization of the discrete vortex and boundary elementmethods is achieved and illustrated in a test case of predicting the forces acting on anoscillating flat plate.The wavemaker, with the fully nonlinear free surface conditions, is used forcalculating the generated wave pattern and wavemaking resistance of a Wigley hull. Theeffects of finite water depth on wavemaking resistance are calculated. The hybridboundary element-discrete vortex method is used for determining the hydrodynamic forcesand moments acting on a yawed Wigley hull.IITable of ContentsAbstractTable of Contents iiiList of Figures vAcknowledgement viiiINTRODUCTION 1Chapter 1 : Review of Relevant Literature 51.1 Boundary element method and free surface flows 51.2 Potential flow around ships 81.3 Interfacial waves generated by a moving body 101.4 Discrete vortex methods 101.5 Vortex-free surface-body interactions 12Chapter 2: Formulation and Calibration 142.1 The wavçmaker problem 142.2 Free surface evolution for higher order elements 182.3 Bow wave generated by a wall-sided wedge 21Chapter 3: Slender Body in Two-Layer Medium 263.1 Boundary conditions 263.2 Numerical procedure 291113.3 Numerical results.31Chapter 4: Hybrid Boundary Element/Discrete Vortex Method 394.1 Discrete vortex method 404.2 Vortex shedding in boundary element method 444.3 Forces on sharp edged cylinders 47Chapter 5 : Applications in Nonlinear Slender Body Ilydrodynaniics 525.1 Wavemakiiig resistance of a slender body 525.2 Wavemaking resistance in shallow water 655.3 Hybrid method in yawed ship calculations 67Chapter 6: Conclusions and Recommendations 816.1 Conclusions 816.1.1 Slender body formulation and implementation 816.1.2 Ships in density stratified water 826.1.3 Hybrid discrete vortex-boundary element method 826.2 Recommendations 84Nomenclature 86Bibliography 88Appendix A: Slender body theory and the wavemaker 93Appendix B: Higher order boundary element method 98ivList of FiguresFigure 2.1: Coordinate system fixed on the body 15Figure 2.2: Boundary definition for wavemaker problem 15Figure 2.3 : Infinite wedge definition sketch 22Figure 2.4: Calculated and measured bow wave amplitude for wedge, = 150 24Figure 2.5 : Calculated and measured wave peak location for wedge, I = 150 25Figure 3.1 : Prolate spheroid and dimensions 27Figure 3.2 : Definition of the computational domains 28Figure 3.3 : Overview of the numerical procedure 30Figure 3.4: Free surface contours for spheroid at F1 = 3.0 34Figure 3.5 : Interfacial displacement for prolate spheroid at F1 = 3.0 35Figure 3.6: Interface contour plot for prolate spheroid at F = 3.0 36Figure 3.7: Location of interfacial crest lines, F1 = 2.95 37Figure 3.8 : Location of interfacial crest lines, F = 3.00 37Figure 3.9 : Decay of first interfacial crest, F = 2.95 37Figure 3.10: Influence of upper layer thickness on first interfacial crest 38Figure 4.1: The physical Z-plane and the transformed .-plane 40Figure 4.2: Vortex sheet rollup at t = 1.00 for a flat plate 43Figure 4.3 : Exclusion of a vortex singularity in domain D 44Figure 4.4: CQmpptational domain and boundary conditions 47Figure 4.5 : Forces on an oscillating flat plate, K 6.6 49VFigure 4.6: Drag coefficients for the flat plate 50Figure 4.7: Inertia coefficients for the flat plate 51Figure 5.1: Wigley hull dimensions and coordinate system 53Figure 5.2 : Wave on side of Wigley hull at FR = 0.267 57Figure 5.3 : Wave on side of Wigley hull at FR = 0.316 58Figure 5.4: Contour plot of wave field for Wigley hull at FR = 0.267 59Figure 5.5 : Measured and calculated wave field (YNU) for Wigley Hull, a = 00 60Figure 5.6 : Free surface displacement for Wigley hull at FR = 0.267 61Figure 5.7 : Pressure distribution on Wigley hull at FR = 0.267 61Figure 5.8 : Wavemaking resistance against YNU results 62Figure 5.9 : Wavemaking resistance against Aanesland (1989) 63Figure 5.10 : Computed and measured sinkage values for Wigley hull 64Figure 5.11 : Wavemaking resistance for Wigley hull in shallow water at FR = 0.3 16 66Figure 5.12: Squat and trim for Wigley hull in shallow water at FR = 0.3 16 66Figure 5.13 : Vortex induced normal velocity on a boundary element 70Figure 5.14 : Wave profile on sides of Wigley hull at FR = 0.267, a = 100 72Figure 5.15 : Wave field for Wigley hull at FR = 0.267, a = 100 73Figure 5.16: Measured wave field (YNU) for Wigley Hull, a = 100 74Figure 5.17: Free surface displacement for Wigley hull at FR = 0.267, a = 10° 75Figure 5.18: Pressure distribution for Wigley hull at FR = 0.267, a = 100 76Figure 5.19 : Forces and moments for yawed Wigley hull 77Figure 5.20: Development of total vortex strength over the hull length 78viFigure 5.21 : Keel vortex evolution for Wigley hull, a = 100 (forward) 79Figure 5.22: Keel vortex evolution for Wigley hull, a = 10° (aft) 80Figure A.1 : Frame of reference and notation 91Figure B.1 : Intersection point and normal velocities 99vuAcknowledgementThis thesis is dedicated to my family, from every member of which I have drawn muchsupport and encouragement. In particular, my wife Julie deserves much praise for herassistance and patience. My children Justin and Sarah are a great source of joy andmotivation. The work reported in this thesis was done under the supervision of ProfessorSander Causal, who has over the years offered much help and guidance. I am also gratefulto Professor Kazuo Suzuki (Yokohama National University, Japan) for providingexperimental data and material on slender body theory. Professor Orner Goren’s (IstanbulTechnical University, Turkey) kind generosity in providing the Dawson programs is muchappreciated. I should also mention the great fellowship I enjoyed with all my past andpresent colleagues. To all these people and others who have in one way or another mademy work easier, I offer my humble thanks.The generous scholarship in computer time on the Fujitsu VPX24O/1Osupercomputer granted to me by the High Performance Computing Centre of Calgary issincerely appreciated.Financial support for this work was provided by the Natural Sciences andEngineering Research Council of Canada.ViiiINTRODUCTIONRobust and reliable methods for determining the hydrodynamic characteristics of a vesselor offshore structure are valuable tools for design engineers and naval architects. Theability to fine tune designs at an early stage means an improvement to the design processand consequently better designs. Apart from that, much cost and time savings accruefrom the reduction of required model testing and the avoidance of remedial work aftercommissioning.Ongoing research over the years, aimed at the development of better and moreefficient numerical tools of practical interest to ship designers, has resulted in numerousinnovative approaches. Despite this, numerical solutions relating to nonlinear ship wavegeneration have proved elusive. A number of computational methods, mostly linearized,were introduced in the 70s. However, as the importance of nonlinearity and viscouseffects become better understood, much work in the last decade or so has attempted toaccount for both.The objective of the present work is to use potential flow theory to study thehydrodynamics of a slender body in calm water, including the effects of nonlinearity, yaw,vortex shedding and fluid density stratification. Numerical implementation was achievedusing the higher order boundary element and discrete vortex methods. Slender bodyapproximations simplified the threc dimensional moving ship problem into a sequence oftwo dimensional wavemaker problems. Working in the cross-flow plane offers manyadvantages from a computational standpoint by allowing the inclusion of vortex shedding,consideration of finite water depth, uneven seabed conditions and density stratification.The schemes developed and reported in this thesis can conveniently be divided into threecategories, which will be described in more detail in the following sections.Nonlinear Slender Body ProcedureThe slender body formulation has been generalized to accomodute oblique motion of aship. This was implemented through a wavemaker approach in which the 3D ship waveproblem is converted into a series of 2D ‘wavemaker’ problems in the plane of stations,marching in small steps starting from the bow section. This approach was used in thestudy of bow wave making for an advancing wall-sided wedge.The algorithm was then extended to include the fully nonlinear free surfaceconditions and applied to the calculation of deep and shallow water wavemakingresistance and the waves generated by an advancing Wigley hull. Vortex shedding effectswere not included in these calculations. Vorticity generated by the round bilges wereassumed unimportant in forward motion problems and hence not accounted for. Thisassumption can be justified considering the success of pure potential flow ship resistancecalculation methods. The nonlinear free surface evolution was achieved using theEulerian-Lagrangian approach. The procedure, when extended to general hull forms canbe an invaluable tool in hull form optimization applications. Suitable ship types for suchapplications are naval vessels, container ships and multi-hull vessels.Ships in Stratified WatersDensity stratification is a common phenomenon occuring both in the oceans and inlandwaters. The variation of the density with water depth is due to changes in the salinity2and/or temperature over a region called a pycnocline. Although the density differencesmay be small, ships moving slowly in the vicinity of such interfaces experience what cameto be known as deadwater effects. When the ship speed approaches that of the fastestgenerated internal waves, the increase in drag over the hull is very large and the availablepower from its engines may not be enough to push the ship beyond that critical speed.This is due to energy being spent in the creation of internal waves below the free surfaceand the disturbance of the flow around the ship due to these internal waves. It is alsoknown that a vessel loses controllability in such situations. In naval applications, fully orpartially submerged vessels could produce undesirable internal wave signatures.The procedure given in this thesis is applicable to a slender hull moving in a twolayer fluid (although it is certainly possible to extend the application to multi-layer flows).The ship draught is taken to be of the order of the, thickness of the upper layer. Two fluiddomains are used, the lower domain with a density slightly higher than the upper foroverall gravitational stability. The classical treatment of two dimensional internal wavesprovides a condition relating to the interface between the two layers. The boundaryelement method was used to obtain the transient solution by solving the lower, then theupper domain sequentially at each time step. This algorithm is expected to be useful forresistance prediction of submerged bodies and for the study of internal wave behaviourwithout the presence of a body.Vortex Shedding EffectsThere are numerous instances in the study of fluid dynamics where vortex shedding forcesare significant and should be taken into account. A good example is in the tumerical3prediction of the roll response of a ship with sharp corners and bilge keels where wavedamping is light and vortex effects are important. However, vortex effects are usually notincluded in potential flow applications to ship dynamics. This is partly due to the complexnature of the fluid-structure interaction which generates vortices from a hull. There is alsothe difficulty of implementing vorticity calculations in wave problems without resorting toviscous flow computations involving the Navier-Stokes equations.A hybridization of the boundary element and discrete vortex methods has beenachieved. This hybrid method allow a simple but effective way to account for the effectsof vortex shedding in potential flow calculations. The vortex generation mechanism istaken to be the keel. Bow and secondary vortices were not included since the complexnature of the flow in both cases cannot be properly modelled using a simple discretevortex approach. The computational advantage in being able to account for vortex effectswhile still working with boundary values is significant.The hybrid method was implemented for the case of a Wigley hull in obliquemotion. Changes in wavemaking resistance and yaw forces and moments due to smallangles of attack, whether intentional or otherwise, can be studied. Such studies are ofpractical interest from an operational standpoint. In addition, forces and moments inyawed ship motion can be used to determine the so called ‘stability derivatives’, importantin the study of ship manouevring. This in turn allows systematic investigations into thedirectional stability of a vessel at the design stage.4Chapter OneReview of Relevant Literature1.1 Boundary element method and free surface flowsThe boundary element method is a much used method for solutions to fluid dynamicsproblems. An advantage of this method is that for potential flow calculations onlyconditions on the boundary need to be specified and calculations are carried out along thediscretized boundary only. This method is well suited for application in boundary valuepotential flow simulations involving relatively complex and possibly time varyinggeometries.Much attention has been focused on finding stable numerical solutions to thenonlinear problem. The satisfaction of the nonlinear free surface conditions at a locationwhich is unknown a priori has been a major stumbling block to the achievement of stablenumerical solutions in nonlinear wave making. Longuet-Higgins and Cokelet (1976) useda mixed Eulerian-Langrangian form of the free surface conditions to trace the temporaldevelopment of high amplitude surface waves. They assumed spatial periodicity, mappingthe domain such that only free surface variables need to be handled in the solutionprocedure. Their formulation was successful in coping with the multi-valuedcharacteristics of an overturning wave.5Faltinsen (1978) applied the boundary element method to the two dimensionalproblem of liquid sloshing in an oscillating container in which the nonlinear free surfaceboundary conditions are treated for consideration of a nodal point constrained to move inthe vertical direction only. Nakayama and Washizu (1981) also reported work on theabove problem. Although not clearly stated, it appears that this work used the kinematicfree surface boundary condition in its linearized form. Cointe (1989), using the mixedEulerian-Lagrangian method, studied wave diffraction from a submerged cylinder.Pawlowski, Baddour and Hookey (1991) discussed the use of linear elements for nonlinearwave simulation. Isaacson and Ng (1993) reported a time domain boundary elementscheme for nonlinear wave radiation from floating cylinders in sinusoidal sway, heave androll motion.In the numerical simulation of nonlinear interactions involving a free surface and afloating body, boundary conditions are to be satisfied at the instantaneous location of theboundaries at any given time step. This requirement, when applied to the fluid-bodyintersection point, results in difficulties arising from the singular behaviour of the fluidvelocities at that point. Consequently, the definition of conditions such as the velocitypotential, the direction and magnitude of the velocities, and the location of the free surfacenode at the fluid-body intersection point has proved to be a challenging task for allresearch into nonlinear wave simulation. Yet, the resolution of this problem is crucial forthe success of any such simulation, both in terms numerical stability and accuracy.Vmje and Brevig (1981) used the body boundary condition at the intersectionpoint. The velocity potential is taken to be the unknown there and emerges as part of the6solution. Un (1984) specified both the velocity potential and stream function value at theintersection point in an algorithm using the Cauchy integral formulation for boundaryelements. Dommermuth and Yue (1987) used a similar technique and a boundaryregridding scheme. The regridding they used replaced the Lagrangian points to give equallength elements on the free surface. The use of double nodes has recently emerged as apopular idea in the fluid-body intersection problem. Grihi and Svendsend (1990), amongothers, used this idea, enforcing compatibility only at the end of each time step.Another major difficulty in numerical calculations of nonlinear wave making is theapplication of the far field condition stemming from the truncation of infinite domains tomanageable dimensions. The lack of proper physical conditions governing the passage ofenergy out of the computational domain has made the radiation condition for full nonlinearproblems difficult to define. The eventual reflection of waves back into the computationaldomain causes errors in the solution for the region of interest. This undesirable situation iscommon to all numerical methods. To cope with this, a variety of ways has been devisedwhich define conditions at a given cut-off boundary which is not excessively distant fromthe source of the disturbance. Obviously, limiting the size of the computational domainimproves computational efficiency.Chan (1975) used a mathematical wave damping device to remove energy of theoutgoing waves. Orlanski (1976), on the other hand, used the Sommerfeld radiationcondition, proposing a finite difference scheme to calculate a depth and time dependentcelerity, which should always be directed out of the domain. The method does not rely onknowledge of the dispersion characteristics of the wave train and as a result may be used7for nonlinear waves. Isaacson and Cheung (1991) obtained satisfactory results for secondorder wave propagation.Yen, Lee and Akai (1977) applied the undisturbed condition at the cut-offboundary which is periodically expanded downstream as the propagating wave approachesit. Apart from the fact that computational efficiency deteriorates with time as the domaingrows, the solutions reportedly became unstable after some time. In a later paper, Yenand Hall (1981) used Orlanski’s method for nonlinear waves together with a proposal toremove the higher frequency error due to the reflected waves having wavelengths of abouttwo grid intervals. This was done using two separate methods: filtering over the wholefree surface using a mathematical smoothing function and a modification of theformulation to include dissipation and dispersion. Dommermuth and Yue (1987) matchedthe inner nonlinear region by a linear outer region and reported this gave a robustalgorithm.More recently, Tanaka and Nakamura (1992) proposed an extension to theOrlanski scheme for application to unsteady and irregular nonlinear waves. Their methodused the nodal point nearest to the free surface-radiation boundary intersection for theprediction of the time dependent phase velocity of the outgoing waves.1.2 Potential flow around shipsDue to the difficulty of obtaining complete analytical solutions to the ship wavemakingproblem, a tremendous amount of work has been done on the development of numericalmethods. Hess and Smith (1967) developed a source panel method for the calculation ofpotential flows around bodies of arbitrary shape. Dawson (1977) devised a method in8which a combination of the double body Rankine source panel method and marchingscheme in the stream-wise direction produced useful information relating to. thehydrodynamics of steady ship motion. Cheng (1989) extended the method for 3D transomstem flows while Tahara (1992) included lifting effects in the double body potential.Gadd’s (1976) iterative numerical procedure for ship resistance prediction was followed upby Causal et al. (1991) for nonlinear wave calculations. The so-called Neumann-Kelvinapproximation is also popular and well reported in the literature, see t’or example the workof Suzuki (1985).Maruo (1982) presented a formulation based on the Neumann-Kelvinapproximation. Slender body approximations result in a simplification with thedistribution of Kelvin sources on the hull surface only. A marching procedure starting atthe bow can be used, with the source strengths at a given section being determined by theship normal velocity at that section and a known quantity involving data evaluatedupstream. Upon completion of the procedure, all the source strengths on the hull surfaceare determined. Quantities of interest are then computed by evaluating the effects all thesesources. Song, Ikehata and Suzuki (1988) applied the method to the Wigley hull,reporting good agreement with measured values. Calisal and Chan (1989) reported on theuse of a simple wavemaker method and slender body theory for the calculation of bowwaves produced by an advancing wall-sided wedge. Wong (1993) pointed out that usingstarting conditions at the bow improves the results from slender body calculations. Grilli(1991) used a similar approach for a cusp-ended ship-like body. More recently, Song andMaruo (1993) described a wavemaker approach to bow impact and deck wetness studies.9Kim and Choi (1993) calculated shallow water effects on hydrodynarnic forces acting onvarious hull forms using a three dimensional iterative procedure involving the Rankinesource panel method.1.3 Interfacial waves generated by a moving bodyOne of the earliest works on ship excited interfacial waves in stratified media was reportedby Ekman (1904), who was able to experimentally confirm that internal waves created bya moving ship cause a large drag on the hull. Sabuncu (1962) calculated the Michell waveresistance in the presence of interfacial waves in a two layer fluid. More recent interest inthe problem was fuelled by attempts to explain the long narrow V-wake left behind a shipmoving in stratified seas as sensed by radar. Tulin and Miloh (1990) used linearizedtheory to obtain a solution in terms of a Green’s function. Wong and Causal (1993b) usedthe cross flow approach to predict the interfacial waves excited by a prolate spheroid andcompared the results to the experimentally measured values of Ma and Tulin (1992).Watson et a!. (1992) reported on full scale experimental work with ships traversinga highly stratified Scottish sea loch. They concluded that interfacial wave amplitudes areonly weakly dependent on ship speed but are more sensitive to changes in fluidstratification.1.4 Discrete vortex methodsFlow separation occurs when the boundary layer on a body surface reaches a sharp edgewhere the radius of curvature of the edge is very much smaller than the boundary layerthickness. In two dimensional flow, the adverse pressure gradient set up results in the10formation of an unstable shear layer which subsequently rolls up into a tight spiral and isthen shed into the flow field as a free vortex sheet. Such a vortex sheet can beapproximated by the introduction of discrete vortices at given time intervals as the flowdevelops.A persistent and yet unresolved problem in such schemes is the irregular roll-up ofthe vortex sheet as well as numerical instability due to the uncharacteristically highvelocities induced when a vortex is very near to another vortex or its image. Theseproblems have been the focus of study in research on vortex methods since the early 70s.The interaction between a pair of vortices in close proximity results in highinduced velocities on each. Chorin (1973) used vortices with a viscous core to modelviscous diffusion as well as to stabilize his numerical procedure. As a result, the maximuminduced velocity of a vortex is finite at a given distance from the vortex position anddecreases to zero as the centre of the vortex is approached. Clements and Maul! (1975),to obtain stable solutions, limited the induced velocities by amalgamating any pair ofvortices that are too close together. Fink and Soh (1974) pioneered a scheme ofrediscretization in which the vortex sheet is rearranged into equidistant positions aftereach time step in the numerical procedure. This same regridding scheme was adopted byDommermuth and Yue (1987) for the treatment of the nonlinear free surface to improvenumerical stability. The above, and many other similar, schemes help in one way oranother to give smooth vortex roll-up as well as extend the computational time span inwhich calculations remain stable, Wong (1990).HGraham (1980) applied the discrete vortex method to calculate the vortex forcesinduced at a sharp edge in oscillatory motion at low Keulegan-Carpenter numbers. Thiswas done by regarding vortex shedding from an infinite wedge as the inner region of flowpast a large but finite body. The underlying assumption in this case is that the body lengthscale is large so that the vortex shed does not affect other parts of the body which are faraway from the edge when compared to the flow length scale. Downie, Bearman andGraham (1988) followed up on this and calculated the vortex damping forces on a rollingrectangular barge. Wong and Causal (1993a) applied a hybrid discrete vortex-boundaryelement method to the determination of the forces acting on a submerged, oscillating flatplate. The results from that study compared well to measured values.Comprehensive general reviews on various vortex methods can be found inLeonard (1980), Lugt (1981) and Sarpkaya (1989).1.5 Vortex-free surface-body interactionsThe study of interactions between vortices and free surfaces, density interfaces andmoving bodies has been motivated by attempts to improve potential flow results forsituations where vortex effects were thought to be important. In particular, bodygenerated vorticity can result in significant influences on the hydrodytiamic forces actingon the body.Bradbury (1986) generated a number of rectangular transverse sections to describea yawed ship-like slender body for the study of hydrodynarnic effect of bilge vortices. Itwas necessary, in that case, to use rectangular blocks where the submerged sharp corners12representing the bilges provide well-defined separation points for the vortex sheddingalgorithm employed. Time marching was used to advance from one cross-flow plane tothe next. At each section, conformal mapping was used as part of an iterative procedurefor the solution. The free surface was assumed rigid for his work.Braathen and Faltinsen (1988) applied a time stepping procedure for twodimensional vortex shedding from rectangular sections in forced harmonic roll oscillation.Tryggvason et al. (1991) studied the collision of a vortex pair and a vortex ring with asharp density interface and free surfaces.13Chapter TwoFormulation and Calibration2.1 The wavemaker problemThis section summarizes the development of a ‘wavemaker’ algorithm for the computationof flow around steadily advancing slender ships. Simplification of the three dimensionalproblem is achieved using the slender body approach given in Appendix A. This reducesthe problem into a series of two dimensional ones in the cross-flow plane, solvedsequentially by marching in time from the bow in the direction of the stern. An effectiveway to visualize this is to view the ship from a fixed yz planc ahead of the bow. As theship advances past this viewplane, different transverse sections are ‘cut’ by this plane suchthat the fixed observer sees the body sections ‘opening up’ until around the midship regionand ‘closing in’ again towards the stern.The frame of reference has the origin on the ship centreline at amidships and on theundisturbed free surface as in Figure 2.1. The angle of incidence of the free streamvelocity, or the yaw angle, is taken to be a. Details concerning the use of slender bodytheory to derive a set of equations describing the wavemaker problem is given in AppendixA and will not be repeated here.14Figure 2.1: Coordinate system fixed on the body.Ideal fluid assumptions, together with slender body approximations, result in thetwo dimensional governing I.aplace equation from Equation (A.15):Pyy + +72 0Figure 2.2: Boundary definition for wavemaker problem.(2.1)APNotationLength LBeamB zDraughtdZFPFree streamvelocity (U)Note: Origin is at undisturbedfree surface level.zFree suace (Sf) Free suace (S///Wall Wall(Sw) /U(YCOSU ± sina)/(1 -t-Y)5/ (Sw)/ // IGoveming Equation : =/// / / // / / /Bottom boundary (S0)15All normal velocities were taken to be positive pointing out of the domain (seeFigure 2.2). The bottom boundary S0 is taken to be impermeable. A ‘mirror image’technique using symmetry to satisfy the bottom condition is often used to reduce thenumber of equations to be considered in the boundary value problem, see for example,Isaacson and Cheung (1991). This can easily be done if collocation points are away fromthe intersection points between the bottom boundary (S0) and the sides (Sw). However,the present work used linear and quadratic elements, that is collocation points are found atall boundary intersections. The additional computational book-keeping resulting from theuse of this technique erodes the savings in computational effort. Moreover, by keepingthe bottom boundary, the formulation can be extended to the more general case of unevenand sloping seabeds. Thus, the bottom boundary was retained for all the calculationsdescribed in this thesis.The impermeable wall S is placed far enough (through numerical experiments)from the ship to ensure that wave reflection does not affect the flow calculations. Therewere two reasons for not imposing the radiation condition at this boundary. Firstlycomputed results are to be compared to values measured in tank testing where flow isalways bounded. The other reason is that it is extremely difficult to implement thenumerical radiation condition for the mixed Lagrangian-Eulerian approach adopted innonlinear cases in the present work.On the body, the normal velocity in the cross-flow plane is due to the rate ofchange of the sections as the marching proceeds. The ‘wavemaker’ displacement is thus16represented by the changing hull sectional shapes along the ship length. Normal velocitieson the hull section (negative into the fluid) at each step is defined using Equation (A.14):a4 -U(Ycosct±sincL)an Ji+YAt the free surface, S, the dynamic and kinematic conditions, to be satisfied at theinstantaneous free surface location, are applicable. Using the conversiona I a .—= UI cosa—+sinct— (2.3)ôt‘ox OyJassuming that U is uniform everywhere in the flow field. From Equations (A.16) and(A.17):atz=rl (2.4)at OyOy z+ +()] + =0 at z = (2.5)The marching procedure from one cross-flow plane to the next can be viewed as atime domain problem where the time step is the amount of time taken to traverse dx, thespatial step along the longitudinal axis.The hydrodynamic forces and moments are evaluated according to the followingintegrals:L12= fcixf p.ndS (2.6)—L/2L/2M= f dxf p(rxn)dS (2.7)-L1217where n is the normal vector at the middle of the element dS and r is taken with referenceto the origin (0,0,0). The integration is carried out over the girth of each consideredcross-section (SB) and, after completion of the sweep, over the ship length usingSimpson’s rule. The pressure p is obtained from the Bernoulli equation:p = P—P0 =— U(cosa+sinu)+_((.)()2÷gz (2.8)ox 2Oy )2.2 Free surface evolution for higher order elementsFor the solution of the cross-flow plane boundary value problem, velocitypotentials are specified at all free surface nodes. The free surface conditions are satisfiedat its exact location for nonlinear calculations. Where linearized conditions were used, thefree surface boundary was fixed at the mean level.In the nonlinear procedure, adopting the Lagrangian-Eulerian approach to tracefree surface particles, the two dimensional dynamic free surface condition in Equation(2.5) can be rewritten as:D1[(O4)(O)2]atz=11 (2.9)using the substantial derivative:(2.10)Dt at OyOy OzOz18The velocities on the free surface can be described by:Di a(2.11)Dt ôzDy a4(2.12)Dt ayFor the linear solution procedure, the kinematic and dynamic free surface conditions arewritten as:= -gq at z = 0 (2.13)atz=O (2.14)öt ôz onThe solution for the time marching procedures require the use of an extrapolatorto find the potential values on the free surface and its location at the next time step. In thepresent work, the first order Adams-Bashforth predictor for time step size ôt was used.(t + ôt) = (t) + [3(t) — 1(t — t)]- (2.15)+ ôt) = q(t) + [3 (t)—‘q (t — ôt)}- (2.16)In this way, values of the velocity potential (4)) and the location of the free surface (‘q) attime (t + ôt) are predicted using Equations (2.15) and (2.16). With Neumann conditionsspecified on all the other boundaries, the problem is then solved using the higher orderboundary method (Appendix B). The normal velocities on the free surface nodes emergeas part of the solution from the boundary element method. From these, time rates ofchange for the potential and the free surface location can be determined using Equations19(2.9), (2.11) and (2.12). For linearized conditions, Equations (2.13) and (2.14) are usedinstead. The values of and ij can then be corrected using the Adams-Moulton corrector,Equations (2.17) and (2.18), and the problem solved a second time before proceeding tothe next time step.p(t + ôt) = p(t) + [p, (t) + p((t + ôt)}. (2.17)rl(t + ôt) = t(t)+ [r(t)+ q(t + ôt)] (2.18)The velocities in the y- and z-directions(p and p) can be determined using thenormal and tangential velocities at each collocation node.=cosO—sinO (2.19)ay as- = -sinO +cosO (2.20)az as anwhere 8 is the angle of the free surface element to the horizontal (anti-clockwise positive).The tangential derivative / aS can be found by taking finite differences, varying linearlyover each quadratic element and being constant over a linear element. Thus, for quadraticelements, the tangential velocity is given by:3P —4Pk÷1 + Pk+2 (2.21)as ASwhere AS is the element length, and k is a nodal index increasing in the anti-clockwisedirection.The singular behaviour at the intersection point S fl SB requires special attention.In the present work, continuity of potential is assumed while it is recognized that the20normal derivatives immediately to either side of the fluid-body intersection node aredifferent. The normal velocity on the body side is imposed for physical compatibility.This means that the normal velocity on the body side should be identical to that of thebody at the intersection point. The normal velocity on the free surface side is linearlyextrapolated from the nearest two free surface nodes. This leaves the potential as the onlyunknown at the intersection point. Details on how this is numerically implemented isgiven in Appendix B.2.3 Bow wave generated by a wall-sided wedgeThe wavemaker algorithm described earlier was used for the prediction of bow wavesgenerated by an infinite (after body) wedge. This simple application provided anopportunity for the calibration of the time stepping procedure. As the midship section isundefined, the frame of reference has to be shifted to the forward perpendicular, seeFigure 2.3. The computational boundary, as shown in Figure 2.3, is similar to that givenin Figure 2.2 except for the presence of the bottom of the wedge section, on which thenormal velocity is zero (ignoring dynamic trim and sinkage effects) and the use of asymmetry x-z plane to save computational effort. It is easily seen that as marchingproceeds, the body is seen as a steadily growing rectangle in the cross-flow plane.For the no yaw condition (a = 0), and taking the internal angle of the wedge to be, the normal velocity on the vertical wall of the wedge section can be expressed, usingEquation (2.2) and noting that Y,, = 0, as:= —uy —u tan(.’i (2.22)on \2J21For this application, the linearized free surface conditions given in Equations(2.13) and (2.14) were satisfied at the undisturbed level. The step size in all cases wasFigure 2.3 : Infinite wedge definition sketchfixed at an x increment of 2.5 mm. Draughts of between 4-16 inches (10.2-40.6 cm) wereconsidered. The water depth and tank width were both taken to be 2 metres. Linearelements were used throughout. Normal velocities for all boundaries except for thevertical boundary of the wedge and the free surface were specified as zero. The Dirichietboundary condition was specified on the free surface through the procedure given in thelast section.The amplitude and location of the first bow wave peak was recorded for allcomputations. These results were compared to measured values for 3 = 150 (Ogilvie,1972) and presented in Figures 2.4 and 2.5. In Figure 2.4, the bow wave amplitudeis in good agreement with experimental values. The analytical results obtained by OgilvieIsectional View]Wedge SectionzFree Surface______lidIPerspective View! —Uidy22(1972) for the bow wave amplitude, non-dirnensionalized using -- /U, showed aconstant value of 1.6, independent of the depth Froude number as well as the draught. Onthe other hand, the same non-dimensionalization of all data points in Figure 2.4 gave arange of 1.9-2.4 for the predicted results while the experimental data has a range 0.5-2.7.The computed bow wave peak location (Figure 2.5) were consistently delayed inall cases. This is probably due to the fact that the marching procedure was started with anundisturbed free surface. Generally, initial conditions in the form of potentials and waveelevations at the first cross-flow plane, that is at the bow section of the ship, can becalculated from Dawson’s (1977) method. This method basically requires the distributionof Rankine sources on the body and free surface, solving for the source strengths byimposing body, free surface and radiation conditions. Evaluation of velocities and freesurface elevations rely on these source strengths. Hence, a closed body is necessary andthe Dawson method was not applicable for the above wedge which is infinite in thelongitudinal direction.232015a)-o10.a)00S..S.• Present method—-——Ogilvie (Expt)• Present method-—-------Ogilvie (Expt)Figure 2.4: Calculated andmeasured bow wave amplitude for wedge, 150.Ioraught = 10.16 cml• Present methodOgilvie (Expt)2015a)010Ea)a-0Draught = 30.48 cml0 100 200 300 400 500Speed of Advance (cm/s)2015a)-o10a)00(2015a)I0Ea)aaIDraught = 20.32 cml• Present method————Ogüvie (Expt)100 200 300 400 500Speed of Advance (cm/s)iDraught = 40.64 cml4:4:...•..•.::...100 200 300 400 500Speed of Advance (cm/s)Th 100 200 300 400 500Speed of Advance (cm/s)24Figure 2.5 Calculated and measured wave peak location for wedge, = 15°.IDraupht =10.16 cml• Present method——Ogilvie (Expt)••Ioraught = 30.48 cml• Present methodOgilvie (Expt).125E010001:0E 750U)0a. 25U)a. 00 100 200 300 400 500Speed of Advance (cm/s)••—. 125E01000 •E 75.oC ..2 50 •/• . /a.25. ••.CD •00 100 200 300 400 500Speed of Advance (cmls)lDrauht = 20.32 cml• Present methodOgilvie (Expt)125E01000E 75C0U)0a. 25U)U)IDraught = 40.64 cml• Present methodOgilvie (Expt)• ..• .• ...• .01250100o•E••..250 ••a. ••U)U) •a.... ..I. .. I.... I....0 100 200 300 400 500Speed of Advance (cm/s)100 200 300 400 500Speed of Advance (cm/s)25Chapter ThreeSlender Body in Two-Layer MediumA ship travelling at very low Froude numbers in stratified water is known to experiencewhat is generally known as dead-water effects. As the ship speed approaches that of thefastest internal waves, the increase in drag over the hull is very large and the availablepower from its engines may not be enough to push the ship beyond that critical speed.The vessel could also lose manouevrability. This is due to energy being expended in thecreation of internal waves below the free surface and the disturbance of the flow aroundthe ship due to these waves.This chapter describes a numerical solution for the time domain calculation ofsurface and interfacial waves generated by a body moving slowly in a gravitationally stabledual density fluid. The wavemaker approach described in Chapter Two and Appendix Awere used. The fluid is assumed to consist of two distinct layers of slightly differentdensities with a discontinuity of the Brunt-Väisälä number at their interface.3.1 Boundary conditionsA moving frame of reference with the origin on the ship centreline at amidships and on theundisturbed free surface as given in Figure 2.1 was adopted. A spheroid with the samedimensions as that in Ma and Tulin (1992) was used for comparison purposes, see Figure3.1. The normal velocity on the body boundary is given by:26=_____onThe prolate spheroid geometry is defined by:y = Y(x,z)= ±[(i_4)2_z2]where B and L are the beam and length of the spheroid respectively.Figure 3.1 : Prolate spheroid and dimensions.(3.1)(3.2)The computational domains, one for the upper layer with density Pi and one forthe lower layer with density P2, have a common boundary at the interface, see Figure 3.2.Such a situation is, of course, a simplification justifiable only if the depth of the pynoclineis small. Consequently, the Brunt-Väisälä number N2 = —gOp / c3z behaves like a ôfunction, being infinite at the interface.Profile FrontLength L — 43.70cm Beam B = 890cm27BodyFigure 3.2: Definition of the computational domains.Boundary conditions have to be defined for the other boundaries in both domains,see Figure 3.2. The bottom is taken to be horizontal, flat and impermeable.‘1’= 0 at z = -h1 (3.3).Normal velocities at ABC and FGH (Figure 3.2) are taken to be zero, simulating a solidwall. Since very low ship speeds are considered, the surface wave generated by the hull isexpected to be small. Hence, the free surface conditions iii Equations (3.4) and (3.5) wereused:÷gi° =0 at z=0 (3.4)(1) (1)I = atz=O (3.5)özwhere i- is the surface wave elevation and g = 9.81 mIs2 is the gravitational constant.The classical treatment of two dimensional internal waves provide a conditionrelating the interface between the two layers. While normal velocities and pressures areF0CBHBottom boundaryA28taken to be unique at the interface, the velocity potentials go through a discontinuityacross that interface. Taking the normal velocity as positive pointing out of a domain, thefollowing two conditions can be written for an upper layer with thickness h2:(2) (1) (2)= atz=—h2 (3.6)at an anap1= p(2) + g2)(p(2) — p(l)) at z = —h2 (3.7)where 2) is the internal wave elevation with respect to the undisturbed level at theinterface. Nodes at the interface belonging to the upper and lower layers coincide andmay be viewed as ‘double nodes’, sharing the same physical location. The normalvelocities at these nodes are assumed identical as given in Equation (3.6). The negativesign for the upper layer is simply due to the convention adopted, that is, normals arepositive pointing out of the domain. Equation (3.7) results from equating pressure at the‘double nodes’. The thickness of the top layer (h2) is taken to be of the order of the shipdraught.3.2 Numerical procedureAn overview of the numerical procedure used for this chapter is given in Figure3.3. The lower and upper domains are solved sequentially in that order. The simulationstarts by considering the lower domain at time ( t + ôt ), which consists of three sideswhere the normal velocity is zero and one side (the interface) where the potential needs tobe predicted. A symmetry plane at y=O was imposed for computational efficiency. Thefirst order Adams-Bashforth predictor given by Equation (2.15) is used to estimate the29potential on the interface. Emerging from the solution are normal velocities at theinterface. On the boundary between B and G (Figure 3.2) these normal velocities are usedas input for the upper domain (see Figure 3.3).Figure 3.3 : Overview of the numerical procedureFor the upper domain, normal velocities are again specified at all nodes except atthe free surface, at which the ABM predictor given in Equation (2.15) is again used. Aftersolving the upper domain, required time and space derivatives at ( t + ôt ) are to becalculated. On the free surface, the time rate of change can be found using Equation30(3.5). Wave elevations at all nodes for C-D and E-F are calculated using Equation (2.18).Then Equation (3.4) can then be used to find pl) on the free surface.At the interface, Equation (3.6) gives the rate of change of the interfacial waveelevation. The internal wave elevations are found using Equation (2.18). For theboundary between B and G belonging to the upper domain, the time derivative of thepotential is found using a three point difference formula in time, as follows:= 3p(1)(t + ót)— 4p’)(t) + p°)(t — ôt)(3 8)2ôtWith pl) at the interface determined, Equation (3.7) is used to determine valuesof p2) This completes the predictor part of the computations. The Adams-Moultoncorrector is then used on potential values on the interface and free surface using Equation(2.17).Using the same equations as before, time and space derivatives are found at theend of the corrector part. The above predictor-corrector pair is repeated for allsubsequent time steps.3.3 Numerical resultsA prolate spheroid of L = 43.70 cm and B = 8.90 cm was used in the calculations. Thedensimetric Froude number, a depth Froude number modified by the density differencebetween the two layers, is defined as:c31F1= )h2(3.9)where U is the advance velocity, op is the difference in density between the two layers andh2 is the thickness of the upper layer. Cases where the value of F1 > 1 are considered‘supersonic’, Tulin and Miloh (1990), a term imported from aerodynamics. Note thatF1 = 1, termed the ‘critical’ densimetric Froude number, is the point where the slope of theship drag curve is very high.The densimetric Froude number is related to the ship Froudenumber (FR ) by the following relationship:FR= F0Ph2 (3.10)Two cases were computed:(1) F1 = 3.00, h2 = 6.70 cm, U = 12.86 cmJs, Op/p = 0.0028.(2) F, = 2.95, h2 = 4.25 cm, U = 11.43 cm/s, Op/p = 0.0036.These cases correspond to two of the experimental cases given in Ma and Tulin (1992),which will be used for all comparisons with results obtained from the present method.Note that Case 1 represents the situation where the draught of the body (d = 4.45 cm) isless than the upper layer thickness while in Case 2, the body extends into the lower layer.For the above cases, the corresponding ship Froude numbers are very small (F =0.055 and 0.062), justifying the use of the linearized free surface condition.A plot of the free surface displacement for Case 1 is given in Figure 3.4. The waveelevation is small but the internal wave signature is noticeable. The interfacial32displacement and wave field contours are plotted in Figures 3.5 and 3.6. The distinctivefeatures of the wave field are very similar to those measured and calculated by otherresearchers. A very deep depression occurs at the interface just below the bow, followedby an upsurge which gave rise to at least five interfacial crests in the far wake.Crest line locations, non-dimensionalized using the upper layer depth, for the twocases are plotted in Figures 3.7 and 3.8. Generally, good agreement with measured crestpatterns were achieved. However, the present method predicts earlier formation of theinterfacial crests. In addition, decay of the first interfacial crest for Case 2 (Case 1 datawas not available) are in excellent agreement with experimental values, see Figure 3.9.The decay of the first interfacial crest can be described using:(3.11)k\h2) h2)The value of the decay parameter a obtained using the present method was 0.96 while Maand Tulin (1992) gave the value of a to be approximately unity.The influence of the upper layer depth on the behaviour of the first interfacial crestheights can be clearly seen in Figure 3.10 to be well predicted by the present method.33x1jC CD I.to C;’+ 1\)C;’to 0CC)CIIIIIIVODOJcnaOII 0C)EC.)>..Figure 3.5 : Interfacial displacement for prolate spheroid at F = 3.0.EC)NEC)><EC.)QSC)c’J C.)c’J wII<a) C) C)> (D03 U35(P .o.zx(cm)0-100-200-300-400-500Levelz(cm)I0.59H0.49G0.40F0.31E0.22D0.19C0.18B0.14A0.1290.108-0.157-0.196-0.245-0.344-0.433-0.522-0.611-0.70Bowabovex=22.85cmy=0.00cmL=45.70cmU=12.86cm/secy(cm)-100-50050100Solid symbols : Present MethodHollow symbols: Ma & Tulin (Expt) 2 = 4.25cmy/h,I 2 3 4 500 50 100 150 200 250 300Figure 3.7: Location of interfacial crest lines, F1 = 2.95.Solid symbols: Present Methodhollow symbols: Ma & Tulin (Expt)y/h225I 2 3 4 520::0 50 100 150 200 250 300Figure 3.8 : Location of interfacial crest lines, F = 3.00.log (i1/h2) j112 = 4.25 cm]-0.50‘ Ma & Tulin (1992)O 75 • Present method-1.00 .-1.25• ,• log (x/h2)1.00 1.25 1.50 1.75 2.00 2.25 2.50Figure 3.9 : Decay of first interfacial crest, F = 2.95.37Symbol hid F, Methodlog(/h2) ‘ 0.61 3.20 Expt0.0 . 0.96 2.95 ExpI• 1.15 3.00 Expt.- 0.96 2.95 Present1.513.00 Present-0.5 -____________________________-1.0 .... •S••....• V•..•-15-2.01.00 1.25 1.50 1.75 200 2.25 2.50log (x/h2)Figure 3.10: Influence of upper layer thickness on first interfacial crest.38Chapter FourHybrid Boundary Element/Discrete Vortex MethodSeparation from a sharp edge results in the formation of a vortex sheet issuing from theedge. This vortex sheet can be modelled by a series of discrete vortices introduced one ata time into the flow field at regular intervals. The motion of each vortex is traced overtime using its convection velocity. If the vortex shedding is restricted to a region close tothe edge, the discrete vortex method can be viewed as the inner region solution in theoverall boundary value problem. This inner region solution has to be matched with theouter potential flow solution. The combination of boundary element and discrete vortexmethods provides this matching and at the same time do not require calculations inside thedomain.It is realized that the boundary element method is compatible with the discretevortex model in that each discrete vortex can be treated as an internal singularity whichcan be handled using an analogue of the Residue Theorem in complex analysis. Significantcomputational advantages result because of the relatively simple approach to handlingseparation at the sharp edges while working only with the boundary values. In thischapter, the hybridization of the discrete vortex and boundary element methods will bedescribed. The hybrid method was applied to the calculation of the forces on a flat plateoscillating in unbounded fluid (Wong 1993a).394.1 Discrete vortex methodThe discrete vortex method is a time stepping procedure which models the shear layerissuing from a sharp edge using discrete vortices (Wong 1990). In the formulation, aninfinite sharp wedge is mapped onto a half plane via a Schwarz-Christoffel transformation:Z=- (4.1)where M is a scaling constant between the physical Z and the transform planes and= 2 — is a parameter dependent on the internal angle iS of the wedge.Physical Plane Transform PlaneFree stream Free streamvelocity (V)/velocity (V)[tvortexatz, [ievortexat rvoexaist,ngthy strength strengthFigure 4.1 : The physical Z-plane and the transformed -pIane.The non-dimensional external flow velocity, v sin(2rrt), wheret = t / T (T being theperiod of oscillation) is in a direction normal to the wedge bisector. The strength of thevortex at z and its corresponding point is given by i. see Figure 4.1. The complexpotential in the -p1ane, with n vortices in the field, is then given by:40F()=iV÷_Lrkln (4.2)k-OApplying Routh’s correction (MauLl, 1979) and nomialising all variables, after Wong(1990), result in an expression for the complex conjugate velocity at a point j in thephysical plane as follows:w.=S.{v+ — _ )—.!L( ‘_— )] (4.3)j k-O,k.j j — + 2t + 2where the last term on the right hand side is due to Routh’s correction. The length scalesused for the normalization were derived in terms of an attached flow velocity scale in thephysical plane. These are (from Wong 1990):L= (1(D)k12lL = (M)1(D)1/(2k_1)L = (M)i(D)1/(2_1)V (4.4)where V0 the amplitude of the oscillatory velocity and D is the maximum distancetraversed (V0T). Assuming that all vortices convect with the flow, the convectionequation in the physical plane is given by:Zk ‘‘k (4.5)A new vortex (called the nascent vortex) is introduced along the external wedge bisectorat the end of each time step. The distance of the nascent vortex from the sharp edge is41taken to be a function of the time step size, the wedge angle, the free stream velocity andthe strength and location of all the other vortices already in the flow field. The expressionfor the initial Location of the nascent vortex, denoted by z0, can be found through thefollowing equations (derived in Wong, 1990):z0 = (KAr)tm (4.6)where m=2). -123t)(4.7)where v is defined in Equation (4.10). The strength of the nascent vortex can beobtained by satisfying the Kutta condition at the edge. The Kutta condition is satisfiedwhen the flow velocity at the shedding edge in the physical plane is finite. Consequently,due to the singularity of the transformation at the shedding edge, the velocity at thecorresponding position in the transform plane is zero. As a result, the followingexpressions can be obtained:dFQ)(4.8)d [ k-O 2rt k )jgiving •ftj VR(;)(4.9)where v = nsin(2ct) — I kRe(k) (4.10)With the strength and location of the nascent vortex determined, the algorithmthen proceeds to the next time step. Note that the location of the first discrete vortex shed42into the flow is defined by Equation (4.6). A vortex decay mechanism is used to accountfor the vortex pairing process described in Graham (1980) and provides a model forvortex diffusion. For this, a simple exponential rule reduces the effective strength of eachvortex to one half of its original value at the end of one time cycle (Equation 4.11). Thevortex decay mechanism is thus an attempt to simulate the loss of total vorticity throughthe cancellation effects in positive and negative strength vortices coming together (inpairing) and viscous influences, both of which are absent in the discrete vortex method.y(t) = lo [1— exp(—0.6932 / t)] (4.11)Although violating conservation of vorticity, this vortex decay model is successfulin maintaining computational stability over a large number of flow cycles. Also, reducingcirculation in the wake limits the magnitude of the vorticity introduced. Consequently, theotherwise over-predicted vortex induced forces are kept at levels which are comparable toexperimental values (as reported in Wong 1990). An example of the calculated vortexroll-up for a flat plate at t=1.O is given in Figure 4.2.Figure 4.2: Vortex sheet rollup at t = 1.00 for a flat plate.Second spiralFlat Plate a‘jo.0.0o0First spiral434.2 Vortex shedding in boundary elenieiit methodIn this section, the inclusion of vortex shedding effects into a boundary element model isdiscussed. We define a domain D bounded in S which includes branch cuts C1 and C2and a circular path S to exclude a vortex singularity located at the point K as in Figure4.3. The circle enclosing point K has a radius E.For the velocity potential (p which satifies the Laplace equation, the divergence andGreen’s theorems can be used to obtain:Figure 4.3 : Exclusion of a vortex singularity in domain D.——G-dS=O (4.12)S on OnThe weighting source function which also satisfies the Laplace equation is taken to beG = ln(r), r being the distance between the point under consideration and any other pointon the boundary.44The integral involving G on SE around the vortex at K vanishes. Thus as S — 0,the line integral around SE, in Equation (4.12), becomes:lim k -—dS (4.13)uflwhere y k is the strength of the vortex at K. The boundary S is discretized into N straightline elements.Putting H = --, the discretized form of the integral equation for the point i onthe boundary with N elements enclosing a domain with M vortices within can be writtenthus:2JHIJ +it = (4.14)j—1 k—N+1 j—1 jwhere i =1 ... N + M, ô is the Kronecker delta function and G and H1 can be foundnumerically (see Appendix B).Note that the summations involving G and H do not include the ‘elements’ aroundvortex singularities (j> N) as the former need not be evaluated since (a / on)1 = 0 fori> N while the latter can be shown to vanish as the radius £ —* 0. It is also noted thatthe branch cuts C1 and C2 do not contribute to Equation (4.12). This result is expectedas the treatment of the vortex singularities is similar to the way the potentials at pointsinside the domain are found. The above problem is overspecifled in that there are N+Mequations in N unknowns and only the first N equations are used.In time domain simulations, the boundary values and vortex strengths are to beprescribed for the next time step (t + bt) in order to march forward in time. The matrices45G and H need to be calculated only once for a fixed boundary. The strengths andpositions of all discrete vortices in the flow at time (t + ö) were found using the discretevortex method described in the previous section. The potential and velocity induced bythe discrete vortices in the flow on any boundary node j is given by:= YkOkj(4.15)k-i 2(-) = Yflj (4.16)vj k-itkjwhere O is the angular position of the node j with reference to the kth vortex. This vortexhas strength k and is at a distance Rk3 from the node j. In Equation (4.16) n is the nomialvector at the node j.Consequently, if we include the vortex induced normal velocities on the body,Equation (4.14) can be rewritten in the following form:=+ () ] — (4.17)for i = 1 ... N. Normal velocities are prescribed on the body. The system of N equationscan be solved for the N unknowns . The dimensional hydrodynamic force per unitdepth acting on Sb at any given instant is taken to be:Fb =p.ndS (4.18)where p is the hydrodynamic pressure (from Bernoulli’s equation) acting on the discreteelement dS with normal vector n. The calculated value of Fb is normalized using pV02d togive a non-dimensional force coefficient C where:46C =--FIn the above, Fb = j23tp(V0 I T)d }F8 for non-dimensional FB.4.3 Forces on sharp edged cylinders(4.19)Although this external flow problem can be solved using the discrete vortex method alone,this example serves to illustrate the time domain implementation of the procedure given inthe last section. Published numerical and experimental results, for example those inGraham (1980), Kudo (1981) amd Keulegan et al. (1958) could be used for comparisonpurposes: not to judge the relative merits of the numerical schemes proposed by otherresearchers but as a test case for the present method.Figure 4.4 : Computational domain and boundary conditions.The simulation of vortex shedding from a finite wedge oscillating in unboundedfluid requires some length scale modifications since the discrete vortex method describedearlier is applicable to an infinite wedge. In other words, the vortex strengths andPn,Discrete vorticeso 0000000 0 Discrete vortices47positions calculated with the discrete vortex method need to be adjusted to account forchanges in Keulegan-Carpenter numbers (KJ. Graham (1980) showed that the ratio ofthe length scales of the infinite and finite wedges is proportional 10 (K)’ . Theresultant non-dimensional vortex strength (y B) and length (L8) scalings to be used forthe boundary element procedure tor finite K and body length (2d) is then given by:Y8 =y(2K) (4.20)kZB=z(2KC) (4.21)It should be noted that the above equation is general in that the vortex sheddingedge angle 0) is embedded in the parameter k. The computational domain arid relevantboundary conditions are given in Figure 4.4. On the body, the non-dimensional nomalvelocity is given by:ap . o— = sin(2tt)cos(—).n (4.22)an 2For this application, all the vorticity shed from a particular edge is placed at a pointvery near to that edge. In other words, we ‘average’ all the vortex positons and theinfluence of all discrete vortices on any given boundary node is calculated through therelative positons of that node and the edge. With that, the value of 0 in Equation (4.15) istaken to be ±(t/2) depending on whether the top or bottom vortex shedding edge (Figure4.4) of the body is considered. In Equation (4.16), R is taken to be the distance from theedge to the node under consideration. This simplification is acceptable for small K sincemuch of the vortex shedding effects are confined to a region near to the edge.48Values of cF for a flat plate oscillating (ô = 0) at K = 6.6 are plotted together withexperimental results obtained by Keulegan et al. (1958) in Figure 4.5. The present methodpredicted C values that are in excellent agreement with experimental results. Drag (Cd)and inertia (C1) coefficients are found by taking Fourier integrals of Morison’s equationover a time cycle and given thus:3 .mn+1)Cd=C sin(2rt)dt (4.23)4nc2K0 j.(u+l)C5= 2 J C cos(2mt)dtwhere n is an integer representing the cycle number.Figure 4.5 : Forces on an oscillating flat plate, iCc = 6.6.(4.24)Calculated drag (Cd) and inertia (C1) coefficients are presented in Figure 4.6 andFigure 4.7. Predicted Cd for the flat plate agree well with numerical and experimentalPresent methodKeulegan et at (expt)Ea)C)a)0C-)a)C.)0LL3.0 4.0Number of Cycles49results from Kudo (1981) and Keulegan et al. (1958) respectively. At K above 3.0,predicted d appears to underestimate experimental values. However, Obasaju’s measureddata, as reported in Graham (1980), showed Cd values lower than those found using thepresent method. Therefore, given the variation of prevailing conditions under whichdifferent researchers obtained their experimental results, the present method predicts dragforces on the normal flat plate very well.15.000Present Method —4—Keulegan et al (sept) o12,500 Kudo (nianencal) A—10.000C)750080) Oo 02 501)Kc0000 • • • I0 5 10 15Figure 4.6 : Drag coefficients for the flat plate.Values of C. obtained using the present method show an upward trend withincreasing K which is similar to the numerical results presented by Kudo (1981). On theother hand the experimental results of both Keulegan and Carpenter (1958) and Obasaju(presented in Graham, 1985) showed an upward trend for C3 up to about K0 7 and fromthere a decreasing Ca, up to K. 12 . Forced symmetric vortex shedding from the twoedges of the plate is thought to be the reason for the poor match for C3. Graham (1985)50reported much better results when vortices are allowed to shed asymmetrically through theintroduction of a small perturbation at the initial stages of the simulation.5.000 Present Method - + -.Keulegan et at (sept) 0Kudo (numericat) A — — —4001) • _r-.-a)3000 —Cl) 02000 80 0o01.000 0Kc0.OOC....0.0 2.5 5.0 7.5 10.0 12.5 15.0Figure 4.7 Inertia coefficients for the flat plate.Finally, it should be mentioned that drag and inertia coefficients become morerealistic through the incorporation of vortex shedding forces. Potential theory alone doesnot give any drag and inertia values for attached flow can be shown to be Ca0 = 1.0 for theflat plate.Use of the same transfonn equation in the discrete vortex model for different edgeangles provide a useful and flexible tool in the investigation of a wide range of flowproblems involving sharp edge vortex shedding. Although the discrete vortex model canbe applied in conjunction with other numerical methods for the given class of problems, itis uncertain whether such an exercise would be more fruitful than working directly withthe Navier-Stokes equations for viscous incompressible flow since the advantage of’working only with elements on the boundary is lost.51Chapter FiveApplications in Nonlinear Slender Body HydrodynamicsThe wavemaker approach opens up the possibility of attaining a nonlinear solution towave generation problems involving slender forms in steady and unsteady motion. Thefully nonlinear free surface conditions, implemented using the Eulerian-Lagrangianapproach, was used for all the work reported in this chapter.In the next section, the utilization of the wavemaker algorithm developed inChapter Two for the prediction of the waves generated and forces acting on a Wigley hullis described. This algorithm was then extended to include the effects of vortex sheddingfrom the keel of the Wigley hull advancing obliquely. The computed results are presentedand assessed by comparison with experimental and numerical results of other researcherswherever possible.5.1 Wavemaking resistance of a slender bodyThe computational boundary, identical to that given in Figure 2.2, is discretized into Nquadratic elements, that is the variation of potential and normal velocity over an elementlength dS is assumed quadratic (see Appendix B). The hull is described by:y = Y(x,z)= (5.1)52a4 —U(Ycosa±sinct)‘/1+Yz2UzyWigley Hullwhere B, L and d are the beam, length and draught of the ship. The normal vectors on thehull, in the slender body approximation, is given by:(_y ,±i,y.(n1,ny,n)=j L) (5.2)where the subscript in Y refers to the partial derivative with respect to that variable.Figure 5.1: Wigley hull dimensions and coordinate system.The normal velocity on the obliquely positioned hull, with yaw angle a, is thengiven by:(5.3)For this section, a Wigley hull advancing without yaw will be considered (a = 0).The nonlinear free surface conditions given in Equations (2.9), (2.11) and (2.12) wereused. Neumann conditions are specified for the rest of the boundaries. At time (t+ôt), theLength 200 cmBeam 20cmDraught 12.5 cm53Adams-Bashforth predictor is used to estimate the velocity potential and elevation of thefree surface. Normal velocities on the free surface are then calculated as part of thesolution. From these, using the procedure described in Chapter Two, the wave elevationis updated and the rate of change of the potential found using Equation (2.9). TheAdams-Moulton corrector is then applied for finding the free surface velocity potential.Solving the problem using the corrected potentials then completes the calculations for thattime step.There are weaknesses associated with applying a wavemaker approach withoutenhancement, as seen in the advancing wedge application in Chapter Two. In the verticaldirection, the intersection between the ship side and the free surface need to be treatedwith care. In addition, there is also the question of numerical accuracy as well as stabilitydue to the high wavemaker acceleration at the start of the time marching for most shiptypes (with the possible exception of cusp-ended ships). The problem is more serious fornon-vertical ship sides. To overcome this, the normal velocity at the ship-waveintersection point is simply extrapolated from values calculated at free surface nodesnearby (see Appendix B for more details on the extrapolation procedure). In addition, theevolution equations for the free surface nodes involve a finite difference approximation tothe tangential velocity component as described in Chapter Two. These gave rise to higherfrequency errors, which over time resulted in stability problems. A rediscretizationscheme was applied to the entire free surface boundary to filter out these high frequencyerrors. The regridding procedure was based on linear interpolation over neighbouringnodes returning equally spaced nodes on the free surface boundary once every time step.54A smoothing scheme identical to that used by Longuet-Higgins and Cokelet (1976) wasapplied after the rediscretization process. The value of 4 at the fluid-body intersectionnode for the latest time step was used as a means of accounting for the variation of thelongitudinal velocity over the body. This modifies the free stream velocity (U) incident onthe hull to give U’ which was then applied to Equation (5.3). The following relationship,based on the value of 4 at the fluid-body intersection, was used:U’ = Ucosa+X(SB flS)[1_.]ap -(U’Y1 ±Usina)g1+y2The principal dimensions of the Wigley hull used in all the calculations in thisChapter are given in Figure 5.1. The impermeable far wall was placed 30 beams awayfrom the body. On the free surface, 60 equal length elements were used. For the freesurface plots, one sweep over the ship length was done in 50 steps. However, for thepressure distribution and forces acting on the hull, extensive numerical experimentation onthe number of steps required (NS) per sweep resulted in the following formulae:NS = INT[108(1.1 — FR X — 4FR)] 0.20 FR 0.25NS = INT[54(3 — 4FR)] 0.25 <FR 0.40 (5.5)More steps were necessary for the lower Froude numbers due to the relativelylarger variations in velocities and pressures in the longitudinal directional. In all cases, thevalue of NS was limited to a range of (70 NS 200). The lower bound on NS was55imposed to ensure adequate resolution in the length-wise pressure distribution while theupper bound limits computation time per sweep over the hull.Initial conditions, that is wave elevation and potential and their rates of changewith respect to x at the start of the time marching were obtained using Dawso&s (1977)method. The Dawson method is an extension of the Hess and Smith (1976) double-bodymethod, handling the linearized free surface condition via a streamwise marchingtechnique.The hydrodynamic forces can be found, from Equations (2.6), using the pressurefrom the Bernoulli equation. Substituting p for the pressure:p = _P[U(coscz + sina) ++ () ) + z] (5.6)L/2F1 = fdxf p.n1dS (5.7)—L/2SHwhere i is the required direction (x, y or z) and 5B represents the girth of the hull sectionunder consideration. The wavemaking resistance and sinkage force are then F1 and Frespectively. The trim moment (corresponding to moments about the y-axis) can be foundusing Equation (2.7):L/2M = fdxf5 p.(nz—nx)dS (5.8)—L12The wave profile along the ship side at Froude numbers F = 0.267 and FR = 0.316are presented in Figures 5.2 and 5.3. For both cases, the results calculated by the pesentmethod are generally in good agreement with experimental and numerical results obtained56by other researchers. The experimental results for FR = 0.267 was obtained from theDepartment of Naval Architecture and Ocean Engineering at the Yokohama NationalUniversity Hayashi and Kunishige, 1988), referred to as YNU hereafter. Measuredvalues for a 6m Wigley huil model reported in Maruo and Song (1990) are used in Figure5.3, for FR = 0.316.Figure 5.2 Wave on side of Wigley hull at F = 0.267.The delay of the first wave peak near the bow may be due to the treatment of thewave-body intersection point for stability reasons. Inaccuracies in the value of modifiedfree stream velocity given in Equation (5.4) may be another cause for the delay. These led2’,-1/L Present method (NonlEnear). Measured (YNU)Maruo & Song (1990)0.040.030.020.01 --0.00 - •-0.01-0.02-0.03FP AP—0.04-1.25 -1.00 -0.75 -0.50 -0.25 -0.00 0.25 0.50 0.75 1.002x/L57to the delay in the crests and troughs over the rest of the hull length. However, since thelag is only about 1-3% of the ship length, this trade-off is deemed acceptable.In contrast, Wong (1993) has shown that use of the simple wavemaker scheme,with linearized free surface conditions and without the enhancements described in thissection, results in a wave profile that totally missed the variations over the middle half ofthe ship length.Figure 5.3 : Wave on side of Wigley hull at FR = 0.316.2r1/L Present method (Nonlinear). Measured (6m model)Maruo & Song (1990)0.040.030.020.01-0.00-0.01-0.02-0.03-0.04-1.25 -1.00 -0.75 -0.50 -0.25 -0.00 025 0.50 0.75 1.002x/L- FP AP58The contour plot of the wave field in Figure 5.4 agrees fairly well with YNUmeasured data (Figure 5.5). A perspective view of the wave field for the Wigley hull atFroude number F = 0.267 is shown in Figure 5.6.Figure 5.4 : Contour plot of wave field for Wigley hull at FR = 0.267.Figure 5.7 shows the distribution of the pressure, non-dimensionalized using0.5pU2 on the hull surface for F = 0.267. The high pressure region in the bow is clearlynoticeable.y (ciii) Level0 5.500 4.50B 3.50A 2.509 1.508 0.507 -0.506 -1.005 -2.004 3003 -4.002 -5.001 -6.000FP 2x/L AP59N.Lu’I-’cr1cr1LulbQ 1Jo\1<L__J I).Figure 5.5 : Measured and calculated wave field (YNU) for Wigley Hull, a = 00.60x1 IJQ VI a’ I = Tj p 0•T1-UQQ U’0%0I—0 OQ CD I0 -U I >The wavernaking resistance, non-dimensionalized by O.5pU2L2, is compared toexperimental and numerical results obtained, elsewhere (Figures 5.8 and 5.9). Valuesobtained using linearized free surface conditions are highly inaccurate (also reported byAllievi, 1993). The general agreement with measured data is considered to acceptable. Itshould also be mentioned that even experimental results from different towing tanksexhibit significant variability depending on the conditions under which the data wasmeasured.Aanesland (1989) showed ‘averaged’ values from towing tank databases for caseswhere the Wigley hull was towed fixed and where the vessel was allowed free to turn andsquat (Figure 5.9). Aanesland’s numerical results were obtained using a three dimensionalpotential flow algorithm, based on Rankine sourcç distribution and strearnwise marchingRJ(.5pU2L)x 1O Numerical (Maruo & Song)Experimental(YNU)0.50________Wave cut analysis ‘Nu)Present method (Nonlinear)0.400.30 ‘ .,.0.20. -, //0.10FR.1... .1.. ..I0.15 0.20 0.25 0.30 0.35 0.40 0.45Figure 5.8 : Wavemaking resistance against YNU results.62schemes similar to that of Dawson. Kelvin sources, distributed in the outer region, werematched to the inner region (in the vicinity of the hull) to yield a radiation condition at theinterface between the inner and outer regions. Calculated results from the present methodare cast among these in Figure 5.9, showing the present results to be within the range‘acceptable’ values.Averagefree (Expt)Averagefixed (Expt)v Aanesland (Numerical)• Present method (Nonlinear)Figure 5.9 : Wavemaking resistance against Aanesland (1989).Sinkage and trim for various Froude numbers are shown in Figure 5.10. Sinkage isnon-dimensionalized using U2 / 2g and trim is given as a percentage of the ship length.These results are compared to averaged measured values and those calculated byAanesland (1989). The non-dimensional sinkage calculated by the present method are inRJ(.5pU2S)x 1O3.02.52.0 ./• ... I• •-..1.5 ,.. -—‘Il.0 V..0.5- .-FRoci0.150 0200 0.250 0.300 0.350 0400 0.45063good agreement with experimental data but tend to underestimate at the higher end of theFroude number range. On the other hand, the non-dimensional trim tends to be overestimated in the Froude number range of 0.27 to 0.35.Sinkage/(U2/2 )0.00-0.02-0.04:-•-. ., .V-0.06 --0.08—010I... .1 I.—..• 0.10 0.15 0.20 0.25 0.30 0.35 0.40Trim/Lx 100 (%)2.00 v1.501.0070.50•0.00 :.FR-0.50 I0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45Figure 5.10 Computed and measured trim and sinkage for Wigley hull.Experimental (Lower Limit)Experimental (Upper Limit)Calcutated (Aanesland)• Present methodFR0.45645.2 Wavemaking resistance in shallow waterThe present algorithm was also used in the calculation of wavemaking resistance,sinkage and trim characteristics of a Wigley hull in water of finite depth. The hull and tankdimensions used were identical to those in the previous section. However, the number ofsteps required (NS) was modified to reflect the fact that higher resolution is necessary forshallower water. Hence, with NS defined in Equation (5.5), the required number of steps(NS*) is given by:NS* = 6 NS) H/d <6.0, 70 NS* 400 (5.9)HidIn Equation (5.9), H and d are the water depth and draught of the shiprespectively. Note that no modification to the value of NS was done for Hid 6.0.Results for Froude number 0.316 are compared to those obtained by Kim and Choi(1993), who used a 3D iterative nonlinear approach, in Figure 5.11. As expected, theincrease in wavemaking resistance near the critical depth Froude number ofFH = U / JiT = 1.0 (corresponding to Hid = 1.60, where H is the water depth) is large.The occurrence of the maximum wavemaking resistance turned out to be at a depthFroude number of slightly above the critical value (at 1.07 where H/d = 1.4). In any case,the two sets of results in Figure 5.11 are generally in good agreement except that thepresent method predicted a slower decline in wavemaking resistance with increase inwater depth. The effect of shallow water on squat and trim are plotted in Figure 5.12.65RJ5pIYS x 1O15Calculated (Kim and Choi, 1993)Present method105.C .1....0.0 2.5 5.0 7.5 10.0H/dFigure 5.11 : Wavemaking resistance for Wigley hull in shallow water at FR = 0.316.Non-dimensionalTrn and Sinkage0.75Sinkage/(LJ2/2 )Tnm/LxlOO%0.500.250.00•. . . .O2 I I0.0 2.5 5.0 7.5 10.0H/dFigure 5.12 : Effect of water depth on squat and trim for Wigley hull at FR = 0.316.665.3 hybrid method in yawed ship calculationsThe steady motion of an advancing ship generates vorticity. The strength of the vorticesshed determines their significance in the computation of the forces and moments acting onthe ship. The influence of vortex shedding is, in general, not very important in the casesinvolving a ship in steady forward motion in calm water. This is exhibited by the usefulresults emerging from pure potential flow methods, for example those from Dawson’s(1977) method. However, for a ship advancing with a yaw angle, vorticity shed from thebilge keels and keel is not necessarily neglible.In the context of the slender body approximation, the translation of the bodysections in the cross-flow plane as the marching proceeds presents a vortex shedding sharpedge at the keel. The influence of such vortex shedding on the forces and moments actingon a yawed Wigley hull and the effects on the free surface wave field will be evaluated inthis section. However, before using the hybrid discrete vortex-boundary element methoddescribed in Chapter Four, it is first necessary to explain the vortex generation mechanismand derivation of the length scales involved in the vortex shedding algorithm.Generally, vortices were taken to be shed from the keel of the yawed hull only.Bilge vortices were not accounted for since separation points cannot be convenientlydefined and the method described in Chapter Four cannot be applied. Also, vorticesoccuring at the bow were not modelled. In addition, secondary vortex shedding wasignored. In other words, the vortex shedding procedure was started at the commencementof the time step marching for the wavemaker algorithm. The introduction of potentialvortices into the flow field results in their influence on the potential and normal velocity at67each boundary node. On the other hand, the assumption that dissipation does not occur inthe overall flow allows the use of Equation (2.8) for the evaluation of pressure on the hullboundary.Before the start of the marching procedure, it was assumed that there are novortices in the cross-flow plane. As the time marching proceeds, the hull transversesection changes in shape and the keel can be viewed to be translating in the horizontaldirection. From the ship reference frame, a lateral velocity of constant maiitude V(given in Equation 5.10) is seen at the keel. A growing spiral vortex hence emits from thekeel as the marching proceeds towards the stern. This spiral can be modelled through theintroduction of discrete vortices, at increments of one per cross-flow plane. It wasnecessary to substitute the oscillatory velocity in the discrete vortex algorithm (ChapterFour) with the constant velocity V. In addition, the ‘wedge’ was taken to be a flat plate.This means that, although the internal angle bounded by the body section at the keel variesalong the x-direction, the vortex shedding edge was taken to be of zero internal angle.This is consistent with the slender body approximation.At each step, the vortex influence on the boundaries were evaluated. These werethen inserted into the modified system of equations (given in Equation 4.17 ) and thesolution of the boundary value problem obtained.The ‘free stream’ velocity for the vortex shedding algorithm, as seen by the keel, istaken to be the translatory velocity of the body section such that:V=Usinct (5.10)68where V here refers to the free stream velocity in the vortex shedding model in ChapterFour.Length and time scales are required to extract the necessary vortex positions andstrengths from the self-similar vortex produced by the discrete vortex method. This canbe done by defining the scaling factors, using Equation (4.4), and defining M = 1., D =I.sin a and V0 = U sin a. Therefore, the length scales for this application becomes:L = (Lalna)M21)= (Lsina)21)= (Lsina)hI(2?_1)Usina (5.11)The non-dimensional time step size for the discrete vortex method is given by:dr=dx(5.12)L cos aFrom the above, the non-dimensional positions and strengths of each vortex in thecross-flow plane can be dimensionalized for the solution of the overall problem in thefollowing manner:YV_YB+Im(ZV)LZz, =—d—Re(Z)L= FL (5.13)where Z is the complex variable Z, = X,, + iY in the vortex physical plane and f is thenon-dimensional strength of any given vortex. In Equation (5.13), yo and d refer to thelateral location of the ship’s keel and the ship draught respectively. The hybrid method69M ,k‘çEkkj‘Pkj-v—i 3tcan then be used, through Equation (4.17), for the solution in the cross-flow plane. Thepotential and normal velocity induced by a vortex on a mid-element node are given by thefollowing pair of equations:(5.14)(oP’_. 1k sin (5.15)\nJk k-l2’kjwhere O and R113 are defined in Figure 5.13, M is the total number of vortices in the fieldand 3 is the difference between the angle of the vortex induced velocity and the elementangle. For nodes at the ends of elements, the angle can be taken as the average of thetwo neighbouring elements. Note that the angle 3 is simply given by t3 - (0 + /2).VRkIVortex0Figure 5.13 : Vortex induced normal velocity on a boundary element.70The far wall was taken to be 30 beams away, with 60 free surface elements oneach side of the body. The number of steps required (NS) for one sweep over the hulllength was determined by Equation (5.5). Yaw angle is takçn positive in the anti-clockwise direction from the negative x-axis, resulting in the bow being at y = -(112) sin a.The keel thus traverses from that position to y = (L/2) sin a at the stem. The yawmoment is given by:1/2M= fdxf p.(nyx_ny)iS (5.16)—1/2BFigure 5.14 compares the calculated free surface profile, with and without vortexshedding, along the side of the Wigley hull at FR = 0.267 with the measured profiles fromHayashi and Kunishige (1988) and calculated ones from Maruo and Song (1990). Theinfluence of the keel vortices on the free surface profile were minimal. Here, the influenceof the longitudinal velocity variations are more important than for the no yaw case.Consequently, the resulting wave profile in Figure 5.14 is less accurate than if there wasno yaw (compare with Figure 5.2). It is noted that the wave profiles were not correctedfor heel, caused both by the potential tiow as well as the vortex shedding, which may besignificant for a ship in yawed motion. The calculated wave field is presented as a contourplot in Figure 5.15. This is in good agreement with measured data from YNU, see Figure5.16.A perspective view of the wave field, with the z-scale exaggerated, is presented inFigure 5.17, while the pressure distribution on the hull surface is given in Figure 5.18.Plots in Figures 5.15, 5.17 and 5.18 are derived from results for the Wigley hull at a = 10071and with vortex shedding from the keel. A typical computer run, with 248 boundaryelements and marching in 100 steps over the hull length, requires approximately 2 hourson a 486 DX2 66 Mhz personal computer, 19 minutes on a HP 715 (75 Mhz) and slightlyunder a minute on the Fujitsu VPX24O/10 supercomputer. Most of the developmentalwork on the programs described in this chapter was carried out on the Fujitsusupercomputer.Present method (With vortex shedding)2ii, L Present method (No vo,tøx shedding)• Measured (YNU)0.08 Maruo&Song (1990)0.06002Windwardsiindwardside-0.04-006 FP AP-0.08 ..I... .-1,25 -1.00 -0.75 -0.50 -0.25 -0.00 0.25 0.50 0.75 1.002xJL2ri/L0.07 Lee side0.05-0.02 .-0.05-0.07 .......................-1.00 -0.75 -0.50 -0.25 -0.00 0.25 0.50 0.75 1.002x/LFigure 5.14 : Wave profile on sides of Wigley hull at FR = 0.267, a = 100.72y (cm)5025_ ‘\‘ Starboardh24 —z._ ‘,‘, ifFP°.25 /1 2 Port1 -025-’ -‘ 1.250.75 -0.5 “ -0.25 / \ 90.25 0.25 0.25 / 0-25 0.5 \ \-0.25-502L-100 -50 0 50 100Figure 5.15 : Wave field for Wigley hull at F = 0.267, a = 10°.73Figure 5.16: Measured wave field (YNU) for Wigley Hull, a = lO74CD cM CD CD C,;0 CD C,;.CD CD z C Iz 0 CD N Cl) C) CD CD x CD CD CD CDLee sidelFigure 5.18 : Pressure distribution for Wigley hull at FR = 0.267, a = 100.lWindward sidelFree surface levelFPFree surface level/APr 0.45-O.15 -fl2lAP1 .490E-80.051 .490E-8FP76Figure 5.19 : Forces and moments for yawed Wigley hull.Calculated values of longitudinal and lateral forces (F1 and F) and the yawmoment (M1) with, and without vortex shedding, for various angles of attack are nondimensionalized and plotted in Figure 5.19. Also plotted in Figure 5.19 are the numericalresults obtained by Maruo and Song (1990), using a Kelvin source based slender bodyalgorithm (linear and without vortex shedding). The present method gave nondimensionalized longitudinal forces that are comparable to those of Maruo and Song(1990) at small angles only, being much larger as the yaw angle increases. On the otherhand, the present method returned values of lateral forces and yaw moments which areF/(.5pU2L)x 101.501.000.500.OC(1 00M21(.5pUL3)x 1020.750.500.250.005 105 10 15 20 26a (degrees)F1(.5pU2L’)x 1021.5015 20 25a (degrees)1.000.50a (degrees)LegendSolid Itnes: Present Method (no 4iortex shedding)Symbols: Present Method (with vonex shedding)Dashed Lines : Maruo & Song (Numerical)77consistent with those of Maruo and Song.. There is, unfortunately, no publishedexperimental data available for assessment of the calculated forces and moments.As would be expected, the influence of vorticity on F is not significant. However,the vortex effects on the lateral forces and yaw moments are obvious from Figure 5.19.Vortex shedding induced forces which augment the lateral forces predicted when novortices were included. This can be clearly seen in Figure 5.19 where the lateral forces arehigher with vortex shedding for all yaw angles, the difference increasing to about 33% at20 degrees leeway. As for the yaw moment, the influence of vorticity increases as theleeway angle increases, the difference being 34% at 20 degrees. Yaw moment is taken tobe positive in the anti-clockwise direction. Due to the higher rate of change of totalcirculation with respect to x forward of amidships (see Figure 5.20), the predicted yawmoment for the vortex shedding cases were higher. In Figure 5.20, it can be seen that thegrowth of the total vortex strength is rapid near the bow and levels off towards the sternregion.Figure 5.20 Development of total vortex strength over the hull length.Fx 102 (m/s)Q = 20°i= 15°0.0u= 10°a =5°-1.0 -0.5 0.0 0.5 1.02L78Figures 5.21 and 5.22 show the keel vortex evolulion at 8 sections over the lengthof the Wigley hull at FR = 0.267, a = 10°. The sequence shows a representation of thevortex roll-up process as the keel traverses laterally.Figure 5.21 : Keel vortex evolution for Wigley hull, a = 100 (forward).The hybrid method just described produced significant vortex effects on the forcesand yaw moment acting on a slender body advancing obliquely, especially at higher leewayangles. Assessment of the accuracy of the computed values by comparing with measureddata cannot be carried out due to the unavailability of the latter. However, the resultinglateral force coefficients are of the same order of magnitude as the experimental data givenby Bradbury (1986) for a simplified ship-like geometry.y (cm) y (cm)-5 0 5 10 15 5 10 15= -0.79] 12x/L = -0.40]12Q’L = :0.601 12x/L = -0.20179Figure 5.22: Keel vortex evolution for Wigley hull, a = 100 (aft).y (cm)12x/L = 0.051I2xIL=0.37112x/L = 0661y (cm)-15 -10 -5 0 5 10 1512x/L = 0.96]80Chapter SixConclusions and Recommendations6.1 ConclusionsThe present study on numerical fluid dynamics has resulted in novel and practicalapproaches which can generally be categorized into the following three areas.6.1.1 Slender body formulation and implementation—‘The slender body formulation has been used to transform the three dimensional steady shipwave problem into a series of two dimensional wavemaker problems. This algorithm wasgeneralized for application to ships advancing at small yaw angles. The method was firstapplied to the study of the bow waves generated by an advancing wedge in homogeneousfluid, then to a spheroid in a dual-density medium, and finally to an advancing Wigley hull.Numerical implementation was achieved through the boundary element method, usinglinear and quadratic elements. This relatively straightforward approach gave results thatcompared very well with published experimental and numerical data.The fully nonlinear free surface conditions were successfully implemented for thecase of a Wigley hull advancing in deep and shallow water in Sections 5.1 and 5.2respectively. The fluid-body intersection singularity problem was avoided using anextrapolation procedure which modifies the influence matrix (Appendix B). Adiscretization scheme was devised to redistribute free surface collocation points. As aresult, the stability of the nonlinear free surface evolution was maintained. A criterion,based on extensive numerical experimentation, was developed for the determination of the81required step size for the marching procedures used. Finally, the calculated wavemakingresistance, sinkage and trim are in good agreement with published measured values.6.1.2 ships in density stratified waterA novel procedure which is able to handle dual (and in principle multi-) density media wasdeveloped and presented in detail. This scheme was successfully utilized in the study ofthe interfacial waves generated by a body moving slowly near and in a density interface.The computed results were in excellent agreement with experimental data. Because of therelatively small influence of the free surface, it appears that the slenderness ratio constraintof the method could be relaxed somewhat. This was apparent in the results given inChapter Three where a prolate spheroid of B/L = 0.20 was used with success. Thisprocedure is general in that it can also be applied to the study of internal wave problemswithout the presence of a moving body.6.1.3 Hybrid discrete vortex-boundary element methodA hybridization of the discrete vortex and boundary element methods was achieved. Thismethod was first tested on the forces acting on an oscillating flat plate, producing resultsthat are in good agreement with other numerical and measured data.The simplicity of the hybrid method accounts for some viscous effects in the flowwithout having to discretize and compute large amounts of data inside the control volume.The inclusion of vortex shedding in the flow while handling only boundary elements in theoverall numerical algorithm represents significant benefits in terms of computationalefficiency.82The hybrid method was successfully used in a nonlinear wavemaker program for aWigley hull travelling in yaw. Forces and yaw moments were calculated and presented inChapter Five. Lack of experimental data did not allow an objective assessment of theforces and moments predicted by the present method. However, these results wereplotted together with published calculated data (from a method using slender bodyapproximations to the Neumann-Kelvin problem). The present method predicted highervalues of longitudinal and lateral forces and yaw moments when vortex shedding wasincluded.The general wavemaker algorithm using nonlinear free surface conditions and thehybrid vortex-boundary element method provides an efficient and fast method for thecomputation of the wavemaking characteristics and resistance of a slender body advancingat small leeway angles. The reduction of the ship problem into a two dimensional onepermits a host of possibilities, which is the real motivating factor behind the developmentof such a scheme. The method can thus be regarded as a base method with a number ofpotential applications in the areas of ship motions and manouevring. Shallow water shipresistance and irregular or sloping seabed effects studies are also considered to be possibleapplications. The algorithms given in Chapter Five may be used to study the effects ofdemi-hull configuration on the wavemaking resistance of multi-hull vessels. In addition,the numerical study of yawed multi-hull hydrodynamics is merely a straightforwardextension of these algorithms.836.2 RecommendationsAlthough the wavemaker algorithm has been shown to be a useful tool in the calculationof hydrodynamic forces and moments on a moving ship, a number of refinements arerecommended for the enhancement of its accuracy and capabilities.The slender body approximation neglects velocity variations in the longitudinaldirection near the bow and the stern. This was accounted for in the present work througha simple parabolic variation of the perturbed x-velocity with depth. Also, even for cuspended ships, the initial step of the marching procedure starts with a finite width transversesection, thus violating the slender body assumption in the local sense. The highlyaccelerated flow at the bow and stern give rise to singularities at the fluid-bodyintersections. This problem was circumvented in the present work through anextrapolation scheme. Consequently, the wave elevation near the hull and wetted surfacearea are sensitive to the length of the free surface elements. This in turn affects theaccuracy of the pressure distribution on the hull surface. In addition, it was observed thatthere are errors in the prediction of the bow wave peak location and lack of recovery ofthe free surface elevation near the stern. Further study on the stability of impulsivelystarted wavemaker flow in relation to the above problems is likely to enhance the accuracyof the method.All slender body calculations reported in this thesis were carried out with the bodyfixed. It is known, from experimental work, that allowing the body the freedom to heel,squat and trim makes a significant difference in the measured resistance (see for exampleAanesland, 1989). In the wavemaker approach, this could probably be achieved by84marching through the length of the body more than once, with heel, trim and sinkagevalues from the previous sweep imposed in the following sweep.Boundary layers created by the ship are likely to alter the effective boundarylocation and hence the normal velocities required in the wavernaker algorithm. Theseeffects were not considered in this thesis.Another area that deserves attention is the determination of the effects of roundbilge vortex shedding. The difficulty in accounting for such effects in the hybrid methodas it stands is due to the fact that separation points are not well defined for roundedbodies. A practical method for the determination of the separation points in generalrounded surfaces would be most useful. This would in turn remove the limitation of thepresent hybrid method to sharp edge vortex shedding.Finally, model testing for the yawed Wigley hull would be extremely useful for thevalidation of the results given in Chapter 5.Extension of the present work for use in general hull forms would necessarily bethe next stage of the research. This would, however, require computer routines for thedetermination of hull geometrical parameters such as slopes and normals at any desiredtransverse section.85NomenclatureAP After perpendicular of the shipB Beam of the shipCa Vortex induced inertia coefficientCd Vortex induced drag coefficientCF Vortex induced force coefficientC Pressure coefficient p/(.5pU2)DIDt Material derivative of a variabled Draught of the shipF1 Densimetric Froude number relating to stratified flow problemFP Forward perpendicular of the shipFR Ship Froude numberG Free space Green function ln(r)g Gravitational constant, 9.81 rn/s2H Water depthh2 Depth of the upper layer in two layer fluidK Keulegan-Carpenter numberL Length of the shipL Length scale in the physical plane in discrete vortex methodI_c Length scale in the transformed plane in discrete vortex methodL. Strength scaling in the physical plane in discrete vortex niethodM Moments about the axes x,y,z (i = 1,2,3)NS Number of steps to march over the length of the hulln Normal vector of a pointp Pressure from the Bernoulli equationR Wavemaking resistance in newtonsr,R Distance between a source point and a field pointS0 Bottom boundary in wavemaker problemFar wall in wavemaker problemS Free surface boundary in wavemaker problemS8 Body boundary in wavemaker problemT Period of oscillatory flow in secondst Dimensional time in secondsU Uniform advance velocity of the shipVo Free stream velocity amplitude in discrete vortex methodw Complex velocity of a discrete vortexx Longitudinal axis in coordinate system fixed on the shipY Offset distance of the body surface from the centreliney Lateral axis in coordinate system fixed on the shipZ Complex variable in physical plane in discrete vortex methodz Vertical axis in coordinate system fixed on the ship86a Yaw angle of the ship in degreesInternal angle of infinite wedge in degreeso Internal angle of wedge in discrete vortex methodE Slenderness parameter B/LOp Difference in density between the upper and lower layersAS,OS Length of a boundary elementt,Ot Time step size in secondsElevation of the free surface from the undisturbed levelPerturbed velocity potentialTotal velocity potential in advancing ship problemWedge angle parameter in discrete vortex methodp Fluid density0 Relative angle between a discrete vortex and a nodal pointNon-dimensional time in discrete vortex methodComplex variable in transformed plane in discrete vortex method87BibliographyAanesland V., “A hybrid model for calculating wave-making resistance”, Proceedings ofthe Fifth International Conference on Numerical Ship Hydrodynamics, Japan, pp.657-666,1989.Allievi A., “On nonlinear free surface potential flow by a Bubnov-Galerkin formulation inspace and a semi-Langrangian semi-implicit scheme in time”, PhD Thesis, MechanicalEngineering Department, University of British Columbia, Vancouver, 1993.Bradbury W.M.S., “A Bilge Vortex Calculation”, Transactions of the Royal institute ofNavalArchitects, pp.189-199, 1986.Braathen A. and Faltinsen O.M., “Application of a vortex tracking method to rolldamping”, Advances in Underwater Technology, Ocean Science and OffshoreEngineering, vol.15, pp.177-193, 1988.Brebbia C.A. and Dominguez J., Boundary Elements: An Introductory Course,Computational Mechanics Publications, Southampton, UK, 1989.Calisal S.M. and Chan J.L.K., “A numerical modelling of ship bow waves”, Journal ofShip Research, vol.33, no.1, pp.21-28, 1989.Calisal S.M., Goren 0. and Okan B., “On an Iterative Solution for Nonlinear WaveCalculations”, Journal ofShip Research, vol. 35, no.1, pp. 944, 1991.Chan R.K.C., “Two-dimensional time-dependent calculations of large-amplitude surfacegravity waves due to a surface disturbance”, Proceedings of the First InternationalConference on Numerical Ship Hydrodynamics, pp. 315-331, 1975.Cheng B.H., “Computation of 3D Transom Stern Flows”, Proceedings of the FifthInternational Conference on Numerical Ship Hydrodynamics, Japan, pp.581-591, 1989.Chorin A.J., “Numerical study of slightly viscous flow”, Journal of Fluid Mechanics,vol.57, pp.787, 1973.Clements R.R. and Maull Di., “The representation of sheets of vorticity by discretevortices”, Progress in Aerospace Science, vol.16, no.2, pp.129-146, Pergamon Press,1975.Cointe R., “Nonlinear simulation of transient free surface flows”, Proceedings of the FifthInternational Conference on Numerical Ship Hydrodynamics, Japan, pp.239-250, 1989.88Dommermuth D.G. and Yue K.P., “Numerical simulations of nonlinear axisymmetricflows with a free surface”,Journal ofFluid Mechanics, vol.178, pp.195-219, 1987.Downie M.J., Bearman P.W. and Graham J.M.R., “Effect of vortex shedding on thecoupled roll response of bodies in waves”, Journal of Fluid Mechanics, vol.189, pp.243-264, 1988.Dawson C.W., “A Practical Computer Method for Solving Ship-Wave Problems”,Proceedings of the Second International Conference on Numerical Ship Hydrodynamics,Berkeley, pp. 30-38, 1977.Ekman, V.W., “On Dead Water: The Norwegian North Polar Expedition 1893-1896”,Christiana, vol.V, 1904.Faltinsen O.M., “A numerical nonlinear method of sloshing in tanks with two dimensionalflow”, Journal ofShip Research, vol.22, no.3, pp.193-202, 1978.Fink P.T. and Soh W.K., “Calculation of vortex sheets in unsteady flow and applicationsin ship hydrodynamics”, Proceedings of the 10th Symposium on Naval Hydrodynamics,Cambridge, Massachussetts,1974.Gadd G.E., “A Method for Computing the Flow and Surface Wave Pattern Around FullForms”, Transactions of the Royal Institute of Naval Architects, pp. 207-2 18, 1976.Graham J.M.R., “The forces on sharp-edged cylinders in oscillatory flow at low KeuleganCarpenter numbers”, Journal of Fluid Mechanics, vol.97, no.1, pp.33 1-346, 1980.Grilli S.T. and Svendsen l.A., “Corner problems and global accuracy in the boundaryelement solution of nonlinear wave flows”, Engineering Analysis with BoundaryElements, vol.7, no.4, pp.178-195, 1990.Grilli S.T., “Wave motion and overturning induced by moving bodies: Application toslender ship wave resistance”, Proceedings of the 1st International Conference onComputational Modelling of Free and Moving Boundary Problems, Southampton, UK,Computational Mechanics Publications, 1991.Hayashi K. and Kunishige Y., “Measurement of wave pattern around a yawed shipmodel”, Graduation thesis, Yokohama National University, Japan, 1988.Hess J.L. and Smith A.M.O., “Calculation of Potential Flow Around Arbitrary Bodies”,Progress in Aeronautical Sciences, vol. 8, pp. 1-138, 1967.Isaacson M. and Cheung K.F., “Second order wave diffraction around two-dimensionalbodies by time domain method”, Applied Ocean Research, vol.13, no.3, 1991.89Isaacson M. and Ng J.Y.T., “Time-domain second-order wave radiation in twodimensions”,Journal ofShip Research, vol.37, no.1, pp.25-33, 1993.Keulegan G.H. and Carpenter L.H., “Forces on cylinders and plates in an oscillating fluid”,J R National Bureau of Standards, vol. 60, no.5, 1958.Kim K. and Chol Y., “A numerical calculation of free surface potential flow field and shipwave resistance in shallow water by fully nonlinear wave theory”, Proceedings of theSecond Japan-Korea Joint Workshop on Ship and Marine Hydrodynamics, Japan, pp.111-120, 1993.Kudo K., “Hydrodynamic forces of the oscillating flat plate”, Proceedings of the 3rdInternational Conference on Numerical Ship Hydrodynamics, Paris, 1981.Leonard A., “Vortex methods for flow simulation”, Journal of Computational Physics,vol. 37, pp.289-335, 1980.Lin W.M., “Nonlinear motion of the free surface near a moving body”, PhD Thesis,Massachusetts Institute of Technology, 1984.Longuet-Higgins M.S. and Cokelet E.D., “The deformation of steep surface waves onwater, A Numerical Method for Computation”, Proceedings of the Royal Society,London, vol.A350, pp.1-26, 1976.Lugt H.J., “Numerical modelling of vortex flows in ship hydrodynamics : A Review”,Proceedings of the 3rd International Conference on Numerical Ship Hydrodynamics, Paris,pp.297-317, 1981.Ma H. and Tulin M.P., “Experimental Study of Ship Internal Waves: The SupersonicCase”, Proceedings of the 11th International Conference on Offshore Mechanics andArctic Engineering, OMAE 92, vol. 1A, pp.51-57, 1992.Maruo H., “New Approach to the theory of slender ships with forward velocity”, Bulletinof the Faculty of Engineering., Yokohama National University, Japan, vol. 31, 1982.Maruo H. and Song W., “Numerical appraisal of the new slender ship formulation insteady motion”, Proceedings of the 18th Symposium on Naval Hydrodynamics, Michigan,pp. 239-257, 1990.Maull D.J., “An introduction to the discrete vortex method”, Proceedings of InternationalUnion of Theoretical and Applied Mechanics, Karlsruhe, 1979.Nakayama T. and Washizu K., “The boundary element method applied to the analysis oftwo dimensional nonlinear sloshing problems”, Journal ofNumerical Methods in Eng, 17,1631-1646, 1981.90Ogilvie T.F., “The wave generated by a fine ship bow”, Report No.127, Department ofNaval Architecture and Marine Engineering, University of Michigan, 1972.Orlanski L., “A simple boundary condition for unbounded hyperbolic flows”, Journal ofComputational Physics, vol.21, pp.251-269, 1976.Pawlowski J.S., Baddour R.E. and Hookey NA., “The development of second generationdirect solution numerical models of large motions of floating bodies in waves”,Proceedings of the First Marine Dynamics Conference, pp.4.3-49, St John’s,Newfoundland, 1991.Sabuncu T., “Some Predictions of the Values of the Wave Resistance and Momentconcerning the Rankine solid under Interfacial Wave Conditions”, Norwegian Ship ModelExperiment Tank, Technical University of Norway, Publication No.65, 1962.Sarpkaya T., Computational methods with vortices - The 1988 Freeman Scholar Lecture”,Journal of Fluids Engineering, vol.111, pp.5-52, 1989.Song W., Ikehata M. and Suzuki K., “Computation of Wave Resistance and Ship WavePattern by the Slender Body Approximation”, Journal of the Kansai Society ofNavalArchitects, Japan, no.209, pp 25-36, 1988 (In Japanese).Song W. and Maruo H., “Bow impact and Deck Wetness: Simulation Based on NonlinearSlender Body Theory”, Proceedings of the 3rd International Offshore and Polar EngConference, Singapore, vol.111, pp.34-38, 1993.Suzuki K., “Boundary Integral Equation Method for the Linear Wave ResistanceProblem”, Proceedings of the 4th International Conference on Numerical ShipHydrodynamics, pp.308-323, 1985.Tahara Y., “A Boundary Element Method for Calculating Free Surface flows around aYawed Ship”, Private communication with Prof K Suzuki of Yokohama NationalUniversity, Japan, 1992.Tanaka Y. and Nakamura T., “Open boundary scheme for nonlinear irregular waves”,Proceedings of the 11th International Conference on Offshore Mechanics and ArcticEngineering, OMAE 92, vol.1A, pp.34.5-351, 1992.Tryggvason G., Unverdi S.O., Song M. and Abdollahi-Alibeik J., “Interaction of vorticeswith a free surface and density interfaces”, Lectures in Applied Mathematics, vol.28,pp.679-698, 1991.Tulin M.P. and Miloh T.,”Ship Internal Waves in a Shallow Thermocline: The SupersonicCase”, Proceedings of the 18th Symposium on Naval Hydrodynamics, pp.567-581, 1990.91Vinje T. and Brevig P., “Nonlinear ship motions”, Proceedings of the 3rd Internationalconference on Numerical Ship Hydrodynamics, Paris, pp.257-266, 1981.Watson G., Chapman R.D. and Apel J.R., “Measurements of the internal wave wake of aship in a highly stratified sea loch”, Journal of Geophysical Research, vol.97, no.C6,pp.9689-97013, 1992.Wong L.H., “A numerical model for vortex shedding from sharp wedges inoscillatory flow”, M.A.Sc. Thesis, Department of Mechanical Engineering, University ofBritish Columbia, Vancouver, 1990.Wong L.H., “Improved slender ship calculations with a wavemaker algorithm”,Proceedings of the 2nd Canadian Marine Dynamics Conference, Vancouver, pp. 65-74,1993.Wong L.H. and Calisal S.M., “A numerical solution for potential flows including theeffects of vortex shedding”, Journal of Offshore Mechanics and Arctic Engineering,Transactions ASME, vol.115, no.2, pp.111-115, 1993(a).Wong L.H. and Calisal S.M., “Waves generated by a ship travelling in stratified water”,Proceedings of the 3rd International Offshore and Polar Engineering, Singapore, pp.514-518, 1993(b).Yen S.M. and Hall D.R., “Implementation of open boundary conditions for nonlinear free-surface wave problems”, Proceedings of the 3rd International Conference on NumericalShip Hydrodynamics, pp.163-177, Paris, 1981.Yen S.M., Lee K.D. and Akai T.J., “Finite element and finite differences solution ofnonlinear free surface wave problems”, Proceedings of the 2nd International Conferenceon Numerical Ship Hydrodynamics, pp.305-318, Berkeley, 1977.92Appendix ASlender body theory and the wavemakerThe following analysis is given to show how the three dimensional problem of a slenderbody advancing at a velocity U is approximated by a series of two dimensional problems.The total potential is given by:P= Uxcosa÷Uysinct+4)(x,y,z) (A.1)where the free stream velocity U is incident on the body at an angle a and 4) is theperturbation potential. Relevant variables are non-dimensionalized (the non-dimensionalform being denoted by a (-), such that:4) — UB4)= BjNotationLength LBeam B[au9ht dFree streamvelocity (UAPNote : Origin is at undisturbedfree surface level.Figure A.1 : Frame of reference and notation.93Y=BYx=Ly = Bz = B (A.2)where r is the wave elevation and Y describes the body boundary. The ship beam andlength are B and L and the hull offset is defined by y = Y(x,z). The spatial coordinates x,y, z are defined in Figure A.i. Since the boundary conditions require information aboutspatial derivatives as well, it is necessary to express these as:a_i a— L aa_i a— B a3(A.3)ôz BazThe fluid is taken to be inviscid and incompressible and flow is assumedirrotational. With these assumptions, the governing Laplace equation is:(A.4)ox2 0y2 0z2where is the velocity potential. From Equations (A.2) and (A.3) and assuming thatderivatives with respect to non-dimensional variables do not change in order of magnitude,the Laplace equation can be written:+ + =0 (A.5)94If the slenderness ratio B / L = e is small (limited to about 0.15), the 3D Laplace equationcan be expressed in the two dimensional form in the plane of the hull sections:+ = 0 (A.6)The 3D impermeable hull boundary is:‘1),n +(1).,n +(I)n =0 (A.7)In non-dimensional form this is rewritten as:—cosaY —-4Y ±(sinct + = 0 (A.8)from which the term of order c2 is eliminated to give:—ecosaY ±(sina+p)+4Y =0 (A.9)The kinematic free surface condition can be written as (in non-dimensional form and theneliminating the term in ):(e cosa+e2 + (na + — = 0(Ecosa)j +(sinct+)j —p = 0 (A.10)The non-dimensional form of the dynamic free surface condition for steady potential flow,eliminating the terms in 2, is given in the development of Equation (A. 11) below.--[i +i i-t_u2]+g9=o+2Ucosa1+4 +2Usina +4j+g=o€cosa +sina ÷‘[E22 +2j+j = 095e cos + sin + ![ + j + =0 (A. ii)To maintain the last term in the dynamic free surface condition, it is necessary to haveO(gB I U) 0(e). Here, the draught (d) at which the body floats is taken to be of theorder of the beam (B). Therefore in this application, the draught Froude number(Fd = U/.J) is taken to beAn expression can be obtained, using the slender body approximation, for thenormal velocity on the body. Equation (A.9) is written in its dimensional form as follows:—UcosaY ±(Usina+p)+Y =0 (A.12)where the normal vectors on the body are given by:—Y +1Y(n,n,n)=Z’ Z (A.13).J1 + + Y2where Y = 0(e) and Y = 0(1). Elimination of the term Y3, the normal velocity(positive pointing out) of the body in the cross-flow plane can be written, using Equations(A.13) and (A.14):+= u(Y cosa ± sin a)on the hull (A.14)an .iJl+Yz2Putting Equations (A.5), (A.10) and (A.11) back into dimensional form and rearranging,the following can be obtained:+ 4 = 0 in the domain A.15)Ucosat +(usina+p)q —• =0 at (A.16)96u(p1 cosa+p na)÷--[p +p]+g’j =0 at z= i (A.17)This set of equations (A.14) to (A.17), together with impermeable bottom and wallboundary conditions, are identical to those used for Ihe numerical modeling of a twodimensional wave tank in the time domain. The only difference here is that the motion ofthe wave paddle is replaced by the geometric change of the ship hull in the longitudinaldirection.97Appendix BHigher order boundary element methodThe following description on the boundary element method applies to two dimensionalpotential flow problems only. We start with Green’s second formula:—GdS=O (B.1)S an onwhere the weighting source function is taken to be:G = --ln( (B.2)2t ‘sr!Both G and i satisfy the governing Laplace equation. The boundary can bediscretized into N number of elements. Collocation points are placed at the ends of linearelements, where the variation of the potential and its normal derivative over the elementsare assumed linear. If the potential and its normal derivative are taken to be varyingquadratically over an element, then the collocation points (nodes) are placed at both theends and at the middle of an element. Substituting the normal derivative of G with H, thepotential at a nodal point on the boundary can be thus be written:=—pHdf (B.3)Equation (13.3) can be rewritten by inserting the term in into the summation sign suchthat:98fdljHujdE’ = (B.4)Note that in Equation (B.4), the term H includes the factor ot 1/2 from Equation (B.3).In the expression for G, the variable r represents the distance between the node i and agiven point on the element j. Dirichiet or Neumann boundary conditions are specified ateach node. All unknowns in Equation (B.4) can then be gathered such that the system ofN equations in N unknowns to be solved for X can be written as (for more details seeBrebbia and Domnguez, 1989):[AJX]=[Bj (B.5)The evaluation of the matrices [A] and [B] involves numerical integration for the Gand H terms in Equation (B 4). For collocation nodes located between two elements, thecontributions to the integrals in G and H come from the two adjacent elements.=——f xiln()dr+----f X2l11()dF2t’ r rH..=x1—-Iln(-ldF + X24ln(ldF(B.6)“ 2t ‘ on L \rJJ 231 F On L \rJJEquations (B.7) and (B.8) for the value of apply to linear and quadratic elementsrespectively.Xi(B.7)X2(11)Xi2 (B.8)X299In the above, each element is parameterized into = —1 to = 1 between the end nodes.In quadratic elements, the mid-element node G and H integrals contain only one termwhere:XL =(‘—X’÷) (B.9)X2 = 0The integrals in Equation (B.6) are evaluated using standard Gaussian andlogarithmic numerical integration quadratures (see for example, Brebbia and Dorninguez,1989).While potential values are taken to be continuous at the collocation nodes, normalvelocities need not be so. This is evident, for example, at the intersection between the freesurface and a body or wall. The normal velocities immediately on the two sides of suchnodes are taken to be different, requiring a split of the integral G at these ‘double nodes’.Thus, the integral is evaluated such that the adjacent free surface element contributes onlyto G for the free surface normal velocity and the adjacent body element contributes only toG for the body side normal velocity.Finally, the procedure for the extrapolation of the normal velocity on the freesurface side of the fluid-body intersection point requires some explanation. Linearvariation of the normal velocity over the three free surface nodes (see Figure B. 1) adjacentto the fluid-body intersection point is assumed. Hence the normal velocity at theintersection point is given by:(4 =2( — (B.10)\iln!A ‘\ônJB \Jc1004,nFigure B.1 : Intersection point and normal velocities.The above eliminates the need to evaluate the normal velocity for the free surface side ofthe node at the intersection point. There, the velocity potential is then taken to be theunknown. Another consequence of the extrapolation is that G integrals at points B and Care modified follows:018 = B +(B.11)Free surface e’ementsNode ANode BNode C101
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Slender ship procedures that include the effects of...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Slender ship procedures that include the effects of yaw, vortex shedding and density stratification Wong, Haw L. 1994
pdf
Page Metadata
Item Metadata
Title | Slender ship procedures that include the effects of yaw, vortex shedding and density stratification |
Creator |
Wong, Haw L. |
Date Issued | 1994 |
Description | The accurate determination of hydrodynamic loads on moving ships is important for hull form design and optimization and structural design purposes. This is especially true at the preliminary design stage during which time quick predictions of the forces and moments acting on a ship advancing steadily with, and without, yaw would be extremely useful. In view of this, simple numerical cross-flow algorithms has been developed. The numerical procedures are based on slender body theory, which is used to convert the three dimensional problem into a series of two dimensional wavemaker problems in the plane of transverse sections, marching in small steps from the bow section towards the stem. Fluid density stratification, vortex shedding, finite water depth and nonlinear free surface effects can be allowed for in the algorithms. A procedure for handling density stratified flow has been developed and successfully used for the calculation of surface and interfacial waves created by a prolate spheroid. Vortex shedding is modelled using the discrete vortex method. A hybridization of the discrete vortex and boundary element methods is achieved and illustrated in a test case of predicting the forces acting on an oscillating flat plate. The wavemaker, with the fully nonlinear free surface conditions, is used for calculating the generated wave pattern and wavemaking resistance of a Wigley hull. The effects of finite water depth on wavemaking resistance are calculated. The hybrid boundary element-discrete vortex method is used for determining the hydrodynamic forces and moments acting on a yawed Wigley hull. |
Extent | 2374742 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-04-15 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0080831 |
URI | http://hdl.handle.net/2429/7135 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mechanical Engineering |
Affiliation |
Applied Science, Faculty of Mechanical Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 1994-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_1994-954108.pdf [ 2.26MB ]
- Metadata
- JSON: 831-1.0080831.json
- JSON-LD: 831-1.0080831-ld.json
- RDF/XML (Pretty): 831-1.0080831-rdf.xml
- RDF/JSON: 831-1.0080831-rdf.json
- Turtle: 831-1.0080831-turtle.txt
- N-Triples: 831-1.0080831-rdf-ntriples.txt
- Original Record: 831-1.0080831-source.json
- Full Text
- 831-1.0080831-fulltext.txt
- Citation
- 831-1.0080831.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0080831/manifest