Slender Ship Procedures That Include the Effects of Yaw, Vortex Shedding and Density Stratification By flaw LWong B.Sc. (Hons), Naval Architecture/Shipbuilding, University of Newcastle Upon Tyne, 1985 M.A.Sc., Mechanical Engineering, University of British Columbia, 1990 A THESIS SUB MITFED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES DEPARTMENT OF MECHANICAL ENGINEERING We accept this thesis as conforming to the required standard October 1994 © Haw L Wong, 1994 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. (Signature) Department of QA The University of British Columbia Vancouver, Canada Date DE-6 (2/88) 0 c%r 99 ‘ . Abstract 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. hi 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. II Table of Contents Abstract Table of Contents iii List of Figures v Acknowledgement viii INTRODUCTION 1 Chapter 1 : Review of Relevant Literature 5 1.1 Boundary element method and free surface flows 5 1.2 Potential flow around ships 8 1.3 Interfacial waves generated by a moving body 10 1.4 Discrete vortex methods 10 1.5 Vortex-free surface-body interactions 12 14 Chapter 2: Formulation and Calibration 2.1 The wavçmaker problem 14 2.2 Free surface evolution for higher order elements 18 2.3 Bow wave generated by a wall-sided wedge 21 Chapter 3: Slender Body in Two-Layer Medium 26 3.1 Boundary conditions 26 3.2 Numerical procedure 29 111 3.3 Numerical results . Chapter 4: Hybrid Boundary Element/Discrete Vortex Method 31 39 4.1 Discrete vortex method 40 4.2 Vortex shedding in boundary element method 44 4.3 Forces on sharp edged cylinders 47 Chapter 5 : Applications in Nonlinear Slender Body Ilydrodynaniics 52 5.1 Wavemakiiig resistance of a slender body 52 5.2 Wavemaking resistance in shallow water 65 5.3 Hybrid method in yawed ship calculations 67 Chapter 6: Conclusions and Recommendations 81 6.1 Conclusions 81 6.1.1 Slender body formulation and implementation 81 6.1.2 Ships in density stratified water 82 6.1.3 Hybrid discrete vortex-boundary element method 82 84 6.2 Recommendations Nomenclature 86 Bibliography 88 Appendix A: Slender body theory and the wavemaker 93 Appendix B: Higher order boundary element method 98 iv List of Figures Figure 2.1: Coordinate system fixed on the body 15 Figure 2.2: Boundary definition for wavemaker problem 15 Figure 2.3 : Infinite wedge definition sketch 22 Figure 2.4: Calculated and measured bow wave amplitude for wedge, Figure 2.5 : Calculated and measured wave peak location for wedge, I = = 150 150 24 25 Figure 3.1 : Prolate spheroid and dimensions 27 Figure 3.2 : Definition of the computational domains 28 Figure 3.3 : Overview of the numerical procedure 30 Figure 3.4: Free surface contours for spheroid at F 1 = 3.0 Figure 3.5 : Interfacial displacement for prolate spheroid at F 1 34 = 3.0 35 Figure 3.6: Interface contour plot for prolate spheroid at F = 3.0 36 Figure 3.7: Location of interfacial crest lines, F 1 = 2.95 37 Figure 3.8 : Location of interfacial crest lines, F 37 Figure 3.9 : Decay of first interfacial crest, F = = 3.00 2.95 37 Figure 3.10: Influence of upper layer thickness on first interfacial crest 38 Figure 4.1: The physical Z-plane and the transformed .-plane 40 1.00 for a flat plate 43 Figure 4.3 : Exclusion of a vortex singularity in domain D 44 Figure 4.4: CQmpptational domain and boundary conditions 47 Figure 4.5 : Forces on an oscillating flat plate, K 49 Figure 4.2: Vortex sheet rollup at t = V 6.6 Figure 4.6: Drag coefficients for the flat plate 50 Figure 4.7: Inertia coefficients for the flat plate 51 Figure 5.1: Wigley hull dimensions and coordinate system 53 Figure 5.2 : Wave on side of Wigley hull at FR = 0.267 57 Figure 5.3 : Wave on side of Wigley hull at FR = 0.316 58 Figure 5.4: Contour plot of wave field for Wigley hull at FR 0.267 = 59 Figure 5.5 : Measured and calculated wave field (YNU) for Wigley Hull, a Figure 5.6 : Free surface displacement for Wigley hull at FR Figure 5.7 : Pressure distribution on Wigley hull at FR = = 00 0.267 = 60 61 0.267 61 Figure 5.8 : Wavemaking resistance against YNU results 62 Figure 5.9 : Wavemaking resistance against Aanesland (1989) 63 Figure 5.10 : Computed and measured sinkage values for Wigley hull 64 Figure 5.11 : Wavemaking resistance for Wigley hull in shallow water at FR Figure 5.12: Squat and trim for Wigley hull in shallow water at FR = = 0.3 16 Figure 5.15 : Wave field for Wigley hull at FR = = 0.267, a 0.267, a = = 70 100 73 = 100 Figure 5.17: Free surface displacement for Wigley hull at FR = 0.267, a = 72 100 Figure 5.16: Measured wave field (YNU) for Wigley Hull, a Figure 5.18: Pressure distribution for Wigley hull at FR 66 66 Figure 5.13 : Vortex induced normal velocity on a boundary element Figure 5.14 : Wave profile on sides of Wigley hull at FR 0.3 16 0.267, a 74 = 100 = 10° 75 76 Figure 5.19 : Forces and moments for yawed Wigley hull 77 Figure 5.20: Development of total vortex strength over the hull length 78 vi Figure 5.21 : Keel vortex evolution for Wigley hull, a = 100 (forward) 79 Figure 5.22: Keel vortex evolution for Wigley hull, a = 10° (aft) 80 Figure A.1 : Frame of reference and notation 91 Figure B.1 : Intersection point and normal velocities 99 vu Acknowledgement This thesis is dedicated to my family, from every member of which I have drawn much support and encouragement. In particular, my wife Julie deserves much praise for her assistance and patience. My children Justin and Sarah are a great source of joy and motivation. The work reported in this thesis was done under the supervision of Professor Sander Causal, who has over the years offered much help and guidance. I am also grateful to Professor Kazuo Suzuki (Yokohama National University, Japan) for providing experimental data and material on slender body theory. Professor Orner Goren’s (Istanbul Technical University, Turkey) kind generosity in providing the Dawson programs is much appreciated. I should also mention the great fellowship I enjoyed with all my past and present colleagues. To all these people and others who have in one way or another made my work easier, I offer my humble thanks. The generous scholarship in computer time on the Fujitsu VPX24O/1O supercomputer granted to me by the High Performance Computing Centre of Calgary is sincerely appreciated. Financial support for this work was provided by the Natural Sciences and Engineering Research Council of Canada. Viii INTRODUCTION Robust and reliable methods for determining the hydrodynamic characteristics of a vessel or offshore structure are valuable tools for design engineers and naval architects. The ability to fine tune designs at an early stage means an improvement to the design process and consequently better designs. Apart from that, much cost and time savings accrue from the reduction of required model testing and the avoidance of remedial work after commissioning. Ongoing research over the years, aimed at the development of better and more efficient numerical tools of practical interest to ship designers, has resulted in numerous innovative approaches. Despite this, numerical solutions relating to nonlinear ship wave generation have proved elusive. A number of computational methods, mostly linearized, were introduced in the 70s. However, as the importance of nonlinearity and viscous effects become better understood, much work in the last decade or so has attempted to account for both. The objective of the present work is to use potential flow theory to study the hydrodynamics of a slender body in calm water, including the effects of nonlinearity, yaw, vortex shedding and fluid density stratification. Numerical implementation was achieved using the higher order boundary element and discrete vortex methods. Slender body approximations simplified the threc dimensional moving ship problem into a sequence of two dimensional wavemaker problems. Working in the cross-flow plane offers many advantages 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 three categories, which will be described in more detail in the following sections. Nonlinear Slender Body Procedure The slender body formulation has been generalized to accomodute oblique motion of a ship. This was implemented through a wavemaker approach in which the 3D ship wave problem 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 the study of bow wave making for an advancing wall-sided wedge. The algorithm was then extended to include the fully nonlinear free surface conditions and applied to the calculation of deep and shallow water wavemaking resistance and the waves generated by an advancing Wigley hull. Vortex shedding effects were not included in these calculations. Vorticity generated by the round bilges were assumed unimportant in forward motion problems and hence not accounted for. This assumption can be justified considering the success of pure potential flow ship resistance calculation methods. The nonlinear free surface evolution was achieved using the Eulerian-Lagrangian approach. The procedure, when extended to general hull forms can be an invaluable tool in hull form optimization applications. Suitable ship types for such applications are naval vessels, container ships and multi-hull vessels. Ships in Stratified Waters Density stratification is a common phenomenon occuring both in the oceans and inland waters. The variation of the density with water depth is due to changes in the salinity 2 and/or temperature over a region called a pycnocline. Although the density differences may be small, ships moving slowly in the vicinity of such interfaces experience what came to be known as deadwater effects. When the ship speed approaches that of the fastest generated internal waves, the increase in drag over the hull is very large and the available power 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 surface and the disturbance of the flow around the ship due to these internal waves. It is also known that a vessel loses controllability in such situations. In naval applications, fully or partially submerged vessels could produce undesirable internal wave signatures. The procedure given in this thesis is applicable to a slender hull moving in a two layer 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 fluid domains are used, the lower domain with a density slightly higher than the upper for overall gravitational stability. The classical treatment of two dimensional internal waves provides a condition relating to the interface between the two layers. The boundary element method was used to obtain the transient solution by solving the lower, then the upper domain sequentially at each time step. This algorithm is expected to be useful for resistance prediction of submerged bodies and for the study of internal wave behaviour without the presence of a body. Vortex Shedding Effects There are numerous instances in the study of fluid dynamics where vortex shedding forces are significant and should be taken into account. A good example is in the tumerical 3 prediction of the roll response of a ship with sharp corners and bilge keels where wave damping is light and vortex effects are important. However, vortex effects are usually not included in potential flow applications to ship dynamics. This is partly due to the complex nature of the fluid-structure interaction which generates vortices from a hull. There is also the difficulty of implementing vorticity calculations in wave problems without resorting to viscous flow computations involving the Navier-Stokes equations. A hybridization of the boundary element and discrete vortex methods has been achieved. This hybrid method allow a simple but effective way to account for the effects of vortex shedding in potential flow calculations. The vortex generation mechanism is taken to be the keel. Bow and secondary vortices were not included since the complex nature of the flow in both cases cannot be properly modelled using a simple discrete vortex approach. The computational advantage in being able to account for vortex effects while still working with boundary values is significant. The hybrid method was implemented for the case of a Wigley hull in oblique motion. Changes in wavemaking resistance and yaw forces and moments due to small angles of attack, whether intentional or otherwise, can be studied. Such studies are of practical interest from an operational standpoint. In addition, forces and moments in yawed ship motion can be used to determine the so called ‘stability derivatives’, important in the study of ship manouevring. This in turn allows systematic investigations into the directional stability of a vessel at the design stage. 4 Chapter One Review of Relevant Literature 1.1 Boundary element method and free surface flows The boundary element method is a much used method for solutions to fluid dynamics problems. An advantage of this method is that for potential flow calculations only conditions on the boundary need to be specified and calculations are carried out along the discretized boundary only. This method is well suited for application in boundary value potential flow simulations involving relatively complex and possibly time varying geometries. Much attention has been focused on finding stable numerical solutions to the nonlinear problem. The satisfaction of the nonlinear free surface conditions at a location which is unknown a priori has been a major stumbling block to the achievement of stable numerical solutions in nonlinear wave making. Longuet-Higgins and Cokelet (1976) used a mixed Eulerian-Langrangian form of the free surface conditions to trace the temporal development of high amplitude surface waves. They assumed spatial periodicity, mapping the domain such that only free surface variables need to be handled in the solution procedure. Their formulation was successful in coping with the multi-valued characteristics of an overturning wave. 5 Faltinsen (1978) applied the boundary element method to the two dimensional problem of liquid sloshing in an oscillating container in which the nonlinear free surface boundary conditions are treated for consideration of a nodal point constrained to move in the vertical direction only. Nakayama and Washizu (1981) also reported work on the above problem. Although not clearly stated, it appears that this work used the kinematic free surface boundary condition in its linearized form. Cointe (1989), using the mixed Eulerian-Lagrangian method, studied wave diffraction from a submerged cylinder. Pawlowski, Baddour and Hookey (1991) discussed the use of linear elements for nonlinear wave simulation. Isaacson and Ng (1993) reported a time domain boundary element scheme for nonlinear wave radiation from floating cylinders in sinusoidal sway, heave and roll motion. In the numerical simulation of nonlinear interactions involving a free surface and a floating body, boundary conditions are to be satisfied at the instantaneous location of the boundaries at any given time step. This requirement, when applied to the fluid-body intersection point, results in difficulties arising from the singular behaviour of the fluid velocities at that point. Consequently, the definition of conditions such as the velocity potential, the direction and magnitude of the velocities, and the location of the free surface node at the fluid-body intersection point has proved to be a challenging task for all research into nonlinear wave simulation. Yet, the resolution of this problem is crucial for the success of any such simulation, both in terms numerical stability and accuracy. Vmje and Brevig (1981) used the body boundary condition at the intersection point. The velocity potential is taken to be the unknown there and emerges as part of the 6 solution. Un (1984) specified both the velocity potential and stream function value at the intersection point in an algorithm using the Cauchy integral formulation for boundary elements. Dommermuth and Yue (1987) used a similar technique and a boundary regridding scheme. The regridding they used replaced the Lagrangian points to give equal length elements on the free surface. The use of double nodes has recently emerged as a popular idea in the fluid-body intersection problem. Grihi and Svendsend (1990), among others, used this idea, enforcing compatibility only at the end of each time step. Another major difficulty in numerical calculations of nonlinear wave making is the application of the far field condition stemming from the truncation of infinite domains to manageable dimensions. The lack of proper physical conditions governing the passage of energy out of the computational domain has made the radiation condition for full nonlinear problems difficult to define. The eventual reflection of waves back into the computational domain causes errors in the solution for the region of interest. This undesirable situation is common to all numerical methods. To cope with this, a variety of ways has been devised which define conditions at a given cut-off boundary which is not excessively distant from the source of the disturbance. Obviously, limiting the size of the computational domain improves computational efficiency. Chan (1975) used a mathematical wave damping device to remove energy of the outgoing waves. Orlanski (1976), on the other hand, used the Sommerfeld radiation condition, proposing a finite difference scheme to calculate a depth and time dependent celerity, which should always be directed out of the domain. The method does not rely on knowledge of the dispersion characteristics of the wave train and as a result may be used 7 for nonlinear waves. Isaacson and Cheung (1991) obtained satisfactory results for second order wave propagation. Yen, Lee and Akai (1977) applied the undisturbed condition at the cut-off boundary which is periodically expanded downstream as the propagating wave approaches it. Apart from the fact that computational efficiency deteriorates with time as the domain grows, the solutions reportedly became unstable after some time. In a later paper, Yen and Hall (1981) used Orlanski’s method for nonlinear waves together with a proposal to remove the higher frequency error due to the reflected waves having wavelengths of about two grid intervals. This was done using two separate methods: filtering over the whole free surface using a mathematical smoothing function and a modification of the formulation to include dissipation and dispersion. Dommermuth and Yue (1987) matched the inner nonlinear region by a linear outer region and reported this gave a robust algorithm. More recently, Tanaka and Nakamura (1992) proposed an extension to the Orlanski scheme for application to unsteady and irregular nonlinear waves. Their method used the nodal point nearest to the free surface-radiation boundary intersection for the prediction of the time dependent phase velocity of the outgoing waves. 1.2 Potential flow around ships Due to the difficulty of obtaining complete analytical solutions to the ship wavemaking problem, a tremendous amount of work has been done on the development of numerical methods. Hess and Smith (1967) developed a source panel method for the calculation of potential flows around bodies of arbitrary shape. Dawson (1977) devised a method in 8 which a combination of the double body Rankine source panel method and marching scheme in the stream-wise direction produced useful information relating to. the hydrodynamics of steady ship motion. Cheng (1989) extended the method for 3D transom stem flows while Tahara (1992) included lifting effects in the double body potential. Gadd’s (1976) iterative numerical procedure for ship resistance prediction was followed up by Causal et al. (1991) for nonlinear wave calculations. The so-called Neumann-Kelvin approximation is also popular and well reported in the literature, see t’or example the work of Suzuki (1985). Maruo approximation. (1982) presented a formulation based on the Neumann-Kelvin Slender body approximations result in a simplification with the distribution of Kelvin sources on the hull surface only. A marching procedure starting at the bow can be used, with the source strengths at a given section being determined by the ship normal velocity at that section and a known quantity involving data evaluated upstream. Upon completion of the procedure, all the source strengths on the hull surface are determined. Quantities of interest are then computed by evaluating the effects all these sources. Song, Ikehata and Suzuki (1988) applied the method to the Wigley hull, reporting good agreement with measured values. Calisal and Chan (1989) reported on the use of a simple wavemaker method and slender body theory for the calculation of bow waves produced by an advancing wall-sided wedge. Wong (1993) pointed out that using starting 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 and Maruo (1993) described a wavemaker approach to bow impact and deck wetness studies. 9 Kim and Choi (1993) calculated shallow water effects on hydrodynarnic forces acting on various hull forms using a three dimensional iterative procedure involving the Rankine source panel method. 1.3 Interfacial waves generated by a moving body One of the earliest works on ship excited interfacial waves in stratified media was reported by Ekman (1904), who was able to experimentally confirm that internal waves created by a moving ship cause a large drag on the hull. Sabuncu (1962) calculated the Michell wave resistance in the presence of interfacial waves in a two layer fluid. More recent interest in the problem was fuelled by attempts to explain the long narrow V-wake left behind a ship moving in stratified seas as sensed by radar. Tulin and Miloh (1990) used linearized theory to obtain a solution in terms of a Green’s function. Wong and Causal (1993b) used the cross flow approach to predict the interfacial waves excited by a prolate spheroid and compared the results to the experimentally measured values of Ma and Tulin (1992). Watson et a!. (1992) reported on full scale experimental work with ships traversing a highly stratified Scottish sea loch. They concluded that interfacial wave amplitudes are only weakly dependent on ship speed but are more sensitive to changes in fluid stratification. 1.4 Discrete vortex methods Flow separation occurs when the boundary layer on a body surface reaches a sharp edge where the radius of curvature of the edge is very much smaller than the boundary layer thickness. In two dimensional flow, the adverse pressure gradient set up results in the 10 formation of an unstable shear layer which subsequently rolls up into a tight spiral and is then shed into the flow field as a free vortex sheet. Such a vortex sheet can be approximated by the introduction of discrete vortices at given time intervals as the flow develops. A persistent and yet unresolved problem in such schemes is the irregular roll-up of the vortex sheet as well as numerical instability due to the uncharacteristically high velocities induced when a vortex is very near to another vortex or its image. These problems 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 high induced velocities on each. Chorin (1973) used vortices with a viscous core to model viscous diffusion as well as to stabilize his numerical procedure. As a result, the maximum induced velocity of a vortex is finite at a given distance from the vortex position and decreases 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 of vortices that are too close together. Fink and Soh (1974) pioneered a scheme of rediscretization in which the vortex sheet is rearranged into equidistant positions after each time step in the numerical procedure. This same regridding scheme was adopted by Dommermuth and Yue (1987) for the treatment of the nonlinear free surface to improve numerical stability. The above, and many other similar, schemes help in one way or another to give smooth vortex roll-up as well as extend the computational time span in which calculations remain stable, Wong (1990). H Graham (1980) applied the discrete vortex method to calculate the vortex forces induced at a sharp edge in oscillatory motion at low Keulegan-Carpenter numbers. This was done by regarding vortex shedding from an infinite wedge as the inner region of flow past a large but finite body. The underlying assumption in this case is that the body length scale is large so that the vortex shed does not affect other parts of the body which are far away from the edge when compared to the flow length scale. Downie, Bearman and Graham (1988) followed up on this and calculated the vortex damping forces on a rolling rectangular barge. Wong and Causal (1993a) applied a hybrid discrete vortex-boundary element method to the determination of the forces acting on a submerged, oscillating flat plate. The results from that study compared well to measured values. Comprehensive general reviews on various vortex methods can be found in Leonard (1980), Lugt (1981) and Sarpkaya (1989). 1.5 Vortex-free surface-body interactions The study of interactions between vortices and free surfaces, density interfaces and moving bodies has been motivated by attempts to improve potential flow results for situations where vortex effects were thought to be important. In particular, body generated vorticity can result in significant influences on the hydrodytiamic forces acting on the body. Bradbury (1986) generated a number of rectangular transverse sections to describe a yawed ship-like slender body for the study of hydrodynarnic effect of bilge vortices. It was necessary, in that case, to use rectangular blocks where the submerged sharp corners 12 representing the bilges provide well-defined separation points for the vortex shedding algorithm employed. Time marching was used to advance from one cross-flow plane to the next. At each section, conformal mapping was used as part of an iterative procedure for the solution. The free surface was assumed rigid for his work. Braathen and Faltinsen (1988) applied a time stepping procedure for two dimensional 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 a sharp density interface and free surfaces. 13 Chapter Two Formulation and Calibration 2.1 The wavemaker problem This section summarizes the development of a ‘wavemaker’ algorithm for the computation of flow around steadily advancing slender ships. Simplification of the three dimensional problem is achieved using the slender body approach given in Appendix A. This reduces the problem into a series of two dimensional ones in the cross-flow plane, solved sequentially by marching in time from the bow in the direction of the stern. An effective way to visualize this is to view the ship from a fixed yz planc ahead of the bow. As the ship advances past this viewplane, different transverse sections are ‘cut’ by this plane such that the fixed observer sees the body sections ‘opening up’ until around the midship region and ‘closing in’ again towards the stern. The frame of reference has the origin on the ship centreline at amidships and on the undisturbed free surface as in Figure 2.1. The angle of incidence of the free stream velocity, or the yaw angle, is taken to be a. Details concerning the use of slender body theory to derive a set of equations describing the wavemaker problem is given in Appendix A and will not be repeated here. 14 Notation Length L BeamB Draughtd z AP ZFP Free stream velocity (U) Note: Origin is at undisturbed free surface level. Figure 2.1: Coordinate system fixed on the body. Ideal fluid assumptions, together with slender body approximations, result in the two dimensional governing I.aplace equation from Equation (A.15): Pyy + +72 0 (2.1) z Free suace (Sf) Free suace (S / Wall (Sw) / / U(YCOSU ± sina)/(1 -t-Y) 5 / / Wall / (Sw) / / IGoveming Equation : = / // / / // / / / Bottom boundary (S ) 0 Figure 2.2: Boundary definition for wavemaker problem. 15 All normal velocities were taken to be positive pointing out of the domain (see Figure 2.2). The bottom boundary S 0 is taken to be impermeable. A ‘mirror image’ technique using symmetry to satisfy the bottom condition is often used to reduce the number 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 from the intersection points between the bottom boundary (S ) and the sides (Sw). However, 0 the present work used linear and quadratic elements, that is collocation points are found at all boundary intersections. The additional computational book-keeping resulting from the use of this technique erodes the savings in computational effort. Moreover, by keeping the bottom boundary, the formulation can be extended to the more general case of uneven and sloping seabeds. Thus, the bottom boundary was retained for all the calculations described 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. There were two reasons for not imposing the radiation condition at this boundary. Firstly computed results are to be compared to values measured in tank testing where flow is always bounded. The other reason is that it is extremely difficult to implement the numerical radiation condition for the mixed Lagrangian-Eulerian approach adopted in nonlinear cases in the present work. On the body, the normal velocity in the cross-flow plane is due to the rate of change of the sections as the marching proceeds. The ‘wavemaker’ displacement is thus 16 represented by the changing hull sectional shapes along the ship length. Normal velocities on the hull section (negative into the fluid) at each step is defined using Equation (A.14): a4 an -U(Ycosct±sincL) Ji+Y At the free surface, S, the dynamic and kinematic conditions, to be satisfied at the instantaneous free surface location, are applicable. Using the conversion a a I —= UI cosa—+sinct— ôt ox OyJ ‘ (2.3) . assuming that U is uniform everywhere in the flow field. From Equations (A.16) and (A.17): OyOy at + z + ()] + =0 atz=rl (2.4) at z (2.5) = The marching procedure from one cross-flow plane to the next can be viewed as a time domain problem where the time step is the amount of time taken to traverse dx, the spatial step along the longitudinal axis. The hydrodynamic forces and moments are evaluated according to the following integrals: L12 = fcixf p.ndS (2.6) —L/2 L/2 M = f dxf p(rxn)dS (2.7) -L12 17 where n is the normal vector at the middle of the element dS and r is taken with reference to the origin (0,0,0). The integration is carried out over the girth of each considered cross-section (SB) and, after completion of the sweep, over the ship length using Simpson’s rule. The pressure p is obtained from the Bernoulli equation: p = 0 P—P =— U(cosa+sinu)+_((.) ox 2Oy ()2 ÷gz ) (2.8) 2.2 Free surface evolution for higher order elements For the solution of the cross-flow plane boundary value problem, velocity potentials are specified at all free surface nodes. The free surface conditions are satisfied at its exact location for nonlinear calculations. Where linearized conditions were used, the free surface boundary was fixed at the mean level. In the nonlinear procedure, adopting the Lagrangian-Eulerian approach to trace free surface particles, the two dimensional dynamic free surface condition in Equation (2.5) can be rewritten as: D1[(O4) ( 2 O)2] atz=11 (2.9) using the substantial derivative: Dt at OyOy (2.10) OzOz 18 The velocities on the free surface can be described by: Di Dt a Dy Dt a4 ay (2.11) ôz (2.12) For the linear solution procedure, the kinematic and dynamic free surface conditions are written as: = öt -gq ôz at z on 0 (2.13) atz=O (2.14) = The solution for the time marching procedures require the use of an extrapolator to find the potential values on the free surface and its location at the next time step. In the present work, the first order Adams-Bashforth predictor for time step size ôt was used. (t + ôt) = + ôt) = (t) + [3(t) q(t) + [3 (t) — — (t 1 ‘q (t — t)]- (2.15) ôt)}- (2.16) (4)) and the location of the free surface (‘q) at — In this way, values of the velocity potential time (t + ôt) are predicted using Equations (2.15) and (2.16). With Neumann conditions specified on all the other boundaries, the problem is then solved using the higher order boundary method (Appendix B). The normal velocities on the free surface nodes emerge as part of the solution from the boundary element method. From these, time rates of change for the potential and the free surface location can be determined using Equations 19 (2.9), (2.11) and (2.12). For linearized conditions, Equations (2.13) and (2.14) are used instead. 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 to the next time step. p(t + ôt) = p(t) + [p, (t) + p((t + ôt)}. (2.17) rl(t ôt) = t(t)+ [r(t)+ q(t (2.18) + + ôt)] The velocities in the y- and z-directions(p and p) can be determined using the normal and tangential velocities at each collocation node. =cosO—sinO as (2.19) -sinO +cosO as an (2.20) ay - az = where 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 linearly over each quadratic element and being constant over a linear element. Thus, for quadratic elements, the tangential velocity is given by: as 3P — 4 P k÷1 AS + Pk+2 (2.21) where AS is the element length, and k is a nodal index increasing in the anti-clockwise direction. 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 the 20 normal derivatives immediately to either side of the fluid-body intersection node are different. 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 the body at the intersection point. The normal velocity on the free surface side is linearly extrapolated from the nearest two free surface nodes. This leaves the potential as the only unknown at the intersection point. Details on how this is numerically implemented is given in Appendix B. 2.3 Bow wave generated by a wall-sided wedge The wavemaker algorithm described earlier was used for the prediction of bow waves generated by an infinite (after body) wedge. This simple application provided an opportunity for the calibration of the time stepping procedure. As the midship section is undefined, the frame of reference has to be shifted to the forward perpendicular, see Figure 2.3. The computational boundary, as shown in Figure 2.3, is similar to that given in Figure 2.2 except for the presence of the bottom of the wedge section, on which the normal velocity is zero (ignoring dynamic trim and sinkage effects) and the use of a symmetry x-z plane to save computational effort. It is easily seen that as marching proceeds, 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, using Equation (2.2) and noting that Y,, on = —uy = 0, as: —u tan(.’i (2.22) \2J 21 For 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 was Isectional View] Wedge Section Free Surface z IPerspective View! —U id y Figure 2.3 : Infinite wedge definition sketch fixed at an x increment of 2.5 mm. Draughts of between 4-16 inches (10.2-40.6 cm) were considered. The water depth and tank width were both taken to be 2 metres. Linear elements were used throughout. Normal velocities for all boundaries except for the vertical boundary of the wedge and the free surface were specified as zero. The Dirichiet boundary condition was specified on the free surface through the procedure given in the last section. The amplitude and location of the first bow wave peak was recorded for all computations. These results were compared to measured values for 3 1972) and presented in Figures 2.4 and 2.5. = 150 (Ogilvie, In Figure 2.4, the bow wave amplitude is in good agreement with experimental values. The analytical results obtained by Ogilvie 22 (1972) for the bow wave amplitude, non-dirnensionalized using -- /U, showed a constant value of 1.6, independent of the depth Froude number as well as the draught. On the other hand, the same non-dimensionalization of all data points in Figure 2.4 gave a range 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 in all cases. This is probably due to the fact that the marching procedure was started with an undisturbed free surface. Generally, initial conditions in the form of potentials and wave elevations at the first cross-flow plane, that is at the bow section of the ship, can be calculated from Dawson’s (1977) method. This method basically requires the distribution of Rankine sources on the body and free surface, solving for the source strengths by imposing body, free surface and radiation conditions. Evaluation of velocities and free surface elevations rely on these source strengths. Hence, a closed body is necessary and the Dawson method was not applicable for the above wedge which is infinite in the longitudinal direction. 23 Ioraught = 10.16 cml 20 • 20 Present method Ogilvie (Expt) 15 15 a) 0 -o 10 a) a) 0 Draught = 30.48 cml • Present method —-——Ogilvie (Expt) a) 10 E a- 0 0 100 200 300 400 0( 500 100 Speed of Advance (cm/s) 20 • Present method ————Ogüvie (Expt) 15 15 4: 100 200 300 • Present method -—-------Ogilvie (Expt) .::... a •.. 0 iDraught = 40.64 cml • a) a 4:... 0 500 I0 E S. a) 400 a) S.. a) -o 10. 300 Speed of Advance (cm/s) IDraught = 20.32 cml 20 200 400 Th 500 100 200 300 400 500 Speed of Advance (cm/s) Speed of Advance (cm/s) Figure 2.4: Calculated and measured bow wave amplitude for wedge, 24 150. E 0 0 1:0 E IDraupht =10.16 cml 125 0 75 0 a. U) a. • C • .2 • 25 0 E o . 0 U) E 0 • Present method ——Ogilvie (Expt) 100 Ioraught —. 200 300 400 500 100 • 75. 50 E 0 0 E 125 0 = 20.32 100 cml 0 a. / / • •. • 0 100 200 300 400 125 Present method Ogilvie (Expt) • o 100 • E 75 .250 25 • • U) U) • 0 . a. . 100 200 300 400 500 • • • ... 0 100 ..I. 200 .. I.... 300 I.... 400 500 Speed of Advance (cm/s) Speed of Advance (cm/s) Figure 2.5 • U) U) ... ••. • a. .. • 500 IDraught = 40.64 cml 0 C 0 U) . Speed of Advance (cmls) Present method Ogilvie (Expt) • . • Speed of Advance (cm/s) lDrauht 30.48 cml = Present method Ogilvie (Expt) . 25 a. CD 100 • • • 0 125 Calculated and measured wave peak location for wedge, 25 = 15°. Chapter Three Slender Body in Two-Layer Medium A ship travelling at very low Froude numbers in stratified water is known to experience what is generally known as dead-water effects. As the ship speed approaches that of the fastest internal waves, the increase in drag over the hull is very large and the available power 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 the creation of internal waves below the free surface and the disturbance of the flow around the ship due to these waves. This chapter describes a numerical solution for the time domain calculation of surface and interfacial waves generated by a body moving slowly in a gravitationally stable dual density fluid. The wavemaker approach described in Chapter Two and Appendix A were used. The fluid is assumed to consist of two distinct layers of slightly different densities with a discontinuity of the Brunt-Väisälä number at their interface. 3.1 Boundary conditions A moving frame of reference with the origin on the ship centreline at amidships and on the undisturbed free surface as given in Figure 2.1 was adopted. A spheroid with the same dimensions as that in Ma and Tulin (1992) was used for comparison purposes, see Figure 3.1. The normal velocity on the body boundary is given by: 26 = (3.1) on The prolate spheroid geometry is defined by: y = Y(x,z)= ±[(i_4 )2 _z2] (3.2) where B and L are the beam and length of the spheroid respectively. Profile Length L — Front 43.70cm Beam B = 890cm Figure 3.1 : Prolate spheroid and dimensions. The computational domains, one for the upper layer with density Pi and one for the 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 pynocline is small. 2 Consequently, the Brunt-Väisälä number N function, being infinite at the interface. 27 = —gOp / c3z behaves like a ô Body F C 0 B H A Bottom boundary Figure 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 = 1 -h (3.3). Normal velocities at ABC and FGH (Figure 3.2) are taken to be zero, simulating a solid wall. Since very low ship speeds are considered, the surface wave generated by the hull is expected to be small. Hence, the free surface conditions iii Equations (3.4) and (3.5) were used: ÷gi° (1) I =0 (1) = öz where i- at z=0 (3.4) atz=O (3.5) is the surface wave elevation and g = 9.81 mIs 2 is the gravitational constant. The classical treatment of two dimensional internal waves provide a condition relating the interface between the two layers. While normal velocities and pressures are 28 taken to be unique at the interface, the velocity potentials go through a discontinuity across that interface. Taking the normal velocity as positive pointing out of a domain, the following two conditions can be written for an upper layer with thickness h : 2 (2) an at 1 ap where = 2) (1) ) 2 p( (2) = an + g2)(p(2) — p(l)) 2 atz=—h (3.6) 2 —h (3.7) at z = is the internal wave elevation with respect to the undisturbed level at the interface. Nodes at the interface belonging to the upper and lower layers coincide and may be viewed as ‘double nodes’, sharing the same physical location. The normal velocities at these nodes are assumed identical as given in Equation (3.6). The negative sign for the upper layer is simply due to the convention adopted, that is, normals are positive pointing out of the domain. Equation (3.7) results from equating pressure at the ) is taken to be of the order of the ship 2 ‘double nodes’. The thickness of the top layer (h draught. 3.2 Numerical procedure An overview of the numerical procedure used for this chapter is given in Figure 3.3. The lower and upper domains are solved sequentially in that order. The simulation starts by considering the lower domain at time ( t + ôt ), which consists of three sides where the normal velocity is zero and one side (the interface) where the potential needs to be predicted. A symmetry plane at y=O was imposed for computational efficiency. The first order Adams-Bashforth predictor given by Equation (2.15) is used to estimate the 29 potential on the interface. Emerging from the solution are normal velocities at the interface. On the boundary between B and G (Figure 3.2) these normal velocities are used as input for the upper domain (see Figure 3.3). Figure 3.3 : Overview of the numerical procedure For the upper domain, normal velocities are again specified at all nodes except at the free surface, at which the ABM predictor given in Equation (2.15) is again used. After solving the upper domain, required time and space derivatives at calculated. On the free surface, the time rate of change 30 ( t + ôt ) are to be can be found using Equation (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 wave elevation. The internal wave elevations are found using Equation (2.18). For the boundary between B and G belonging to the upper domain, the time derivative of the potential is found using a three point difference formula in time, as follows: = )(t 1 3p( + ót)— 4p’)(t) + p°)(t — ôt) (3 8) 2ôt With of ) 2 p pl) at the interface determined, Equation (3.7) is used to determine values This completes the predictor part of the computations. The Adams-Moulton corrector 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 the end of the corrector part. The above predictor-corrector pair is repeated for all subsequent time steps. 3.3 Numerical results A prolate spheroid of L = 43.70 cm and B = 8.90 cm was used in the calculations. The densimetric Froude number, a depth Froude number modified by the density difference between the two layers, is defined as: c 31 1 F (3.9) 2 )h = where U is the advance velocity, op is the difference in density between the two layers and 2 is the thickness of the upper layer. Cases where the value of F h 1 > 1 are considered ‘supersonic’, Tulin and Miloh (1990), a term imported from aerodynamics. Note that 1 F = 1, termed the ‘critical’ densimetric Froude number, is the point where the slope of the ship drag curve is very high.The densimetric Froude number is related to the ship Froude number (FR ) by the following relationship: FR = F0Ph2 (3.10) Two cases were computed: (1) 1 F = 3.00, h 2 = 6.70 cm, U = 12.86 cmJs, Op/p = 0.0028. (2) F, = 2.95, h 2 = 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) is less 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 wave elevation is small but the internal wave signature is noticeable. 32 The interfacial displacement and wave field contours are plotted in Figures 3.5 and 3.6. The distinctive features of the wave field are very similar to those measured and calculated by other researchers. A very deep depression occurs at the interface just below the bow, followed by 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 two cases are plotted in Figures 3.7 and 3.8. Generally, good agreement with measured crest patterns were achieved. However, the present method predicts earlier formation of the interfacial crests. In addition, decay of the first interfacial crest for Case 2 (Case 1 data was not available) are in excellent agreement with experimental values, see Figure 3.9. The decay of the first interfacial crest can be described using: ) 2 k\h (3.11) ) 2 h The value of the decay parameter a obtained using the present method was 0.96 while Ma and 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 crest heights can be clearly seen in Figure 3.10 to be well predicted by the present method. 33 I. CD C 1j C) + 1\) C;’ C to C;’ to 0 C C) x II II II 0 V ODOJ cnaO II E C) N E C.) >.. E C) >< E C.) Q S C) c’J c’J II < a) > 0 3 C.) w C) C) (D U Figure 3.5 : Interfacial displacement for prolate spheroid at F 35 = 3.0. z .o. (P 0 -500 -400 -300 -200 -100 x(cm) -100 -50 0 50 100 = U 12.86 cm/sec 45.70cm y (cm) = L Bow above x = 22.85 cm y= 0.00cm 3 2 1 6 5 4 9 8 7 -0.15 -0.19 -0.24 -0.34 -0.43 -0.52 -0.61 -0.70 0.22 0.19 0.18 0.14 0.12 0.10 E D C B A 0.59 0.49 0.40 0.31 I H G F Level z (cm) Solid symbols : Present Method Hollow symbols: Ma & Tulin (Expt) y/h, I 2 3 2 = 4.25cm 4 5 0 0 50 100 150 200 250 300 Figure 3.7: Location of interfacial crest lines, F 1 = 2.95. Solid symbols: Present Method y/h 2 25 hollow symbols: Ma & Tulin (Expt) I : 2 3 4 5 20 50 0 100 150 200 250 300 Figure 3.8 : Location of interfacial crest lines, F /h 1 (i ) log 2 -0.50 = 3.00. j112 = 4.25 cm] ‘ • O 75 -1.00 Ma & Tulin (1992) Present method . -1.25 • 1.00 1.25 1.50 1.75 2.00 , • 2.25 log (x/h ) 2 2.50 Figure 3.9 : Decay of first interfacial crest, F 37 = 2.95. Symbol ) 2 log(/h ‘ 0.0 . • .- -0.5 Method hid F, 0.61 0.96 1.15 0.96 3.20 2.95 3.00 2.95 3.00 Expt ExpI Expt Present Present - -1.0 .... • -2.0 1.00 1.25 1.50 •.... • V • .. -15 • S 1.75 200 2.25 2.50 ) 2 log (x/h Figure 3.10: Influence of upper layer thickness on first interfacial crest. 38 Chapter Four Hybrid Boundary Element/Discrete Vortex Method Separation from a sharp edge results in the formation of a vortex sheet issuing from the edge. This vortex sheet can be modelled by a series of discrete vortices introduced one at a time into the flow field at regular intervals. The motion of each vortex is traced over time using its convection velocity. If the vortex shedding is restricted to a region close to the edge, the discrete vortex method can be viewed as the inner region solution in the overall boundary value problem. This inner region solution has to be matched with the outer potential flow solution. The combination of boundary element and discrete vortex methods provides this matching and at the same time do not require calculations inside the domain. It is realized that the boundary element method is compatible with the discrete vortex model in that each discrete vortex can be treated as an internal singularity which can be handled using an analogue of the Residue Theorem in complex analysis. Significant computational advantages result because of the relatively simple approach to handling separation at the sharp edges while working only with the boundary values. In this chapter, the hybridization of the discrete vortex and boundary element methods will be described. The hybrid method was applied to the calculation of the forces on a flat plate oscillating in unbounded fluid (Wong 1993a). 39 4.1 Discrete vortex method The discrete vortex method is a time stepping procedure which models the shear layer issuing from a sharp edge using discrete vortices (Wong 1990). In the formulation, an infinite 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 = 2 — planes and is a parameter dependent on the internal angle iS of the wedge. Transform Plane Physical Plane Free stream velocity (V) [tvortexatz, st,ngthy / [ievortexat strength Free stream velocity (V) rvoexai strength Figure 4.1 : The physical Z-plane and the transformed -pIane. The non-dimensional external flow velocity, v sin(2rrt), wheret = t / T (T being the period of oscillation) is in a direction normal to the wedge bisector. The strength of the vortex at z and its corresponding point is given by i. see Figure 4.1. The complex potential in the -p1ane, with n vortices in the field, is then given by: 40 F()=iV÷_Lrkln (4.2) k-O Applying 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 the physical plane as follows: w.=S.{v+ — k-O,k.j j j _ + — )—.!L( 2t ‘_ 2 — + )] (4.3) where the last term on the right hand side is due to Routh’s correction. The length scales used for the normalization were derived in terms of an attached flow velocity scale in the physical plane. These are (from Wong 1990): = (1 l 2 ( 1 (D)k L L = (M)1(D)1/(2k_1) L = (M)i(D)1/(2_1)V (4.4) where V 0 the amplitude of the oscillatory velocity and D is the maximum distance traversed (V T). 0 Assuming that all vortices convect with the flow, the convection equation in the physical plane is given by: (4.5) Zk ‘‘k A new vortex (called the nascent vortex) is introduced along the external wedge bisector at the end of each time step. The distance of the nascent vortex from the sharp edge is 41 taken to be a function of the time step size, the wedge angle, the free stream velocity and the strength and location of all the other vortices already in the flow field. The expression for the initial Location of the nascent vortex, denoted by z , can be found through the 0 following equations (derived in Wong, 1990): 0 z = where m= tm (KAr) (4.6) 2). -1 (4.7) 3 23t) where v is defined in Equation (4.10). The strength of the nascent vortex can be obtained by satisfying the Kutta condition at the edge. The Kutta condition is satisfied when 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 the corresponding position in the transform plane is zero. As a result, the following expressions can be obtained: dFQ) [ d giving •ftj where v V = k-O 2rt k (4.8) )j (4.9) R(;) nsin(2ct) — Ik Re(k) (4.10) With the strength and location of the nascent vortex determined, the algorithm then proceeds to the next time step. Note that the location of the first discrete vortex shed 42 into the flow is defined by Equation (4.6). A vortex decay mechanism is used to account for the vortex pairing process described in Graham (1980) and provides a model for vortex diffusion. For this, a simple exponential rule reduces the effective strength of each vortex to one half of its original value at the end of one time cycle (Equation 4.11). The vortex decay mechanism is thus an attempt to simulate the loss of total vorticity through the cancellation effects in positive and negative strength vortices coming together (in pairing) 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 successful in maintaining computational stability over a large number of flow cycles. Also, reducing circulation in the wake limits the magnitude of the vorticity introduced. Consequently, the otherwise over-predicted vortex induced forces are kept at levels which are comparable to experimental values (as reported in Wong 1990). An example of the calculated vortex roll-up for a flat plate at t=1.O is given in Figure 4.2. 0. 0 Second spiral o Flat Plate 0 First spiral a ‘jo. Figure 4.2: Vortex sheet rollup at t 43 = 1.00 for a flat plate. 4.2 Vortex shedding in boundary elenieiit method In this section, the inclusion of vortex shedding effects into a boundary element model is discussed. We define a domain D bounded in S which includes branch cuts C 1 and C 2 and a circular path S to exclude a vortex singularity located at the point K as in Figure 4.3. The circle enclosing point K has a radius E. For the velocity potential (p which satifies the Laplace equation, the divergence and Green’s theorems can be used to obtain: Figure 4.3 : Exclusion of a vortex singularity in domain D. S (4.12) ——G-dS=O On on The weighting source function which also satisfies the Laplace equation is taken to be G = ln(r), r being the distance between the point under consideration and any other point on the boundary. 44 The 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) ufl where y k is the strength of the vortex at K. The boundary S is discretized into N straight line elements. Putting H --, = the discretized form of the integral equation for the point i on the boundary with N elements enclosing a domain with M vortices within can be written thus: 2JHIJ +it where i =1 ... N + (4.14) = k—N+1 j—1 j—1 j 1 can be found M, ô is the Kronecker delta function and G and H numerically (see Appendix B). Note that the summations involving G and H do not include the ‘elements’ around vortex singularities (j> N) as the former need not be evaluated since (a / on) 1 i> N while the latter can be shown to vanish as the radius £ —* = 0 for 0. It is also noted that the branch cuts C 1 and C 2 do not contribute to Equation (4.12). This result is expected as the treatment of the vortex singularities is similar to the way the potentials at points inside the domain are found. The above problem is overspecifled in that there are N+M equations in N unknowns and only the first N equations are used. In time domain simulations, the boundary values and vortex strengths are to be prescribed for the next time step (t + bt) in order to march forward in time. The matrices 45 G and H need to be calculated only once for a fixed boundary. positions of all discrete vortices in the flow at time (t + The strengths and ö) were found using the discrete vortex method described in the previous section. The potential and velocity induced by the discrete vortices in the flow on any boundary node j is given by: = YkOkj (-) vj (4.15) 2 k-i Yflj = (4.16) k-itkj where O is the angular position of the node j with reference to the kth vortex. This vortex has strength k and is at a distance Rk 3 from the node j. In Equation (4.16) n is the nomial vector 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: = + for i = 1 ... () ] — (4.17) N. Normal velocities are prescribed on the body. The system of N equations can be solved for the N unknowns . The dimensional hydrodynamic force per unit depth acting on Sb at any given instant is taken to be: (4.18) Fb =p.ndS where p is the hydrodynamic pressure (from Bernoulli’s equation) acting on the discrete d to 2 0 element dS with normal vector n. The calculated value of Fb is normalized using pV give a non-dimensional force coefficient C where: 46 C =--F In the above, Fb = (4.19) 0 I T)d }F j23tp(V 8 for non-dimensional FB. 4.3 Forces on sharp edged cylinders 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 in the last section. Published numerical and experimental results, for example those in Graham (1980), Kudo (1981) amd Keulegan et al. (1958) could be used for comparison purposes: not to judge the relative merits of the numerical schemes proposed by other researchers but as a test case for the present method. Discrete vortices o 0 00 Pn, 00 0 0 0 Discrete vortices Figure 4.4 : Computational domain and boundary conditions. The simulation of vortex shedding from a finite wedge oscillating in unbounded fluid requires some length scale modifications since the discrete vortex method described In other words, the vortex strengths and earlier is applicable to an infinite wedge. 47 positions calculated with the discrete vortex method need to be adjusted to account for changes in Keulegan-Carpenter numbers (KJ. Graham (1980) showed that the ratio of the length scales of the infinite and finite wedges is proportional 10 0 (K)’ . The resultant non-dimensional vortex strength (y B) and length (L ) scalings to be used for 8 the boundary element procedure tor finite K and body length (2d) is then given by: 8 =y(2K) Y (4.20) k (4.21) ZB=z(2KC) It should be noted that the above equation is general in that the vortex shedding edge angle 0) is embedded in the parameter k. The computational domain arid relevant boundary conditions are given in Figure 4.4. On the body, the non-dimensional nomal velocity is given by: ap — an o sin(2tt)cos(—).n 2 . = (4.22) For this application, all the vorticity shed from a particular edge is placed at a point very near to that edge. In other words, we ‘average’ all the vortex positons and the influence of all discrete vortices on any given boundary node is calculated through the relative positons of that node and the edge. With that, the value of 0 in Equation (4.15) is taken to be ±(t/2) depending on whether the top or bottom vortex shedding edge (Figure 4.4) of the body is considered. In Equation (4.16), R is taken to be the distance from the edge to the node under consideration. This simplification is acceptable for small K since much of the vortex shedding effects are confined to a region near to the edge. 48 Values of cF for a flat plate oscillating (ô = 0) at K = 6.6 are plotted together with experimental results obtained by Keulegan et al. (1958) in Figure 4.5. The present method predicted C values that are in excellent agreement with experimental results. Drag (Cd) and inertia (C ) coefficients are found by taking Fourier integrals of Morison’s equation 1 over a time cycle and given thus: Cd 5 C 3 .mn+1) = 4nc 0 2K = 2 C sin(2rt)dt j.(u+l) J (4.23) C cos(2mt)dt (4.24) where n is an integer representing the cycle number. Present method Keulegan et at (expt) E a) C) a) 0 C-) a) C.) 0 LL 3.0 4.0 Number of Cycles Figure 4.5 : Forces on an oscillating flat plate, iCc Calculated drag (Cd) = 6.6. and inertia (C ) coefficients are presented in Figure 4.6 and 1 Figure 4.7. Predicted Cd for the flat plate agree well with numerical and experimental 49 results from Kudo (1981) and Keulegan et al. (1958) respectively. At K above 3.0, predicted d appears to underestimate experimental values. However, Obasaju’s measured data, as reported in Graham (1980), showed Cd values lower than those found using the present method. Therefore, given the variation of prevailing conditions under which different researchers obtained their experimental results, the present method predicts drag forces on the normal flat plate very well. 15.000 Present Method Keulegan et al (sept) Kudo (nianencal) 12,500 — —4— o A 10.000 C) 8 7500 0) 0 Oo 2 501) 0000 • 0 • • 10 5 Kc I 15 Figure 4.6 : Drag coefficients for the flat plate. Values of C. obtained using the present method show an upward trend with increasing K which is similar to the numerical results presented by Kudo (1981). On the other hand the experimental results of both Keulegan and Carpenter (1958) and Obasaju 0 3 up to about K (presented in Graham, 1985) showed an upward trend for C there a decreasing Ca, up to K. 12 . 7 and from Forced symmetric vortex shedding from the two . Graham (1985) 3 edges of the plate is thought to be the reason for the poor match for C 50 reported much better results when vortices are allowed to shed asymmetrically through the introduction of a small perturbation at the initial stages of the simulation. 5.000 4001) Present Method Keulegan et at (sept) - + -. 0 A Kudo (numericat) — — — • a) _r -.- 3000 — Cl) 8 2000 0 0 0 o0 0 1.000 Kc 0.OOC.... 0.0 2.5 5.0 7.5 10.0 12.5 15.0 Figure 4.7 Inertia coefficients for the flat plate. Finally, it should be mentioned that drag and inertia coefficients become more realistic through the incorporation of vortex shedding forces. Potential theory alone does not give any drag and inertia values for attached flow can be shown to be Ca 0 = 1.0 for the flat plate. Use of the same transfonn equation in the discrete vortex model for different edge angles provide a useful and flexible tool in the investigation of a wide range of flow problems involving sharp edge vortex shedding. Although the discrete vortex model can be applied in conjunction with other numerical methods for the given class of problems, it is uncertain whether such an exercise would be more fruitful than working directly with the Navier-Stokes equations for viscous incompressible flow since the advantage of’ working only with elements on the boundary is lost. 51 Chapter Five Applications in Nonlinear Slender Body Hydrodynamics The wavemaker approach opens up the possibility of attaining a nonlinear solution to wave generation problems involving slender forms in steady and unsteady motion. The fully nonlinear free surface conditions, implemented using the Eulerian-Lagrangian approach, was used for all the work reported in this chapter. In the next section, the utilization of the wavemaker algorithm developed in Chapter Two for the prediction of the waves generated and forces acting on a Wigley hull is described. This algorithm was then extended to include the effects of vortex shedding from the keel of the Wigley hull advancing obliquely. The computed results are presented and assessed by comparison with experimental and numerical results of other researchers wherever possible. 5.1 Wavemaking resistance of a slender body The computational boundary, identical to that given in Figure 2.2, is discretized into N quadratic elements, that is the variation of potential and normal velocity over an element length dS is assumed quadratic (see Appendix B). The hull is described by: y = Y(x,z)= (5.1) 52 where B, L and d are the beam, length and draught of the ship. The normal vectors on the hull, in the slender body approximation, is given by: (_y ,±i,y. L) ny,n)=j ,(n 1 (5.2) where the subscript in Y refers to the partial derivative with respect to that variable. U z y Wigley Hull Length 200 cm Beam 20cm Draught 12.5 cm Figure 5.1: Wigley hull dimensions and coordinate system. The normal velocity on the obliquely positioned hull, with yaw angle a, is then given by: a4 —U(Ycosa±sinct) (5.3) 2 ‘/1+Yz 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) were used. Neumann conditions are specified for the rest of the boundaries. At time (t+ôt), the 53 Adams-Bashforth predictor is used to estimate the velocity potential and elevation of the free surface. Normal velocities on the free surface are then calculated as part of the solution. From these, using the procedure described in Chapter Two, the wave elevation is updated and the rate of change of the potential found using Equation (2.9). The Adams-Moulton corrector is then applied for finding the free surface velocity potential. Solving the problem using the corrected potentials then completes the calculations for that time step. There are weaknesses associated with applying a wavemaker approach without enhancement, as seen in the advancing wedge application in Chapter Two. In the vertical direction, the intersection between the ship side and the free surface need to be treated with care. In addition, there is also the question of numerical accuracy as well as stability due to the high wavemaker acceleration at the start of the time marching for most ship types (with the possible exception of cusp-ended ships). The problem is more serious for non-vertical ship sides. To overcome this, the normal velocity at the ship-wave intersection point is simply extrapolated from values calculated at free surface nodes nearby (see Appendix B for more details on the extrapolation procedure). In addition, the evolution equations for the free surface nodes involve a finite difference approximation to the tangential velocity component as described in Chapter Two. These gave rise to higher frequency errors, which over time resulted in stability problems. A rediscretization scheme was applied to the entire free surface boundary to filter out these high frequency errors. The regridding procedure was based on linear interpolation over neighbouring nodes returning equally spaced nodes on the free surface boundary once every time step. 54 A smoothing scheme identical to that used by Longuet-Higgins and Cokelet (1976) was applied after the rediscretization process. The value of 4 at the fluid-body intersection node for the latest time step was used as a means of accounting for the variation of the longitudinal velocity over the body. This modifies the free stream velocity (U) incident on the hull to give U’ which was then applied to Equation (5.3). The following relationship, based on the value of U’ = ap 4 at the fluid-body intersection, was used: Ucosa+X(SB flS)[1_.] -(U’Y ±Usina) 1 g1+y2 The principal dimensions of the Wigley hull used in all the calculations in this Chapter are given in Figure 5.1. The impermeable far wall was placed 30 beams away from the body. On the free surface, 60 equal length elements were used. For the free surface plots, one sweep over the ship length was done in 50 steps. However, for the pressure distribution and forces acting on the hull, extensive numerical experimentation on the number of steps required (NS) per sweep resulted in the following formulae: X NS = INT[108(1.1 — NS = INT[54(3 4FR)] — FR — 4FR)] 0.20 FR 0.25 <FR 0.25 0.40 (5.5) More steps were necessary for the lower Froude numbers due to the relatively larger variations in velocities and pressures in the longitudinal directional. In all cases, the value of NS was limited to a range of (70 NS 55 200). The lower bound on NS was imposed to ensure adequate resolution in the length-wise pressure distribution while the upper bound limits computation time per sweep over the hull. Initial conditions, that is wave elevation and potential and their rates of change with 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-body method, handling the linearized free surface condition via a streamwise marching technique. The hydrodynamic forces can be found, from Equations (2.6), using the pressure from the Bernoulli equation. Substituting p for the pressure: p 1 F = _P[U(coscz + sin a) + + () ) + (5.6) z] L/2 (5.7) dS 1 fdxf p.n = —L/2 SH where i is the required direction (x, y or z) and B 5 represents the girth of the hull section 1 and F under consideration. The wavemaking resistance and sinkage force are then F respectively. The trim moment (corresponding to moments about the y-axis) can be found using Equation (2.7): M L/2 = (5.8) 5 p.(nz—nx)dS fdxf —L12 The wave profile along the ship side at Froude numbers F = 0.267 and FR = 0.316 are presented in Figures 5.2 and 5.3. For both cases, the results calculated by the pesent method are generally in good agreement with experimental and numerical results obtained 56 by other researchers. The experimental results for FR = 0.267 was obtained from the Department of Naval Architecture and Ocean Engineering at the Yokohama National University Hayashi and Kunishige, 1988), referred to as YNU hereafter. Measured values for a 6m Wigley huil model reported in Maruo and Song (1990) are used in Figure 5.3, for FR = 0.316. /L 1 2’,- . 0.04 Present method (NonlEnear) Measured (YNU) Maruo & Song (1990) 0.03 0.02 0.01 -0.00 - • - -0.01 -0.02 -0.03 —0.04 AP FP -1.25 -1.00 -0.75 -0.50 -0.25 -0.00 0.25 0.75 0.50 1.00 2x/L 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 the wave-body intersection point for stability reasons. Inaccuracies in the value of modified free stream velocity given in Equation (5.4) may be another cause for the delay. These led 57 to the delay in the crests and troughs over the rest of the hull length. However, since the lag 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 this section, results in a wave profile that totally missed the variations over the middle half of the ship length. /L 1 2r . 0.04 Present method (Nonlinear) Measured (6m model) Maruo & Song (1990) 0.03 0.02 0.01 -0.00 -0.01 -0.02 -0.03 AP FP - -0.04 -1.25 -1.00 -0.75 -0.50 -0.25 -0.00 025 0.50 0.75 1.00 2x/L Figure 5.3 : Wave on side of Wigley hull at FR 58 = 0.316. The contour plot of the wave field in Figure 5.4 agrees fairly well with YNU measured data (Figure 5.5). A perspective view of the wave field for the Wigley hull at Froude number F = 0.267 is shown in Figure 5.6. y (ciii) Level 0 0 B A 9 8 7 6 5 4 3 2 1 5.50 4.50 3.50 2.50 1.50 0.50 -0.50 -1.00 -2.00 300 -4.00 -5.00 -6.00 0 FP AP 2x/L 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 using 2 on the hull surface for F 0.5pU = 0.267. The high pressure region in the bow is clearly noticeable. 59 N. Lu’ I-’ cr1 cr1 Lul bQ 1 Jo \1< L __J I). Figure 5.5 : Measured and calculated wave field (YNU) for Wigley Hull, a 60 = 00. I— 0% OQ CD I 0 0 U’ QQ -U I 0 > •T1 -U x1 0 p Tj = I VI a’ IJQ The wavernaking resistance, non-dimensionalized by O.5pU 2L , is compared to 2 experimental and numerical results obtained, elsewhere (Figures 5.8 and 5.9). Values obtained using linearized free surface conditions are highly inaccurate (also reported by Allievi, 1993). The general agreement with measured data is considered to acceptable. It should also be mentioned that even experimental results from different towing tanks exhibit significant variability depending on the conditions under which the data was measured. Aanesland (1989) showed ‘averaged’ values from towing tank databases for cases where the Wigley hull was towed fixed and where the vessel was allowed free to turn and squat (Figure 5.9). Aanesland’s numerical results were obtained using a three dimensional potential flow algorithm, based on Rankine sourcç distribution and strearnwise marching ) x 1O L 2 RJ(.5pU Numerical (Maruo & Song) (YNU) Wave cut analysis ‘Nu) Present method (Nonlinear) 0.50 0.40 0.30 ‘ . ,. 0.20 . -, / / 0.10 FR 0.15 0.20 .1... .1.. 0.25 0.30 ..I 0.35 0.40 0.45 Figure 5.8 : Wavemaking resistance against YNU results. 62 schemes similar to that of Dawson. Kelvin sources, distributed in the outer region, were matched to the inner region (in the vicinity of the hull) to yield a radiation condition at the interface between the inner and outer regions. Calculated results from the present method are cast among these in Figure 5.9, showing the present results to be within the range ‘acceptable’ values. v free (Expt) fixed (Expt) Aanesland (Numerical) • Present method (Nonlinear) S) x 1O 2 RJ(.5pU 3.0 2.5 2.0 ,.. 1.5 • • • -.. . . ./ I -— ‘I l.0 V.. 0.5 oci 0.150 - FR .- 0200 0.250 0.300 0.350 0400 0.450 Figure 5.9 : Wavemaking resistance against Aanesland (1989). Sinkage and trim for various Froude numbers are shown in Figure 5.10. Sinkage is non-dimensionalized using U 2 / 2g and trim is given as a percentage of the ship length. These results are compared to averaged measured values and those calculated by Aanesland (1989). The non-dimensional sinkage calculated by the present method are in 63 good agreement with experimental data but tend to underestimate at the higher end of the Froude number range. On the other hand, the non-dimensional trim tends to be over estimated in the Froude number range of 0.27 to 0.35. Experimental (Lower Limit) Experimental (Upper Limit) Calcutated (Aanesland) Sinkage/(U / 2 2g) Present method • 0.00 -0.02 -0.04 ., :-•-. . V -0.06 - -0.08 —010 • 0.10 FR 0.15 I... 0.20 .1 0.25 0.30 0.35 I.—.. 0.40 0.45 Trim/Lx 100 (%) 2.00 v 1.50 1.00 7 0.50 • 0.00 -0.50 0.10 :. FR I 0.15 0.20 0.25 0.30 0.35 0.40 0.45 Figure 5.10 Computed and measured trim and sinkage for Wigley hull. 64 5.2 Wavemaking resistance in shallow water The 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 tank dimensions used were identical to those in the previous section. However, the number of steps required (NS) was modified to reflect the fact that higher resolution is necessary for shallower water. Hence, with NS defined in Equation (5.5), the required number of steps (NS*) is given by: NS* = 6 NS) Hid H/d <6.0, 70 NS* 400 (5.9) In Equation (5.9), H and d are the water depth and draught of the ship respectively. 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, the increase in wavemaking resistance near the critical depth Froude number of FH = 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 depth Froude 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 the present method predicted a slower decline in wavemaking resistance with increase in water depth. The effect of shallow water on squat and trim are plotted in Figure 5.12. 65 RJ5pIYS x 1O 15 Calculated (Kim and Choi, 1993) Present method 10 5 . .1.... C 0.0 2.5 5.0 7.5 10.0 H/d Figure 5.11 : Wavemaking resistance for Wigley hull in shallow water at FR = 0.316. = 0.316. Non-dimensional Trn and Sinkage 0.75 /2g) 2 Sinkage/(LJ Tnm/LxlOO% 0.50 0.25 0.00 O2 0.0 •. . 2.5 . . I 5.0 I 7.5 10.0 H/d Figure 5.12 : Effect of water depth on squat and trim for Wigley hull at FR 66 5.3 hybrid method in yawed ship calculations The steady motion of an advancing ship generates vorticity. The strength of the vortices shed determines their significance in the computation of the forces and moments acting on the ship. The influence of vortex shedding is, in general, not very important in the cases involving a ship in steady forward motion in calm water. This is exhibited by the useful results 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 the bilge keels and keel is not necessarily neglible. In the context of the slender body approximation, the translation of the body sections in the cross-flow plane as the marching proceeds presents a vortex shedding sharp edge at the keel. The influence of such vortex shedding on the forces and moments acting on a yawed Wigley hull and the effects on the free surface wave field will be evaluated in this section. However, before using the hybrid discrete vortex-boundary element method described in Chapter Four, it is first necessary to explain the vortex generation mechanism and 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 conveniently defined and the method described in Chapter Four cannot be applied. occuring at the bow were not modelled. Also, vortices In addition, secondary vortex shedding was ignored. In other words, the vortex shedding procedure was started at the commencement of the time step marching for the wavemaker algorithm. The introduction of potential vortices into the flow field results in their influence on the potential and normal velocity at 67 each boundary node. On the other hand, the assumption that dissipation does not occur in the overall flow allows the use of Equation (2.8) for the evaluation of pressure on the hull boundary. Before the start of the marching procedure, it was assumed that there are no vortices in the cross-flow plane. As the time marching proceeds, the hull transverse section changes in shape and the keel can be viewed to be translating in the horizontal direction. 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 the keel as the marching proceeds towards the stern. This spiral can be modelled through the introduction of discrete vortices, at increments of one per cross-flow plane. It was necessary to substitute the oscillatory velocity in the discrete vortex algorithm (Chapter Four) 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 varies along 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 were then inserted into the modified system of equations (given in Equation 4.17 ) and the solution of the boundary value problem obtained. The ‘free stream’ velocity for the vortex shedding algorithm, as seen by the keel, is taken to be the translatory velocity of the body section such that: (5.10) V=Usinct 68 where V here refers to the free stream velocity in the vortex shedding model in Chapter Four. Length and time scales are required to extract the necessary vortex positions and strengths from the self-similar vortex produced by the discrete vortex method. This can be done by defining the scaling factors, using Equation (4.4), and defining M = 1., D = I.sin a and V 0 = 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 a From the above, the non-dimensional positions and strengths of each vortex in the cross-flow plane can be dimensionalized for the solution of the overall problem in the following manner: YV_YB+Im(ZV)LZ z, =—d—Re(Z)L = (5.13) FL where Z is the complex variable Z, = X,, + iY in the vortex physical plane and f is the non-dimensional strength of any given vortex. In Equation (5.13), yo and d refer to the lateral location of the ship’s keel and the ship draught respectively. The hybrid method 69 can then be used, through Equation (4.17), for the solution in the cross-flow plane. The potential and normal velocity induced by a vortex on a mid-element node are given by the following pair of equations: M , ‘çEkkj k ‘Pkj - v—i (oP’ \nJk (5.14) 3t _. 1k ’kj 2 k-l sin (5.15) 113 are defined in Figure 5.13, M is the total number of vortices in the field where O and R and 3 is the difference between the angle of the vortex induced velocity and the element angle. For nodes at the ends of elements, the angle can be taken as the average of the two neighbouring elements. Note that the angle 3 is simply given by t3 (0 + /2). - V RkI 0 Vortex Figure 5.13 : Vortex induced normal velocity on a boundary element. 70 The far wall was taken to be 30 beams away, with 60 free surface elements on each side of the body. The number of steps required (NS) for one sweep over the hull length 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 The keel thus traverses from that position to y = = -(112) sin a. (L/2) sin a at the stem. The yaw moment is given by: 1/2 M= fdxf p.(nyx_ny)iS —1/2 (5.16) B Figure 5.14 compares the calculated free surface profile, with and without vortex shedding, along the side of the Wigley hull at FR = 0.267 with the measured profiles from Hayashi and Kunishige (1988) and calculated ones from Maruo and Song (1990). The influence of the keel vortices on the free surface profile were minimal. Here, the influence of 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 was no yaw (compare with Figure 5.2). It is noted that the wave profiles were not corrected for heel, caused both by the potential tiow as well as the vortex shedding, which may be significant for a ship in yawed motion. The calculated wave field is presented as a contour plot in Figure 5.15. This is in good agreement with measured data from YNU, see Figure 5.16. A perspective view of the wave field, with the z-scale exaggerated, is presented in Figure 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 71 = 100 and with vortex shedding from the keel. A typical computer run, with 248 boundary elements and marching in 100 steps over the hull length, requires approximately 2 hours on a 486 DX2 66 Mhz personal computer, 19 minutes on a HP 715 (75 Mhz) and slightly under a minute on the Fujitsu VPX24O/10 supercomputer. Most of the developmental work on the programs described in this chapter was carried out on the Fujitsu supercomputer. Present method (With vortex shedding) 2ii, L Present method (No vo,tøx shedding) • 0.08 0.06 Measured (YNU) Maruo&Song (1990) Windwardsiindwardside 002 -0.04 -006 FP -0.08 ..I... -1,25 -1.00 AP . -0.75 -0.50 -0.25 -0.00 0.25 0.50 0.75 1.00 2xJL 2ri/L Lee side 0.07 0.05 -0.02 . -0.05 -0.07 ....................... -1.00 -0.75 -0.50 -0.25 -0.00 0.25 0.50 1.00 0.75 2x/L Figure 5.14 : Wave profile on sides of Wigley hull at FR 72 = 0.267, a = 100. y (cm) 50 25 FP° —z._ h24 .25 1 0.75 -25 0.5 ‘, -‘ 0.25 \ 0.25 -0.5 “ -0.25 \ / / ‘\‘ ‘, if 2 /1 -025-’ 0.25 _ 1.25 Starboard Port \ 90 -0.25 -50 -100 -50 50 0 Figure 5.15 : Wave field for Wigley hull at F = 0.267, a 73 2L 100 = 10°. Figure 5.16: Measured wave field (YNU) for Wigley Hull, a 74 = lO C CD z CD I C,;. 0 CD C,; CD CD CD cM CD CD CD CD x CD CD Cl) C) N 0 CD z lWindward sidel Free surface level FP AP Lee sidel Free surface level 1 .490E-8 0.05 / r 1 .490E-8 0.45 -O.15 -fl2l FP AP Figure 5.18 : Pressure distribution for Wigley hull at FR 76 = 0.267, a = 100. 1.50 ) x 10 L 2 F/(.5pU 1 00 1(.5pU x 102 M ) 3 L 2 0.75 1.00 0.50 0.50 0.OC( 0.25 5 10 15 0.00 20 26 a (degrees) 5 10 15 20 25 a (degrees) L’) x 102 2 F1(.5pU 1.50 Legend 1.00 Solid Itnes: Present Method (no 4iortex shedding) Symbols: Present Method (with vonex shedding) Dashed Lines : Maruo & Song (Numerical) 0.50 a (degrees) Figure 5.19 : Forces and moments for yawed Wigley hull. Calculated values of longitudinal and lateral forces (F 1 and F) and the yaw moment (M ) with, and without vortex shedding, for various angles of attack are non 1 dimensionalized and plotted in Figure 5.19. Also plotted in Figure 5.19 are the numerical results obtained by Maruo and Song (1990), using a Kelvin source based slender body algorithm (linear and without vortex shedding). The present method gave non dimensionalized 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 other hand, the present method returned values of lateral forces and yaw moments which are 77 consistent with those of Maruo and Song.. There is, unfortunately, no published experimental 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 no vortices were included. This can be clearly seen in Figure 5.19 where the lateral forces are higher with vortex shedding for all yaw angles, the difference increasing to about 33% at 20 degrees leeway. As for the yaw moment, the influence of vorticity increases as the leeway angle increases, the difference being 34% at 20 degrees. Yaw moment is taken to be positive in the anti-clockwise direction. Due to the higher rate of change of total circulation with respect to x forward of amidships (see Figure 5.20), the predicted yaw moment for the vortex shedding cases were higher. In Figure 5.20, it can be seen that the growth of the total vortex strength is rapid near the bow and levels off towards the stern region. Fx 102 /s) 2 (m Q = 20° i= 15° u= 10° a =5° 0.0 -1.0 -0.5 0.0 0.5 1.0 2L Figure 5.20 Development of total vortex strength over the hull length. 78 Figures 5.21 and 5.22 show the keel vortex evolulion at 8 sections over the length of the Wigley hull at FR = 0.267, a = 10°. The sequence shows a representation of the vortex roll-up process as the keel traverses laterally. 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.201 = Figure 5.21 : Keel vortex evolution for Wigley hull, a = 100 (forward). The hybrid method just described produced significant vortex effects on the forces and yaw moment acting on a slender body advancing obliquely, especially at higher leeway angles. Assessment of the accuracy of the computed values by comparing with measured data cannot be carried out due to the unavailability of the latter. However, the resulting lateral force coefficients are of the same order of magnitude as the experimental data given by Bradbury (1986) for a simplified ship-like geometry. 79 y (cm) 12x/L = 0.051 12x/L = 0661 y (cm) -15 I2xIL=0.371 -10 -5 5 0 10 15 12x/L = 0.96] Figure 5.22: Keel vortex evolution for Wigley hull, a 80 = 100 (aft). Chapter Six Conclusions and Recommendations 6.1 Conclusions The present study on numerical fluid dynamics has resulted in novel and practical approaches 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 ship wave problem into a series of two dimensional wavemaker problems. This algorithm was generalized for application to ships advancing at small yaw angles. The method was first applied to the study of the bow waves generated by an advancing wedge in homogeneous fluid, 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, using linear and quadratic elements. This relatively straightforward approach gave results that compared very well with published experimental and numerical data. The fully nonlinear free surface conditions were successfully implemented for the case of a Wigley hull advancing in deep and shallow water in Sections 5.1 and 5.2 respectively. The fluid-body intersection singularity problem was avoided using an extrapolation procedure which modifies the influence matrix (Appendix B). A discretization scheme was devised to redistribute free surface collocation points. As a result, the stability of the nonlinear free surface evolution was maintained. A criterion, based on extensive numerical experimentation, was developed for the determination of the 81 required step size for the marching procedures used. Finally, the calculated wavemaking resistance, sinkage and trim are in good agreement with published measured values. 6.1.2 ships in density stratified water A novel procedure which is able to handle dual (and in principle multi-) density media was developed and presented in detail. This scheme was successfully utilized in the study of the 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 the relatively small influence of the free surface, it appears that the slenderness ratio constraint of the method could be relaxed somewhat. This was apparent in the results given in Chapter Three where a prolate spheroid of B/L = 0.20 was used with success. This procedure is general in that it can also be applied to the study of internal wave problems without the presence of a moving body. 6.1.3 Hybrid discrete vortex-boundary element method A hybridization of the discrete vortex and boundary element methods was achieved. This method was first tested on the forces acting on an oscillating flat plate, producing results that are in good agreement with other numerical and measured data. The simplicity of the hybrid method accounts for some viscous effects in the flow without 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 the overall numerical algorithm represents significant benefits in terms of computational efficiency. 82 The hybrid method was successfully used in a nonlinear wavemaker program for a Wigley hull travelling in yaw. Forces and yaw moments were calculated and presented in Chapter Five. Lack of experimental data did not allow an objective assessment of the forces and moments predicted by the present method. However, these results were plotted together with published calculated data (from a method using slender body approximations to the Neumann-Kelvin problem). The present method predicted higher values of longitudinal and lateral forces and yaw moments when vortex shedding was included. The general wavemaker algorithm using nonlinear free surface conditions and the hybrid vortex-boundary element method provides an efficient and fast method for the computation of the wavemaking characteristics and resistance of a slender body advancing at small leeway angles. The reduction of the ship problem into a two dimensional one permits a host of possibilities, which is the real motivating factor behind the development of such a scheme. The method can thus be regarded as a base method with a number of potential applications in the areas of ship motions and manouevring. Shallow water ship resistance and irregular or sloping seabed effects studies are also considered to be possible applications. The algorithms given in Chapter Five may be used to study the effects of demi-hull configuration on the wavemaking resistance of multi-hull vessels. In addition, the numerical study of yawed multi-hull hydrodynamics is merely extension of these algorithms. 83 a straightforward 6.2 Recommendations Although the wavemaker algorithm has been shown to be a useful tool in the calculation of hydrodynamic forces and moments on a moving ship, a number of refinements are recommended for the enhancement of its accuracy and capabilities. The slender body approximation neglects velocity variations in the longitudinal direction near the bow and the stern. This was accounted for in the present work through a simple parabolic variation of the perturbed x-velocity with depth. Also, even for cusp ended ships, the initial step of the marching procedure starts with a finite width transverse section, thus violating the slender body assumption in the local sense. The highly accelerated flow at the bow and stern give rise to singularities at the fluid-body intersections. This problem was circumvented in the present work through an extrapolation scheme. Consequently, the wave elevation near the hull and wetted surface area are sensitive to the length of the free surface elements. This in turn affects the accuracy of the pressure distribution on the hull surface. In addition, it was observed that there are errors in the prediction of the bow wave peak location and lack of recovery of the free surface elevation near the stern. Further study on the stability of impulsively started wavemaker flow in relation to the above problems is likely to enhance the accuracy of the method. All slender body calculations reported in this thesis were carried out with the body fixed. 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 example Aanesland, 1989). In the wavemaker approach, this could probably be achieved by 84 marching through the length of the body more than once, with heel, trim and sinkage values from the previous sweep imposed in the following sweep. Boundary layers created by the ship are likely to alter the effective boundary location and hence the normal velocities required in the wavernaker algorithm. These effects were not considered in this thesis. Another area that deserves attention is the determination of the effects of round bilge vortex shedding. The difficulty in accounting for such effects in the hybrid method as it stands is due to the fact that separation points are not well defined for rounded bodies. A practical method for the determination of the separation points in general rounded surfaces would be most useful. This would in turn remove the limitation of the present hybrid method to sharp edge vortex shedding. Finally, model testing for the yawed Wigley hull would be extremely useful for the validation of the results given in Chapter 5. Extension of the present work for use in general hull forms would necessarily be the next stage of the research. This would, however, require computer routines for the determination of hull geometrical parameters such as slopes and normals at any desired transverse section. 85 Nomenclature AP B Ca Cd CF C DIDt d 1 F FP FR G g H 2 h K L L I_c L. M NS n p R r,R 0 S S 8 S T t U Vo w x Y y Z z After perpendicular of the ship Beam of the ship Vortex induced inertia coefficient Vortex induced drag coefficient Vortex induced force coefficient Pressure coefficient p/(.5pU ) 2 Material derivative of a variable Draught of the ship Densimetric Froude number relating to stratified flow problem Forward perpendicular of the ship Ship Froude number Free space Green function ln(r) Gravitational constant, 9.81 rn/s 2 Water depth Depth of the upper layer in two layer fluid Keulegan-Carpenter number Length of the ship Length scale in the physical plane in discrete vortex method Length scale in the transformed plane in discrete vortex method Strength scaling in the physical plane in discrete vortex niethod Moments about the axes x,y,z (i = 1,2,3) Number of steps to march over the length of the hull Normal vector of a point Pressure from the Bernoulli equation Wavemaking resistance in newtons Distance between a source point and a field point Bottom boundary in wavemaker problem Far wall in wavemaker problem Free surface boundary in wavemaker problem Body boundary in wavemaker problem Period of oscillatory flow in seconds Dimensional time in seconds Uniform advance velocity of the ship Free stream velocity amplitude in discrete vortex method Complex velocity of a discrete vortex Longitudinal axis in coordinate system fixed on the ship Offset distance of the body surface from the centreline Lateral axis in coordinate system fixed on the ship Complex variable in physical plane in discrete vortex method Vertical axis in coordinate system fixed on the ship 86 a o E Op AS,OS t,Ot p 0 Yaw angle of the ship in degrees Internal angle of infinite wedge in degrees Internal angle of wedge in discrete vortex method Slenderness parameter B/L Difference in density between the upper and lower layers Length of a boundary element Time step size in seconds Elevation of the free surface from the undisturbed level Perturbed velocity potential Total velocity potential in advancing ship problem Wedge angle parameter in discrete vortex method Fluid density Relative angle between a discrete vortex and a nodal point Non-dimensional time in discrete vortex method Complex variable in transformed plane in discrete vortex method 87 Bibliography Aanesland V., “A hybrid model for calculating wave-making resistance”, Proceedings of the 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 in space and a semi-Langrangian semi-implicit scheme in time”, PhD Thesis, Mechanical Engineering Department, University of British Columbia, Vancouver, 1993. Bradbury W.M.S., “A Bilge Vortex Calculation”, Transactions of the Royal institute of NavalArchitects, pp.189-199, 1986. Braathen A. and Faltinsen O.M., “Application of a vortex tracking method to roll damping”, Advances in Underwater Technology, Ocean Science and Offshore Engineering, 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 of Ship Research, vol.33, no.1, pp.21-28, 1989. Calisal S.M., Goren 0. and Okan B., “On an Iterative Solution for Nonlinear Wave Calculations”, Journal ofShip Research, vol. 35, no.1, pp. 944, 1991. Chan R.K.C., “Two-dimensional time-dependent calculations of large-amplitude surface gravity waves due to a surface disturbance”, Proceedings of the First International Conference on Numerical Ship Hydrodynamics, pp. 315-331, 1975. Cheng B.H., “Computation of 3D Transom Stern Flows”, Proceedings of the Fifth International 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 discrete vortices”, 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 Fifth International Conference on Numerical Ship Hydrodynamics, Japan, pp.239-250, 1989. 88 Dommermuth D.G. and Yue K.P., “Numerical simulations of nonlinear axisymmetric flows with a free surface”,Journal of Fluid Mechanics, vol.178, pp.195-219, 1987. Downie M.J., Bearman P.W. and Graham J.M.R., “Effect of vortex shedding on the coupled roll response of bodies in waves”, Journal of Fluid Mechanics, vol.189, pp.243264, 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 dimensional flow”, Journal ofShip Research, vol.22, no.3, pp. 1 93-202, 1978. Fink P.T. and Soh W.K., “Calculation of vortex sheets in unsteady flow and applications in 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 Full Forms”, 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 Keulegan Carpenter 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 boundary element solution of nonlinear wave flows”, Engineering Analysis with Boundary Elements, vol.7, no.4, pp.178-195, 1990. Grilli S.T., “Wave motion and overturning induced by moving bodies: Application to slender ship wave resistance”, Proceedings of the 1st International Conference on Computational 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 ship model”, 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-dimensional bodies by time domain method”, Applied Ocean Research, vol.13, no.3, 1991. 89 Isaacson M. and Ng J.Y.T., “Time-domain second-order wave radiation in two dime nsions”,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 ship wave resistance in shallow water by fully nonlinear wave theory”, Proceedings of the Second Japan-Korea Joint Workshop on Ship and Marine Hydrodynamics, Japan, pp.111120, 1993. Kudo K., “Hydrodynamic forces of the oscillating flat plate”, Proceedings of the 3rd International 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 on water, 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 Supersonic Case”, Proceedings of the 11th International Conference on Offshore Mechanics and Arctic Engineering, OMAE 92, vol. 1A, pp. 5 1-57, 1992. Maruo H., “New Approach to the theory of slender ships with forward velocity”, Bulletin of the Faculty of Engineering., Yokohama National University, Japan, vol. 31, 1982. Maruo H. and Song W., “Numerical appraisal of the new slender ship formulation in steady 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 International Union of Theoretical and Applied Mechanics, Karlsruhe, 1979. Nakayama T. and Washizu K., “The boundary element method applied to the analysis of two dimensional nonlinear sloshing problems”, Journal of Numerical Methods in Eng, 17, 1631-1646, 1981. 90 Ogilvie T.F., “The wave generated by a fine ship bow”, Report No.127, Department of Naval Architecture and Marine Engineering, University of Michigan, 1972. Orlanski L., “A simple boundary condition for unbounded hyperbolic flows”, Journal of Computational Physics, vol.21, pp. 25 1-269, 1976. Pawlowski J.S., Baddour R.E. and Hookey NA., “The development of second generation direct 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 Moment concerning the Rankine solid under Interfacial Wave Conditions”, Norwegian Ship Model Experiment 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 Wave Pattern by the Slender Body Approximation”, Journal of the Kansai Society ofNaval Architects, Japan, no.209, pp 25-36, 1988 (In Japanese). Song W. and Maruo H., “Bow impact and Deck Wetness: Simulation Based on Nonlinear Slender Body Theory”, Proceedings of the 3rd International Offshore and Polar Eng Conference, Singapore, vol.111, pp.34-38, 1993. Suzuki K., “Boundary Integral Equation Method for the Linear Wave Resistance Problem”, Proceedings of the 4th International Conference on Numerical Ship Hydrodynamics, pp.308-323, 1985. Tahara Y., “A Boundary Element Method for Calculating Free Surface flows around a Yawed Ship”, Private communication with Prof K Suzuki of Yokohama National University, Japan, 1992. Tanaka Y. and Nakamura T., “Open boundary scheme for nonlinear irregular waves”, Proceedings of the 11th International Conference on Offshore Mechanics and Arctic Engineering, OMAE 92, vol.1A, pp.34.5-351, 1992. Tryggvason G., Unverdi S.O., Song M. and Abdollahi-Alibeik J., “Interaction of vortices with 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 Supersonic Case”, Proceedings of the 18th Symposium on Naval Hydrodynamics, pp.567-581, 1990. 91 Vinje T. and Brevig P., “Nonlinear ship motions”, Proceedings of the 3rd International conference 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 a ship 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 in oscillatory flow”, M.A.Sc. Thesis, Department of Mechanical Engineering, University of British 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 the effects 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.514518, 1993(b). Yen S.M. and Hall D.R., “Implementation of open boundary conditions for nonlinear freesurface wave problems”, Proceedings of the 3rd International Conference on Numerical Ship Hydrodynamics, pp.163-177, Paris, 1981. Yen S.M., Lee K.D. and Akai T.J., “Finite element and finite differences solution of nonlinear free surface wave problems”, Proceedings of the 2nd International Conference on Numerical Ship Hydrodynamics, pp.305-318, Berkeley, 1977. 92 Appendix A Slender body theory and the wavemaker Notation Length L Beam B [au9ht d AP Free stream velocity (U Note : Origin is at undisturbed free surface level. Figure A.1 : Frame of reference and notation. The following analysis is given to show how the three dimensional problem of a slender body advancing at a velocity U is approximated by a series of two dimensional problems. The total potential is given by: (A.1) P= Uxcosa÷Uysinct+4)(x,y,z) where the free stream velocity U is incident on the body at an angle a and 4) is the perturbation potential. Relevant variables are non-dimensionalized (the non-dimensional form being denoted by a 4) — = (-), such that: UB4) Bj 93 Y=BY x=L y = B z = B (A.2) where r is the wave elevation and Y describes the body boundary. The ship beam and length 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 about spatial derivatives as well, it is necessary to express these as: a_i L — a a a_i a B a3 — ôz (A.3) Baz The fluid is taken to be inviscid and incompressible and flow is assumed irrotational. With these assumptions, the governing Laplace equation is: 2 ox where 2 0y (A.4) 2 0z is the velocity potential. From Equations (A.2) and (A.3) and assuming that derivatives with respect to non-dimensional variables do not change in order of magnitude, the Laplace equation can be written: + + (A.5) =0 94 If the slenderness ratio B / L = e is small (limited to about 0.15), the 3D Laplace equation can be expressed in the two dimensional form in the plane of the hull sections: + = 0 (A.6) The 3D impermeable hull boundary is: 1 +(1).,n +(I)n =0 ‘1),n (A.7) In non-dimensional form this is rewritten as: —cosaY —-4Y ±(sinct + = 0 (A.8) from which the term of order c 2 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 then eliminating the term in ): (e cosa +e2 (Ecosa)j + (na + +(sinct+)j —p — = = 0 (A.10) 0 The non-dimensional form of the dynamic free surface condition for steady potential flow, eliminating the terms in 2, --[i +i i-t is given in the development of Equation (A. 11) below. ]+g9=o 2 _u 1 +2Ucosa €cosa +sina +4 +2Usina +4j+g=o ÷‘[E22 +2j+j 95 = 0 e cos + sin ![ + j+ + =0 (A. ii) To maintain the last term in the dynamic free surface condition, it is necessary to have 0(e). Here, the draught (d) at which the body floats is taken to be of the O(gB I U) order of the beam (B). (Fd = Therefore in this application, the draught Froude number U/.J) is taken to be An expression can be obtained, using the slender body approximation, for the normal 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: (n,n,n)= where Y = —YZ’ +1YZ .J1 + 0(e) and Y + = (A.13) 2 Y 0(1). , the normal velocity 3 Elimination of the term Y (positive pointing out) of the body in the cross-flow plane can be written, using Equations (A.13) and (A.14): = + an u(Y cosa ± sin a) 2 .iJl+Yz on the hull (A.14) Putting Equations (A.5), (A.10) and (A.11) back into dimensional form and rearranging, the following can be obtained: + 4 = 0 Ucosat +(usina+p)q —• =0 96 in the domain A.15) at (A.16) 1 cosa+p na)÷--[p u(p +p]+g’j =0 at z= i (A.17) This set of equations (A.14) to (A.17), together with impermeable bottom and wall boundary conditions, are identical to those used for Ihe numerical modeling of a two dimensional wave tank in the time domain. The only difference here is that the motion of the wave paddle is replaced by the geometric change of the ship hull in the longitudinal direction. 97 Appendix B Higher order boundary element method The following description on the boundary element method applies to two dimensional potential flow problems only. We start with Green’s second formula: (B.1) —GdS=O S on an where the weighting source function is taken to be: G = (B.2) --ln( 2t ‘sr! Both G and i satisfy the governing Laplace equation. The boundary can be discretized into N number of elements. Collocation points are placed at the ends of linear elements, where the variation of the potential and its normal derivative over the elements are assumed linear. If the potential and its normal derivative are taken to be varying quadratically over an element, then the collocation points (nodes) are placed at both the ends and at the middle of an element. Substituting the normal derivative of G with H, the potential 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 that: 98 into the summation sign such fdljHujdE’ (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 a given point on the element j. Dirichiet or Neumann boundary conditions are specified at each node. All unknowns in Equation (B.4) can then be gathered such that the system of N equations in N unknowns to be solved for X can be written as (for more details see Brebbia and Domnguez, 1989): (B.5) [AJX]=[Bj The evaluation of the matrices [A] and [B] involves numerical integration for the G and H terms in Equation (B 4). For collocation nodes located between two elements, the contributions to the integrals in G and H come from the two adjacent elements. xiln()dr+----f =——f r 2t’ H..“ = 2t ‘ x —-Iln(-ldF 1 on L \rJJ + lX 2 11()dF r 231 F Equations (B.7) and (B.8) for the value of ln(ldF X2 4 On L \rJJ (B.6) apply to linear and quadratic elements respectively. Xi (B.7) X2(11) Xi 2 (B.8) X2 99 In 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 term where: XL =(‘—X’÷) X2 = (B.9) 0 The integrals in Equation (B.6) are evaluated using standard Gaussian and logarithmic numerical integration quadratures (see for example, Brebbia and Dorninguez, 1989). While potential values are taken to be continuous at the collocation nodes, normal velocities need not be so. This is evident, for example, at the intersection between the free surface and a body or wall. The normal velocities immediately on the two sides of such nodes 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 only to G for the free surface normal velocity and the adjacent body element contributes only to G for the body side normal velocity. Finally, the procedure for the extrapolation of the normal velocity on the free surface side of the fluid-body intersection point requires some explanation. Linear variation of the normal velocity over the three free surface nodes (see Figure B. 1) adjacent to the fluid-body intersection point is assumed. Hence the normal velocity at the intersection point is given by: (4 \iln!A = 2( ‘\ônJB (B.10) — \Jc 100 4,n Free surface e’ements Node A Node B Node C Figure B.1 : Intersection point and normal velocities. The above eliminates the need to evaluate the normal velocity for the free surface side of the node at the intersection point. There, the velocity potential is then taken to be the unknown. Another consequence of the extrapolation is that G integrals at points B and C are modified follows: 018 = B + (B.11) 101
- 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 |
File Format | 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 |
Graduation Date | 1994-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | 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:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0080831/manifest